Python | 降水概率密度分布图

写在前面

最近终于把之前的工作收尾了，就等老外再看一遍了。开始一个新的项目，从复现老外的文章开始。

主要包括：

• 1、降水概率密度分布
• 2、根据个别日降水率计算年平均降水量和热带平均降水量

• 3、从2的结果中计算累计降水

降水概率密度分分布

第一个较为简单，通过设置不同的bins，统计每个bins中降水的分布即可。

BIN_EDGES_MM_DAY = (
0.0, 0.1, 1, 2, 3, 4, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50,
60, 70, 80, 90, 100, 125, 150, 175, 200, 225, 250, 275, 300,

)
pdf_density, bin_edges = np.histogram(ds, bins=BIN_EDGES_MM_DAY, density=True)
plt.plot(BIN_EDGES_MM_DAY[:-1], pdf_density, drawstyle='steps-post')
plt.xscale('log')
plt.yscale('log')
plt.xlabel('Precipitation (mm/day)')
plt.ylabel('Probability Density')
plt.title('Precipitation Probability Density Function')
plt.grid(True, which='both', linestyle='--', linewidth=0.5)
plt.show()

每个降水区间对总降水的贡献占比 = 每个区间的总降水 / 总降水

本质上是为了，基于之前分的不同降水区间，去计算每个降水区间贡献了多少总降水：降水率落在 τ 到 τ+δτ 这一区间内时，对全年热带平均降水的贡献，这可以回答以下问题：

• 总降水主要来自哪个强度区间？
• 小雨、中雨、暴雨，谁对总降水贡献最大？

上述公式的分子A（lon,lat,time）类似一个掩膜的作用，如果某个格点在某一天的降水数值满足τ 到 τ+δτ，例如[0,1],那么将这个格点mask为1，否则为0，目的就是为了只保留该降水区间内降水格点，其他都置为0.这是为了方便统计这一降水区间[0，1]对热带平均降水的贡献。

在此基础上，需要进一步计算每个降水区间实际占据了多少面积，所以这里需要计算实际的经纬度网格的格点之间的球面积，考虑到lat不同范围进行cos加权。

defcal_area(dset, x_step=None, y_step=None):
"""Calculate the horizontal area of each lat-lon grid box.
    dx=RcosϕΔλ
    dy=RΔϕ
    area=dx*dy=R^2cosϕΔλΔϕ
    delta x = x_step (degree) x pi/180 x Earth radius
    delta y = y_step (degree) x pi/180 x Earth radius
    """
if x_step isNone:
        x_step = float(np.abs(np.diff(dset.lon.values)).mean())
if y_step isNone:
        y_step = float(np.abs(np.diff(dset.lat.values)).mean())

    earthr = 6.371e6
    lon2d, lat2d = np.meshgrid(dset.lon.values, dset.lat.values)
    area = xr.DataArray(
        data=(x_step / 180 * np.pi) * earthr * np.cos(np.deg2rad(lat2d)) * (y_step / 180 * np.pi) * earthr,
        dims=["lat", "lon"],
        coords={"lat": dset.lat.values, "lon": dset.lon.values},
    )
return area

area = cal_area(ds)
precip_values = ds.values.reshape(-1) # flattern to 1D array
area_values = np.broadcast_to(area.values, ds.shape).reshape(-1) # flattern
area_total, _ = np.histogram(precip_values, bins=BIN_EDGES_MM_DAY, weights=area_values) 

再次调用histogram函数，将面积权重加进去，得到每个降水区间占了多少面积，预期统计结果大概就是：

• 0–1 mm/day 的降水，对应总面积是多少
• 1–5 mm/day 的降水，对应总面积是多少
• 5–10 mm/day 的降水，对应总面积是多少 以下结果很显而易见，小雨部分占据了最多的面积

进一步的，但是我们需要思考的是，小雨占据了最大的空间面积，他对于总降水的贡献就是最重要的吗？

• 计算每个降水区间对应的总降水量贡献，计算落进不同降水区间以后加多少权重，比如得到：
• 20–30 mm/day 的区域一共下了多少雨？

precip_total, _ = np.histogram(
        precip_values,
        bins=   BIN_EDGES_MM_DAY,
        weights=precip_values * area_values,
    )     
plt.plot(BIN_EDGES_MM_DAY[:-1], precip_total, drawstyle='steps-post')
plt.xscale('log')
plt.xlabel('Precipitation (mm/day)')
plt.ylabel('Probability Density')
plt.title('Precipitation Probability Density Function')
plt.grid(True, which='both', linestyle='--', linewidth=0.5)
plt.show()

• 贡献更大的反而是5-10mm/day左右的降水，对于总降水贡献最大。

进一步计算，第 i 个降水区间对区域平均降水的贡献。

total_area = area_values.sum() # 计算总面积
bin_mean_contribution = precip_total / total_area 
plt.plot(BIN_EDGES_MM_DAY[:-1], bin_mean_contribution, drawstyle='steps-post')
plt.xscale('log')
plt.xlabel('Precipitation (mm/day)')
plt.ylabel('Contribution to Total Precipitation')
plt.title('Contribution of Each Precipitation Bin to Total Precipitation')
plt.grid(True, which='both', linestyle='--', linewidth=0.5)
plt.show()

进一步计算：从小雨到大雨逐步累加时，各降水强度对总降水的累计贡献。cumulative_precip = np.cumsum(bin_mean_contribution)

绘图比较：