求一个从全球数据提取中国区域数据进行降水-时间变化分析的脚本

顺利毕业！

小白想从一个全球降水数据（16年）提取出中国区域的数据，然后进行中国区域的格点平均，做一个降水随时间变化的图，看了很多论坛的帖子，不知道该用用mask函数还是用shapefile_mask_data函数，菜鸡琢磨捣鼓了快两个星期仍然毫无进展，看了很多还是不知道该从何下手，跪求大佬提供一个可参考的脚本，救救孩子！！！！

dhl-521

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import xarray as xrimport shapely.geometry as sgeom
import shapefile
import math
from shapely.prepared import prep

shp_file = r'F:\data\XXXXXXXXX.shp'  ####掩膜的shp文件
path_file = r'F:\XXXXXXXXX.nc'  ####需要读取的数据
#########----用shp创建掩膜文件----#####
def polygon_to_mask(polygon, x, y):
    '''生成落入多边形的点的掩膜数组.'''
    x = np.atleast_1d(x)
    y = np.atleast_1d(y)
    if x.shape != y.shape:
        raise ValueError('x和y的形状不匹配')
    prepared = prep(polygon)

    def recursion(x, y):
        '''递归判断坐标为x和y的点集是否落入多边形中.'''
        xmin, xmax = x.min(), x.max()
        ymin, ymax = y.min(), y.max()
        xflag = math.isclose(xmin, xmax)
        yflag = math.isclose(ymin, ymax)
        mask = np.zeros(x.shape, dtype=bool)

        # 散点重合为单点的情况.
        if xflag and yflag:
            point = sgeom.Point(xmin, ymin)
            if prepared.contains(point):
                mask[:] = True
            else:
                mask[:] = False
            return mask

        xmid = (xmin + xmax) / 2
        ymid = (ymin + ymax) / 2

        # 散点落在水平和垂直直线上的情况.
        if xflag or yflag:
            line = sgeom.LineString([(xmin, ymin), (xmax, ymax)])
            if prepared.contains(line):
                mask[:] = True
            elif prepared.intersects(line):
                if xflag:
                    m1 = (y >= ymin) & (y <= ymid)
                    m2 = (y >= ymid) & (y <= ymax)
                if yflag:
                    m1 = (x >= xmin) & (x <= xmid)
                    m2 = (x >= xmid) & (x <= xmax)
                if m1.any(): mask[m1] = recursion(x[m1], y[m1])
                if m2.any(): mask[m2] = recursion(x[m2], y[m2])
            else:
                mask[:] = False
            return mask

        # 散点可以张成矩形的情况.
        box = sgeom.box(xmin, ymin, xmax, ymax)
        if prepared.contains(box):
            mask[:] = True
        elif prepared.intersects(box):
            m1 = (x >= xmid) & (x <= xmax) & (y >= ymid) & (y <= ymax)
            m2 = (x >= xmin) & (x <= xmid) & (y >= ymid) & (y <= ymax)
            m3 = (x >= xmin) & (x <= xmid) & (y >= ymin) & (y <= ymid)
            m4 = (x >= xmid) & (x <= xmax) & (y >= ymin) & (y <= ymid)
            if m1.any(): mask[m1] = recursion(x[m1], y[m1])
            if m2.any(): mask[m2] = recursion(x[m2], y[m2])
            if m3.any(): mask[m3] = recursion(x[m3], y[m3])
            if m4.any(): mask[m4] = recursion(x[m4], y[m4])
        else:
            mask[:] = False

        return mask

    return recursion(x, y)

##########-----读取shp----######
with shapefile.Reader(shp_file) as reader:
    china = sgeom.shape(reader.shape(0))
#########

data = xr.open_dataset(path_file)   ###读取数据
pre_mean = data.pre                     ####读取降水数据
lon = pre_mean.lon                       ####读取经纬度
lat = pre_mean.lat                        #####读取经纬度
X,Y = np.meshgrid(lon,lat)      
mask_check = polygon_to_mask(qinghai, X, Y)      ####返回mask
masked_check = np.ma.masked_equal(mask_check, False)  ####
mask_mean = pre_mean*masked_check    ####掩膜
mean_pre = np.nanmean(mask_mean)  ####求区域平均值


原贴地址：   https://mp.weixin.qq.com/s/srmwRdlRCcqNXRI8j7f5ZA

顺利毕业！

dhl-521 发表于 2022-11-18 09:14
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

感谢！！！！ 但是我需要的NCL脚本呜呜

小余

尝试做一个只有中国区域的nc文件（0，1），然后对数据进行对应的处理就可以了。利用这个文件去和你现有的数据相乘即可得到中国区域的数据，需要注意经纬度需要一致

来自苹果APP客户端

顺利毕业！

小余 发表于 2022-11-19 15:05
尝试做一个只有中国区域的nc文件（0，1），然后对数据进行对应的处理就可以了。利用这个文件去和你现有的数 ...

有相应的脚本可以参考吗