gma 教程 | 气候气象 | 计算标准化降水蒸散指数（SPEI）

Posted 2023-04-04 洛的地理研学

目标

【基于 Excel 降水和蒸散数据计算 SPEI】
【基于 GTiff 栅格降水和蒸散数据计算 SPEI】

环境

系统： Window 10+ (X64)
Python 版本： 3.8.8 +
gma 版本： 1.0.10 +

gma 安装和详细功能帮助见：地理与气象分析库

函数

gma.climet.SPEI(PRE, PET, Axis = None, Scale = 1, Periodicity = 12)

---

功能：【标准化降水蒸散指数】。基于 Log-Logistic 分布计算标准化降水蒸散指数。

参数：

 PRE: array。降水量（mm）。

 PET: array。潜在蒸散量（mm）。

可选参数：

  Axis = int。计算轴。如果不设置（None），多维数据会将所有数据展开到一维计算。

  Scale = int。时间尺度。默认为 1。例如：1月、3月或其他。

  Periodicity = int。周期。默认为 12。例如：月数据可以以 12 为周期。

*注意：Scale、Periodicity 基于计算轴！

返回：array。

参考文献：

 Vicente-Serrano S M , S Beguería, JI López-Moreno. A Multiscalar Drought Index Sensitive to Global Warming: The Standardized Precipitation Evapotranspiration Index[J]. Journal of Climate, 2010, 23(7):1696-1718.

---

示例

数据

• 洛阳市辖区内某点 1981-2020 年月降水和蒸散数据（共 480 个月）
• 洛阳市辖区1981-2020 年月降水和蒸散空间栅格数据（GTiff 格式）（降水和蒸散各一个文件，每个文件含 1981-2020 对应月数据，共 480 个波段），由全国数据裁剪合并获取。全国范围蒸散（ET0）（基于彭曼 Penman-Monteith 法）获取地址：https://download.csdn.net/download/weixin_42155937/85575273

基于 Excel 表数据的 SPEI 计算

import gma
import pandas as pd
# 读取数据
Data = pd.read_excel('PRE_ET0.xlsx')

# 看一下 Excel 文件前 5 行内容
Data.head()

PRE	ET0
15.2	53.4
7.2	66.9
34.5	114.4
22.8	144.1
3.7	208.5

# 提取 PRE 和 ET0 用于 SPEI 运算
PRE = Data['PRE'].values
ET0 = Data['ET0'].values
# 分别计算1个月、3个月、6个月、12个月、24个月、60个月尺度的 SPEI 数据
SPEI1 = gma.climet.SPEI(PRE, ET0)
SPEI3 = gma.climet.SPEI(PRE, ET0, Scale = 3)
SPEI6 = gma.climet.SPEI(PRE, ET0, Scale = 6)
SPEI12 = gma.climet.SPEI(PRE, ET0, Scale = 12)
SPEI24 = gma.climet.SPEI(PRE, ET0, Scale = 24)
SPEI60 = gma.climet.SPEI(PRE, ET0, Scale = 60)

绘制计算结果

import matplotlib.pyplot as plt
PAR = 'font.sans-serif': 'Times New Roman',
       'axes.unicode_minus': False,
      
plt.rcParams.update(PAR)

## 标记一下不同时间尺度结果变量的不同
S = [1,3,6,12,24,60]
## 准备横坐标（年份）的标签
X = range(len(PRE))
Date = gma.osf.DateSeries('198101','202101',DateDelta='M', Format='%Y%m').strftime('%Y-%m')

## 循环绘制 6 个尺度的 SPEI 结果
plt.figure(figsize = (18, 14), dpi = 300)
for i in range(6):
    ax = plt.subplot(4, 2, i + 1) 
    ### 绘制数据
    ax.plot(X, eval(f'SPEIS[i]'), linewidth = 0.8, c = 'gray')
    ### 添加图例
    ax.legend([f'SPEIS[i]'],frameon = False)
    ### 添加横坐标标签
    ax.set_xticks(X[::72], Date[::72], rotation = 0)
    ### 定义横纵坐标显示范围
    ax.set_xlim(-12)
    ax.set_ylim(-3.8, 3.8)
    ### 绘制干（-1）湿（1）分界线
    plt.axhline(y = -1,ls=(0,(6,6)), c="r", linewidth = 0.4)
    plt.axhline(y = 1,ls=(0,(6,6)), c="b", linewidth = 0.4)
    ### 绘制其他网格
    ax.grid(True, linestyle = (0,(6,6)), linewidth = 0.3)
### 修改子图边距
plt.subplots_adjust(wspace = 0.04, hspace = 0.18)
plt.show()

绘图结果如下：

将结果保存到文件

OUT = pd.DataFrame([SPEI1, SPEI3, SPEI6, SPEI12, SPEI24, SPEI60],
                   index = ['SPEI1','SPEI3','SPEI6','SPEI12','SPEI24','SPEI60']).T
OUT.to_excel('SPEI.xlsx', index = False)

基于 GTiff 栅格数据的 SPEI 计算

import gma
## 读取数据集
PERSet = gma.Open(r'D:\\BaiduNetdiskDownload\\PRE_Luoyang_1981-2020.tif')
ET0Set = gma.Open(r'D:\\BaiduNetdiskDownload\\ET0_Luoyang_1981-2020.tif')
PRE = PERSet.ToArray()
ET0 = ET0Set.ToArray()
### 读取的数据为三维数据（波段，行，列），第一维为时间序列（月数据）。因此按照轴 0 来计算
SPEI1 = gma.climet.SPEI(PRE, ET0, Axis = 0)
SPEI3 = gma.climet.SPEI(PRE, ET0, Axis = 0, Scale = 3)
SPEI6 = gma.climet.SPEI(PRE, ET0, Axis = 0, Scale = 6)
SPEI12 = gma.climet.SPEI(PRE, ET0, Axis = 0, Scale = 12)
SPEI24 = gma.climet.SPEI(PRE, ET0, Axis = 0, Scale = 24)
SPEI60 = gma.climet.SPEI(PRE, ET0, Axis = 0, Scale = 60)

绘制最后一个月（2020年12月）计算结果

绘图方法请参考：基于 Python 的地理空间绘图指南

存储计算结果

for i in S:
	# 保存所有结果中的非全 nan 波段。即：去除时间尺度累积时序列前无效的波段。
    gma.rasp.WriteRaster(fr'.\\1981-2020_SPEIi.tif', 
                         eval(f'SPEIi[i-1:]'), 
                         Projection = PRESet.Projection,
                         Transform = PRESet.GeoTransform, 
                         DataType = 'Float32', 
                         NoData = np.nan)  

反馈与讨论

微信号：Luo_Suppe

相关链接

地理与气象分析库

全国多年平均降水量空间分布数据1951-2021年气温分布数据蒸散量数据蒸发量数据太阳辐射数据日照数据风速数据

引言

地理遥感生态网结合1951-2021年各地区地面气象监测站数据，应用气候数据空间插值软件Anusplin预测全国年降雨量分布数据成果。得出全国各个省市自治区降雨分布图，全国各省市自治区降雨量数据产品是地理遥感生态网推出的气象气候类数据产品之一，产品数据精度高，质量良好。

 数据下载链接：数据下载链接

正文

栅格数据源

数据名称	年降雨量空间分布数据
数据类型	栅格
数据格式	TIFF
时间序列	1951年-2021年
空间分辨率	1km/500m/30m
投影坐标系	Krasovsky_1940_Albers
覆盖范围	全域

数据简介

    年平均降雨量，是指某地多年降雨量总和除以年数得到的均值，或某地多个观测点测得的年降雨量均值。年平均降雨量是一地气候的重要衡量指标之一。

    降雨数据对预测和防范洪涝灾害、对水文水资源系统的预报预测、农业作物生产及生态环境变化都有着极其重要的意义，为农业水利等相关部门工作及制定相应的政策措施提供决策支持。

    地理遥感生态网结合2000-2010年各地区地面气象监测站数据，应用气候数据空间插值软件Anusplin预测全国年降雨量分布数据成果。得出全国各个省市自治区降雨分布图，全国各省市自治区降雨量数据产品是地理国情监测云平台推出的气象气候类数据产品之一，产品数据精度高，质量良好。

    辽宁省降雨量分布特点为从沿海向内陆递减，主要表现为东西差异，降雨量在480-1452mm之间。东部以抚顺、本溪和丹东为降雨量多的地区，均在1150mm以上，降雨量较低的地区为朝阳市和阜新市，降雨量小于750mm。

辽宁省2010年降雨量分布

数据来源引用：地理遥感生态网www.gisrs.cn

同时，地理遥感生态网www.gisrs.cn上分享了很多地理遥感领域的科学数据（土地利用数据、npp净初级生产力数据数据、NDVI数据、气象数据（降雨量、气温、蒸散量、辐射、湿度、日照时数、风速、水汽压数据）、径流量数据、夜间灯光数据、统计年鉴、道路网、POI兴趣点数据、GDP分布、人口密度分布、三级流域矢量边界、地质灾害分布数据、土壤类型、土壤质地、土壤有机质、土壤PH值、土壤质地、土壤侵蚀、植被类型、自然保护区分布、建筑轮廓分布等等地理数据，以及关于gis、遥感从方面的操作教程）。