栅格尺度的SPEI采用python，主要是参照
https://climate-indices.readthedocs.io/en/latest/
该网站详细介绍了计算SPEI以及其他气候指数的过程，不懂的同学就先下载例子进行试验。
第一步配置环境变量
该python编译器是采用的Anaconda3，如果没有安装的童靴先安装，这里安装以及环境变量的配置步骤就省略了。
接下来就是SPEI运行环境的配置了
conda create -n indices_env python=3.6
conda activate indices_env
pip install climate-indices
conda install -c conda-forge nco
image.png
在安装的Anaconda3的这里运行代码，每一次新打开运行之前都需要执行 conda activate indices_env
第二步 数据处理
该程序计算的原始数据格式为nc5，并且数据（以降水为例）为维度（第一维），经度（第二维），时间（第三维），如果数据不一样的话要进行处理。包括单位都要与例子的数据一致
接下来以CMIP5数据为例进行数据处理,以下是处理了多个nc文件，一般情况下只有一个nc文件，我这里的数据是有多个
import os
from netCDF4 import Dataset
import netCDF4 as nc
import numpy as np
import pandas as pd
import datetime as dt
import calendar
input = r"E:\Paper\paper5\01-data\SPEI--RCP\OriginalData\Cal_SPEI_Th\pr\historical"
output = r"E:\Paper\paper5\01-data\SPEI--RCP\OriginalDataDownScallingChangeUnit"
def ReadData(filePath):
with nc.Dataset(filePath) as file:
file.set_auto_mask(False)
variables = {x: file[x][()] for x in file.variables}
return variables
def WriteData(inputVariables,outputData,outputFilePath,variableName,variableUnits):
newDataFile = nc.Dataset(outputFilePath, 'w', format='NETCDF4')
#define dimensions
long=newDataFile.createDimension('lon',size=length_lon)
lati=newDataFile.createDimension('lat',size=length_lat)
times=newDataFile.createDimension('time',size=length_time)
#define variables
lon=newDataFile.createVariable('lon','f8',dimensions='lon')
lat=newDataFile.createVariable('lat','f8',dimensions='lat')
time=newDataFile.createVariable('time','S10',dimensions='time')
var=newDataFile.createVariable(variableName,'f8',dimensions=('lat','lon','time'))
#add data to variables
lon[:]=inputVariables['lon']
lat[:]=inputVariables['lat']
#time=inputVariables['time']
var[...]=outputData
timeRange=pd.date_range(dt.datetime(1850,1,1,0),dt.datetime(2005,12,1,0),freq='MS')
for i in range(timeRange.shape[0]):
time[i]=timeRange[i].strftime('%Y-%m-%d %H:%M')
year = int(time[i][0:4])
month = int(time[i][5:7])
day_temp = calendar.monthrange(year,month)
day = day_temp[1]
#add attributes
#global attributes
newDataFile.times=time.shape[0]
newDataFile.start_time=time[0]
newDataFile.end_time = time[-1]
##variables attributes
###lon
lon.units = "degrees_east"
lon.long_name = "longitude"
lon.standard_name = "longitude"
lon.axis = "X"
# lon.valid_min = 0.25
# lon.valid_max = 359.75
###lat
lat.units = "degrees_north"
lat.long_name = "latitude"
lat.standard_name = "latitude"
lat.axis = "Y"
# lat.valid_min = -89.75
# lat.valid_max = 89.75
###lwe
var.units = variableUnits
var.grid_mapping = "WGS84"
var.coordinates = "lat lon time"
##close file
newDataFile.close()
def getDay():
day=np.zeros(length_time)
timeRange = pd.date_range(dt.datetime(1850, 1, 1, 0), dt.datetime(2005, 12, 1, 0), freq='MS')
for i in range(timeRange.shape[0]):
year = int(timeRange[i].strftime('%Y-%m-%d %H:%M')[0:4])
month = int(timeRange[i].strftime('%Y-%m-%d %H:%M')[5:7])
day_temp = calendar.monthrange(year,month)
day[i] = day_temp[1]
return day
def datachange(filepath):
variables = ReadData(filepath)
var_data=variables['pr']
a=np.swapaxes(var_data,2,0)
b=np.swapaxes(a,0,1)
day=getDay()
b = b * day * 24 * 36 * 100 # 这里代表修改了原始数据的降水单位
WriteData(variables, b, fileOutPath, 'pr', 'millimeter')
filenames = os.listdir(input)
for i in range (len(filenames)):
filepath = input + "\\" + filenames[i]
fileOutPath = output + "\\" + filenames[i]
data = Dataset(filepath)
# all_vars = data.variables.keys() # 获取所有变量名称
# all_vars_info = data.variables.items() # 查看每一个变量的信息
var2 = 'lat'
var_info2 = data.variables[var2]
length_lat = len(list(var_info2))
# print(var_info2)
var3 = 'lon'
var_info3 = data.variables[var3]
length_lon = len(list(var_info3))
# print(var_info3)
var4 = 'time'
var_info4 = data.variables[var4]
length_time = len(list(var_info4))
datachange(filepath)
计算PET
上述数据处理好之后，就可以计算计算PET了
PET通过这里的程序进行计算，只提供了两种计算方式，Thornthwaite和Hargreaves
以Thornthwaite为例
process_climate_indices --index pet --periodicity monthly --netcdf_temp E:/Paper/paper5/01-data/SPEI--RCP/OriginalDataDownScallingChangeUnit/tas/tas_Amon_CanESM2_rcp85_r1i1p1_200601-210012.nc --var_name_temp tas --output_file_base E:/Paper/paper5/01-data/SPEI--RCP/output_PET/CanESM2_rcp85 --multiprocessing all_but_one
计算SPEI
process_climate_indices --index spei --periodicity monthly --netcdf_precip E:/Paper/paper5/01-data/SPEI--RCP/OriginalDataDownScallingChangeUnit/pr/pr_Amon_CanESM2_rcp85_r1i1p1_200601-210012.nc --var_name_precip pre --netcdf_pet E:/Paper/paper5/01-data/SPEI--RCP/output_PET/CanESM2_rcp85_pet_thornthwaite.nc --var_name_pet pet_thornthwaite --output_file_base E:/Paper/paper5/01-data/SPEI--RCP/output_SPEI/CanESM2_rcp85 --scales 1 3 6 12 --calibration_start_year 2006 --calibration_end_year 2100 --multiprocessing all
image.png
上面的计算步骤需要修改路径和变量名称，见图片上说明
批量处理
上面的计算过程只是针对单个的文件，但是如果有多个nc文件需要计算SPEI，就可以采用以下程序，将以下程序复制保存成bat文件，然后将bat文件拖进Anaconda Prompt里运行即可
echo off
setlocal enabledelayedexpansion
for %%i in (E:\Paper\paper5\01-data\SPEI--RCP\OriginalDataDownScallingChangeUnit\tas\*.nc) do (
set file=%%~ni
echo %%~ni
echo !file:~3!
set outputfile=E:\Paper\paper5\01-data\SPEI--RCP\output_PET1\pet!file:~3!
echo !outputfile!
process_climate_indices --index pet --periodicity monthly --netcdf_temp %%i --var_name_temp tas --output_file_base !outputfile! --multiprocessing all_but_one
echo off
setlocal enabledelayedexpansion
for %%i in (E:\Paper\paper5\01-data\SPEI--RCP\OriginalDataDownScallingChangeUnit\pr\*.nc) do (
set prfile=%%i
set prfilename=%%~ni
set petfile=E:\Paper\paper5\01-data\SPEI--RCP\output_PET1\pet!prfilename:~2!_pet_thornthwaite.nc
set outputfile=E:\Paper\paper5\01-data\SPEI--RCP\output_SPEI1\spei!prfilename:~2!
echo !prfile!
echo !petfile!
echo !outputfile!
process_climate_indices --index spei --periodicity monthly --netcdf_precip !prfile! --var_name_precip pre --netcdf_pet !petfile! --var_name_pet pet_thornthwaite --output_file_base !outputfile! --scales 1 3 6 12 --calibration_start_year 2006 --calibration_end_year 2100 --multiprocessing all
)
