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

---

前言

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

1. 概述

2.1 目的

• 针对栅格温度和降水数据，调用Climate Indices库进行栅格SPEI的计算

2.2 说明

• 此示例计算SPEI过程中，运用的是桑斯维特算法计算潜在蒸散发。
• Climate Indices库计算SPEI的基本说明请参考基于Python的Climate Indices库计算SPEI(标准化降水蒸散发指数)02—站点SPEI计算
• 如果有潜在蒸散发（PET）的数据，则可以直接计算SPEI。
• 测试数据和程序可以点击基于Python的Climate Indices库实现栅格SPEI计算下载。

2. 版本

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

3. 微信公众号GISRSGeography

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

---

一、数据

1. 输入数据

• 本次示例数据是海河流域1971-2099年逐月的栅格温度和降水数据，数据存储为.tif格式。
• 该.tif数据是逐月数据组成的多波段数据，每一个气象要素的波段数是（2099-1971+1）*12=1548
  

2. 输出数据

• 本程序的输出数据是逐月存储的特定时间尺度SPEI的计算结果，本示例计算的是3个月尺度的SPEI，输出结果组织形式如下：
  

![输出数据](https://img-blog.csdnimg.cn/direct/5b8737ab438d4fb19b5f3e4d5d15e134.png

二、程序代码

# -*- coding: utf-8 -*-
"""
1. 程序目的(1) 基于Python的Climate Indices库在栅格尺度上计算海河流域的SPEI2. 版本2.1 2024年1月7日  星期日  Version 13. 数据根路径：'E:\05_Study\05_DroughtIndex\02_SPEI\'3.1 输入数据'01_Data\02_TifData\'路径下的数据3.1 输出数据'03_Result\02_TifSPEI\'路径下的数据4. 参考资料"""# %% 包的导入import numpy as npfrom osgeo import gdal,osrfrom climate_indices import indices
from climate_indices import compute# %% 函数预定义# %% 获取tif数据基本信息
def tif_infor(tif_infor_file: str) -> dict:'''(1) 函数功能-------------获取已知的.tif文件的基本信息，包括数据的BandCount、XSize，YSize，GeoTransform和Projection(2) 参数-------------tif_infor_file : str, tif数据的完整路径(3) 返回值-------------tif_infor_dict : dict, 包括BandCount、XSize，YSize，GeoTransform和Projection的字典'''# 获取数据集tif_ds = gdal.Open(tif_infor_file)# 投影信息tif_proj = tif_ds.GetProjection()# 仿射矩阵信息geo_transform = tif_ds.GetGeoTransform()# 波段数目band_count = tif_ds.RasterCount# 行数roww = tif_ds.RasterYSize# 列数coll = tif_ds.RasterXSize# 存储tif文件的基本信息tif_infor_dict = {'tif_Band':band_count,'tif_Row':roww,'tif_Col':coll,'tif_GeoTrans':geo_transform,'tif_Proj':tif_proj}return tif_infor_dict# %% tif数据读取函数
def tif_read(tif_file: str,) -> np.array:'''(1)函数功能------------读取tif格式的栅格数据(2)参数------------tif_file : str, 存储tif数据的绝对路径(3)返回值------------tiff_data : np.array, 数据读取结果   '''tif_infor_ds = gdal.Open(tif_file)im_width = tif_infor_ds.RasterXSize  # 栅格矩阵的列数im_height = tif_infor_ds.RasterYSize # 栅格矩阵的行数im_bands = tif_infor_ds.RasterCount  # 栅格矩阵的层数if im_bands > 2: # 读取的.tif数据是三维数组tiff_data = np.full((im_bands,im_height,im_width),fill_value=np.NaN,dtype='float32')for band in range(0,im_bands):# 需要注意的是GetRasterBand()的第1个波段标记为1，而不是0，因此idx=band+1tif_infor_band = tif_infor_ds.GetRasterBand(band+1) tif_data_band = tif_infor_band.ReadAsArray() # 将波段数据读取为数组tiff_data[band,:,:] = tif_data_bandelse: # 读取的.tif数据是二维数组tiff_data = np.full((im_height,im_width),fill_value=np.NaN,dtype='float32')tif_infor_band = tif_infor_ds.GetRasterBand(1) # 获取波段tif_data_band = tif_infor_band.ReadAsArray() # 将波段数据读取为数组tiff_data[:,:] = tif_data_bandreturn tiff_data# %% 纬度矩阵创建函数
def lat_matrix_get(tif_infor_dic: dict) -> np.array:'''(1)函数功能------------创建纬度矩阵，辅助PET的计算(2)参数------------tif_infor_dic: dict, 包括tif数据基本信息的字典(3)返回值------------lat_matrix: np.array, 纬度矩阵'''tif_geotrans = tifinfor_dict['tif_GeoTrans']lat_left_left = tif_geotrans[3] # 左上角像元的左上角纬度raster_res = tif_geotrans[1] # 空间分辨率# 左上角像元的中心纬度lat_left_center = lat_left_left - raster_res/2# 创建存储纬度矩阵的空数组roww = tif_infor_dic['tif_Row']coll = tif_infor_dic['tif_Col']lat_matrix = np.full((roww,coll),fill_value=np.NaN,dtype='float32')# lat_matrix填充for ii in np.arange(roww):lat_matrix[ii,:] = lat_left_center - raster_res*iireturn lat_matrix# %% 基于栅格数据的特定时间尺度的spei计算函数定义
def raster_spei_cal(tas_dataa: np.array,pre_dataa: np.array,lat_matrixx: np.array,styr: int,edyr: int,scale_aim: int) -> np.array:"""(1) 功能：1) 实现基于栅格数据的特定时间尺度的SPEI的计算---(2) 输入参数1) tas_dataa: np.array，温度数据，月值，三维数组2) pre_dataa: np.array, 降水数据，月值，三维数组3) lat_matrixx: np.array，纬度矩阵，二维数组3) styr: int, 开始年4) edyr: int, 结束年5) scale_aim: int, spei的时间尺度---(3) 返回值1) spei_array: np.array, 特定时间尺度的SPEI，三维数组注：SPEI计算结果中的-99表示无效值"""# 创建存储计算结果的三维数组# 创建列名spei_aim_data = np.full(pre_data.shape,fill_value=np.NaN,dtype='float32')# 逐栅格计算SPEIroww = pre_dataa.shape[1]coll = pre_dataa.shape[2]for roww_temp in range(0,roww):for coll_temp in range(0,coll):tas_temp = tas_dataa[:,roww_temp,coll_temp]pre_temp = pre_dataa[:,roww_temp,coll_temp]lat_temp = lat_matrixx[roww_temp,coll_temp]# pet计算-桑斯维特方法pet_data = indices.pet(temperature_celsius = tas_temp,latitude_degrees = lat_temp,data_start_year = styr)# spei计算spei_temp = indices.spei(precips_mm=pre_temp,pet_mm = pet_data,scale = scale_aim,distribution = indices.Distribution.gamma, # 选择gamma分布拟合periodicity = compute.Periodicity.monthly,data_start_year = styr,calibration_year_initial =styr,calibration_year_final = edyr)spei_aim_data[:,roww_temp,coll_temp] = spei_tempspei_aim_data[np.isnan(spei_aim_data)] = -99 # nan转换为-99#print(roww_temp)return spei_aim_data# %% 结果写出
def TifWrite(out_file: str,tif_data: np.array,geo_transform: tuple,noDataValue : int) -> None:'''(1) 函数功能：----------将目标数据以.tif格式写出(2) 输入参数----------out_file : str，.tif文件的完整输出路径tif_data : np.array，待写出为.tif文件的数据geo_transform : tuple，地理坐标转换信息noDataValue : tuple，设置为NaN的值(3) 返回值-------无'''shape_size = np.shape(tif_data)gtiff_driver = gdal.GetDriverByName('Gtiff')if len(shape_size) == 3: # 如果是三维数组【存在多个波段】band_num = tif_data.shape[0] # 三维数组层数，即波段数x_size = tif_data.shape[2] # 列数y_size = tif_data.shape[1] # 行数# tiff文件数据集gtiff_frame = gtiff_driver.Create(out_file,x_size,y_size,band_num,gdal.GDT_Float32)# tiff文件数据集中存储不同的波段# 注意：gtiff_frame中波段索引从1开始，tif_data[numpy数组]中索引从0开始for ii_band in range(1,band_num+1):band_data = tif_data[ii_band-1,:,:] # numpy数组的索引从0开始out_band = gtiff_frame.GetRasterBand(ii_band) # gtiff_frame中波段索引从1开始out_band.SetNoDataValue(noDataValue)out_band.WriteArray(band_data)else: # 如果是二维数组【只有一个波段】# 启动geotiff数据驱动band_num = 1x_size = tif_data.shape[1] # 列数y_size = tif_data.shape[0] # 行数# tiff文件数据集gtiff_frame = gtiff_driver.Create(out_file,x_size,y_size,band_num,gdal.GDT_Float32)# 数据写入，此时只有一个波段的数据需要写入gtiff_frame.GetRasterBand(1).SetNoDataValue(noDataValue)gtiff_frame.GetRasterBand(1).WriteArray(tif_data)# 影像显示范围设置gtiff_frame.SetGeoTransform(geo_transform)#gtiff_frame.SetProjection(img_proj)# 地理坐标系统设置srs = osr.SpatialReference()srs.ImportFromEPSG(4326) # WGS1984地理坐标系gtiff_frame.SetProjection(srs.ExportToWkt())# 从内存写入磁盘gtiff_frame.FlushCache()gtiff_frame = None # 关闭指针    # %%
if __name__ == '__main__':# %% 路径处理和基本变量定义rootdir = r'E:\05_Study\05_DroughtIndex\02_SPEI'inpath = rootdir + '\\01_Data\\02_TifData'outpath = rootdir + '\\03_Result\\02_TifSPEI'styr = 1971edyr = 2099stmt = 1edmt = 12aim_scale = 3if aim_scale < 10:scale_aim_str = '0' + str(aim_scale)else:scale_aim_str = str(aim_scale)no_data = -99   # 温度和降水数据完整路径tas_file = inpath + '\\HaiRiver_Hadgem2_Tas_RCP26_1971-2099.tif'pre_file = inpath + '\\HaiRiver_Hadgem2_Pre_RCP26_1971-2099.tif'# %% 获取tif数据基本信息tifinfor_dict = tif_infor(tas_file)# %% 温度和降水数据读取tas_data = tif_read(tas_file)pre_data = tif_read(pre_file)# %% 获取纬度矩阵lat_matrix_aim = lat_matrix_get(tifinfor_dict)# %% SPEI计算spei_raster = raster_spei_cal(tas_data,pre_data,lat_matrix_aim,styr,edyr,aim_scale)# %% 逐月写出for year in range(styr,edyr+1):for month in range(stmt,edmt+1):aim_file_temp = outpath + '\\' + 'HaiHe_matrix_' + 'spei' + scale_aim_str + '_' + str(year*100+month)+'.tif'aim_data_temp = spei_raster[(month-1+(year-styr)*12),:,:]TifWrite(aim_file_temp,aim_data_temp,tifinfor_dict['tif_GeoTrans'],no_data)print(year,month)# %%print('Finished.')