使用Python从逐月降水量数据中获取某位置的全年降水量

1 数据下载

可以在中国科技资源共享网中进行下载1km分辨率逐月降水量数据集。

2 数据说明

该数据为.nc格式，里面有lat、lon、time和pre四个数据块，其中lat、lon是[H, W]形状的矩阵，每一个像素代表该像素对应坐标的经纬度；time有12个，分别为1-12月；pre是对应的降水数据块，形状为[12, H, W]，代表某个像素12个月的降雨量。

3 数据处理

我们需要的是获取一些经纬度点的12个月降雨量，因此可以使用Python的netCDF4进行读取，然后进行处理，代码如下：

import numpy as np
import netCDF4 as nc


class WaterReader:
    def __init__(self, nc_path):
        self.src = nc.Dataset(nc_path)
        self.data = np.asarray(self.src["pre"])
        self.rang = self._get_range()
        self.lon_step = (self.rang[1] - self.rang[0]) / len(self.src["lon"])
        self.lat_step = (self.rang[3] - self.rang[2]) / len(self.src["lat"])

    def _get_range(self):
        return (
            np.min(self.src["lon"]),
            np.max(self.src["lon"]),
            np.min(self.src["lat"]),
            np.max(self.src["lat"]))
    
    def _get_rowcol_by_lonlat(self, lon, lat):
        lon_min, lon_max, lat_min, lat_max = self.rang
        if lon < lon_min or lon > lon_max or lat < lat_min or lat > lat_max:
            raise ValueError("lon or lat is out of range")
        lon_index = int((lon - lon_min) / self.lon_step)
        lat_index = int((lat_max - lat) / self.lat_step)
        return lon_index, lat_index
    
    def get_water(self, lon, lat):
        lon_index, lat_index = self._get_rowcol_by_lonlat(lon, lat)
        return self.data[:, lat_index, lon_index]

if __name__ == "__main__":
    wr = WaterReader("pre_2021.nc")
    print(wr.get_water(104.4214641, 33.731619))

其中关键在于通过_get_rowcol_by_lonlat将输入的经纬度对应找到图上最近的像元，然后通过该像元的行列号找到降水量。这里还可以插找距离最近的四个像元，通过线性插值的方法得到该点的12个月降水量。

如果需要批量处理的话，写成循环就好了。