python由站点数据插值到网格数据方法对比

见贤思齐

文章目录1、前言2、结果对比2.1原始散点站位图2.2griddata插值2.3krige插值2.4RBF插值2.5IDW插值3、总结4、其它5、国家站能见度插值结果对比1、前言气象海洋中空间数据类型有站点数据、格点数据。站点数据空间分布不连续，不利于进行时空分析；有时需要将站点数据插值到网格中。本文基于一个散点数据对比了griddata、IDW、krige、rbf的插值结果。interp2d插值，地理数据中的分辨率转换2、结果对比2.1原始散点站位图原始散点站位如下：2.2griddata插值fromscipy.interpolateimportgriddatalat=df_filter['Lat']lon=df_filter['Lon']VIS_Min=df_filter['VIS_Min']olon=np.linspace(lon.min(),lon.max(),40)olat=np.linspace(lat.min(),lat.max(),40)olon,olat=np.meshgrid(olon,olat)VIS_Min_grd=griddata((lon,lat),VIS_Min,(olon,olat),method='linear')plt.figure(figsize=(8,6))plt.pcolormesh(olon,olat,VIS_Min_grd)plt.scatter(lon,lat,s=50,c=VIS_Min,cmap='viridis',edgecolors='k')1234567891011可视化代码importosimportnumpyasnpimportmatplotlib.pyplotaspltimportmatplotlib.colorsasmcolorsimportmatplotlibasmplimportcartopy.crsasccrsimportcartopy.io.shapereaderasshpreaderfrommpl_toolkits.axes_grid1.inset_locatorimportinset_axesfrommatplotlib.colorsimportLinearSegmentedColormapimportregionmaskimportgeopandasasgpdimportwarningswarnings.filterwarnings("ignore")importmatplotlib.tickerastickerextents=[115.6,117.2,39.7,40.8]crs=ccrs.PlateCarree()fig=plt.figure(figsize=(12,6))ax=fig.add_subplot(111,projection=crs)ax.set_extent(extents,crs)prov_reader=shpreader.Reader(r'.\bou2_4p.shp')nineline_reader=shpreader.Reader(r'.\九段线.shp')ax.add_geometries(prov_reader.geometries(),crs,edgecolor='gray',facecolor='none',lw=1)ax.add_geometries(nineline_reader.geometries(),crs,edgecolor='gray',facecolor='none',lw=1)ax.set_xticks(np.arange(extents[0],extents[1]+0.3,0.3))ax.set_yticks(np.arange(extents[2],extents[3]+0.2,0.2))plt.tick_params(axis='both',which='major',labelsize=18)norm=mcolors.Normalize(vmin=0,vmax=np.round(VIS_Min.max()))cmap_jet=plt.cm.jetcmap_gray=plt.cm.graygray=cmap_gray(np.linspace(0,1,9))colors=np.vstack((gray[8],cmap_jet(np.linspace(0,1,9))))cmaps=LinearSegmentedColormap.from_list('Combined',colors,N=64)map=plt.pcolormesh(olon,olat,VIS_Min_grd,cmap=cmaps)plt.scatter(lon,lat,s=50,c=VIS_Min,cmap=cmaps,edgecolors='k',norm=norm)cb_ax=inset_axes(ax,width="3%",height="100%",loc='lowerleft',bbox_to_anchor=(1.01,0.,1,1),bbox_transform=ax.transAxes,borderpad=0)cbar=plt.colorbar(mpl.cm.ScalarMappable(norm=norm,cmap=cmaps),cax=cb_ax,ax=map)cbar.ax.set_title('能见度',size=18,fontname='SimHei',pad=10)cbar.ax.tick_params(labelsize=15)ax.set_xlabel('Lon(°E)',fontsize=18)ax.set_ylabel('Lat(°N)',fontsize=18)ax.set_title('原始散点站位',size=18,fontname='SimHei',pad=10)123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051可视化结果2.3krige插值插值核心代码frompykrige.okimportOrdinaryKriginglat=df_filter['Lat']lon=df_filter['Lon']VIS_Min=df_filter['VIS_Min']olon=np.linspace(lon.min(),lon.max(),40)olat=np.linspace(lat.min(),lat.max(),40)Krin=OrdinaryKriging(lon,lat,VIS_Min,variogram_model='gaussian',coordinates_type="geographic",nlags=1,)data2,ssl=Krin.execute('grid',olon,olat)olon,olat=np.meshgrid(olon,olat)plt.figure(figsize=(8,6))plt.pcolormesh(olon,olat,data2)plt.scatter(lon,lat,s=50,c=VIS_Min,cmap='viridis',edgecolors='k')12345678910111213141516插值结果2.4RBF插值插值核心代码fromscipy.interpolateimportRbflat=df_filter['Lat']lon=df_filter['Lon']VIS_Min=df_filter['VIS_Min']olon=np.linspace(lon.min(),lon.max(),40)olat=np.linspace(lat.min(),lat.max(),40)olon,olat=np.meshgrid(olon,olat)func1=Rbf(lon,lat,VIS_Min,function='linear')VIS_Min_rbf=func1(olon,olat)plt.figure(figsize=(8,6))plt.pcolormesh(olon,olat,VIS_Min_rbf)plt.scatter(lon,lat,s=50,c=VIS_Min,cmap='viridis',edgecolors='k')12345678910111213插值结果2.5IDW插值插值核心代码frommathimportradians,sin,cos,sqrt,asinimportnumpyasnp#degreetokm(Haversinemethod)defhaversine(lon1,lat1,lon2,lat2):R=6372.8#EarthradiusinkilometersdLon=radians(lon2-lon1)dLat=radians(lat2-lat1)lat1=radians(lat1)lat2=radians(lat2)a=sin(dLat/2)**2+cos(lat1)*cos(lat2)*sin(dLon/2)**2c=2*asin(sqrt(a))d=R*creturnddefIDW(x,y,z,xi,yi):lstxyzi=[]forpinrange(len(xi)):lstdist=[]forsinrange(len(x)):d=(haversine(x[s],y[s],xi[p],yi[p]))lstdist.append(d)sumsup=list((1/np.power(lstdist,2)))suminf=np.sum(sumsup)sumsup=np.sum(np.array(sumsup)*np.array(z))u=sumsup/suminfxyzi=[xi[p],yi[p],u]lstxyzi.append(xyzi)return(lstxyzi)lat=df_filter['Lat']lon=df_filter['Lon']VIS_Min=df_filter['VIS_Min']olon=np.linspace(lon.min(),lon.max(),40)olat=np.linspace(lat.min(),lat.max(),40)olon,olat=np.meshgrid(olon,olat)VIS_Min_idw=IDW(lon,lat,VIS_Min,olon.flatten(),olat.flatten())VIS_Min_idw=np.array(VIS_Min_idw)[:,2].reshape(olon.shape)plt.figure(figsize=(8,6))plt.pcolormesh(olon,olat,VIS_Min_idw)plt.scatter(lon,lat,s=50,c=VIS_Min,cmap='viridis',edgecolors='k')123456789101112131415161718192021222324252627282930313233343536373839404142矩阵运算importnumpyasnpdefhaversine(lon1,lat1,lon2,lat2):R=6372.8#EarthradiusinkilometersdLon=np.radians(lon2-lon1)dLat=np.radians(lat2-lat1)lat1=np.radians(lat1)lat2=np.radians(lat2)a=np.sin(dLat/2)**2+np.cos(lat1)*np.cos(lat2)*np.sin(dLon/2)**2c=2*np.arcsin(np.sqrt(a))d=R*creturnddefIDW(x,y,z,xi,yi):xi,yi=np.meshgrid(xi,yi)#Creategridxi=xi.ravel()yi=yi.ravel()dists=haversine(np.tile(x,(len(xi),1)).T,np.tile(y,(len(yi),1)).T,xi,yi)weights=1/dists**2weights[np.isinf(weights)]=1e12#Handledivisionbyzeroweighted_sum=np.sum(weights*z[:,None],axis=0)sum_of_weights=np.sum(weights,axis=0)u=weighted_sum/sum_of_weightsresult=np.column_stack((xi,yi,u))returnresult#示例数据x=np.array([0,1,2,3,4])y=np.array([0,1,2,3,4])z=np.array([1,2,3,4,5])xi=np.linspace(0,4,10)yi=np.linspace(0,4,10)#进行IDW插值interpolated_data=IDW(x,y,z,xi,yi)print(interpolated_data[:,2])123456789101112131415161718192021222324252627282930313233343536373839404142插值结果3、总结插值范围上，除了griddata以外的其它插值方法，都可以插值到站点以外的区域；插值效果上，从本案例的结果上看，我认为反距离权重IDW结果最好，之后为RBF，其次是克里金。4、其它在数据插值过程中，发现了MetPy库下的inverse_distance_to_grid插值函数。用于执行基于反距离权重（InverseDistanceWeighting,IDW）的插值算法，常用于气象学和其他领域的空间数据分析。这个函数允许你将一系列离散观测点的数据插值到一个网格上，从而生成连续的空间场。使用上述数据，调整了相关参数，做了插值，结果并不理想。frommetpy.interpolateimportinverse_distance_to_griddata=inverse_distance_to_grid(lon,lat,VIS_Min,olon,olat,r=0.5)plt.figure(figsize=(8,6))plt.pcolormesh(olon,olat,data)plt.scatter(lon,lat,s=50,c=VIS_Min,cmap='viridis',edgecolors='k')123455、国家站能见度插值结果对比基于全国2000多各国家站逐小时能见度结果进行插值结果对比。对比方法包括griddata、rbf、IDW插值效率：griddata>rbf>IDW插值效果：rbf>griddata>IDW