气象绘图加强版（十二）——白化杂谈

气海无涯 2021-08-31

5197

本节提要：对绘图的白化部分做一个总结。

什么是白化？我在一年前也是头一次接触到这个词语，其实就是将你不需要的部分的等值线、等值线填色、风场、流场等挖去。目前气象领域流行的是花式利用地图shp文件进行操作，达到白化的目的。

我现在能够看到和掌握的白化方法有以下几种：

①气象家园平流层的萝卜首发、晋陵小生优化的maskout模块白化功能；

②公众号好奇心Log推送的salem库包下的白化；

③公众号DataCharm推送的geopandas库包下的clip白化；

④公众号DataCharm推送的Fiona、shapely库包赋nan值白化；

⑤气象家园Masterpiece提出的shp转path判别法。

这几种白化方法中，最知名的应该是maskout这个方法，而与后面三种不同的是，这个方法是真正的裁剪，因为按照其使用流程：

ac=ax.contourf(...)
clip=maskout.shp2clip(ac, ax,r'E:\dijishi\cn_province.shp' ,420000)

可以看出，这种方法是先画出等值线，然后裁剪。其他四种都是先判别是不是在指定shp文件内部，然后再画，但是geopandas.clip的办法与Masterpiece的办法会改变数据的维度，导致无法还原为2D数组，不能用在contourf绘制等值线图上，但是可以用在scatter上，利用c与cmap两个关键字参数进行颜色映射，达到类似的效果。兹简单谈谈各种白化的使用方法。

一、maskout白化

这个方法其实很多地方都使用过，比如气象学家公众号上面，还有很多同志也自主开发了类似方法。但是由于气象家园的这个是影响最大的，也是比较完善的，我单独提出来。这里只简单讲解一下白化的原理：

平流层的萝卜的程序与气象家园上的程序都是使用set_clip_path功能来进行白化，进入官网查询本命令，可知set_clip_path主要对collection类作白化。我们常用的等值线contour与等值线填色contourf查询其结构，可知其返回一个quadcontourset：

于是我们进一步查询该命令，可以知道等值线填色图确实是linecollection或pathcollection的集合，所以可以用set_clip_path来裁剪等值线图：

首先展示在普通子图中对等值线进行裁剪，首先需要生成一个边界path，我们将其命名为boundary，然后对集合中的每一个collection进行裁剪：

#如何简单理解白化
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.path as mpath
plt.rcParams['font.sans-serif']=['KaiTi']#英文新罗马字体
plt.rcParams['axes.unicode_minus']=False
fig=plt.figure(figsize=(3,1.5),dpi=500)
vertices=[(-5,-5),(-5,5),(5,5),(5,-5),(-5,-5)]
boundary=mpath.Path(vertices)
x=np.arange(-10,10,1)
y=np.arange(-10,10,1)
x,y=np.meshgrid(x,y)
z=x**2+y**2
ax1=plt.subplot(121)
ax2=plt.subplot(122)
ac1=ax1.contourf(x,y,z)
ac2=ax2.contourf(x,y,z)
for i in ac1.collections:
    i.set_clip_path(boundary,transform=ax1.transData)
for i in vertices:
    ax1.text(i[0],i[1],i,fontsize=4)
for ax in [ax1,ax2]:
    ax.tick_params(direction='in',
                   labelsize=2,
                   top=True,
                   right=True,
                   width=0.3,
                   length=2)
ax1.set_title('白化等值线后',fontsize=5)
ax2.set_title('白化等值线前',fontsize=5)

当然，这只是构造了一个最简单的矩形框来对等值线进行裁剪。我们平时对地图白化，一般不用这种粗糙的方式。出于精细化的考虑，我们一般用shp文件的地图作精细化裁剪。而裁剪是需要path（路径）的，这个路径从哪里来，我们肯定不能徒手去构造。在上一期提到的三维地图的添加时我们就提到了各种绘图集合图形之间是可以相互转化的，我们还是得使用那个geos_to_path的命令，将shp里的集合形状geometry转化成path，然后用set_clip_path命令来裁剪。这也是平流层的萝卜编辑maskout文件时的方法。不过大佬为了丰富白化功能，增添了一个region参数，使得我们能够凭关键词筛选要白化的省份，并且不使用cartopy中的geos_to_path命令，而是使用的shapefile库包来读取边界。

这里我们只简单讲解基层原理：

import numpy as np
import cartopy.crs as ccrs 
import cartopy.feature as cf
import matplotlib.pyplot as plt
import cartopy.io.shapereader as shpreader
from cartopy.mpl.ticker import LatitudeFormatter,LongitudeFormatter
from matplotlib.path import Path
from cartopy.mpl.patch import geos_to_path
plt.rcParams['font.sans-serif']=['KaiTi']
shp_path=r'E:\My Documents\桌面\利川市地图\利川.shp'
shp_data=shpreader.Reader(shp_path)
fig=plt.figure(figsize=(3,2),dpi=500)
ax1=plt.subplot(121,projection=ccrs.PlateCarree())
ax2=plt.subplot(122,projection=ccrs.PlateCarree())
for i,ax in enumerate([ax1,ax2]):
    ax.add_geometries(shp_data.geometries(),crs=ccrs.PlateCarree(),edgecolor='k',facecolor='none',lw=0.5)
    ax.set_extent([108.3,109.35,29.7,30.7],crs=ccrs.PlateCarree())
    ax.set_xticks(np.arange(108.3,109.35,0.2))
    ax.set_yticks(np.arange(29.7,30.7,0.2))
    ax.xaxis.set_major_formatter(LongitudeFormatter())
    ax.yaxis.set_major_formatter(LatitudeFormatter())
    ax.tick_params(direction='in',labelsize=3,top=True,right=True,length=2,width=0.5)
    if i==0:
        ax.set_title('未白化',fontsize=6)
    else:
        ax.set_title('白化后',fontsize=6)
########定义绘图数据######################
x=np.arange(108.3,109.35,0.01)
y=np.arange(29.7,30.7,0.01)
X,Y=np.meshgrid(x,y)
Z=(X-108)**2+(Y-29)**2
#######循环画图#########################
for i,ax in enumerate([ax1,ax2]):
    if i==0:
        ax.contourf(X,Y,Z)
    else:
        ac=ax.contourf(X,Y,Z)
#######获取path#######################
records=shp_data.records()
for record in records:
    path=Path.make_compound_path(*geos_to_path([record.geometry]))
#######白化###########################
for collection in ac.collections:
    collection.set_clip_path(path,transform=ax2.transData)

同样的等值线也可以用这种方法白化：

上述只是通过普通子图的白化推广到地图子图的白化，所以默认的PlateCarree投影。这种方法的过程都是如下：

通过地图库包获取当前shp文件信息→将geometry转化为path→绘制等值线→使用得到的path对等值线的collection进行裁剪。另外，此命令对大部分绘图后生成的为collection的都可以白化，比如scatter绘图后生成pathcollections。

maskout库包就是对上述流程的优化和功能增强，特别是优化后的maskout库包还可以裁剪兰勃脱、麦卡托等其他投影。当前白化原理简略即如此。

具体如何使用maskout程序，可以参考我在一年前刚学习时候的推文Python气象绘图教程特刊（一）。

这里简单回复很多人提过的问题——怎么引入maskout。新手苦恼的地方就是import maskout时报错，查找不到这个库包。其实我们知道，各个库包也就是打包起来的后缀为.py的程序而已。如果看过气象家园的文章，可以知道作者提供了一个名为maskout.py的Python程序。但是不知道怎么引入。这里我提供了一个懒办法——添加指定的读取路径。比如我将下载下来的maskout.py文件放在我的桌面上。然后在引入时添加如下一句：

import sys
sys.path.append(r"C:\Users\lenovo\Desktop")
import maskout

这样，平台会在桌面上找maskout这个库包。

你也可以将这个程序放置到当前anaconda下库包安装位置，也是可以找到的。

这个白化方法，有一个问题就是每更换一次shp文件，就必须重新查找record并对maskout中的相关部分修改。而且是先画后裁剪，并不能筛选指定地区的数据。

二、salem白化

这个方法是先筛选数据，将处于指定shp内部的数据点筛选出来后，再画图。具体使用方法请参考上面提到的公众号原文。这个数据裁剪不改变数组的维度，所以后期可以用在contourf上。

三、geopandas.clip白化

本方法我是第一次在DataCharm公号上看到的，具体使用的就是geopandas自带的clip功能。推文中说道，该方法适合plotnine库包。我试验了一下，确实如此，由于构建DataFrame时需要将经纬度即数据扁平化，裁剪之后的数据难以还原为原本的维度和形状，所以不太适合用于matplotlib的contourf可视化，但是我们可以用scatter来展示其原理。scatter不要求数组为2D形状。

前面步骤与其他方法一致，先将站点数据使用径向基函数插值为格点数据：

from scipy.interpolate import Rbf#引入径向基函数
import pandas as pd
import numpy as np
filename=r'C:\Users\lenovo\Desktop\example.xlsx'
df=pd.read_excel(filename)
lon=df['lon']#站点经度
lat=df['lat']#站点纬度
data=df['data']#站点数据
olon=np.linspace(108.3,109.35,500)#格点经度
olat=np.linspace(29.7,30.7,500)#格点纬度
olon,olat=np.meshgrid(olon,olat)
func=Rbf(lon,lat,data,function='linear')
data_new=func(olon,olat)
olon
olat
data_new #插值后的格点数据

至此，我们得到三个插值完成后的格点数据olon,olat,data_new。随后进入下一步，构造用于裁剪的GeoDataFrame。在这里插一句，geopandas.clip这个函数，只能用来裁剪固定格式的point,line,polygon:

如此说来，我们用上述哪种方式裁剪最好呢？由于我们得到的是格点资料，说以最容易构造的就是points。首先我们将olon,olat,data_new构造为pandas里的DataFrame:

df_new=pd.DataFrame(dict(lon=olon.flatten(),lat=olat.flatten()))
df_new['data']=data_new.flatten()
df_new

然后利用经纬度构造GeoDataFrame，这里用crs='EPSG:4326'的原因是保证df_geo与后面a.crs一致，就是cartopy里的PlateCarree():

df_geo=gpd.GeoDataFrame(df_new,
                        geometry=gpd.points_from_xy(df_new['lon'],df_new['lat']),
                        crs='EPSG:4326')
df_geo

然后读取shp文件：

a=gpd.read_file(r'E:\map\利川.shp',encoding='UTF-8')
a

此时还可以查验a的crs：

a.crs

随后利用a为范围，裁剪df_new

df_clip=gpd.clip(df_geo,a)
df_clip

看统计信息，250000行减少为102345行数据。失去的那一部分就是被裁剪的。接下来我们使用scatter函数来绘制图形：

ac1=ax.scatter(olon,
               olat,
               s=0.1,
               marker='s',#s代表正方形marker
               c=data_new,
               cmap='Spectral_r')
ac=ax2.scatter(df_clip['lon'],
               df_clip['lat'],
               s=0.1,
               marker='s',
               c=df_clip['data'],
               cmap='Spectral_r')
fig.colorbar(ac,ax=[ax1,ax2])

能这么圆滑是因为在前面网格化时我们划分了500格点，如果使用粗格点就能看出scatter效果与geopandas.clip的裁剪痕迹，反而有点pcolormesh的味道：

olon=np.linspace(108.3,109.35,40)
olat=np.linspace(29.7,30.7,40)

四、Fiona、shapely库包赋nan值

这又是DataCharm大佬讲述的一种方法，不过在原文中，其使用的是basemap，我们更换为cartopy。这种方法可以保住数据的维度为2D结构，所以既可以用scatter也可以用contourf绘图。不过在网格比较粗的情况下会出现裁剪劣化，不该白化的地方掉白变多。

由于shp边界之外的地方变为nan值，所以会变白。这里和其他几种不一样，其他几种直接将点减掉了，所以其他几种不能回复原来的网格维度，而这里可以回复原来的网格维度。

ac=ax.scatter(olon,olat,s=9,marker='s',c=data_new,cmap='Spectral_r')
print(data_new.shape)
df_grid=pd.DataFrame(dict(lon=olon.flatten(),lat=olat.flatten()))
df_grid['data']=data_new.flatten()
lichuan_shp=fiona.open(r'E:\家园\利川市地图\利川.shp',encoding='utf-8')
a=gpd.read_file(r'E:\家园\利川市地图\利川.shp',encoding='utf-8')
pol=lichuan_shp.next()
poly_data=pol['geometry']['coordinates'][0]
shp_polygeon=Polygon(poly_data)
mask_value=[value if shp_polygeon.contains(Point(lon,lat))==True\
                  else np.nan for lon,lat,value in zip(df_grid['lon'],df_grid['lat'],df_grid['data'])]
ax.scatter(olon,olat,s=9,marker='s',c=mask_value,cmap='Spectral_r')

ac=ax.contourf(olon,olat,data_new,cmap='Spectral_r')
print(data_new.shape)
df_grid=pd.DataFrame(dict(lon=olon.flatten(),lat=olat.flatten()))
df_grid['data']=data_new.flatten()
lichuan_shp=fiona.open(r'E:\家园\利川市地图\利川.shp',encoding='utf-8')
a=gpd.read_file(r'E:\家园\利川市地图\利川.shp',encoding='utf-8')
pol=lichuan_shp.next()
poly_data=pol['geometry']['coordinates'][0]
shp_polygeon=Polygon(poly_data)
mask_value=[value if shp_polygeon.contains(Point(lon,lat))==True\
                  else np.nan for lon,lat,value in zip(df_grid['lon'],df_grid['lat'],df_grid['data'])]
df_grid['mask_value']=mask_value
mask_grid=df_grid['mask_value'].values.reshape(40,40)
ax2.contourf(olon,olat,mask_grid,cmap='Spectral_r')

五、气象家园Masterpiece提出的方法

这种方法其实十分类似geopandas的方法，最终生成的mask数组都被降维，难以恢复到原来的维度，如果直接contourf会报出错误

但是我们可以用scatter映射的方法达到类似的结果。

import matplotlib.pyplot as plt
import numpy as np
import cartopy.crs as ccrs
import cartopy.io.shapereader as shpreader
import shapefile
from matplotlib.path import Path
from cartopy.mpl.ticker import LatitudeFormatter,LongitudeFormatter
plt.rcParams['font.sans-serif']=['KaiTi']
shp_path=r'E:\My Documents\桌面\利川市地图\利川.shp'
shp_data=shpreader.Reader(shp_path)
fig=plt.figure(figsize=(3,2),dpi=600)
ax1=plt.subplot(121,projection=ccrs.PlateCarree())
ax2=plt.subplot(122,projection=ccrs.PlateCarree())
for i,ax in enumerate([ax1,ax2]):
    ax.add_geometries(shp_data.geometries(),crs=ccrs.PlateCarree(),edgecolor='k',facecolor='none',lw=0.5)
    ax.set_extent([108.3,109.35,29.7,30.7],crs=ccrs.PlateCarree())
    ax.set_xticks(np.arange(108.3,109.35,0.2))
    ax.set_yticks(np.arange(29.7,30.7,0.2))
    ax.xaxis.set_major_formatter(LongitudeFormatter())
    ax.yaxis.set_major_formatter(LatitudeFormatter())
    if i==0:
        ax.set_title('未白化',fontsize=6)
        ax.tick_params(direction='in',labelsize=3,top=True,right=True,length=2,width=0.5)
    else:
        ax.set_title('白化后',fontsize=6)
        ax.tick_params(direction='in',labelsize=3,top=True,right=True,length=2,width=0.5)
########定义绘图数据######################
x=np.arange(108,109.5,0.05)
y=np.arange(29.7,30.7,0.05)
lon_num=len(x)
lat_num=len(y)
lons,lats=np.meshgrid(x,y)
Z=(lons-108)**2+(lats-29)**2
#######################################
sf=shapefile.Reader(shp_path)
Shapes=sf.shapes()
boundary=Shapes[0].points
p= Path(boundary)
grib=np.stack((lons,lats),axis=-1) 
mask=p.contains_points(grib.reshape(lat_num*lon_num,2)) 
mask=mask.reshape(lat_num,lon_num)
ax2.scatter(lons[mask],lats[mask],s=12,c=Z[mask],cmap='viridis',marker='s')
ax1.scatter(lons,lats,s=12,c=Z,cmap='viridis',marker='s')

总结

单纯从白化的角度出发，最好用的当然是maskout这个成熟的程序，经过各方面的大佬的打磨，已经实现了对散点scatter、等值线contour、等值线填色contourf、风矢杆barb等常见气象绘图形式的裁剪。

而salem库包则是裁剪最为简便的，而且裁剪之后的数据不会改变维度和形状。

geopandas裁剪由于自身函数属性的限制，对点状数据的裁剪效果最好。

fiona和shapely的方法可以原汁原味的返回原来网格。

Masterpiece的方法就是shpaefile和matplotlib结合的clip方法，与geopandas.clip类似，破坏原来网格。后四种可以用来筛选数据。

maskout:可以用于contour、contourf、barbs、scatter的裁剪，不能筛选数据

salem：nc格式的再分析资料用起来非常好

geopandas.clip：一定局限性，可以筛选数据

fiona和shapely：性能稳定，但是比起maskout来复杂许多

path:一定局限性，可以筛选数据

气象绘图加强版（十二）——白化杂谈

评论