日分辨率的干旱指标SPEI计算，干旱事件识别

置顶探寻TUT

已于 2024-09-15 12:09:49 修改

阅读量3.2k

点赞数 47

文章标签：人工智能 python 开发语言

于 2024-09-15 11:57:06 首次发布

本文链接：https://blog.csdn.net/weixin_44021581/article/details/142212588

版权

主要任务

SPEI是最常用的干旱指标，考虑了降雨和潜在蒸散发的水平衡状况，并通过不同时间尺度上累积水平衡状况反映不同时间长度的干旱情况，具体来说3个月尺度的SPEI反映农业、土壤干旱，6个月尺度的SPEI反映水文干旱。

在现有能搜到的SPEI计算代码中，R库SPEI包可以计算月分辨率的SPEI，python库的Climate_indices包也是计算月分辨率的SPEI，没有公开的代码计算日分辨率的SPEI。考虑计算日分辨率的SPEI是因为，月分辨率的SPEI不能捕捉持续仅几周的短期干旱事件，不能精确捕捉草地生产力变化情况。为此，一些文章提到了构建日分辨率的SPEI指标，如Wang et al.,2015 (https://doi.org/10.1002/joc.4244), 李军（https://doi.org/10.5194/hess-25-1587-2021），以捕捉短期干旱事件。由于具体的代码文章的作者没用公布，本文主要目的是介绍如何计算日分辨率的SPEI，帮助广大计算日分辨率SPEI指标。

得到SPEI后，通常可以分析一个地区的干湿趋势，同时也可以基于游程理论提取干旱事件（Redirecting）。这部分代码写起来也简单，但如果是初学者，不了解python可能会比较费劲。如果是matlab的朋友，可以去CHATGPT把这部分python代码转成matlab代码就行。

代码

完整的代码在我附上的资源连接：https://download.csdn.net/download/weixin_44021581/89755906。

下面简单介绍每一部分代码的作用：

【0】test_data文件夹，input为读者需要准备的气象数据，包括：降雨、温度、日照时长、风速、压强、湿度，output为后续程序计算出来的结果，包括净辐射、潜在蒸散发、SPEI。

【1】test_run_cal_spei _share.py：调用test_data文件夹中的气象数据，依次计算净辐射、潜在蒸散发、SPEI。供读者自定义的基础参数包括：文件根目录、纬度、海拔、气象数据的开始结束日期。与SPEI相关的参数包括：1. 计算潜在蒸散发的方式；2. 需要的分辨率，即日、8day、月三种；3. 不同的累计水平衡的时间，即不同的时间尺度；4.SPEI的拟合分布及拟合分布的参数的优化方法。

###################### global variables #####################
root_path = r'D:\0_pro_Data\Clim\SPEI\SPEI_code_share'  # 根目录
lat = 46.7  # 纬度
zalt = 1000 # 海拔
start_date = '1982-01-01' # test data第一个值的日期
end_date = '2019-12-31'

PET_type = 'PM' # 提供 PM公式和 HG公式两种方式计算日蒸散发
periodicity = 'daily'  # 提供不同分辨率 or '8day' or 'monthly'

if periodicity == 'daily':
    AccumulatedTime = 90  # 提供不同的累积水平衡的时间，农业干旱一般是3个月，could also be 30/60/90/120 days
elif periodicity == '8day':
    AccumulatedTime = 12  # 提供不同的累积水平衡的时间，农业干旱一般是3个月，could also be 4*8days, 12*8days days
elif periodicity == 'monthly':
    AccumulatedTime = 3   # 提供不同的累积水平衡的时间，农业干旱一般是3个月，could also be 3/6/9/12 months
SPEI_distribution_type = 'glo' # 提供两种SPEI的拟合分布：'glo' 广义逻辑分布；'gev'：广义极值分布
SPEI_param_optim_type = 'ub-pwm' # 提供两种SPEI拟合分布的最优参数优化方法： 'ub-pwm'和 'pb-pwm'

【2】utilsPY.py：包括数据预处理、净辐射、潜在蒸散发、SPEI具体计算四个部分的代码。

【3】draw_spei.py：后处理，绘制三个分辨率的SPEI的图像。

结果

1. 不同分辨率的SPEI结果的对比

不同分辨率的SPEI值基本相似

具体以2016年的干旱年为例，三种分辨率得到的SPEI轨迹相似，日SPEI和8daySPEI基本重合，Monthly SPEI的干旱发展略早一点，与月尺度拟合的数据少、误差大有关。

2. 不同潜在蒸散发计算结果的对比

两种公式得到的PET基本相近，R2=0.75，主要误差在于，HG公式未考虑风速的影响，往往在春夏风速大会低估潜在蒸散发的值。

尽管PET计算的幅值有差异，但是由于SPEI是计算前90天的累积水平衡，所以单独几天的低估并不会太影响SPEI的取值。2016年6月和2016年10月，HG公式下的SPEI更大，是因为在春秋季节风速的影响没有考虑，ET被低估，导致SPEI表征的水分匮缺程度是被低估的，和上面的ET误差原因分析一致。考虑到这点，在气象数据充足的区域，建议使用PM公式计算潜在蒸散发和SPEI。

3. 极端事件识别代码

    def pick_drought_events(SPEI, threshold):
        '''
        :param SPEI: SPEI time series
        :param threshold: SPEI threshold for extreme drought
        :return: a list of drought events
        '''
        drought_month = []
        for i, val in enumerate(SPEI):
            if val < threshold:
                drought_month.append(i)
            else:
                drought_month.append(-999999)
        events = []
        event_i = []
        for ii in drought_month:
            if ii > -99:
                event_i.append(ii)
            else:
                if len(event_i) > 0:
                    events.append(event_i)
                    event_i = []
                else:
                    event_i = []
        return events