Python对数据白化

雨落森林

比如，当画东亚的温度等值线的图的时候，如果只想画中国区域的，那么这个过程叫做白化。用Python做白化一直都是个稍微麻烦一点的棘手事。目前已经有很多人提出了白化方法，使用最为广泛的是maskout方法，可以参见这个帖子[经验总结] Python完美白化，然后有很多人提出了各种补充，比如[经验总结] python白化去除区域外标签,去除区域外风矢量，还有这个[经验总结] 如何更优雅地绘制/白化中国的地图，可以更加方便好用地白化地图，很多时候也完全够用。
但是有些时候，比如，如果我想要对某个省，计算区域平均，这个就很不好办了，因为上面那些白化的原理是把画出来的图，裁掉区域外的，它是对图动手，不改变数据。所以我个人更倾向于对数据动手的白化，而且目前我更喜欢制作地图掩码的nc文件，这个文件是一个只有1和nan的数组，使用的时候直接乘以数据就可以了，在区域内的部分乘以的是1，所以数据没变，在区域外的部分乘以的是nan，就被mask掉了。然后后期画图或者计算的时候直接读取这个文件，用最临近插值调整到自己的研究区域和自己的分辨率，然后乘一下插值后的掩码就完事了。包括如果要计算区域的平均，直接np.nanmean就可以了。

首先，对于一个shp文件，我们先手动将其打开，然后画出来，看看它大概的位置，然后用稍微大一点范围（能覆盖掉那个区域）的格点，用每个具体的经纬度来判断是不是在形状内，这个方法参见[url=http://bbs.06climate.com/关于Python做地图白化的又一思路%20http://bbs.06climate.com/forum.p ... &fromuid=120047][经验总结] Python白化的又一思路[/url]，然后是把这个格点矩阵根据True或者False转换成1和nan，然后创建nc数据并保存，方便以后调用。


代码如下，注意看代码红色部分的代表的范围：
main.py是制作掩码文件的
from matplotlib.patches import Polygon
import matplotlib.pyplot as plt
import xarray as xr
import numpy as np
from shapefile import Reader
def make_mask(shpfiles, lon, lat, encoding='gbk', radius=1e-4, mesh=True):
    if isinstance(shpfiles, str):
        shpfiles = [shpfiles]
    X, Y = np.meshgrid(lon, lat) if mesh else (lon, lat)
    shp = X.shape
    points = np.stack([X.flatten(), Y.flatten()], 1)
    masks = []
    for file in shpfiles:
        f = Reader(file, encoding=encoding)
        pts = np.array(f.shape().points)
        p = Polygon(pts)
        idx = p.contains_points(points, radius=radius)
        masks.append(~idx.reshape(shp))
    return np.stack(masks).all(0)
   
def shp2da(shpfiles, lon, lat, **kwargs):
    mask = make_mask(shpfiles, lon, lat, **kwargs)
    da = xr.DataArray(np.where(mask, np.nan, 1), coords=dict(lon=lon, lat=lat),
                      dims=("lat", "lon"),name="mask")
    return da
   
if __name__ == "__main__":
    files = ["云南", "广东", "广西", "海南", "贵州"]
    files = [f"./shaps/{file}.shp" for file in files]
    lon, lat = np.arange(96, 118.1, 0.1), np.arange(18, 30.01, 0.1)
    import oyl
    m = oyl.map([96, 118, 0.1], [18, 30, 0.1], coast=False)
    for file in files:
        m.add_shape(file, method=1, edgecolor='red')
    mask = shp2da(files, lon, lat)
    mask.to_netcdf("./mask.nc")
    m.contourf(mask.data)
    m.show()
   
test.py是画图测试用的
from oyl import np, plt, xr, map
ds = xr.open_dataset("./mask.nc")
x, y = np.arange(90, 120.1, 0.15), np.arange(10, 40.1, 0.15)
data = np.arange(len(x)*len(y)).reshape(len(y), -1)
m = map([90, 120, 0.2], [10, 40, 0.2])
m.figure(1, [4.8, 6])
m.subplot(2,1,1)
m.contourf(x, y, data)
m.subplot(2,1,2)
m.contourf(x, y, data*ds.mask.interp(lon=x, lat=y, method='nearest'))
m.show()


见测试图：

yqcxm

很好的思路，沙发先

dhl-521

https://mp.weixin.qq.com/s/srmwRdlRCcqNXRI8j7f5ZA

雨落森林

dhl-521 发表于 2022-11-16 09:47
https://mp.weixin.qq.com/s/srmwRdlRCcqNXRI8j7f5ZA

确实，这篇文章也很不错

灭火器

提前做好mask并用最邻近插值调整形状确实方便。
另外推荐

1. import shapefile
   
2. import shapely.geometry as sgeom
   
3. from shapely.vectorized import contains
   
4. 
   
5. with shapefile.Reader('xxx.shp') as reader:
   
6.     geom = sgeom.shape(reader.shape(0))
   
7. 
   
8. X, Y = np.meshgrid(lon, lat)
   
9. mask = contains(geom, X, Y)

复制代码

相比matplotlib.patches速度更快，并且能处理带洞的shp文件

雨落森林

灭火器 发表于 2022-11-17 17:55
提前做好mask并用最邻近插值调整形状确实方便。
另外推荐
相比matplotlib.patches速度更快，并且能处理带 ...

好的，感谢指教

Pearl王川

顶呱呱，

雨落森林

Pearl王川 发表于 2022-11-19 17:23
顶呱呱，

{:eb511:}{:eb511:}