1.引言

?????????標(biāo)準(zhǔn)化降水指數(shù)（SPI）是一個廣泛使用的指數(shù)，用于描述一系列時間尺度上的氣象干旱的特征。但是經(jīng)過研究發(fā)現(xiàn)，目前的處理方法基本都是單點(diǎn)進(jìn)行計(jì)算，缺少多點(diǎn)（大區(qū)域）的批量計(jì)算過程。因此本博客從氣象數(shù)據(jù)下載，處理成NC格式文件以及依靠climate indices庫完成多點(diǎn)的SPI指數(shù)計(jì)算，并可以在ARCGIS中利用反距離權(quán)重生成指數(shù)SPI柵格數(shù)據(jù)。本文有博主整理的完成Python代碼資源以及樣例數(shù)據(jù)。

2.SPI原理概述

????????SPI計(jì)算原理是將某時間尺度(如1、3、6、12個月等)降水量的連續(xù)時間序列(最好是長期記錄，一般最少30年)看作服從某種概率密度函數(shù)分布(如gamma分布)，然后推導(dǎo)出相應(yīng)的累積概率函數(shù)，再通過累積概率函數(shù)轉(zhuǎn)換為標(biāo)準(zhǔn)正態(tài)分布。轉(zhuǎn)換之后，某時間尺度樣本的SPI為：該樣本降水量的累積概率所對應(yīng)的標(biāo)準(zhǔn)正態(tài)分布的x軸的值。

????????例如：以3個月為時間尺度，使用1981-2010年30年的降水?dāng)?shù)據(jù)。計(jì)算2010年1月的SPI值。因?yàn)闀r間尺度是3個月，所以2010年1月的累積降水量被定義為2009年11月-2010年1月期間的總降水量，記作P；使用的時間序列為往年同期的降水量數(shù)據(jù)，即各年11月-1月的降水。首先按照原理，將時間序列數(shù)據(jù)假設(shè)為滿足gamma分布g(x)，然后推導(dǎo)其累積概率函數(shù)H(x)，再轉(zhuǎn)換為標(biāo)準(zhǔn)正態(tài)分布。然后，查找P對應(yīng)的累積概率H(P)，然后查找與H(P)相同累積概率的標(biāo)準(zhǔn)正態(tài)分布所對應(yīng)的x軸的值，即為SPI。示意圖如下，左圖為累積概率H(x)，右圖為轉(zhuǎn)換之后的標(biāo)準(zhǔn)正態(tài)分布。

圖1?從虛線伽馬分布到標(biāo)準(zhǔn)正態(tài)分布的等概率變換的例子

詳細(xì)的數(shù)學(xué)原理請參考博客

3.技術(shù)方法

3.1 氣象數(shù)據(jù)下載

????????首先，我們需要下載30年某一區(qū)域的降水?dāng)?shù)據(jù)。因?yàn)殚L時間序列數(shù)據(jù)一般來說不太好獲取，博主使用的是NASA POWER提供的氣象數(shù)據(jù)，是免費(fèi)的。

網(wǎng)站：NASA POWER | Data Access Viewer

• 入下圖所示，我們選擇東北的區(qū)域的數(shù)據(jù)，時間范圍為：1992-2022。因?yàn)檫x擇的是月度數(shù)據(jù)，網(wǎng)站只更新到1992-2020，后面2021與2022只能下載日度的數(shù)據(jù)進(jìn)行累積后處理。（文尾提供Python代碼與樣例氣象數(shù)據(jù)，1_merge_daily_precipitation.py）

圖2 NASA POWER VIEWER 網(wǎng)站頁面

• 下圖為博主下載的區(qū)域氣象數(shù)據(jù)，格式為csv,打開可以看到。數(shù)據(jù)其實(shí)是按照經(jīng)度緯度各網(wǎng)點(diǎn)排列的，因?yàn)镹ASA提供的是0.5°×0.5°的分辨率數(shù)據(jù)，所以整個吉林省為20×12個點(diǎn)。列JAN-DEC為每個月的降水?dāng)?shù)據(jù)。

• 需要注意的是，因?yàn)?021與2022為日度數(shù)據(jù)所以需要進(jìn)行自己累積計(jì)算，算完后合并到1992-2020中即可。

j2022=r"raw_data\jilinday2022.csv"
    weather_table = pd.read_csv(j2022)
    IS_ALL_Month=False
    if IS_ALL_Month:
        out_dic={"PARAMETER":[],"YEAR":[],"LAT":[],"LON":[],"1":[],"2":[],"3":[],"4":[],"5":[],"6":[],"7":[],
            "8":[],"9":[],"10":[],"11":[],"12":[]}
    else:
        out_dic={"PARAMETER":[],"YEAR":[],"LAT":[],"LON":[],"1":[],"2":[],"3":[],"4":[],"5":[],"6":[],"7":[],
            "8":[]}
    
    weather_table=weather_table.set_index(['LAT', 'LON','MO'])    
    for y in (weather_table.index.get_level_values(0).unique()):#緯度方向
        data_y=weather_table[(weather_table.index.get_level_values(0) == y)]
        for x in (data_y.index.get_level_values(1).unique()):#經(jīng)度方向
            data_yx=data_y[(data_y.index.get_level_values(1) == x)]
            for month in (data_yx.index.get_level_values(2).unique()):
                data_month = data_yx[(data_yx.index.get_level_values(2) == month)]
                out_dic[str(month)]=out_dic[str(month)]+[np.sum(data_month["PRECTOTCORR"].values)]
            #循環(huán)日累加
            out_dic["PARAMETER"]=out_dic["PARAMETER"]+["PRECTOTCORR_SUM"]
            out_dic["YEAR"]=out_dic["YEAR"]+[2022]
            out_dic["LAT"]=out_dic["LAT"]+[y]
            out_dic["LON"]=out_dic["LON"]+[x]
    df_out=pd.DataFrame(out_dic)
    df_out.to_csv(r"process_data\jilin2022.csv",index=False)

3.2 安裝 climate indices庫

? ? ? ?本博客在計(jì)算SPI指數(shù)時候需要安裝climate indices庫。

????????climate indices?是由James Adams利用Python開發(fā)的一個計(jì)算各種氣象指數(shù)，包括SPEI, SPI, PET, PDSI, scPDSI的庫，可以使用pip install 來安裝該庫。

????????在cmd 中啟動python，導(dǎo)入climate_indices庫，如果沒提示錯誤，則表示安裝成功。

????????同時，在python\Scripts 文件夾中會生成這個process_climate_indices.exe,后面批處理主要依靠這個exe。

注意：如果按照報錯，可以使用文尾資源中的lib\climate_indices-py3.8.whl文件按照，這個已經(jīng)讓地理所老師重新編譯，適用于window10-11，python3.8環(huán)境。?

3.3 轉(zhuǎn)換NC氣象格式數(shù)據(jù)

? ? ? ? 因?yàn)橄胍褂胏limate indices，輸入數(shù)據(jù)必須為nc格式的數(shù)據(jù)，因此我們必須要將處理得到的xlsx降水?dāng)?shù)據(jù)生成nc文件（2_write_ncfile.py）。

#---netcdf foramt---#
    f_w = nc.Dataset(outpath,'w',format = 'NETCDF4')
    f_w.createDimension('time',times)   
    f_w.createDimension('lat',y_size)   
    f_w.createDimension('lon',x_size)
    ##創(chuàng)建變量。參數(shù)依次為：‘變量名稱’，‘?dāng)?shù)據(jù)類型’，‘基礎(chǔ)維度信息’
    time=f_w.createVariable('time',"S19",('time')) 
    for i in range(times):
        #time[i] = data_serise[i].strftime('%Y-%m-%d')
        time[i] = data_serise[i].strftime('%Y-%m-%d %H:%M:%S')
    time.units = 'times since {0:s}'.format(time[0])
    time.standard_name = 'Time'
    time.axis = 'T' 
    f_w.createVariable('lat',np.float32,('lat'))  
    f_w.createVariable('lon',np.float32,('lon'))
    #t=np.linspace(0,times-1,times,dtype=int)
    lon=np.linspace(min_x,max_x,x_size,dtype=float)
    lat=np.linspace(min_y,max_y,y_size,dtype=float)
    #寫入變量time的數(shù)據(jù)。維度必須與定義的一致。
    #f_w.variables['time'][:]=data_serise#np.array(list_data)#data_serise
    f_w.variables['lon'][:]=lon
    f_w.variables['lat'][:]=lat
    #xarray_data_r=np.transpose(xarray_data,(3,2,0,1))
    #f_w.createVariable( "prcp", np.float32, ('time','lat','lon'),fill_value=fill_value)
    f_w.createVariable( "prcp", np.float32, ('lat','lon',"time"),fill_value=fill_value)
    #f_w.variables["prcp"][:]=xarray_data_r[0]
    f_w.variables["prcp"][:]=xarray_data[:,:,:,0]
    f_w.variables["prcp"].units="millimeter"
    f_w.close

?3.4 運(yùn)行批處理文件

? ? ? ? 看到這一步，我先恭喜大家，數(shù)據(jù)獲取與處理確實(shí)不易。下一步就可以運(yùn)行我們的批處理文件了（3_cal_SPIindex.py）。下面對其中的一些參數(shù)進(jìn)行解釋：

3.5 轉(zhuǎn)成ArcGIS可以導(dǎo)入的格式

????????運(yùn)行完3.4后，可以生成后綴為nc_spi_gamma_03.nc與nc_spi_pearson_03.nc格式的結(jié)果文件，運(yùn)行4_SPI2xlsx.py文件，可以轉(zhuǎn)成excel文件，導(dǎo)入arcgis結(jié)果如下：

4.結(jié)果展示

4.1 SPI 點(diǎn)位結(jié)果展示

? ? ? ? 如下圖所示，這是通過我們批處理計(jì)算得到的吉林省多點(diǎn)的SPI點(diǎn)狀數(shù)據(jù)：

4.2 生成整個區(qū)域面狀結(jié)果

? ? ? ? 這里我們需要在4.1的基礎(chǔ)上，結(jié)合ArcGIS工具中的Spatial?Analyst tools /Interpolation/IDW插值得到插值柵格影像。

????????其中分辨率那一欄，可以根據(jù)需求自己設(shè)置，需要注意的是這里的分辨率單位是°，不是米。插值結(jié)果如下：

4.3 寫在最后

? ? ? ? 完整資源和樣例數(shù)據(jù)地址如下https://download.csdn.net/download/u010329292/88291585：需要的小伙伴可以供下載學(xué)習(xí):)文章來源地址http://www.zghlxwxcb.cn/news/detail-821797.html

到了這里，關(guān)于基于Python的大區(qū)域SPI標(biāo)準(zhǔn)降水指數(shù)自動批量化處理的文章就介紹完了。如果您還想了解更多內(nèi)容，請?jiān)谟疑辖撬阉鱐OY模板網(wǎng)以前的文章或繼續(xù)瀏覽下面的相關(guān)文章，希望大家以后多多支持TOY模板網(wǎng)！