基于Sen+MK的多站点不同季节和年尺度的SPEI趋势分析.md

再大的风浪，不过只短暂喧哗。

---

前言

1. 概述

2.1 问题情景

• 假设目前我们有若干个站点的春、夏、秋、冬和四季的spei，需要利用Sen+MK方法计算各站点四季和年尺度的SPEI的趋势。

2.2 说明

• 测试数据和程序可以点击基于Sen+MK的多站点不同季节和年尺度的SPEI趋势分析下载。
• 请尊重知识产权，切勿复制本程序和示例数据到闲鱼、淘宝等平台，一经发现，必将追究相关责任。

2. 版本

2.1 天津，2024年5月4日，Version1

3. 微信公众号GISRSGeography

• 欢迎大家关注公众号 GISRSGeography，谢谢！。
  

---

一、数据

1. 输入数据

• 存储在同一个excel文件中的三个站点的四季和年尺度的SPEI。
  

2. 输出数据

• 存储在一个excel文件中的三个站点的四季和年的SPEI的SEN+MK的趋势分析结果。
  

二、程序代码

# -*- coding: utf-8 -*-
"""
1. 程序目的
   (1) 基于SEN+MK分析不同站点四季和年的干旱趋势

2. 2024年5月4日 Version 1

3. 数据
   3.1 输入数据
       '...\JJJ_SensonYearSPEI.xlsx'
       
   3.2 输出数据
       '...\JJJ_SensonYearSPEI_SEN.xlsx'
   
"""
# %% 包的导入
# 包的导入
import numpy as np
import pandas as pd

import pymannkendall as mk

# %% 四季和年spei数据读取
def spei_get(
        inpath: str
        ) -> list:
    """
    (1) 功能：
        1) 读取四季和年的SPEI
    ----------
    
    (2) 输入参数
        1) inpath: str
           文件存储绝对路径
    ----------
      
    (3) 输出参数
        1) spei_data: DataFrame
           spei数据
  
    """
    spei_data = pd.read_excel(inpath)
    
    return spei_data

# %%   
def spei_trend(
        speidata: pd.DataFrame,
        styr: int,
        edyr: int
        )  -> pd.DataFrame:
    """
    (1) 功能
        1) 提取目标时段的spei
        2) 计算不同时段的spei的sen趋势，并利用mk进行显著性检验
    ----------
    
    (2) 输入参数
    speidata : pd.DataFrame
               存储着目标站点四季和年尺度的SPEI
    -------
    
    (3) 输出数据
    spei_senslope : pd.DataFrame
                    存储着目标站点四季和年尺度的SPEI的Sen+MK的趋势分析结果

    """
    # 目标时段数据选取
    speidata_aim = speidata.loc[(speidata['Year']>=styr) & (speidata['Year']<=edyr),:]
    
    # 特定站点数据提取并调用pymannkendall进行趋势计算和显著性检验
    seasonyear = ['Spring','Summer','Autumn','Winter','YearSPEI']
    sen_p = ['_SenSlope','_PValue']
    columns = ['Sta_ID','Period']
    for seasonyear_temp in seasonyear:
        for senp_temp in sen_p:
           columns.append(seasonyear_temp+senp_temp) 
    
    spei_senslope = pd.DataFrame([],columns=columns)
    sta_all = np.unique(speidata_aim.Sta_ID)
    ii = 0
    for sta_temp in sta_all:
        # 站点数据提取
        spei_sta = speidata_aim.loc[speidata_aim['Sta_ID']==sta_temp,:] 
        
        # Sen趋势度计算和显著性检验
        for seasonyear_temp in seasonyear:
            spei_sta_seasonyear = spei_sta[seasonyear_temp]
            
            senmkresult = mk.original_test(spei_sta_seasonyear)
            spei_senslope.loc[ii,seasonyear_temp+'_SenSlope']=senmkresult.slope
            spei_senslope.loc[ii,seasonyear_temp+'_PValue']=senmkresult.p
            #print(senmkresult)
        
        spei_senslope.loc[ii,'Sta_ID'] = sta_temp
        
        ii = ii + 1
        
        print(ii)
            
    spei_senslope['Period'] = str(styr) + '-' + str(edyr)
    
   
    return spei_senslope

                     
# %%
if __name__ == '__main__':
    
# %% 路径处理和变量预定义
    infile = r'...\JJJ_SensonYearSPEI.xlsx'     # 需要修改为自己的路径
    outfile = r'...\JJJ_SensonYearSPEI_SEN.xlsx' # 需要修改为自己的路径
    
    styr=1980;
    edyr=2020;
       
# %% spei数据读取
    spei_data  = spei_get(infile)
    
# %% Sen+MK
    spei_sen = spei_trend(spei_data,styr,edyr)
    
# %% 结果写出
    spei_sen.to_excel(outfile,index=False,sheet_name='JJJ')
    
    print('Finished.')

三、参考资料

1. pymannkendall 1.4.3