python |
您所在的位置:网站首页 › cartoee › python |
使用python对海洋气象数据做显著性检验,并绘制空间pattern
选择数据集: 1 SST (Daily Sea Surface Temperature) NOAA High-resolution Blended Analysis daily分辨率:2.5时间:2010http://www.esrl.noaa.gov/psd/data/gridded/data.noaa.oisst.v… 2 OLR (Outgoing Longwave Radiation)- NOAA Interpolated OLRdaily分辨率:2.5时间覆盖范围:1974-2013http://www.esrl.noaa.gov/psd/data/gridded/data.interp_OLR.html使用编程工具: python主要使用函数: stats.linregress() 函数说明点这里内容: 对2010年的SST和OLR数据分布进行显著性检验,并绘制空间分布先各自对数据计算一年的趋势以及相关,再计算两个数据之间的相关绘制空间pattern具体过程: 1 导入库2 读取数据3 计算相关和趋势4 绘图5 保存数据过程还是比较清晰的,直接附上结果, ps(标题时间打错了、横轴的label也搞错了,懒得重画了,代码中应该没问题了) :
最好,还是附上完整的代码吧,也没啥好保留的: # -*- coding: utf-8 -*- """ Created on Sat Oct 22 20:50:09 2022 @author: Administrator """ import cartopy.feature as cfeature import numpy as np import xarray as xr from cartopy.mpl.ticker import LongitudeFormatter,LatitudeFormatter import cmaps import matplotlib.pyplot as plt import cartopy.crs as ccrs import matplotlib.ticker as mticker from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER from scipy import stats ############################################################################### olr_path = r'J:/olr.day.mean.nc' sst_path = r'J:/sst.intep.nc' da1 = xr.open_dataset(sst_path).sortby('lat') da2 = xr.open_dataset(olr_path).sortby('lat') sst = da1.sst.sel(lat=slice(-30,30),lon=slice(100,200)) olr = da2.olr.sel(lat=slice(-30,30),lon=slice(100,200), time=slice('2010','2010')) ############### calculate ################################################ trend = np.zeros((sst.lat.shape[0],sst.lon.shape[0])) p_value = np.zeros((sst.lat.shape[0],sst.lon.shape[0])) for i in range (0,sst.lat.shape[0]): for j in range (0,sst.lon.shape[0]): trend[i,j], intercept, r_value, p_value[i,j], std_err=stats.linregress(np.arange(1,366),sst[:,i,j]) ############################################################################## ################ plot ####################################################### lon = sst.lon.data lat = sst.lat.data ############################################################################## box = [100,200,-20,20] xstep,ystep = 20,10 proj = ccrs.PlateCarree(central_longitude=180) plt.rcParams['font.family'] = 'Times New Roman', ############################################################################## fig = plt.figure(figsize=(8,7),dpi=200) fig.tight_layout() ax = fig.add_axes([0.1,0.2,0.8,0.7],projection = proj) ax.set_extent(box,crs=ccrs.PlateCarree()) ax.coastlines('50m') ax.set_xticks(np.arange(box[0],box[1]+xstep, xstep),crs=ccrs.PlateCarree()) ax.set_yticks(np.arange(box[2], box[3]+1, ystep),crs=ccrs.PlateCarree()) lon_formatter = LongitudeFormatter(zero_direction_label=False)#True/False lat_formatter = LatitudeFormatter() ax.xaxis.set_major_formatter(lon_formatter) ax.yaxis.set_major_formatter(lat_formatter) ax.spines[['right','left','top','bottom']].set_linewidth(1.1) ax.spines[['right','left','top','bottom']].set_visible(True) ax.set_xlabel('Lon',fontsize=14) ax.set_title('Significance Test',fontsize=16,pad=8,loc='right') ax.set_title('this is title',fontsize=16,pad=8,loc='left') ax.tick_params( which='both',direction='in', width=0.7, pad=5, labelsize=14, bottom=True, left=True, right=True, top=True) c = ax.contourf(lon,lat,trend, # levels=np.linspace(-0.008,0.008,17), # levels=np.linspace(-0.32,0.32,17), #levels=np.linspace(-0.012,0.012,17), extend = 'both', transform=ccrs.PlateCarree(), cmap=cmaps.BlueWhiteOrangeRed) c1b = ax.contourf(lon,lat, p_value, [np.nanmin(p_value),0.05,np.nanmax(p_value)], hatches=['.', None], colors="none", transform=ccrs.PlateCarree()) cb=plt.colorbar(c, shrink=0.85, pad=0.15, orientation='horizontal', aspect=25, ) cb.ax.tick_params(labelsize=10,which='both',direction='in',) plt.show() # save picture # fig.savefig(r'D:\SST_OLR.png',format='png',dpi=500) ############################################################################## # xtick = np.arange(100, 200, 20) # ytick = np.arange(-30,31, 10) # gl=ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True, # xlocs=[120,140,160,180,], # ylocs=[-30,-20,-10,0,10,20,30,], # x_inline=False, # y_inline=False, # linewidth=1.5, color='gray', alpha=0.5, linestyle='--', # ) # gl.xlabels_top = False # gl.ylabels_right = False # gl.xlines = True # plt.xlabel('Lon',fontsize=15) # plt.ylabel('Lat',fontsize=15) # gl.xformatter = LONGITUDE_FORMATTER # gl.yformatter = LATITUDE_FORMATTER |
CopyRight 2018-2019 办公设备维修网 版权所有 豫ICP备15022753号-3 |