用了氣象家園大佬的shp2clip函數(shù)繁扎，他是基于basemap寫的育勺，我稍微修改了一下
直接上白化函數(shù)
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匹配的唯一標(biāo)識(shí)符线定，record[]中必有一項(xiàng)是對(duì)應(yīng)的野建。
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
傳統(tǒng)的方法是根據(jù)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文件咋整嗽测？ 自己造一個(gè)唄。
其實(shí)肿孵，基于cnmaps使得這個(gè)很容易實(shí)現(xiàn)
以廣東省為例, 以geopandas的GeoDataFrame格式返回唠粥，可以向engine參數(shù)傳入'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')
好了, 省級(jí)地圖做好了 停做。
由于shp2clip函數(shù)需要讀取地形文件,它是按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蛉腌，然后畫圖實(shí)驗(yàn)一下
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
大功告成
