基于GDAL庫的溫度推算模型的研究與實現(xiàn)

劉文毫，衛(wèi)建國,2，張春梅，楊 豫，劉兆宇

(1.北方民族大學 計算機科學與工程學院，寧夏 銀川 750021；2.中國氣象局旱區(qū)特色農(nóng)業(yè)氣象災害監(jiān)測預警與風險管理重點實驗室，寧夏 銀川 750002；3.寧夏大學 農(nóng)學院，寧夏 銀川 750021)

0 引 言

為對中寧太陽梁試驗區(qū)農(nóng)作物霜凍情況進行研究，需要對研究區(qū)域的溫度分布狀況進行推算。目前對特定區(qū)域的溫度分布狀況進行推算研究的成果較多，Doug M. Smith等[1]通過建立的模型對全球的溫度進行了推算，Y.Radhika等[2]采用支持向量機技術對區(qū)域的大氣溫度進行了預測。國內(nèi)學者在溫度推算方面也有很深入的研究。例如，唐圣鈞等[3]對貴州省南北兩個區(qū)域，結(jié)合多元回歸等方法，分區(qū)建立了基于緯度、海拔和坡向以及坡度因子的氣溫降水模型，在研究區(qū)域進行氣候要素的小網(wǎng)格推算，結(jié)果表明緯度因子對溫度和降水的貢獻最大。衛(wèi)建國[4]根據(jù)溫度分布和熱量理論，在250 m的地理網(wǎng)格資料中，采用多點推算方法，綜合改進距離加權平均法，對運算中的多層推算結(jié)果進行權值訂正，結(jié)果表明經(jīng)過多點推算的溫度圖像平滑而富有變化，從根本上消除了單點推算的臺階現(xiàn)象。劉靜等[5]利用線性內(nèi)插確定界限溫度的初終日期，利用逐旬氣候資料進行小網(wǎng)格訂正，從而得出了紅寺堡等地區(qū)的溫度分布。但是受限于自然條件及其他方面的限制，這些對溫度推算的研究都是采用傳統(tǒng)實測數(shù)據(jù)對地區(qū)的溫度進行分布，而沒有采用地形數(shù)據(jù)結(jié)合經(jīng)緯度信息等資料對區(qū)域的溫度進行推算，在實現(xiàn)方式上大多是使用傳統(tǒng)的C/C++編程技術進行實現(xiàn)，沒有采用最新的Python開發(fā)技術進行實現(xiàn)，因此處理效率相對較低，無法應對基于大批量數(shù)據(jù)的溫度推算工作?；诖爽F(xiàn)狀，文中基于GDAL開源庫，以無人機航測5 m分辨率DEM數(shù)據(jù)及實測的經(jīng)緯度數(shù)據(jù)為數(shù)據(jù)源，實現(xiàn)了基于小區(qū)域溫度推算模型的溫度分布推算。

1 試驗區(qū)概況

太陽梁位于衛(wèi)寧平原之上的中寧縣渠口農(nóng)場，南與中寧縣棗園鄉(xiāng)接壤，以白馬梁為界；北與青銅峽與中寧交界線的胡子溝為界；西倚內(nèi)蒙古阿拉善左旗頭道湖，東靠包蘭鐵路。試驗區(qū)地形復雜多變，由土石山地和緩坡丘陵組成，絕大部分為緩坡丘陵，局地為微坡地。

試驗區(qū)處于中溫帶干旱區(qū)，大陸性季風氣候特征明顯，光照充足、熱量豐富、蒸發(fā)強烈、晝夜溫差大，氣候干旱少雨，春暖快、夏熱短、秋涼早、冬寒長。主要氣象災害有干旱、沙塵暴、霜凍、冰雹、暴雨、干熱風等。

由于研究區(qū)域氣候特征獨特，農(nóng)業(yè)災害類型多樣，因此該試驗區(qū)適合農(nóng)業(yè)氣象災害的研究。又因試驗區(qū)沒有詳細的溫度分布數(shù)據(jù)，因此，需要建立一套合理的溫度推算模型，推算該試驗區(qū)溫度分布狀況。

太陽梁桃園試驗區(qū)布局如圖1所示。

2 數(shù)據(jù)來源及溫度推算模型

項目采用太陽梁試驗區(qū)無人機航測的5 m分辨率細網(wǎng)格DEM數(shù)據(jù)作為數(shù)據(jù)源。通常情況下，某地區(qū)的平均溫度分布情況與該地區(qū)大氣環(huán)境、海拔高度、經(jīng)緯度及小氣候環(huán)境有關[6]。不同季節(jié)太陽輻射[7]和大氣環(huán)流背景是導致氣溫季節(jié)變化的主要原因，同時，坡度坡向也會影響該地區(qū)的輻射分布情況[8]。該項目采用前人[3]對貴州全區(qū)氣候要素的推算模型為基礎，并結(jié)合太陽梁試驗區(qū)的具體氣候特征，設定每項氣候要素的權重系數(shù)，使其能夠更好地實現(xiàn)對太陽梁試驗區(qū)的溫度推算。具體推算模型如下：

圖1 太陽梁桃園試驗區(qū)布局 (左側(cè)及底部陰影區(qū)為研究區(qū)域)

設基準站的平均氣溫為S，則預測點的溫度值可表示為：

S=f(λ,φ,h,χ,ω)+T

(1)

其中，S為某預測點氣候要素值；λ為預測點與基準點的經(jīng)度差值；φ為預測點與基準點的緯度差值；h為預測點與基準點的海拔差值；χ為預測點與基準點的坡向差值；ω為預測點與基準點的坡度差值。

基于氣象科研人員對模型中的每個因子進行分析評估得到各個因子的系數(shù)，同時結(jié)合研究區(qū)域2018年3月的基準點實測的平均溫度為14.8℃，確定最終的溫度推算模型如下：

T=[1λ+1φ+0.03h+(-0.01)χ+0.01ω]*

0.05+14.8

(2)

3 關鍵技術

3.1 GDAL地理空間數(shù)據(jù)抽象庫

GDAL[9]全稱為地理空間數(shù)據(jù)抽象庫(geospatial data abstraction library)，是眾多操作地理空間柵格數(shù)據(jù)的開源庫之一。它利用抽象數(shù)據(jù)模型來表達所支持的各種文件格式，并使用一系列命令行工具來進行數(shù)據(jù)的轉(zhuǎn)換和處理。GDAL包含兩個庫：用于操作地理空間柵格數(shù)據(jù)[10]的GDAL和用于操作地理空間矢量數(shù)據(jù)的OGR，且其支持的數(shù)據(jù)格式也相當廣泛，包括GeoTiff(.tiff)、Erdas Image(.img)、(.grd)等155種柵格數(shù)據(jù)，以及(.shp)、(.gml)、(.GeoPackage)等95種矢量數(shù)據(jù)。GDAL庫還提供了一系列算法的接口[11]，如矢量柵格化、柵格矢量化、遙感數(shù)據(jù)的空間糾正等，并對這些算法提供了可以運行的文件，方便用戶的使用。利用GDAL和OGR這兩個開源庫，可實現(xiàn)對柵格和矢量數(shù)據(jù)的操作及處理。該項目采用開源GDAL庫讀取了.grd格式的原始DEM文件并轉(zhuǎn)換成通用GeoTiff文件，基于GeoTiff文件提取溫度推算所需坡度、坡向等因子并寫入到柵格文件中。通過調(diào)用GDAL庫操作柵格數(shù)據(jù)的具體流程如圖2所示。

圖2 GDAL庫操作柵格數(shù)據(jù)流程

3.2 Numpy庫及Matplotlib庫

Numpy庫是Python進行科學計算的基礎模塊，它是一個提供多維數(shù)據(jù)對象的Python庫，包含了多種衍生對象以及一系列為快速計算數(shù)組而生的例程[12]。Numpy庫最核心的部分是ndarray對象，通過ndarray對象可以實現(xiàn)對柵格數(shù)據(jù)的有效存儲。同時，Numpy庫具有強大的矩陣運算能力，可以實現(xiàn)矩陣求逆、求特征值、求特征向量等操作。Numpy科學運算庫為溫度推算工作的實現(xiàn)提供了支撐。

Matplotlib庫[13]是Python編程語言及Numpy擴展包的可視化操作界面，其提供了面向?qū)ο蟮腁PI用于通用GUI工具包。Matplotlib庫類似于MATLAB中的Pyplot模塊?；贛atplotlib庫可將溫度推算中的每個因子以圖像的形式直觀顯示，并自動生成圖例，生成的圖例不僅方便檢查數(shù)據(jù)的準確性，而且可以幫助靈活地調(diào)整權重系數(shù)，從而實現(xiàn)對溫度分布狀況更加準確的推算。

4 項目實現(xiàn)流程

項目首先采用無人機航拍的細網(wǎng)格DEM數(shù)據(jù)生成海拔格點數(shù)據(jù)，生成該地區(qū)的坡度格點數(shù)據(jù)、坡向格點數(shù)據(jù)。依據(jù)實測的研究區(qū)域經(jīng)緯度范圍，生成包含該區(qū)域每個格點的經(jīng)度緯度值的格點數(shù)據(jù)。由于生成的格點數(shù)據(jù)結(jié)果為標準矩形矩陣，因此需要使用DEM數(shù)據(jù)生成二值矩陣將格點數(shù)據(jù)裁剪為研究區(qū)形狀。獲取了溫度推算所需5個因子的格點數(shù)據(jù)以后，基于Python的ConfigParser模塊，讀取了每個因子的權重系數(shù)并基于Numpy庫實現(xiàn)模型的運算操作。溫度推算技術路線圖如圖3所示。

圖3 溫度推算技術路線圖

4.1 溫度推算因子的提取

(1)海拔因子。

由于DEM數(shù)據(jù)[14-15]已經(jīng)包含研究區(qū)域的海拔數(shù)據(jù)，因此只需對DEM數(shù)據(jù)進行格式轉(zhuǎn)換，將原始.grd格式文件轉(zhuǎn)換為模型推算所需GeoTiff格式。具體步驟如下：通過GDAL.open()函數(shù)打開.grd格式的DEM文件，讀取為Numpy庫支持的ndarray格式，然后通過GDAL庫提供的WriteArray()函數(shù)將數(shù)據(jù)寫入為GeoTiff格式。

(2)坡度坡向因子。

由于GDAL庫包含生成坡度和坡向的函數(shù)，因此該項目通過調(diào)用GDAL庫的函數(shù)生成了包含坡度坡向兩個因子的文件。坡度坡向因子生成步驟如圖4所示。

圖4 坡度坡向因子生成步驟

生成坡度坡向因子需調(diào)用subprocess模塊，subprocess模塊可以通過命令行的方式執(zhí)行GDAL命令，這是執(zhí)行生成坡度坡向因子命令的前提。引用完成subprocess后開始編寫生成坡度坡向因子的命令，需調(diào)用GDAL庫中gdaldem模塊，然后指定如下參數(shù)：指定需要生成的因子(slope、aspect)；指定輸出圖像的壓縮方式，該項目采用LZW壓縮方式；指定其他參數(shù)，如-s指定計算比例，-alg指定計算算法。命令編寫完成后，便可通過subprocess.call()函數(shù)將命令執(zhí)行，分別生成包含研究區(qū)域坡度坡向因子的柵格文件。

(3)經(jīng)度緯度因子。

由于無人機航測的DEM數(shù)據(jù)中包含研究區(qū)的經(jīng)緯度范圍，因此使用GDAL庫中的GetGeoTransform()函數(shù)即可獲取研究區(qū)域的經(jīng)緯度范圍，但獲取到的坐標是投影坐標，還需要編寫函數(shù)將投影坐標轉(zhuǎn)換為經(jīng)緯度坐標，函數(shù)主要用到的是TransformPoint()方法。最終得到研究區(qū)域的經(jīng)緯度范圍：東經(jīng)105°46′～105°50′；北緯37°38′～37°43′?；诮?jīng)度緯度的區(qū)間，使用Numpy的linspace()函數(shù)分別生成該地區(qū)的經(jīng)度緯度分布矩陣。

4.2 溫度推算模型的實現(xiàn)

溫度推算模型的實現(xiàn)需要通過GDAL庫將各個因子數(shù)據(jù)的GeoTiff文件進行讀取，轉(zhuǎn)化為可以進行運算的矩陣。主要操作如下：

(1)疊加因子權重系數(shù)：由于每個因子對該區(qū)域的溫度影響作用不同，因此需要對每個因子疊加不同權重系數(shù)；

(2)疊加總系數(shù)：每個因子疊加完權重系數(shù)后，推算出的溫度區(qū)間可能與實際溫度分布區(qū)間不符，因此需要設置總系數(shù)對溫度分布區(qū)間進行調(diào)整；

(3)疊加基準站溫度T：前兩項運算得出的是目標點與基準站溫度的差值，因此需要與基準站溫度進行疊加，從而實現(xiàn)對整個區(qū)域的溫度進行推算。

以上三種對矩陣的操作皆可通過Numpy庫進行實現(xiàn)，這里需要注意的是，在執(zhí)行完每一步操作以后，都要使用Matplotlib庫中的imshow方法將運算完成的矩陣以圖形化展示，這樣可以保證在運算出現(xiàn)異常值時及時發(fā)現(xiàn)，同時可以對每一步運算結(jié)果進行驗證，提升溫度推算工作的魯棒性。

基于用戶體驗需求，工具在使用時需要隱藏具體的實現(xiàn)細節(jié)，同時又需要將權重系數(shù)配置接口開放給用戶。因此引入了Python的ConfigParser模塊，通過此模塊可以實現(xiàn)讀取配置文件中的權重系數(shù)實現(xiàn)溫度的推算，因此用戶僅需在配置文件中配置好每個因子的權重系數(shù)即可生成研究區(qū)域的溫度分布圖。

4.3 溫度分布圖的生成

基于Python技術對提取到的5個因子通過溫度推算模型進行運算以后，最終生成包含每個網(wǎng)格點溫度值的GeoTiff文件。雖然基于python的Matplotlib模塊可以將溫度推算中的每個因子進行圖形化展示，方便檢查數(shù)據(jù)的準確性，但是在實際應用中發(fā)現(xiàn)Matplotlib模塊在繪制最終正果溫度分布圖時存在很多不足的地方，如無法設置圖例樣式，無法添加指南針比例尺等。因此決定采用ArcMap制圖工具[16-17]制作最終的正果溫度分布圖。ArcMap是一款強大的制圖工具，可以方便地添加圖例、比例尺等組件，同時可以靈活地根據(jù)數(shù)值范圍設置不同的配色。

5 實驗結(jié)果

最終生成的正果溫度分布如圖5所示，由圖中不同標注的區(qū)域可以看到，研究區(qū)域的西側(cè)區(qū)域溫度較低，而研究區(qū)域北部區(qū)域及東南部區(qū)域溫度偏高，其中溫度最高的區(qū)域為東北部區(qū)域，溫度達到了15.3 ℃以上。

圖5 正果溫度分布

為了對項目建立的溫度推算模型進行驗證，分別在太陽梁試驗區(qū)選取三個觀測點，對觀測點的溫度進行逐時監(jiān)測，統(tǒng)計了2018年3月1日至31日的每小時溫度。具體觀測位置如圖6所示。

圖6 監(jiān)測點位置分布

監(jiān)測點A位于試驗區(qū)的北部區(qū)域，監(jiān)測點B位于試驗區(qū)的西南部區(qū)域，監(jiān)測點C位于試驗區(qū)的東南部區(qū)域。通過統(tǒng)計得出了26日至30日的平均溫度：觀測點A的平均溫度為15.1℃，觀測點B的平均溫度為15.13℃，觀測點C的平均溫度為14.9℃。經(jīng)過與模型推算結(jié)果對比，得出推算結(jié)果與地區(qū)實際溫度分布狀況比較吻合，達到了預期的目標。同時，溫度推算設定的五個因子中，海拔對溫度的推算影響最大，其次是坡度和坡向?？梢娫谛^(qū)域網(wǎng)格溫度推算中，經(jīng)度和緯度的影響權重較小，基本上可以忽略不計，這一點與大區(qū)域的網(wǎng)格點溫度推算不同。

6 結(jié)束語

基于對溫度推算模型及開源柵格處理GDAL庫的研究，設計并實現(xiàn)了一種基于DEM數(shù)據(jù)的小網(wǎng)格溫度推算模型，并基于開源柵格處理庫GDAL庫及科學運算庫Numpy庫進行實現(xiàn)。首先介紹了研究區(qū)概況及模型的設計，其次介紹了實現(xiàn)模型時用到的關鍵技術，并詳細介紹了溫度推算模型的實現(xiàn)步驟，最后對溫度推算結(jié)果進行了分析與驗證?；陂_源GDAL庫實現(xiàn)的溫度推算模型相較于傳統(tǒng)的實現(xiàn)方式具有可擴展性強且運算速度快等特點。該成果同樣可推廣應用于其他模型的推算工作，如基于該模型中的坡度、坡向、海拔因子，可以應用于小區(qū)域風場分布情況的推算，這也是該項目接下來的研究方向。