用了气象家园大佬的shp2clip函数，他是基于basemap写的，我稍微修改了一下
直接上白化函数

import shapefile
from matplotlib.path import Path
from matplotlib.patches import PathPatch
from shapely.geometry import Point as ShapelyPoint
from shapely.geometry import Polygon as ShapelyPolygon
from collections import Iterable
def shp2clip(originfig,ax,shpfile,region, clabel = None, vcplot = None):
    print('your mask region=',region)
    sf = shapefile.Reader(shpfile)
    vertices = []
    codes = []
    for shape_rec in sf.shapeRecords():
        ####这里需要找到和region匹配的唯一标识符，record[]中必有一项是对应的。
         if shape_rec.record[3] in region:   #####
        #if shape_rec.record[1] in region:  #####
            pts = shape_rec.shape.points
            prt = list(shape_rec.shape.parts) + [len(pts)]
            for i in range(len(prt) - 1):
                for j in range(prt[i], prt[i+1]):
                    vertices.append((pts[j][0], pts[j][1]))
                codes += [Path.MOVETO]
                codes += [Path.LINETO] * (prt[i+1] - prt[i] -2)
                codes += [Path.CLOSEPOLY]
            clip = Path(vertices, codes)
            clip = PathPatch(clip, transform=ax.transData)

    if vcplot:
        if isinstance(originfig,Iterable):
            for ivec in originfig:
                ivec.set_clip_path(clip)
        else:
            originfig.set_clip_path(clip)
    else:
        for contour in originfig.collections:
            contour.set_clip_path(clip)

    if  clabel:
        clip_map_shapely = ShapelyPolygon(vertices)
        for text_object in clabel:
            if not clip_map_shapely.contains(ShapelyPoint(text_object.get_position())):
                text_object.set_visible(False)    

    return clip

这里利用cnmaps来白化，cnmaps是真的香，简单暴力，没有的直接pip install cnmaps
传统的方法是根据shape文件读取里面的地理信息，匹配，裁剪

import numpy
import netCDF4 as nc
import sys, os
import cartopy.crs as crs
import cartopy.crs as ccrs
from cnmaps import get_adm_maps, draw_maps
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter

box   = [95, 135, 15, 45]
xstep = 15    #x axis step
ystep = 15    #y axis step
fig=plt.figure(num=1,figsize=(8.5,7))   ###????.fig
axlevel1 = [0.1, 0.1, 0.65, 0.8]

#using cnmap
draw_maps(get_adm_maps(level='国')) 
draw_maps(get_adm_maps(level='省'),linewidth=0.3,color='gray') 
plt.axis(box)
ax = plt.gca()

c1=plt.contourf(lon,lat,isop,cmap= plt.cm.get_cmap('jet') )

#below is maskwhite
shpfile='/home/wangnan/data/baihua/country1.shp'
shp2clip(c1, ax, shpfile,'China',clabel=None,vcplot=None)
#ax.add_geometries(shapereader.Reader(shpfile).geometries(), crs=cart_proj, facecolor='none',edgecolor='k',linewidth=2.5)

#axis setting
ax.set_xticks(np.arange(box[0], box[1] + xstep, xstep), )
ax.set_yticks(np.arange(20, 40 + ystep, ystep), )
ax.set_xlim(box[0], box[1])
ax.set_ylim(box[2], box[3])
lon_formatter = LongitudeFormatter(zero_direction_label=False)
lat_formatter = LatitudeFormatter()
ax.xaxis.set_major_formatter(lon_formatter)
ax.yaxis.set_major_formatter(lat_formatter)
ax.tick_params(labelsize=22)

image.png

那么问题来了，我要是没有地图shape文件咋整？ 自己造一个呗。
其实，基于cnmaps使得这个很容易实现
以广东省为例, 以geopandas的GeoDataFrame格式返回，可以向engine参数传入'goepandas'。

from cnmaps import get_adm_maps,get_adm_names
prov=get_adm_maps(level='省',engine='geopandas')
prov

image.png

该shape文件可以直接保存

prov.to_file('/home/wangnan/data/shapefile2/wnProvince.shp',encoding='utf-8')

好了, 省级地图做好了 。
由于shp2clip函数需要读取地形文件,它是按region来匹配的，所以我要看看我的region在第几列

shpfile2= '/home/wangnan/data/shapefile2/wnProvince.shp'  
sf = shapefile.Reader(shpfile2)
for shape_rec in sf.shapeRecords():
    print (shape_rec.record)

image.png

在第2列，所以要把shape2clip的if判断中的if shape_rec.record[3] in region换成1，然后画图实验一下

import numpy
import netCDF4 as nc
import sys, os,cmaps
import cartopy.crs as crs
import cartopy.crs as ccrs
from cnmaps import get_adm_maps, draw_maps
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
import matplotlib.colors as colors
box   = [109, 118, 19, 26]
xstep = 3    #x axis step
ystep = 3    #y axis step
fig=plt.figure(num=1,figsize=(8.5,7))   ###????.fig

draw_maps(get_adm_maps(province='广东省')) 
plt.axis(box)
ax = plt.gca()

cmap_color = cmaps.WhiteBlueGreenYellowRed   #WhiteGreen
c1=ax.contourf(lon,lat,terp,cmap=cmap_color )

#below is maskwhite
shpfile2= '/home/wangnan/data/shapefile2/wnProvince.shp'
region=['广东省']
shp2clip(c1, ax, shpfile2,region,clabel=None,vcplot=None)

#axis setting
ax.set_xlim(box[0], box[1])
ax.set_ylim(box[2], box[3])
lon_formatter = LongitudeFormatter(zero_direction_label=False)
lat_formatter = LatitudeFormatter()
ax.xaxis.set_major_formatter(lon_formatter)
ax.yaxis.set_major_formatter(lat_formatter)
ax.tick_params(labelsize=22)

image.png

大功告成