Python对风场进行显著性检验(文末免费包邮送书)

import xarray as xrimport numpy as npimport datetime as dtimport cartopy.crs as ccrsimport cartopy.feature as cfeatureimport cartopy.mpl.ticker as ctickerimport matplotlib.pyplot as pltfrom scipy.stats import linregressfrom cartopy.io.shapereader import Readerfrom matplotlib import rcParamsconfig = {"font.family":'Times New Roman',"font.size": 16,"mathtext.fontset":'stix'}rcParams.update(config)f_t = xr.open_dataset('./air.mon.mean.nc')# t_nea = np.array(f_t.air.loc[f_t.time.dt.month.isin([12,1,2])].loc['1979-12-01':'2020-03-01',850,50:30,110:130].mean('lat').mean('lon')).reshape(41,3).mean((1))t = np.array(f_t.air.loc[f_t.time.dt.month.isin([12,1,2])].loc['1979-12-01':'2020-03-01',850,50:30,110:130]).mean((1,2)).reshape(41,3).mean((1))f_u = xr.open_dataset('./uwnd.mon.mean.nc')u = np.array(f_u.uwnd.loc[f_u.time.dt.month.isin([12,1,2])].loc['1979-12-01':'2020-03-01',500]).reshape(41,3,73,144).mean((1))f_v = xr.open_dataset('./vwnd.mon.mean.nc')v = np.array(f_v.vwnd.loc[f_v.time.dt.month.isin([12,1,2])].loc['1979-12-01':'2020-03-01',500]).reshape(41,3,73,144).mean((1))lat = f_u.latlon = f_u.lonprint(t.shape,u.shape)su,pu = np.zeros((73,144)),np.zeros((73,144))sv,pv = np.zeros((73,144)),np.zeros((73,144))for i in range(len(lat)):for j in range(len(lon)):        su[i,j],_,_, pu[i,j],_  = linregress(t,u[:,i,j])        sv[i,j],_,_, pv[i,j],_  = linregress(t,v[:,i,j])fig = plt.figure(1, figsize=[16,12])proj=ccrs.PlateCarree()ax = plt.subplot(1, 1, 1, projection=ccrs.PlateCarree())extent = [0,180,0,90]leftlon, rightlon, lowerlat, upperlat = (0,180,0,90)img_extent = [leftlon, rightlon, lowerlat, upperlat]ax.set_extent(extent, crs=ccrs.PlateCarree())ax.stock_img()ax.add_feature(cfeature.COASTLINE.with_scale('50m'),linewidth=0.5,zorder=2,color='k')# 添加海岸线ax.add_feature(cfeature.LAKES.with_scale('50m'))ax.add_feature(cfeature.RIVERS.with_scale('50m'))ax.add_feature(cfeature.OCEAN.with_scale('50m'))ax.add_feature(cfeature.LAND.with_scale('50m'))ax.add_geometries(Reader(r'./henan1.shp').geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor='k',linewidth=1.2,zorder=1)ax.add_geometries(Reader(r'./river1.shp').geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor='b',linewidth=0.2)ax.add_geometries(Reader(r'./china1.shp').geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor='k',linewidth=1.2)ax.add_geometries(Reader(r'./china2.shp').geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor='k',linewidth=0.8)ax.add_geometries(Reader(r'./ne_50m_lakes.shp').geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor='k',linewidth=0.2)ax.add_geometries(Reader(r'./1级河流.shp').geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor='RoyalBlue',linewidth=0.4)ax.add_geometries(Reader(r'./2级河流.shp').geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor='DodgerBlue',linewidth=0.3)ax.add_geometries(Reader(r'./3级河流.shp').geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor='DeepSkyBlue',linewidth=0.2)ax.add_geometries(Reader(r'./4级河流.shp').geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor='SkyBlue',linewidth=0.15)ax.add_geometries(Reader(r'./5级河流.shp').geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor='LightSkyBlue',linewidth=0.05)ax.add_geometries(Reader(r'./主要湖泊.shp').geometries(),ccrs.PlateCarree(),edgecolor='none',linewidth=0,facecolor='#BEE8FF')ax.add_geometries(Reader(r'./country1.dbf').geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor='k',linewidth=1.2)ax.add_geometries(Reader(r'./bou2_4l.dbf').geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor='r',linewidth=1.2)ax.set_xticks(np.arange(leftlon,rightlon+30,30), crs=ccrs.PlateCarree())ax.set_yticks(np.arange(lowerlat,upperlat+30,30), crs=ccrs.PlateCarree())lon_formatter = cticker.LongitudeFormatter()lat_formatter = cticker.LatitudeFormatter()ax.xaxis.set_major_formatter(lon_formatter)ax.yaxis.set_major_formatter(lat_formatter)ax.quiver(lon[::2],lat[::2],su[::2,::2],sv[::2,::2],transform=ccrs.PlateCarree(),scale=60,cmap=plt.cm.jet)plt.savefig('./plot295.1.png',dpi=600,bbox_inches='tight',pad_inches=0)plt.show()

su_tmp  = su.copy()sv_tmp  = sv.copy()su_tmp[pu>0.1]=np.nansu_tmp[pv>0.1]=np.nansv_tmp[pu>0.1]=np.nansv_tmp[pv>0.1]=np.nanax.quiver(lon[::2],lat[::2],su_tmp[::2,::2],sv_tmp[::2,::2],transform=ccrs.PlateCarree(),scale=60,cmap='gist_rainbow')plt.savefig('./plot295.2.png',dpi=600,bbox_inches='tight',pad_inches=0)plt.show()