基于Python的Climate Indices库计算SPEI(标准化降水蒸散发指数)03—单站点不同时间尺度的SPEI计算

忘记过去就等于背叛自己。

---

前言

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

1. 概述

2.1 目的

• 针对某一站点，调用Climate Indices库进行不同时间尺度的SPEI的计算

2.2 说明

• 此示例计算SPEI过程中，运用的是桑斯维特算法计算潜在蒸散发。
• Climate Indices库计算SPEI的基本说明请参考基于Python的Climate Indices库计算SPEI(标准化降水蒸散发指数)02—站点SPEI计算

2. 版本

2.1 山东青岛，2021年9月23日，Version1

3. 微信公众号GISRSGeography

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

---

一、数据

1. 输入数据

• 基于桑斯维特算法计算SPEI时需要纬度、平均温和降水三要素，其中纬度和温度用于计算蒸散发，同时Climate Indices的库计算SPEI时需要数据的开始年和结束年用于校正，因此本示例的数据（青岛站）组织形式如下：
  
• 数据可以通过公众号GISRSGeography回复“SPEI青岛”获取
• spei.exe软件可以通过公众号GISRSGeography回复“spei.exe”获取

2. 输出数据

• 本程序的输出数据是以.csv格式存储的不同时间尺度SPEI的计算结果，输出结果组织形式如下：
  

3. 设计思路

• 此程序的设计思路核心是修改climate_indices.spei()方法中的scale（即时间尺度） 参数的值，计算获得不同时间尺度的spei并写出。
  

二、程序示例

1. 程序代码

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

2. 山东青岛，2021年9月23日

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

from climate_indices import indices
from climate_indices import compute  # 计算SPEI的包

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



# %% 不同时间尺度spei计算函数定义
def dif_scale_spei(
        sta_id: int,
        lat: float,
        tas_data: np.array,
        pre_data: np.array,
        styr: int,
        edyr: int,
        scale: tuple
        ) -> pd.DataFrame:
    
    """
    (1) 功能：
        1) 计算单一站点不同时间尺度的spei
    ---
    
    (2) 输入参数
        1) sta_id: int, 站点号
        2) lat: float, 站点纬度
        3) tas_data: np.array, 平均温序列，月值
        4) pre_data: np.array, 降水序列，月值
        5) styr: int, 开始年
        6) edyr: int, 结束年
        7) scale: tuple, spei的时间尺度
    ---
    
    (3) 返回值
        1) spei_df: pd.DataFrame, 不同时间尺度的SPEI
           列名：站点号，年，月，不同时间尺度的SPEI
           
        注：SPEI计算结果中的-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
      # 空的数据框的创建
    spei_df = pd.DataFrame(data=[],columns=col_names)
    
    # 潜在蒸散发计算-桑斯维特方法
    pet_data = indices.pet(temperature_celsius=tas_data,
                      latitude_degrees=lat,
                      data_start_year=styr)
    
    # 不同时间尺度的spei的计算
    for scale_temp in scale:
        spei = indices.spei(precips_mm=pre_data,
                    pet_mm=pet_data,
                    scale=scale_temp, # 时间尺度
                    distribution=indices.Distribution.gamma,
                    periodicity=compute.Periodicity.monthly,
                    data_start_year=styr,
                    calibration_year_initial=styr,
                    calibration_year_final=edyr,
                    )
        
        spei[np.isnan(spei)] = -99 # nan转换为-99
        if scale_temp < 10:
            spei_df['Scale_0'+str(scale_temp)]=spei
        else:
            spei_df['Scale_'+str(scale_temp)]=spei
    
    # 站点信息补充
    spei_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)
    
    spei_df['Year'] = year_list   
    spei_df['Month'] = month_list      
      
    return spei_df
    

# %%
if __name__ == '__main__':
    
# 路径处理和变量预定义
    inpath = r'D\01_Data\\'
    infile = inpath + '54857ClimMonthData.csv'
    outpath = r'D\03_Result\\'
    
# 数据读取    
    sta_id,styr,edyr,lat,tas_mean,pre = read_oridata(infile)
    
# 不同时间尺度的SPEI的计算
    spei_dif_scale = dif_scale_spei(sta_id,lat,tas_mean,pre,styr,edyr,scale=(1,3,6,12,9))
    
# 计算结果的写出
    spei_dif_scale.to_csv(outpath+'spei_dif_scale.csv',index=False)    
    
    print('Finished.')

2. 输出结果

• 输出结果为特定站点不同时间尺度的SPEI