SPEI干旱指数计算与可视化

发表于 分类于 阅读次数： Disqus：

标准化降水和蒸散量指数(SPEI)在农业、气象、水文乃至保险等领域有着广泛应用，本文简单介绍了几种主流的干旱指数，并以SPEI为例介绍了数据的获取、处理和展示过程。

干旱指数

干旱作为一种极端气候事件和自然灾害，一直以来是造成农业损失和环境破坏的重要原因，对经济、农业、水资源、旅游乃至生态系统等不同领域造成巨大影响。

“冰冻三尺非一日之寒”，与一般的自然灾害不同，干旱往往是在长时间的降水短缺之后，才会展现出显著的影响和巨大的破坏，因此很难确定一次干旱事件的起始、结束和严重程度。此外，干旱从范围上呈区域性特征，每个地区都有特定的气候特征。发生在北美大平原的干旱与那些发生在巴西东北部、南部非洲、西欧、澳大利亚东部、或中国华北平原的干旱是不同的。上述地区中每个地区之间在降水量、降水季节性和降水形式上的差异很大。

干旱对不同的用户有不同的含义，这些用户可能是水管理者、农业生产者、水力发电厂经营者、以及野生动物生物学工作者。即使在行业内部对干旱也有诸多不同的观点，因为影响可能明显不同。通常以气象、农业和水文为类别将干旱进行分类，但各类之间的强度、持续时间和空间分布有所不同。

所以，如何客观、准确地量化干旱事件的强度、幅度、持续时间和空间范围等特征，已经成为全球关注的难题，世界各国都在此领域投入了大量精力。

其中，干旱指数已经逐步成为衡量干旱事件严重程度的重要工具和普遍方法。但是，干旱定义中的主观性导致难以建立一个具有普适性的干旱指数，各类干旱指数不断推陈出新。目前，大多数干旱监测系统和分析研究主要基于以下两种指数衍变：

1.帕尔默干旱强度指数

帕尔默干旱强度指数（the Palmer Drought Severity Index, PDSI)，也称帕尔默干旱指数（Palmer Drought Index, PDI），是一个基于土壤水平衡方程的指数，作为作为利用不局限于降水资料确定干旱的初步尝试之一，PDSI在干旱指数的发展史上具有里程碑的意义，自上世纪60年代提出以来在全球各地广泛使用，美国国家海洋和大气管理局每周都会发布Palmer地图。在很长一段时间内，PDSI是唯一的业务干旱指数，后来许多干旱指数都是针对PDSI的不足进行改进而诞生的。Palmer指数在判定长期干旱——比如长达数个月的干旱——方面被证明是非常有效的，但在判断持续数周的干旱上稍差。此外，PDSI要求连续完整的资料并存在季节性和滞后性的缺点。

2.标准化降水指数

标准化降水指数（the Standardised Precipitation Index, SPI）是一种基于降水概率的方法，它可以将观测获得的任意时间尺度的降水概率转化为一个指数。20世纪90年代，McKee等人将土壤湿度、地下水、积雪、河流径流、水库蓄水等在内的可用水资源纳入考虑，清楚地阐述了干旱多标量这一基本特征。

由于从水资源输入到真正可用的时间差别很大，因此干旱持续的时间尺度至关重要，并且需要将水文上的、环境上的、农业上的以及其他意义上的干旱加以区分。所以，干旱指数必须与特定的时间尺度相关联，以便于监测和管理不同的可用水资源。由于SPI可以利用任何地点的历史降水记录得出可在任何数量时间尺度（1个月至48个月或更长）计算的降水概率，并且在不同时间、空间之间具有可比性，SPI被广泛用于不同水资源干旱的监控。
在不同时间尺度上计算SPI的能力有利于多种用途。根据上述的干旱影响，3个月以下的SPI值或许可用于基本干旱监测，而6个月以下的SPI值可用于监测农业影响，12个月以上的SPI值适用于水文影响。

2009年，WMO推荐SPI作为主要的气象干旱指数，各国可将其用于监测和跟踪旱情。

标准化降水和蒸散量指数

基于降水的干旱指数（包括SPI在内）依赖于两个基本假设：
1. 降水的变化远远高于其他变量，如温度和潜在蒸散发（potential evapotranspiration, PET）；
2. 其他变量在时间上是平稳的（stationary），即没有时间上的变化趋势。
在这两个前提下，其他变量的重要性可以忽略不计，但是，研究者警告不要系统性地忽略温度对干旱程度影响的重要性。实证研究表明，气温上升对干旱的严重程度具有显著影响，最近的研究分析了变暖引起的干旱胁迫的作用，该研究分析了干旱对净初级生产和树木死亡率的影响（Williams et al。，2011;Martínez-Villalta） et al。，2008; McGuire et al。，2010; Linares and Camarero，2011）。在这种情况下，标准化降水和蒸散量指数从众多干旱指数中脱颖而出。

标准化降水和蒸散量指数（the Standardized precipitation evapotranspiration index, SPEI）,顾名思义，是对标准化降水止水SPI的扩展，旨在考虑降水和潜在蒸散量（potential evapotranspiration, PET）来确定干旱。 SPEI以SPI为基础，但增加了温度组成部分，使得该指数可通过基本的水平衡计算，来考虑温度对干旱发展的影响。 因此，与SPI不同，SPEI说明了较高温度导致的干旱持续时间和幅度的增加。SPEI有强度范围，按照这个范围来计算正负值，以确定干湿事件。它可以计算短至1个月长至48个月或更长的时间尺度。SPEI的输入包括月降水和温度资料，输出适用于所有气候体系，其结果是类似的，因为它们都是标准化的。由于使用了温度资料，在审视不同未来情景下模式输出中的气候变化影响时，SPEI是一个较为理想的指数。

SPEI相对于考虑PET对干旱严重程度的影响的其他广泛使用的干旱指数的一个关键优势是其多标量特征使得能够在气候变化的背景下识别不同的干旱类型和影响。气候变化模型预测，21世纪气温将显着增加。作为后果，预计未来干旱条件将急剧发展，由于蒸发蒸腾导致需水量增加。

在本文中，我们将采用简单的Thornthwaite方法估计潜在蒸散量（PET），然后通过拟合Gamma分布得到SPEI指数。所需要的数据可以在http://www.cru.uea.ac.uk/data/ 上找到。

noaa-Zus94oboIsM-unsplash

准备数据

我们从https://crudata.uea.ac.uk/cru/data/hrg/ 上下载数据，crt_ts的最新版本为4.02，计算SPEI需要降水和温度，分别在pre和tmp两个文件夹下面，数据是按照10年为一个文件组织的，例如“cru_ts4.02.1981.1990.pre.dat.nc.gz”表示1981~1990年的降水数据，文件格式为nc的压缩包。

假设我们下载了1961~2017年6个压缩包，在本地解压缩以后，发现6个文件大约1.32GB，为了方便使用，可以首先将6个文件拼成一个，可以自己编写Python代码，也可以用现成的工具，如NCO和CDO。但是这些数据似乎都主要针对UNIX/LINUX环境，Windows下面用起来并不方便。NCO的使用命令（我用cmd不行，用的bash）为：

1
ncra -Y ncrcat *.tmp.dat.nc cru_ts4.02.1961.2017.tmp.dat.nc

准备好数据，就进入Python环境下进行分析。

1
2
3
4
5
import numpy as np
from timeit import default_timer as timer
import datetime             
from netCDF4 import Dataset
from netcdftime import utime

加载数据

加载降水数据

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
pre_file = r'D:\ClimateData\PRCP\cru_ts4.02.1961.2017.pre.dat.nc'
pre_fh = Dataset(pre_file, mode='r') # file handle, open in read only mode
pre_fh.set_auto_mask(False)
lons = pre_fh.variables['lon'][:]
lats = pre_fh.variables['lat'][:]
nctime = pre_fh.variables['time'][:]
t_unit = pre_fh.variables['time'].units
prcp = pre_fh.variables['pre'][:]
#  undef = 9.96921E36
undef = pre_fh.variables['pre'].missing_value

try :
    t_cal = pre_fh.variables['time'].calendar
except AttributeError : # Attribute doesn't exist
    t_cal = u"gregorian" # or standard

pre_fh.close() # close the file

mask = (prcp==undef)
prcp = np.ma.MaskedArray(prcp, mask=mask)
prcp = prcp.astype(np.float64)

nt,nlat,nlon = prcp.shape

加载温度数据

1
2
3
4
5
6
7
8
tmp_file = r'D:\ClimateData\TEMP\cru_ts4.02.1961.2017.tmp.dat.nc'
tmp_fh = Dataset(tmp_file, mode='r') # file handle, open in read only mode
tmp_fh.set_auto_mask(False)
tmp = tmp_fh.variables['tmp'][:]
tmp_fh.close()

tmp = np.ma.MaskedArray(tmp, mask=mask)
tmp = tmp.astype(np.float64)
1
2
3
4
data_utime   = utime(t_unit, calendar=t_cal)
datevar = data_utime.num2date(nctime)

datevar[0:5]
array([datetime.datetime(1961, 1, 16, 0, 0),
       datetime.datetime(1961, 2, 15, 0, 0),
       datetime.datetime(1961, 3, 16, 0, 0),
       datetime.datetime(1961, 4, 16, 0, 0),
       datetime.datetime(1961, 5, 16, 0, 0)], dtype=object)

SPEI计算与呈现

计算过程

我们通过调用在NOAA干旱指数基础上开放的climate_indices(https://github.com/monocongo/climate_indices)库计算SPEI，SPEI需要用到PET，这个也可以通过调用climate_indices完成。

虽然这些指数的计算都并非十分复杂，但是如果针对全球的网格进行计算还是需要耗费一定的时间。由于这个for循环结构十分简单，完全可以并行化，遗憾的是在python中没有找到类似MATLAB的parfor简单的做法，搞了很久也没有弄明白（多进程遇到了数据共享的问题），最后只好作罢。（climate_indices实现了多进程的脚本，但是并不想用）

1
2
from climate_indices import indices
from climate_indices import compute
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
spei = np.zeros(tmp.shape)
spei[:,:,:] = np.nan
pet = np.zeros(tmp.shape)
pet[:,:,:] = np.nan
nrun = 24

s = timer()
for ilat in np.arange(nlat):
    lat = lats[ilat]
    for ilon in np.arange(nlon):
        one_pr = prcp[:,ilat,ilon]
        one_tmp = tmp[:,ilat,ilon]

        if(not np.ma.is_masked(one_pr)):
            one_pet = indices.pet(temperature_celsius=one_tmp,
                                 latitude_degrees=lat,
                                 data_start_year=1961)
            pet[:,ilat,ilon] = one_pet
            spei[:,ilat,ilon] = indices.spei(precips_mm=one_pr,
                                             pet_mm=one_pet,
                                             scale=nrun,
                                             distribution=indices.Distribution.gamma,
                                             periodicity=compute.Periodicity.monthly,
                                             data_start_year=1961,
                                             calibration_year_initial=1961,
                                             calibration_year_final=2016,
                                            )
pet = np.ma.MaskedArray(pet, mask=mask)
spei = np.ma.MaskedArray(spei, mask=mask)
e = timer()
print(e - s)      
1288.7284692

我用的是1961~2016年期间共56年的数据来计算全球SPEI（网格大小为 0.5∘×0.5∘0.5^\circ\times0.5^\circ ，从上面的计时可以发现，大概用了20多分钟。此外，计算产生的数据也比较大。为了便于下次使用，可以将计算的结果存入到文件当中。

1
2
pet.dump('pet')
spei.dump('spei')

全球SPEI数据可视化

1
spei = np.load('spei')
1
2
3
4
5
6
import matplotlib.path as mpath
import matplotlib.pyplot as plt
from matplotlib.offsetbox import AnchoredText

import cartopy.crs as ccrs
import cartopy.feature as cfeature

利用cartopy可以方便地实现地理数据的可视化，具体如下所示。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
fig = plt.figure(figsize=(30,24))

ax = fig.add_subplot(2, 1, 1, projection=ccrs.PlateCarree())
ax.set_extent([min(lons), max(lons), min(lats), max(lats)], crs=ccrs.PlateCarree())
# Put a background image on for nice sea rendering.
# ax.stock_img()


SOURCE = 'Natural Earth'
LICENSE = 'public domain'
nrun=24
ax.add_feature(cfeature.LAND)
ax.add_feature(cfeature.COASTLINE)

# Plot the first one
# ax = fig.add_subplot(211)
idx = 498
x, y = np.meshgrid(lons, lats)
# cs = m.contourf(x, y,, clevs, cmap=plt.cm.RdBu)
cs = ax.contourf(x, y, spei[idx,:,:],
            transform=ccrs.PlateCarree(),
            cmap=plt.cm.RdBu)
cb = fig.colorbar(cs)
plt.title('SPEI-' + str(nrun) + ' In '+ datetime.date.strftime(datevar[idx], "%m/%Y"), fontsize=16)

# Add a text annotation for the license information to the
# the bottom right corner.
text = AnchoredText(r'$\mathcircled{{c}}$ {}; license: {}'
                    ''.format(SOURCE, LICENSE),
                    loc=4, prop={'size': 12}, frameon=True)
ax.add_artist(text)


#  plot the second one
ax = fig.add_subplot(2, 1, 2, projection=ccrs.PlateCarree())
ax.set_extent([min(lons), max(lons), min(lats), max(lats)], crs=ccrs.PlateCarree())
# Put a background image on for nice sea rendering.
# ax.stock_img()


SOURCE = 'Natural Earth'
LICENSE = 'public domain'
nrun=24
ax.add_feature(cfeature.LAND)
ax.add_feature(cfeature.COASTLINE)

idx = 600
x, y = np.meshgrid(lons, lats)
# cs = m.contourf(x, y,, clevs, cmap=plt.cm.RdBu)
cs = ax.contourf(x, y, spei[idx,:,:],
            transform=ccrs.PlateCarree(),
            cmap=plt.cm.RdBu)
cb = fig.colorbar(cs)

plt.title('SPEI-' + str(nrun) + ' In '+ datetime.date.strftime(datevar[idx], "%m/%Y"), fontsize=16)

# Add a text annotation for the license information to the
# the bottom right corner.
text = AnchoredText(r'$\mathcircled{{c}}$ {}; license: {}'
                    ''.format(SOURCE, LICENSE),
                    loc=4, prop={'size': 12}, frameon=True)

plt.show()

output_16_0

参考文献

  1. Palmer, W.C., 1965: Meteorological Drought. Research Paper No. 45, US Weather Bureau, Washington, DC
  2. http://spei.csic.es/home.html
  3. https://library.wmo.int/pmb_ged/wmo_1090_zh.pdf
  4. http://www.droughtmanagement.info/literature/WMO-GWP_Drought-Indices_ch_2016.pdf
  5. Guttman, N.B., 1999: Accepting the Standardized Precipitation Index: a calculation algorithm. Journal of the American Water Resources Association, 35:311–322,
    doi:10.1111/j.1752-1688.1999.tb03592.x
  6. Vicente-Serrano, S.M., S. Begueria and J.I. Lopez-Moreno, 2010: A multi-scalar drought index sensitive to global warming: the Standardized Precipitation Evapotranspiration Index. Journal of Climate, 23:1696–1718
  7. https://climatedataguide.ucar.edu/climate-data/standardized-precipitation-evapotranspiration-index-spei
  8. http://www.cru.uea.ac.uk/data/
  9. https://www.linkedin.com/pulse/standardized-precipitation-evapotranspiration-index-spei-chonghua-yin/?published=t
  10. https://github.com/monocongo/climate_indices