热闹的尽头是孤寂，在虚浮的欢闹中保持自己，纷繁世间，可报期望者不过二三。

---

前言

  此系列博文的目的是基于Python的Climate Indices库计算标准化降水指数（SPI）。

1. 概述

2.1 目的

• 针对不同站点数据，调用Climate Indices库进行不同时间尺度的SPI的计算

2.2 说明

• 测试数据和程序可以点击基于Pyton的Climate Indices库计算不同站点、不同时间尺度的SPI下载。

2. 版本

2.1 天津，2024年1月18日，Version1

3. 微信公众号GISRSGeography

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

---

一、数据

1. 输入数据

• 本示例的输入数据是北京、天津和石家庄三个站点的逐月气象数据。
  

2. 输出数据

• 本示例的输出数据是逐站点存储的不同时间尺度的SPI的计算结果。
  

二、程序代码

# -*- coding: utf-8 -*-
"""
1. 程序目的
   (1) 调用climate indices库计算不同站点、不同时间尺度的spi

2. 2024年3月26日

"""
# %% 包的导入
# 包的导入
import os
import numpy as np
import pandas as pd

from climate_indices import indices
from climate_indices import compute  # 包含计算SPI的包

# %% 获取某一路径下特定格式的文件
def get_file(
        inpath: str,
        type_str: str
        ) -> list:
    """
    (1) 功能：
        1) 获取某一路径下特定格式的文件
    ----------
    
    (2) 输入参数
        1) inpath: str
           文件存储绝对路径
        2) type_str: str
           文件后缀名
    ----------
      
    (3) 输出参数
        1) file_list: list
           特定文件格式的数据
           
    """
    file_items = os.listdir(inpath)
    
    file_list = []
    file_path_list = []
    for ii in file_items:
        if ii.endswith(type_str):
            file_list.append(ii)
            file_path_list.append(inpath+ii)
    return file_list,file_path_list

# %% 数据读取函数定义
def read_oridata(
        filepath: str
        ):
    
    """
    (1) 功能：
        读取用于计算SPI的气象数据
        
        注：此函数读取的是测试数据，也可以读取按照相同方式存储的气象数据
    ---
    
    (2) 输入参数：
        1) filepath: str,输入数存储路径
    ---
    
    (3) 输出参数
        1) sta_id: int,站点号
        2) styr: int, 开始年
        3) edyr: int, 结束年
        4) pre: np.array, 降水
       
    """
    # 数据读取
    climdata = pd.read_csv(filepath)
    
    # 站点号
    sta_id = climdata.Sta_ID[0]
    # 开始年和结束年
    styr = climdata.Year[0]
    edyr = max(climdata.Year)
    # 降水
    pre = np.asarray(climdata.Pre_2020)
    
    return sta_id,styr,edyr,pre
    

# %% 不同时间尺度spi计算函数定义
def dif_scale_spi(
        sta_id: int,
        pre_data: np.array,
        styr: int,
        edyr: int,
        scale: tuple
        ) -> pd.DataFrame:
    
    """
    (1) 功能：
        1) 计算单一站点不同时间尺度的spi
    ---
    
    (2) 输入参数
        1) sta_id: int, 站点号
        2) pre_data: np.array, 降水序列，月值
        3) styr: int, 开始年
        4) edyr: int, 结束年
        5) scale: tuple, spei的时间尺度
    ---
    
    (3) 返回值
        1) spi_df: pd.DataFrame, 不同时间尺度的SPI
           列名：站点号，年，月，不同时间尺度的SPI
           
        注：SPI计算结果中的-99表示无效值
    
    """
    
    # 创建存储计算结果的Dataframe(站点号，年，月，不同时间尺度的SPEI)
      # 创建列名
    scale = sorted(scale)
    scale_list = []
    for scale_temp in scale:
        if scale_temp < 10:
            scale_list.append('Scale_0'+str(scale_temp))
        else:
            scale_list.append('Scale_'+str(scale_temp))
    col_names = ['Sta_ID','Year','Month'] + scale_list
      # 空的数据框的创建
    spi_df = pd.DataFrame(data=[],columns=col_names)
    
    # 不同时间尺度的spi的计算
    for scale_temp in scale:
        spi = indices.spi(values=pre_data,
                    scale=scale_temp, # 时间尺度
                    distribution=indices.Distribution.gamma,
                    periodicity=compute.Periodicity.monthly,
                    data_start_year=styr,
                    calibration_year_initial=styr,
                    calibration_year_final=edyr
                    )
        
        spi[np.isnan(spi)] = -99 # nan转换为-99
        if scale_temp < 10:
            spi_df['Scale_0'+str(scale_temp)]=spi
        else:
            spi_df['Scale_'+str(scale_temp)]=spi
    
    # 站点信息补充
    spi_df['Sta_ID'] = sta_id
    
    # 生成年月序列
    year_list = []
    month_list = []
    
    for year in range(styr,edyr+1):
        for month in range(1,13):
            year_list.append(year)
            month_list.append(month)
    
    spi_df['Year'] = year_list   
    spi_df['Month'] = month_list      
      
    return spi_df
    

# %%
if __name__ == '__main__':
    
# 路径处理和变量预定义
    rootdir = r'E:\05_Study\05_DroughtIndex\03_SPI'
    inpath = rootdir + '\\02_Data\\'
    
    outpath = rootdir + '\\04_Result\\'
    type_str = '.csv'
    
# 文件名称获取
    file_list,file_list_path  = get_file(inpath,type_str)
    
# 逐站点计算SPEI
    ii = 0    
    for file_path_temp in file_list_path:
        sta_id,styr,edyr,pre = read_oridata(file_path_temp)
        
        # %% climate_indices库遇到降水为0时，对数据的处理会出现问题，使结果为NaN，为防止0值影响，将降水0值设置为0.001
        pre[pre==0] = 0.001 
        
        # %% 不同时间尺度SPI的计算
        spi_dif_scale = dif_scale_spi(sta_id,pre,styr,edyr,(1,3,6,9,12))
        
        # %% 结果的写出
        outfile = str(int(sta_id)) + '_SPI.xlsx'
        outfile_path = outpath + outfile
        spi_dif_scale.to_excel(outfile_path,index=False,sheet_name=str(sta_id))
        
        print(str(sta_id)+' has been finished.')
 
    print('Finished.')