常见地图白化方法（一）

import numpy as np import matplotlib.cm as cm import matplotlib.pyplot as plt from matplotlib.path import Path from matplotlib.patches import PathPatch  delta = 0.025x = y = np.arange(-3.0, 3.0, delta) X, Y = np.meshgrid(x, y)   Z1 = np.exp(-X**2 - Y**2) Z2 = np.exp(-(X - 1)**2 - (Y - 1)**2) Z = (Z1 - Z2) * 2   path = Path([[0, 1], [1, 0], [0, -1],             [-1, 0], [0, 1]]) patch = PathPatch(path, facecolor ='none')    fig, ax = plt.subplots() ax.add_patch(patch)    im = ax.imshow(Z,                interpolation ='bilinear',                 cmap = cm.gray,                origin ='lower',                 extent =[-3, 3, -3, 3],                clip_path = patch,                 clip_on = True)   im.set_clip_path(patch)   fig.suptitle('matplotlib.axes.Axes.set_clip_path()function Example\n\n', fontweight ="bold")   plt.show() 

fig, ax = plt.subplots() ax.imshow(Z,     interpolation ='bilinear',      cmap = cm.gray,      origin ='lower',      extent =[-3, 3, -3, 3]) plt.show() 

读取era5数据import xarray as xrnc = xr.open_dataset('/home/mw/input/1107125177/2023110720.nc')
数据处理data= nc.z[0,2,:,:]lon=data.longitudelat=data.latitude
导入库与生成路径import cartopy.io.shapereader as shpreaderfrom cartopy.mpl.patch import geos_to_pathimport matplotlib.pyplot as pltfrom matplotlib.path import Pathimport cartopy.crs as ccrsimport cmaps# 读取shp文件shp_path = '/home/mw/input/china1656/china_map/china_map/China_Province_2022.shp'shp = shpreader.Reader(shp_path)geo_list = list(shp.geometries())proj = ccrs.PlateCarree()poly = geo_list[17]path = Path.make_compound_path(*geos_to_path(poly))
绘图部分In [16]:fig = plt.figure(figsize=(12, 8))ax = fig.add_subplot(projection=ccrs.PlateCarree())
cf = ax.contourf(lon, lat, data/98 ,levels=np.arange(520,600,4),transform=proj, cmap=cmaps.radar)cr = ax.contour(lon,lat, data/98, colors='k',levels=np.arange(520,600,4))ax.clabel(cr, inline=1, fontsize=10, fmt="%i")ax.add_geometries(shp.geometries(), proj, facecolor='none', edgecolor='k') ax.set_extent([100,120,20,40])
# 白化中国以外的区域for col in cf.collections:    col.set_clip_path(path, transform=ax.transData)# 色标 cb = plt.colorbar(cf, orientation='vertical', shrink=0.8)plt.show()

geo_list[17]

shp = gpd.read_file("/home/mw/input/china1656/china_map/china_map/China_Province_2022.shp")shp

import osimport numpy as npimport geopandas as gpdimport rioxarrayimport xarray as xrfrom shapely.geometry import mappingimport matplotlib.pyplot as plt
# 读取中国省份地图数据shp = shpreader.Reader(shp_path)shpgpd = gpd.read_file(shp_path)
# 设置数据投影data.rio.write_crs("epsg:4326", inplace=True)
# 进行裁剪并绘图cliped = data.rio.clip(shpgpd.geometry.apply(mapping), shpgpd.crs, drop=False)
fig = plt.figure(figsize=(12, 8))ax = fig.add_subplot(projection=ccrs.PlateCarree())cf = ax.contourf(lon, lat, cliped/98,levels=np.arange(520,600,4), transform=proj, cmap=cmaps.radar)cr = ax.contour(lon,lat, cliped/98,levels=np.arange(520,600,4), colors='k')ax.clabel(cr, inline=1, fontsize=10, fmt="%i")ax.add_geometries(shp.geometries(), proj, facecolor='none', edgecolor='k') ax.set_extent([100,120,20,40])cb = plt.colorbar(cf, orientation='vertical', shrink=0.8)plt.show()

gd = shpgpd[shpgpd['省'] == '广东省']gd.to_file('/home/mw/project/gd.geojson', driver='GeoJSON')


# 读取中国省份地图数据js_path = '/home/mw/project/gd.geojson'gpdjs = gpd.read_file(js_path)
# 设置数据投影data.rio.write_crs("epsg:4326", inplace=True)
# 进行裁剪并绘图cliped = data.rio.clip(gpdjs.geometry.apply(mapping), gpdjs.crs, drop=False)
fig = plt.figure(figsize=(12, 8))ax = fig.add_subplot(projection=ccrs.PlateCarree())cf = ax.contourf(lon, lat, cliped/98,levels=np.arange(520,600,4), transform=proj, cmap=cmaps.radar)cr = ax.contour(lon,lat, cliped/98,levels=np.arange(520,600,4), colors='k')ax.clabel(cr, inline=1, fontsize=10, fmt="%i")ax.add_geometries(shp.geometries(), proj, facecolor='none', edgecolor='k') ax.set_extent([100,120,20,40])plt.show()