Python遥感可视化 — Basemap将地面观测站点进行空间插值可视化

欢迎关注博主的微信公众号：“智能遥感”。

该公众号将为您奉上Python地学分析、爬虫、数据分析、Web开发、机器学习、深度学习等热门源代码。

本人的GitHub代码资料主页（持续更新中，多给Star，多Fork）：

/xbr

CSDN也在同步更新：

/XBR_

“本节主要内容是将地面观测站点数据在地图上显示，并且对站点数据进行空间插值，从而可以将点数据转化为栅格数据，栅格数据具有较好的空间连续性。”

今天的遥感之美封面图—白雪下的麦肯齐河系统。乍一看，有一种山回路转不见君，雪上空留马行处的景象。麦肯齐河系统是加拿大最大的流域，也是世界上第十大流域。这条河从加拿大落基山脉的哥伦比亚冰原到达北冰洋，全长4,200公里。

11月7日由Landsat 8上OLI传感器拍摄的图像，该图显示了Mackenzie River Delta的一部分 - 包括沿着河道东部的这条冰冻高速公路。由白色冰雪覆盖的水道在绿色松树覆盖的土地上脱颖而出。

PM2.5可视化

在平时遥感数据可视化时，等经纬度投影是我们最长使用的一种方法，也是最简单的投影。等经纬度投影又称等距圆柱投影，是假想球面与圆筒面相切于赤道，赤道为没有变形的线。经纬线网格，同一般正轴圆柱投影，经纬线投影成两组相互垂直的平行直线。其特性是：保持经距和纬距相等，经纬线成正方形网格；沿经线方向无长度变形；角度和面积等变形线与纬线平行，变形值由赤道向高纬逐渐增大。该投影适合于低纬地区制图。

首先来看一下Basemap官网上的参考代码，这是最简单的显示全球地图的方式（图1）：大陆廓线（未添加国家行政边界）以及陆地（橘黄色）、海洋（浅绿色）、内陆湖泊（蓝色）的颜色填充，让人更加直观地了解各大洲的分布情况。

代码实现：

#_*_coding:utf-8_*___author__='xbr'__date__='/1/1114:09'frommpl_toolkits.basemapimportBasemapimportmatplotlib.pyplotaspltmaps=Basemap(projection='cyl',lat_0=0,lon_0=0)#首先给地球涂上蓝色maps.drawmapboundary(fill_color='aqua')#再给大陆涂上橘黄色,给江河湖泊涂上大海一样的颜色maps.fillcontinents(color='coral',lake_color='blue')maps.drawcoastlines()plt.show()#plt.savefig('test.png')

结果图：

图1等经纬度投影下的世界地图

在上述基础之上，我们将地面观测站点数据叠加到底图上，来分析数据的时空变化。下面以01月01日中国地区PM2.5观测数据为例，图2展示全国PM2.5站点空间分布，我们发现中国东部地区观测站点数量加多，并且点相对较大，点的大小对应PM2.5数值的高低，说明东部地区这一天的大气污染较西部地区较为严重。

图2全国PM2.5站点分布，点越大说明PM2.5浓度越高

但是上图不能够很好地反应PM2.5的空间变化，从视觉效果来看，颜色较为单一。因此，我们给每个站点添加颜色，颜色越红，表示PM2.5浓度越大，也就是空气质量越差（见下图）。点的颜色和大小可以同时表示数值大小，更加增强了数据的可视化效果。

图3全国PM2.5站点分布，颜色越红且点越大说明PM2.5浓度越高

代码实现：

#_*_coding:utf-8_*___author__='xbr'__date__='/1/1114:49'frommpl_toolkits.basemapimportBasemapimportmatplotlib.pyplotaspltimportpandasaspdimportnumpyasnp#用来正常显示中文plt.rcParams['font.sans-serif']=['SimHei']#用来正常显示负号plt.rcParams['axes.unicode_minus']=False#获取PM2.5数据df=pd.read_excel(r'D:\data\0101PM25-CHINA.xlsx')#剔除无效值NANdf=df.dropna(axis=0,how='any')lat=np.array(df["lat"][:])#获取维度之维度值lon=np.array(df["lon"][:])#获取经度值PM25=np.array(df["PM25"][:],dtype=float)#画图fig=plt.figure(figsize=(16,9))plt.rc('font',size=15,weight='bold')ax=fig.add_subplot(111)#添加标题,PM2.5下标设置plt.title(u'01月01日中国地区$\mathrm{PM}_{2.5}$质量浓度分布',size=25,weight='bold')#创建底图,等经纬度投影mp=Basemap(llcrnrlon=73.,llcrnrlat=17.,urcrnrlon=135.,urcrnrlat=55,projection='cyl',resolution='h')#添加海岸线mp.drawcoastlines()#添加国家行政边界mp.drawcountries()#设置colorbar中颜色间隔个数levels=np.linspace(0,np.max(PM25),20)#cf=mp.scatter(lon,lat,PM25,marker='o',color='r')#设置颜色表示数值大小cf=mp.scatter(lon,lat,PM25,c=PM25,cmap='jet',alpha=0.75)#设置上下标以及单位的希腊字母cbar=mp.colorbar(cf,location='right',format='%d',size=0.3,ticks=np.linspace(0,np.max(PM25),10),label='$\mathrm{PM}_{2.5}$($\mu$g/$\mathrm{m}^{3}$)')plt.show()

地面站点数据观测精度相对较高，但是不能够覆盖所有的空间范围。为了研究大范围连续空间的特征变化，通常需要对站点数据进行空间插值，例如克里金插值、样条插值、反距离加权。下面以三种插值为例（最邻近插值、线性插值和三次样条插值）。如果喜欢克里金插值的小伙伴，可以pykrige模块去实现，本节不作讲解。

图4最邻近插值后，全国PM2.5空间分布，颜色越红说明PM2.5浓度越高

图5线性插值后，全国PM2.5空间分布，颜色越红说明PM2.5浓度越高

图6三次样条插值后，全国PM2.5空间分布，颜色越红说明PM2.5浓度越高

说明：三次样条插值后的效果其实分布也是连续的，与图5很相似，但是在Basemap上显示就出现破碎现象，目前还不知道什么原因。

代码实现：

#_*_coding:utf-8_*___author__='xbr'__date__='/1/1114:49'frommpl_toolkits.basemapimportBasemapimportmatplotlib.pyplotaspltimportpandasaspdimportnumpyasnpfromscipy.interpolateimportgriddataimportcmaps#用来正常显示中文plt.rcParams['font.sans-serif']=['SimHei']#用来正常显示负号plt.rcParams['axes.unicode_minus']=False#获取PM2.5数据df=pd.read_excel(r'D:\data\0101PM25-CHINA.xlsx')#剔除无效值NANdf=df.dropna(axis=0,how='any')#获取纬度lat=np.array(df["lat"][:])#获取经度lon=np.array(df["lon"][:])#获取PM2.5PM25=np.array(df["PM25"][:])#创建格网grid_x,grid_y=np.mgrid[73.56:135.04,18.2:53.56]#插值方法：'nearest','linear','cubic'grid_z=griddata((lon,lat),PM25,(grid_x,grid_y),method='cubic')#画图fig=plt.figure(figsize=(16,9))plt.rc('font',size=15,weight='bold')ax=fig.add_subplot(111)#添加标题,PM2.5下标设置plt.title(u'01月01日中国地区$\mathrm{PM}_{2.5}$质量浓度分布',size=25,weight='bold')#创建底图,等经纬度投影mp=Basemap(llcrnrlon=73.,llcrnrlat=17.,urcrnrlon=135.,urcrnrlat=55,projection='cyl',resolution='h')#添加海岸线mp.drawcoastlines()#添加国家行政边界mp.drawcountries()#设置colorbar中颜色间隔个数levels=np.linspace(0,np.max(PM25),20)cf=mp.contourf(grid_x,grid_y,grid_z,levels=levels,cmap=cmaps.GMT_panoply)cbar=mp.colorbar(cf,location='right',format='%d',size=0.3,ticks=np.linspace(0,np.max(PM25),10),label='$\mathrm{PM}_{2.5}$($\mu$g/$\mathrm{m}^{3}$)')plt.show()