请选择 进入手机版 | 继续访问电脑版
爱气象,爱气象家园! 

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

新浪微博登陆

只需一步, 快速开始

搜索
查看: 4566|回复: 7

[经验总结] Python对数据白化

[复制链接]

新浪微博达人勋

发表于 2022-11-13 19:14:53 | 显示全部楼层 |阅读模式

登录后查看更多精彩内容~

您需要 登录 才可以下载或查看,没有帐号?立即注册 新浪微博登陆

x
本帖最后由 雨落森林 于 2022-11-13 19:18 编辑

比如,当画东亚的温度等值线的图的时候,如果只想画中国区域的,那么这个过程叫做白化。用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()


见测试图:


im.jpg

评分

参与人数 3金钱 +27 贡献 +5 收起 理由
xfz147256 + 2 很给力!
Pearl王川 + 20 + 5 很给力!
bearok999 + 5 很给力!

查看全部评分

密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2022-11-13 20:19:17 | 显示全部楼层
很好的思路,沙发先
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2022-11-16 09:47:37 | 显示全部楼层
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2022-11-17 11:51:42 | 显示全部楼层
dhl-521 发表于 2022-11-16 09:47
https://mp.weixin.qq.com/s/srmwRdlRCcqNXRI8j7f5ZA

确实,这篇文章也很不错
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2022-11-17 17:55:54 | 显示全部楼层
提前做好mask并用最邻近插值调整形状确实方便。
另外推荐
  1. import shapefile
  2. import shapely.geometry as sgeom
  3. from shapely.vectorized import contains

  4. with shapefile.Reader('xxx.shp') as reader:
  5.     geom = sgeom.shape(reader.shape(0))

  6. X, Y = np.meshgrid(lon, lat)
  7. mask = contains(geom, X, Y)
复制代码
相比matplotlib.patches速度更快,并且能处理带洞的shp文件
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2022-11-18 14:05:40 | 显示全部楼层
灭火器 发表于 2022-11-17 17:55
提前做好mask并用最邻近插值调整形状确实方便。
另外推荐
相比matplotlib.patches速度更快,并且能处理带 ...

好的,感谢指教
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2022-11-19 17:23:17 | 显示全部楼层
顶呱呱,
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2022-11-19 23:36:30 | 显示全部楼层

{:eb511:}{:eb511:}
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册 新浪微博登陆

本版积分规则

Copyright ©2011-2014 bbs.06climate.com All Rights Reserved.  Powered by Discuz! (京ICP-10201084)

本站信息均由会员发表,不代表气象家园立场,禁止在本站发表与国家法律相抵触言论

快速回复 返回顶部 返回列表