常见地图白化方法（一）
本帖最后由 气ython风雨 于 2024-3-3 13:59 编辑前言
地图白化是一种绘制地图的技术，它可以实现对感兴趣区域以外的数据进行遮盖或填充白色的效果，从而突出显示目标区域的特征。
地图白化的原理是利用 shapefile 文件中的多边形坐标来创建一个剪切路径，然后将这个路径应用到 matplotlib 的绘图对象上，使得只有路径内的数据可见，路径外的数据被隐藏或覆盖。
气象家园的五星上将兰溪之水说过，其实所谓的“白化”一般就两种途径：①mask数据；②mask图形
方法一 set_clip_path方法
matplotlib 官方函数示例
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([, , ,
[-1, 0], ])
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()
代码的目的是生成一个二维网格图，并在其中添加一个路径（path）作为剪裁（clip）区域，只显示路径内部的部分。
原图如下
fig, ax = plt.subplots()
ax.imshow(Z,
interpolation ='bilinear',
cmap = cm.gray,
origin ='lower',
extent =[-3, 3, -3, 3])
plt.show()
属于mask图片类型
读取era5数据
import xarray as xr
nc = xr.open_dataset('/home/mw/input/1107125177/2023110720.nc')
数据处理
data= nc.z
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
path = Path.make_compound_path(*geos_to_path(poly))
绘图部分
In :
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()
# 白化中国以外的区域
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是什么，根本看不懂
哎呀，别急，马上给你print出来
geo_list
你问我为什么知道广东省的序号是17，那当然是用geopandas看的
shp = gpd.read_file("/home/mw/input/china1656/china_map/china_map/China_Province_2022.shp")
shp
这时这个方法就显示出弊端，读者的shp千变万化，没安装geopandas又对自己的shp列表不熟，则要一一排查geolist
方法二 rioxarray的rio.clip
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()
cb = plt.colorbar(cf, orientation='vertical', shrink=0.8)
plt.show()
以上就是比较典型的mask数据类，此类锯齿较明显
拓展：一个省怎么画：使用裁剪后的行政区json进行绘图
gd = 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()
plt.show()
方法一参考云台书使：气象绘图加强版（三十）——白化补叙
方法二参考：rioxarray官方文档
文章来源于微信公众号：气ython风雨
页:
[1]