本节提要：如何对数据进行地形缺测。
今天摸鱼大佬新录了一期关于地形缺测的内容，那我也顺手补充一下。地形缺测，在大气科学方向，最常用的有两个方向，一个是平面图中对青藏高原的数据进行覆盖；一个是剖面图中进行覆盖。青藏高原这个很容易理解，基本气压下限就是500百帕，常用的700,850层次一般都是补插出来的。剖面图覆盖一般是天气学常用，比如雅安暴雨分析，就经常要补地形，这是因为常用的天气学再分析资料，一般是完整的从1000-1hPa全层次，但是由于地形限制，一般下限在850左右。这个时候，就需要引入绘图缺测。
在之前，学习白化我们就知道，有提前筛选数据与绘图白化两种方法，这里，由于地形的不规则性，我们最好不要使用筛选数据的方法，一般用绘图白化。
平面图的地形缺测，有两种办法，第一种是利用shp文件，将shp叠放到绘制的图片上层，起到遮挡的效果，另一种，是利用绘制等值线填色图生成的pathcollection叠加起遮挡效果。第一种方法比较简便，第二种方法则更精确。比如雅鲁藏布大峡谷，是一个非常深的槽，应该能有700-850层次，但是若使用shp文件直接遮盖，则显然会被缺测。这个是比较简单的，首先将一个850温度场绘制出来。这里，依然使用ncep再分析的气温资料。file=r'E:\aaaa\air.mon.mean.nc'ds=xr.open_dataset(file)lon=ds['lon']lat=ds['lat']time=ds['time']air=ds['air']air=air.loc['1948-01-01',850,:,:]ac=ax.contourf(lon,lat,air, levels=np.arange(-34,38,4), cmap='Spectral_r', transform=ccrs.PlateCarree())cb=fig.colorbar(ac,ax=ax,shrink=0.8,pad=0.05)cb.ax.tick_params(labelsize=6)cb.set_ticks(np.arange(-34,38,4))cb.ax.set_title('气温（℃）',fontsize=8)
然后叠加shp文件。add_geometries的绘图zorder是比contourf高的，所以一般会直接覆盖，但不绝对，若没有正确覆盖，则需要修改层次。这里仅做介绍，后续大家可以改参数来使图片更美观。shp=r'E:\tibetan\tibetan.shp'mapshp=shpreader.Reader(shp)ax.add_geometries(mapshp.geometries(), ccrs.PlateCarree(), edgecolor='r', facecolor='lightgrey', hatch='////', lw=0.5,ls='-')
b.contourf返回值叠加添加
这个方法比较复杂，不建议大家使用。但是精度比添加shp的方式高。主要办法是，将contourf绘制的等高线图转化，使遮盖阴影与高程相关性较高。
首先，利用nc格式的全球高程数据绘制等值填色图，这个数据在气象家园有分享。那个链接肯定比我分享百度网盘快。
file=r'E:\aaaa\world_geo.nc'ds=xr.open_dataset(file)x=ds['x']x=x.loc[69.95:101.95]y=ds['y']y=y.loc[25.95:40.95]z=ds['z'].loc[25.95:40.95,69.95:101.95]ax.contourf(x,y,z,levels=[0,4000,9000],transform=ccrs.PlateCarree())
这里为了看出青藏高原，我设置了两个等次，下面叠加的时候不用这样，只取4000以上就行了。
然后设置地形遮盖为高zorder，气温为低zorder，就可以起到遮蔽效果。
file=r'E:\aaaa\world_geo.nc'ds=xr.open_dataset(file)x=ds['x']x=x.loc[69.95:101.95]y=ds['y']y=y.loc[25.95:40.95]z=ds['z'].loc[25.95:40.95,69.95:101.95]ax.contourf(x,y,z,levels=[4000,9000], cmap='Greys', transform=ccrs.PlateCarree(), zorder=5)file=r'E:\aaaa\analyse\air.mon.mean.nc'ds=xr.open_dataset(file)lon=ds['lon']lat=ds['lat']time=ds['time']air=ds['air']air=air.loc['1948-01-01',850,:,:]ac=ax.contourf(lon,lat,air,
levels=np.arange(-34,38,4), cmap='Spectral_r', transform=ccrs.PlateCarree())cb=fig.colorbar(ac,ax=ax,shrink=0.8,pad=0.05)cb.ax.tick_params(labelsize=6)cb.set_ticks(np.arange(-34,38,4))cb.ax.set_title('气温（℃）',fontsize=8)
当然，这都是看使用者的意愿进行遮盖。
二、剖面图缺测
这个摸鱼大佬已经讲过了，我以前也写过一点。这里以上次剖面图的一张图为代表进行简单介绍。
这里可以看出，沿北纬30度的经高图，有一个诡异的直上直下的白图区，中间又突然出现一个上下贯通的水汽富区。不妨查看地形图：
我怀疑左侧诡异白条是青藏高原主体造成的，中间的水汽柱是雅鲁藏布峡谷，右侧白条是横断山脉。
按照摸鱼大佬的方法，我们首先将y轴变为指数型分布的高度轴。因为大气主要质量集中在近地面，越往上空气越稀薄，所以从上往下，气压层密度应该是越来越大的。
主要使用的语句是：
ax.set_yscale('symlog')ax.set_yticks([1000,925,850,700,500,300,200,100])ax.set_yticklabels([1000,925,850,700,500,300,200,100])ax.set_ylim(100,1000)ax.invert_yaxis()
除了摸鱼大佬的那种方法外，我们还可以使用寄生轴系统，新生成一个y轴，作为高度km轴。然后大致按照100hPa约在16公里处设置范围，并提取地图数据绘制。由于我们是推算平均海平面与100hPa距离的方法，不如摸鱼的公式计算精确，但气压层叠地形本来误差就大，只是一个大概的说明即可。
ax_h=ax.twinx()ax_h.set_ylim(0,16)ax_h.set_ylabel('Height(km)',fontsize=12)file=r'E:\aaaa\world_geo.nc'ds=xr.open_dataset(file)x=ds['x'].loc[60.95:140.95]z=ds['z'].loc[30.016666,60.95:140.95]ax_h.plot(x,z*0.001,c='k',lw=1)
还是能看出，我们之前的推测基本正确，青藏高原沿着诡异白区贴合分布。然后使用fill_between命令将地形上色。
ax_h.fill_between(x,y1=z*0.001,y2=0,where=(z>0),facecolor='grey')
