前言
这几天在画图，需要用到白化功能，一直在搜集各位大佬的资料学习（如下）。结合自己的使用经历，给大家简单分享一下亲身体验。
【自学气象人】公众号获取
气象绘图——白化杂谈
cnmaps，你值得拥有！
python绘图 | salem一招解决所有可视化中的掩膜(Mask)问题
安装
maskout最简单
我提供一个如下的maskout.py代码（这份代码由于经过多位大佬们的完善，具体出自谁手我已经不太清楚了，反正感谢大佬们辛苦开发），大家使用前直接运行一下下面的代码或者import maskout即可。
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.abc import Sized, Container, Iterable
def shp2clip(originfig,ax,region_shpfile,clabel=None,vcplot=None):
sf = shapefile.Reader(region_shpfile)
for shape_rec in sf.shapeRecords():
vertices = []
codes = []
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安装难度正常
cnmpas是Clarmy吱声公众号的大佬开发的，非常好用！安装方法如下：
conda install -c conda-forge cnmaps
唯一一个小小的问题，就是在使用过程中必须保持shapely=1.8.5，不然部分函数在使用时候会出现报错（如下）。
遇到这种报错直接运行如下命令再重启内核即可解决。
pip install shapely＝＝1.8.5
salem安装稍有难度
salem是github上fmaussion大佬开发的，也非常好用，安装方法如下：
conda install -c conda-forge salem
但安装好后，第一次import salem的时候，会要求联网下载salem-sample-data。而这个东西比较难下载，很少直接下成功。 推荐自己直接去fmaussion的github主页上download下来，然后把压缩包重新命名后放在~/.salem_cache文件夹下面即可。
下载后放在~/.salem_cache文件夹下。 
注意：请注意下载下来解压后文件夹的命名。具体如何命名，可以根据import salem下面的报错，找到这个~/anaconda3/envs/cnmaps/lib/python3.10/site-packages/salem/utils.py文件，加入print(ofile)看看具体的命名。
这是报错。 
这是根据报错找到文件加入print(ofile)。 
这是修改后再次运行import salem，会告诉我们应该修改的具体名字。
使用
先给个亲身体验后的总结：
-
cnmpas自带了中国的各级省市县shp，使用起来非常方便，不需要我们额外去下载shp文件了。而且支持的功能也非常强大，不仅支持等值线、填色图、风场图的裁剪等，还支持掩膜选取数据。但就是目前不支持导入自己下载的shp文件，所以对于类似青藏高原这种不在默认shp里面的地方，就没办法处理了。
-
maskout使用方便，需要自己导入shp文件，但不支持掩膜选取数据,且速度稍慢。也就是说是完全的裁剪：先画好图，然后根据shp文件从完整的图中扣出自己感兴趣的区域。注意，用maskout时候无法调整中心经度。因为裁剪的位置在整张图片中是固定的，但调整中心经度后图片偏移了，而裁剪位置没变，故裁剪结果就不对了。
-
salem安装略有点麻烦，但功能依然非常强大，需要自己导入shp文件，且支持掩膜选取数据，还可以对地图投影做转换，对WRF的后处理也有一定的支持。
综上，如果不想麻烦，只需要裁剪自己感兴趣的区域，直接上maskout就行，但注意不要使用central_longitude等使图片偏移的参数；如果处理省市县级数据，建议使用cnmaps；salem可以作为cnmasp的补充，处理那些不在默认shp内的区域，比如青藏高原。
具体使用方法可以看本文最前面的几篇推送，介绍比较详细，下面放几个代码示例：
原始数据
import xarray as xr
import numpy as np
import cartopy.crs as ccrs
import cartopy.feature as cfeat
from cartopy.mpl.ticker import LongitudeFormatter,LatitudeFormatter
import matplotlib.pyplot as plt
ds = xr.open_dataset('analysis/data.nc')
fig, ax = plt.subplots(1, 1, figsize=(16, 12), dpi=200, subplot_kw={'projection': ccrs.PlateCarree(central_longitude=180)})
fig.subplots_adjust(hspace=0.3, wspace=0.3)
ax.coastlines(lw=0.7)
clevs = np.arange(-10, 10, 1)
im = ax.contourf(ds.lon, ds.lat, ds['t2m'],
clevs, cmap='coolwarm', extend='both', transform=ccrs.PlateCarree())
ax.set_xticks([ 60,120,180,240,300], crs=ccrs.PlateCarree())
ax.set_yticks([ -90,-60, -30, 0, 30, 60,90], crs=ccrs.PlateCarree())
lon_formatter = LongitudeFormatter(zero_direction_label=False)
lat_formatter = LatitudeFormatter()
ax.xaxis.set_major_formatter(lon_formatter)
ax.yaxis.set_major_formatter(lat_formatter)
plt.show()
maskout选取青藏高原
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.abc import Sized, Container, Iterable
def shp2clip(originfig,ax,region_shpfile,clabel=None,vcplot=None):
sf = shapefile.Reader(region_shpfile)
for shape_rec in sf.shapeRecords():
vertices = []
codes = []
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
#处理画图----------------------------------------------------
import xarray as xr
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.mpl.ticker import LongitudeFormatter,LatitudeFormatter
import numpy as np
TP_shp = '/work/home/zhxia/my_tool/TP_shp/TPBoundary_new_2021/TPBoundary_new(2021).shp'
ds = xr.open_dataset('analysis/Pangu_ERA5_fct360_winter.nc')
lon = ds.lon.values
lat = ds.lat.values
lons, lats = np.meshgrid(lon, lat)
fig, ax = plt.subplots(1, 1, figsize=(16, 12), dpi=200, subplot_kw={'projection': ccrs.PlateCarree()}) #不能使用central_longitude=180，否则裁剪的位置不对
fig.subplots_adjust(hspace=0.3, wspace=0.3)
ax.coastlines(lw=0.7)
clevs = np.arange(-10, 10, 1)
im = ax.contourf(lons, lats, ds['t2m'],
clevs, cmap='coolwarm', extend='both', transform=ccrs.PlateCarree())
ax.set_xticks([ -120, -60, 0, 60, 120], crs=ccrs.PlateCarree()) #不使用central_longitude=180后，这里相应修改一下方便显示
ax.set_yticks([ -90,-60, -30, 0, 30, 60,90], crs=ccrs.PlateCarree())
lon_formatter = LongitudeFormatter(zero_direction_label=False)
lat_formatter = LatitudeFormatter()
ax.xaxis.set_major_formatter(lon_formatter)
ax.yaxis.set_major_formatter(lat_formatter)
clip1=shp2clip(im,ax,TP_shp) #这一步裁剪
plt.show()
salem选取青藏高原
import xarray as xr
import numpy as np
from xarray.backends import NetCDF4DataStore
import salem
from datetime import datetime
import cartopy.crs as ccrs
import cartopy.feature as cfeat
from cartopy.mpl.ticker import LongitudeFormatter,LatitudeFormatter
from cartopy.io.shapereader import Reader, natural_earth
import matplotlib.pyplot as plt
import matplotlib.ticker as mticker
import geopandas
ds = xr.open_dataset('data.nc')
TP_shp = '/work/home/zhxia/my_tool/TP_shp/TPBoundary_new_2021/TPBoundary_new(2021).shp'
shp = geopandas.read_file(TP_shp)
t = ds['t2m'].salem.roi(shape=shp)
fig, ax = plt.subplots(1, 1, figsize=(16, 12), dpi=200, subplot_kw={'projection': ccrs.PlateCarree(central_longitude=180)})
fig.subplots_adjust(hspace=0.3, wspace=0.3)
ax.coastlines(lw=0.7)
clevs = np.arange(-10, 10, 1)
im = ax.contourf(ds.lon, ds.lat, t,
clevs, cmap='coolwarm', extend='both', transform=ccrs.PlateCarree())
ax.set_xticks([ 60,120,180,240,300], crs=ccrs.PlateCarree())
ax.set_yticks([ -90,-60, -30, 0, 30, 60,90], crs=ccrs.PlateCarree())
lon_formatter = LongitudeFormatter(zero_direction_label=False)
lat_formatter = LatitudeFormatter()
ax.xaxis.set_major_formatter(lon_formatter)
ax.yaxis.set_major_formatter(lat_formatter)
plt.show()
cnmpas选取中国区域
通过裁剪图片的方法
import numpy as np
from cnmaps import get_adm_maps
import matplotlib.pyplot as plt
map_polygon = get_adm_maps(country='中华人民共和国', record='first', only_polygon=True)
ds = xr.open_dataset('data.nc')
fig, ax = plt.subplots(1, 1, figsize=(16, 12), dpi=200, subplot_kw={'projection': ccrs.PlateCarree(central_longitude=180)})
fig.subplots_adjust(hspace=0.3, wspace=0.3)
ax.coastlines(lw=0.7)
clevs = np.arange(-10, 10, 1)
im = ax.contourf(ds.lon, ds.lat, ds['t2m'],
clevs, cmap='coolwarm', extend='both', transform=ccrs.PlateCarree())
ax.set_xticks([ 60,120,180,240,300], crs=ccrs.PlateCarree())
ax.set_yticks([ -90,-60, -30, 0, 30, 60,90], crs=ccrs.PlateCarree())
lon_formatter = LongitudeFormatter(zero_direction_label=False)
lat_formatter = LatitudeFormatter()
ax.xaxis.set_major_formatter(lon_formatter)
ax.yaxis.set_major_formatter(lat_formatter)
clip_contours_by_map(im, map_polygon)
draw_map(map_polygon, color='k', linewidth=1)
ax.set_extent(map_polygon.get_extent(buffer=1))
plt.show()
通过掩膜数据的方法
import numpy as np
from cnmaps import get_adm_maps
import matplotlib.pyplot as plt
ds = xr.open_dataset('data.nc')
lon = ds.lon.values
lat = ds.lat.values
lons, lats = np.meshgrid(lon, lat)
map_polygon = get_adm_maps(country='中华人民共和国', record='first', only_polygon=True)
maskout_data = map_polygon.maskout(lons, lats, ds['t2m'])
fig, ax = plt.subplots(1, 1, figsize=(16, 12), dpi=200, subplot_kw={'projection': ccrs.PlateCarree(central_longitude=180)})
fig.subplots_adjust(hspace=0.3, wspace=0.3)
ax.coastlines(lw=0.7)
clevs = np.arange(-10, 10, 1)
im = ax.contourf(lons, lats, maskout_data,
clevs, cmap='coolwarm', extend='both', transform=ccrs.PlateCarree())
ax.set_xticks([ 60,120,180,240,300], crs=ccrs.PlateCarree())
ax.set_yticks([ -90,-60, -30, 0, 30, 60,90], crs=ccrs.PlateCarree())
lon_formatter = LongitudeFormatter(zero_direction_label=False)
lat_formatter = LatitudeFormatter()
ax.xaxis.set_major_formatter(lon_formatter)
ax.yaxis.set_major_formatter(lat_formatter)
clip_contours_by_map(im, map_polygon)
draw_map(map_polygon, color='k', linewidth=1)
ax.set_extent(map_polygon.get_extent(buffer=1))
plt.show()
联系我加入交流群
本文由 mdnice 多平台发布
2147
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
