关于Python做地图白化的又一思路

雨落森林

首先，遇到地图白化问题，家园有帖子提出了解决办法，另外，“云台书使”公众号还汇总了常见的几种做法。比如maskout白化，但是没有源文件很麻烦，要说抓住核心语句，对填色或者等值线的collections做set_clip_path，但是我回头一看，它是用cartopy.io的部分读取shp文件的。有些以国标编码的shp不能用cartopy读取（因为不知道怎么设置编码方式），所以我使用shapefile来读shp。裁剪的核心函数set_clip_path需要的path，用Path(坐标点)是不能满足的，所以这种方法失败了。而且，还有一个问题，由于set_clip_path把一个shp外部屏蔽了，如果我要画多个shp，就会互相屏蔽，最后只剩下一个shp有效白化了，除非把很多shp拼接成一个。所以说，这种方法不好。另外两种方法要用到salem和geopandas，对于我这种懒人，懒得装，也懒得学。（而且以前在windows试着安装了geopandas，报错，装不成）。所以乍一看好像方法挺多，但实用的好像不多。今天我发现一个很重要的函数Polygon.contains_points,它可以判断点是否在区域内。So, 问题解决了。我完全可以暴力点，只用我自己的方法去做。不需要它cartopy的添加feature，我直接自己用Ploygon来画shp的地图，而且用它来生成我需要的填色的区域。我试着用2017年的月平均位势高度来画，白化青海、安徽和江苏三个省。废话不多少，看代码，看注释。另外，oyl库的map是一个画地理地图的库，本质是基于cartopy的，各位观众可以不必理会，只要知道它相当于跳过了定义底图就可，不影响阅读。


from oyl import map
import xarray as xr
import numpy as np
import matplotlib.pyplot as plt
from shapefile import Reader
import cartopy.crs as ccrs
x, y = np.arange(70,140.1,0.25),  np.arange(0,60.1,0.25)
##读取shp文件里面的点的坐标
xy1 = np.array(Reader("./shape/青海.shp",encoding='gbk').shape().points)
p1 = plt.Polygon(xy1,fill=False)
xy2 = np.array(Reader("./shape/江苏.shp",encoding='gbk').shape().points)
p2 = plt.Polygon(xy2,fill=False)
xy3 = np.array(Reader("./shape/安徽.shp",encoding='gbk').shape().points)
p3 = plt.Polygon(xy3,fill=False)
#oyl.map定义底图
m = map([70,140,0.25],[0,60,0.25])
m.subplot()
ax = m.ax
##读取nc数据
z=xr.open_dataset("./2017z.nc").z.data[0,::-1]/98
##填色
plt.colorbar(ax.contourf(x,y,z,cmap=plt.cm.gist_rainbow,levels=np.arange(560,593),extend='both',transform=ccrs.PlateCarree()))
##添加shp
ax.add_patch(plt.Polygon(xy1,transform=ccrs.PlateCarree(),fill=False))
ax.add_patch(plt.Polygon(xy2,transform=ccrs.PlateCarree(),fill=False))
ax.add_patch(plt.Polygon(xy3,transform=ccrs.PlateCarree(),fill=False))
plt.show()

##定义填色的数据的经纬度
xys = np.stack(np.meshgrid(x,y)).reshape(2,-1).T
##开始筛选满足区域里面的点
idx1 = p1.contains_points(xys,radius=0.1)
idx2 = p2.contains_points(xys,radius=0.1)
idx3 = p3.contains_points(xys,radius=0.1)
idx = (idx1|idx2|idx3).reshape(241,281)
##把不满足区域的点变成缺失值
z[~idx] = np.nan
##以下开始是重新再画一遍图，只不过因为有缺失值，所以相当于做了裁剪
m = map([70,140,0.25],[0,60,0.25])
m.subplot()
ax = m.ax
plt.colorbar(ax.contourf(x,y,z,cmap=plt.cm.gist_rainbow,levels=np.arange(560,593),extend='both',transform=ccrs.PlateCarree()))
ax.add_patch(plt.Polygon(xy1,transform=ccrs.PlateCarree(),fill=False))
ax.add_patch(plt.Polygon(xy2,transform=ccrs.PlateCarree(),fill=False))
ax.add_patch(plt.Polygon(xy3,transform=ccrs.PlateCarree(),fill=False))
plt.show()

雨落森林

oyl现在已经在pypi上发布了，可以pip install oyl
欢迎感兴趣的朋友参与体验

yqcxm

oyl库里哪里有呢？

雨落森林

yqcxm 发表于 2021-11-1 11:29
oyl库里哪里有呢？

我自己写的库，不用在意，看代码整体思路就行

笺疏

见教。不过这种方法适用于细网格数据（0.5×0.5或更细），遇上粗网格就会出现大面积挖白，遇上复杂的曲折shp也会出现大面积挖白。不过这种方法确实更明白些。

雨落森林

笺疏 发表于 2021-11-4 19:03
见教。不过这种方法适用于细网格数据（0.5×0.5或更细），遇上粗网格就会出现大面积挖白，遇上复杂的曲折sh ...

确实我现在也发现了{:cry:}{:cry:}

love2100

学习一下

vision2761

阳霖牛逼哦

雨落森林

vision2761 发表于 2021-11-22 15:46
阳霖牛逼哦

{:eb511:}{:eb511:}

夢之夜

阳霖np，准备用白化函数，结果在气象家园看到了你的白化函数