之前处理 NetCDF 一般采用netCDF4包进行处理，近期处理降雨数据时发现，可采用xarray包进行处理。

  • netCDF4：是一个直接操作 NetCDF 文件的库，它提供了读取和创建 NetCDF 文件的基本接口。NetCDF（Network Common Data Form）是一种用于存储和分发科学数据的文件格式，特别适合于存储多维数组（如时间序列、空间图像和多维空间数据）。netCDF4库利用了 HDF5 的一些特性，如压缩、分块和高级数据模型，以提高存储和访问效率。
  • xarray ：是基于pandas的数据结构构建的，专门用于处理标签化的多维数组数据。它提供了一个简单易用的数据模型和分析工具，可以非常方便地处理 NetCDF 数据。xarray的数据模型是围绕着Dataset和DataArray两个核心数据结构设计的，它们可以存储原始数据和坐标信息，并提供了丰富的数据操作方法，如选择、分组、合并等。

根据经纬度范围裁剪 nc 文件

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
import os
import xarray as xr
import numpy as np

dir = 'data/'
# 定义经纬度范围，对nc数据进行裁剪
# 108.5706773°东 29.8471430°北, 108.7320040°东 29.9436699°北
lon_range = (108.50, 108.74)
lat_range = (29.95, 29.80)

fileList = os.listdir(dir)

print('Files in directory: ', len(fileList))

# %% 读取nc文件
for file in fileList:
    nc = xr.open_dataset(dir + file, decode_times=False, engine='netcdf4')
    print(nc)
    # 选择经纬度范围, longitude, latitude
    nc_new = nc.sel(longitude=slice(*lon_range), latitude=slice(*lat_range))
    nc_new.to_netcdf('nc_crop/'+ file)
    print('Processed file: ', file)

根据计算极端降雨量

依据前 30 年的年降雨数据，按 Gumbel 分布估算 10 年一遇、20 年一遇、50 年一遇和 100 年一遇的极端降雨工况下的降雨分布，并生成 nc 文件。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
#%% 导入包

import netCDF4 as nc
import numpy as np
import os
import time
import multiprocessing as mp

def Gumble_func(data,N,length):
    #  data应为一列数据
    #  N为多少年一遇
    #  length为所用数据年数(从大到小排序后希望获得前几个数据)
    #  Recurrence_period_value为N年一遇值
    # data = pd.read_excel(file_path)
    # data = np.array(data)
    # N = 100
    # length = 100
    if np.isnan(data).any():
        return np.nan
    datasort = np.sort(data, axis=0)[::-1] # 从大到小排列
    data_select_sorts = datasort[0:length - 1]
    mu = np.mean(data_select_sorts)
    sum_mm = 0
    for data_select_sort in data_select_sorts:
        sum_mm = sum_mm + (data_select_sort - mu) ** 2
    sigma = np.sqrt(sum_mm / (len(data_select_sorts) - 1))  # 标准差
    alpha = np.pi / (np.sqrt(6) * sigma)
    u = mu - 0.57721 / alpha
    for i in np.arange(1, 99999, 0.1):  # 降水量设置为1~99999mm
        if 1 / N >= 1 - np.exp(-np.exp(-alpha * (i - u))):  # N年一遇
            Recurrence_period_value = i  # Recurrence_period_value即为N年一遇降水值
            # print(str(N)+'年一遇降水值：',Recurrence_period_value,'毫米')
            break
    return Recurrence_period_value

def cal_N_year_appear_single(YEARLY_RAINFALL, i, j, N=10):
    rainfall_history = YEARLY_RAINFALL[:,i,j].flatten()
    if np.isnan(rainfall_history).any():
        return (i,j,np.nan)
    return (i,j,Gumble_func(rainfall_history, N, rainfall_history.shape[0]))

def cal_N_year_appear(YEARLY_RAINFALL, N=10):
    yearly_rainfall_N = np.empty(YEARLY_RAINFALL.shape[1:])
    # 采用多线程加速计算

    # 找到所有非nan的位置
    idx = np.where(~np.isnan(YEARLY_RAINFALL[0,:,:]))
    print('Number of non-nan elements: ', len(idx[0]))
    with mp.Pool(processes=6) as pool:
        results = pool.starmap(cal_N_year_appear_single, [(YEARLY_RAINFALL, i, j, N) for i, j in zip(idx[0], idx[1])])
        # wait for all processes to finish
        # pool.join()
        pool.close()
    for i, j, result in results:
        yearly_rainfall_N[i,j] = result
    return yearly_rainfall_N

def main():
    dir = 'nc_crop/'

    fileList = os.listdir(dir)

    print('Files in directory: ', len(fileList))

    #%% 读取nc文件

    yearly_rainfall_sum = []

    src = None

    for file in fileList[:3]:
        print('Processing file: ', file)
        data = nc.Dataset(dir + file)
        print('Variables in file: ', data.variables)
        print('Dimensions in file: ', data.dimensions)
        print('Shape of variable: ', data.variables['time'].shape)
        print('Data type of variable: ', data.variables['time'].dtype)
        print('Values of variable: ', data.variables['time'][:])
        # 读取pre数据
        pre = np.array(data.variables['pre'][:])
        print('Shape of pre: ', pre.shape)
        # 从文件名中提取年份
        year = int(file.split('_')[-1].split('.')[0])
        print('Year: ', year)
        # 根据年份计算每个月的天数, 需要考虑闰年
        monthly_days = [31, 28 if year % 4 != 0 else 29, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31]
        # 针对time维度, 与monthly_days进行加权求和
        pre_sum = np.dot(pre.transpose((1,2,0)), monthly_days)
        # pre_sum = np.sum(pre, axis=0)
        print('Shape of pre_sum: ', pre_sum.shape)
        yearly_rainfall_sum.append(pre_sum)
        print('-----------------------------------')
        src = data

    yearly_rainfall_sum = np.array(yearly_rainfall_sum)

    # %% 后处理

    # 将年降水量求平均值，写入文件
    print('Shape of yearly_rainfall_sum: ', yearly_rainfall_sum.shape)
    yearly_rainfall_sum_mean = np.mean(yearly_rainfall_sum, axis=0)
    print('Shape of yearly_rainfall_sum_mean: ', yearly_rainfall_sum_mean.shape)
    # 计算N年一遇降水量, N=10, 20, 50, 100

    # 计算10年一遇降雨量
    print('开始计算10年一遇降雨量...')
    yearly_rainfall_10 = cal_N_year_appear(yearly_rainfall_sum, 10)
    print('重现期为10年的平均降水量: ', yearly_rainfall_10.mean())
    # 计算20年一遇降雨量
    print('开始计算20年一遇降雨量...')
    yearly_rainfall_20 = cal_N_year_appear(yearly_rainfall_sum, 20)
    print('重现期为20年的平均降水量: ', yearly_rainfall_20.mean())
    # 计算50年一遇降雨量
    print('开始计算50年一遇降雨量...')
    yearly_rainfall_50 = cal_N_year_appear(yearly_rainfall_sum, 50)
    print('重现期为50年的平均降水量: ', yearly_rainfall_50.mean())
    # 计算100年一遇降雨量
    print('开始计算100年一遇降雨量...')
    yearly_rainfall_100 = cal_N_year_appear(yearly_rainfall_sum, 100)
    print('重现期为100年的平均降水量: ', yearly_rainfall_100.mean())

    # %% 写入文件
    # 写入文件
    out = nc.Dataset('yearly_rainfall_%s.nc'%time.time(), 'w', format='NETCDF4')
    # 维度与原始数据一致
    out.createDimension('latitude', src.dimensions['latitude'].size, )
    out.createDimension('longitude', src.dimensions['longitude'].size)
    # 设置经纬度变量
    lat=out.createVariable('latitude', 'f4', ('latitude',))
    out.variables['latitude'][:] = src.variables['latitude'][:]
    long=out.createVariable('longitude', 'f4', ('longitude',))
    out.variables['longitude'][:] = src.variables['longitude'][:]
    lat.units = 'degreesN'
    long.units = 'degreesE'
    # 创建变量
    out.createVariable('yearly_rainfall_sum', 'f8', ('longitude','latitude'))
    out.variables['yearly_rainfall_sum'][:] = yearly_rainfall_sum_mean
    out.createVariable('yearly_rainfall_10', 'f8', ('longitude','latitude'))
    out.variables['yearly_rainfall_10'][:] = yearly_rainfall_10
    out.createVariable('yearly_rainfall_20', 'f8', ('longitude','latitude'))
    out.variables['yearly_rainfall_20'][:] = yearly_rainfall_20
    out.createVariable('yearly_rainfall_50', 'f8', ('longitude','latitude'))
    out.variables['yearly_rainfall_50'][:] = yearly_rainfall_50
    out.createVariable('yearly_rainfall_100', 'f8', ('longitude','latitude'))
    out.variables['yearly_rainfall_100'][:] = yearly_rainfall_100
    out.close()
    print('File written successfully.')

if __name__ == '__main__':
    print('Start processing...')
    main()
    print('Processing finished.')