文章目录
前言
??此系列博文的目的是基于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青岛”获取
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
