空间可视化！超实用的pcolormesh方法汇总！

geoallan 2022-02-16

展开全文

本节提要：关于色码图pcolormesh的入门内容。

色码图pcolormesh是除开contourf的一种比较常用的二维色彩图，经常在格点资料中使用。其数据需求与contourf非常类似。

一、pcolormesh基础参数

由于pcolormesh经常对二维数据进行可视化，所以其参数非常类似于等值线填色图，主要需求为三个参数：x,y,z。x,y表示数据横纵坐标数据，在空间场一般为经度和纬度。这里我们以再分析资料的相对湿度场为例：

ds=xr.open_dataset(r'E:\aaaa\datanc\2020.7.15\fnl_20200715_18_00.nc')lat=ds['lat_0'].loc[50:20]lon=ds['lon_0'].loc[95:125]RH=ds['RH_P0_L100_GLL0'].loc[85000,50:20,95:125]

以上表示提取北纬20-50，东经95-125，850百帕相对湿度场。然后绘图：

ap=ax.pcolormesh(lon,lat,RH,transform=ccrs.PlateCarree(),cmap='viridis_r')

二、norm与cmap参数

在上面这张基础的色码图图上，观察其色条，可知pcolormesh在默认情况下使用的是均匀映射的norm，这显然是不能满足我们的需求的。pcolormesh没有contourf的levels参数，所以只能使用更改norm参数的方式。norm的使用参考该章节：气象绘图加强版（二十八）—cmap、cbar。

colorlevel=np.arange(0,110,10)cmap=mpl.cm.Bluesnorm=mcolors.BoundaryNorm(colorlevel,cmap.N)ap=ax.pcolormesh(lon,lat,RH,transform=ccrs.PlateCarree(),cmap=cmap,norm=norm)

三、vmin、vmax

按照官网文档的说法，这两个参数不可以与norm同时指定。该参数可以指定绘制图片的上下限，例如相对湿度70%以下进行钝化：

ap=ax.pcolormesh(lon,lat,RH,transform=ccrs.PlateCarree(),cmap='Blues',vmin=75)

这个指定的缺点是，实际上小于70的地方还是进行了上色过程。当我们更换色条时，会更加明显：

这与contourf的levels参数不同，若levels限制了范围，则范围之外不会上色：

四、使用缺测方法凸显

这里我们使用官网用过的缺测方法，替换相对湿度低于70%的值为缺测。

RH=ds['RH_P0_L100_GLL0'].loc[85000,50:20,95:125]colorlevel=np.arange(70,110,10)cmap=mpl.cm.Bluesnorm=mcolors.BoundaryNorm(colorlevel,cmap.N)RH_mask=np.ma.masked_where(RH<70,RH)ap=ax.pcolormesh(lon,lat,RH_mask,transform=ccrs.PlateCarree(),cmap='viridis_r',norm=norm)

可以看出，低于70%相对湿度的地方完全没有上色。

五、常见图

pcolormesh与地图结合主要展示格点资料，与统计结合主要是相关系数图或排列时序图。

比如下面这张图：

该图就是比较明显色码风格。scatter与pcolormesh都可以实现这样的样式，而数字则属于另外text上去的。我们同样以再分析资料计算水汽通量散度场来绘制类似风格的图片。

首先加载库包与绘制底图：

import cartopy.crs as ccrs import cartopy.feature as cfimport cartopy.io.shapereader as shpreaderimport matplotlib.pyplot as pltimport matplotlib as mplimport numpy as npfrom cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTERimport cartopy.feature as cfeatureimport xarray as xrimport netCDF4 as ncfrom metpy.units import unitsimport metpy.calc as mpcalcfrom scipy import interpolateimport warningsimport matplotlib.ticker as mtickerwarnings.filterwarnings('ignore')plt.rcParams['font.sans-serif']=['FangSong']plt.rcParams['axes.unicode_minus']=False proj = ccrs.PlateCarree()fig=plt.figure(figsize=(4,4),dpi=300)ax=fig.add_axes([0,0,1,1],projection=ccrs.LambertConformal(central_longitude=110, central_latitude=30))ax.add_feature(cf.COASTLINE.with_scale('50m'),lw=0.5)ax.add_geometries(shpreader.Reader(r'E:\hubei\hubei.shp').geometries(),lw=0.5,crs=proj,edgecolor='k',facecolor='none')gl=ax.gridlines(draw_labels=True,linewidth=0.5,linestyle='--',color='grey',alpha=0.75)gl.xlocator=mticker.FixedLocator(np.arange(100,130,5))gl.ylocator=mticker.FixedLocator(np.arange(20,40,5))gl.top_labels=Falsegl.right_labels=Falsegl.bottom_labels=Truegl.rotate_labels=Nonegl.xformatter=LONGITUDE_FORMATTERgl.yformatter=LATITUDE_FORMATTERgl.x_inline=Falsegl.y_inline=Falsegl.xlabel_style={'size':7}gl.ylabel_style={'size':7}ax.set_extent([100,120,20,40],crs=ccrs.PlateCarree())

然后计算水汽通量散度并绘图：

def div_calculate(path): ds = xr.open_dataset(path) T=ds['TMP_P0_L100_GLL0'][35][:][:] rh=ds['RH_P0_L100_GLL0'][35][:][:] u=ds['UGRD_P0_L100_GLL0'][35][:][:] v=ds['VGRD_P0_L100_GLL0'][35][:][:] es=6.112*np.power(np.e,(17.67*(T-273.15)/(T-35.86))) e=(rh/100*es) q=rh*(0.622*es/(850-0.378*e))/100 qu=q*u/9.8 qv=q*v/9.8 lat=ds['lat_0'] lon=ds['lon_0'] dx,dy=mpcalc.lat_lon_grid_deltas(lon,lat) div=mpcalc.divergence(u=qu[:], v=qv[:], dx=dx, dy=dy) return lon,lat,div,u,vimport matplotlib.colors as mcolorscmap=mpl.cm.viridiscolorlevel=np.arange(-15,1,1)norm=mcolors.BoundaryNorm(colorlevel,cmap.N)lon,lat,div,u,v=div_calculate(r'E:\aaaa\datanc\2021.6.18\fnl_20210617_12_00.nc')lon=lon.loc[95:125]lat=lat.loc[35:25]div=div.loc[35:25,95:125]div_mask=np.ma.masked_where(div>0,div)ax.pcolormesh(lon,lat,div_mask*100000000,cmap=cmap,norm=norm,transform=proj)lon,lat=np.meshgrid(lon,lat)for i in range(div.shape[0]): for j in range(div.shape[1]):

        ax.text(lon[i,j]-0.5,lat[i,j]-0.1,r'${}$'.format(int(div[i,j].values*100000000)),zorder=5,fontsize=8,transform=proj,color='w')

另外，pcolormesh还可以用来绘制热力图或累积时序图：

上面这张图在横轴上表现为时间序列，在纵轴上表现为累积序列。要实现这样的图，最主要的就是构造二维数据。编造一个工厂车间的生成绩效表如下：

读取数据后，重构为二维数据：

file=r'C:\Users\lenovo\Desktop\1-12.xlsx'data=pd.read_excel(file)time=data['时间']mix_data=np.array([data['一车间'],data['二车间'],data['三车间'],                   data['四车间'],data['五车间'],data['六车间']])

由于pandas中的标准时间序列可以被matplotlib识别，所以我们可以直接绘图：

fig=plt.figure(figsize=(6,4),dpi=300)ax=fig.add_axes([0,0,1,1])colorlevel=np.arange(0,21,1)cmap=mpl.cm.viridis_rnorm=mcolors.BoundaryNorm(colorlevel,cmap.N)ap=ax.pcolormesh(time,range(mix_data.shape[0]),mix_data,shading='auto',edgecolor='w',cmap=cmap,norm=norm,lw=3)divider=make_axes_locatable(ax)cax=divider.new_horizontal(size='5%', pad=0.05,axes_class=plt.Axes)fig.add_axes(cax)cb=fig.colorbar(ap,cax=cax)cb.ax.tick_params(labelsize=10)cb.ax.set_ylabel('生产件数',fontsize=10)ax.set_yticks([0,1,2,3,4,5])ax.set_yticklabels(['一车间','二车间','三车间','四车间','五车间','六车间'])ax.tick_params(axis='x',labelsize=7)ax.set_xlabel('时间',fontsize=13)ax.set_ylabel('生产车间',fontsize=13)