常见地图白化方法（一）

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.025
x = 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 xr
nc = xr.open_dataset('/home/mw/input/1107125177/2023110720.nc')

数据处理
data= nc.z[0,2,:,:]
lon=data.longitude
lat=data.latitude

导入库与生成路径
import cartopy.io.shapereader as shpreader
from cartopy.mpl.patch import geos_to_path
import matplotlib.pyplot as plt
from matplotlib.path import Path
import cartopy.crs as ccrs
import 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 os
import numpy as np
import geopandas as gpd
import rioxarray
import xarray as xr
from shapely.geometry import mapping
import 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()