格點(diǎn)插值計(jì)算面雨量方法的Python編程實(shí)現(xiàn)

李勝

關(guān)鍵詞：面雨量；格點(diǎn)；插值；射線(xiàn)法

中圖分類(lèi)號(hào)：TP311.1 文獻(xiàn)標(biāo)識(shí)碼：A

文章編號(hào)：1009-3044（2022）36-0039-04

1 前言

近若干年來(lái)，隨著全球氣候變暖，極端天氣氣候事件呈現(xiàn)增多增強(qiáng)趨勢(shì)。根據(jù)《WMO天氣、氣候和水極端事件造成的死亡和經(jīng)濟(jì)損失圖集（1970—2019）》數(shù)據(jù)顯示，50年間，全球共發(fā)生1.1萬(wàn)多起與天氣、氣候和水相關(guān)的災(zāi)害，共造成了200多萬(wàn)人死亡和3.64萬(wàn)億美元的經(jīng)濟(jì)損失。中國(guó)是世界上受氣象災(zāi)害影響最嚴(yán)重的國(guó)家之一，氣象災(zāi)害所造成的損失占到了自然災(zāi)害損失的70%以上。因極端天氣，特別是當(dāng)暴雨、短時(shí)強(qiáng)降水等惡劣天氣發(fā)生時(shí)，由于降雨過(guò)于集中很快造成大的地表徑流，在城市或平原地區(qū)可能因排水不及時(shí)而產(chǎn)生大量的積水，致使土地、房屋等漬水、受淹而造成雨澇災(zāi)害。在山區(qū)，由于受地形阻擋作用，常常會(huì)形成繞流和爬流等，從而易誘發(fā)山體滑坡和泥石流等次生災(zāi)害。因此了解降雨區(qū)域內(nèi)降雨總體情況十分必要。

面雨量從定義上講是指計(jì)算統(tǒng)計(jì)區(qū)域內(nèi)單位面積上降水量的平均值，能比較客觀(guān)真實(shí)地反映某個(gè)區(qū)域的整體降水情況。面雨量是水利部門(mén)開(kāi)展洪水預(yù)報(bào)非常重要的基礎(chǔ)數(shù)據(jù)，開(kāi)展面雨量計(jì)算能更好地為各級(jí)政府組織防汛抗洪以及水庫(kù)涵閘調(diào)度、攔蓄泄洪等重大決策提供依據(jù)，在氣象與水文水利部門(mén)的合作中發(fā)揮著十分重要的作用，是服務(wù)地方經(jīng)濟(jì)建設(shè)和防災(zāi)減災(zāi)工作的一個(gè)重要手段。

2 面雨量簡(jiǎn)介

氣象學(xué)中常用的降雨量，是指一個(gè)氣象觀(guān)測(cè)站點(diǎn)上測(cè)得的可代表觀(guān)測(cè)站點(diǎn)周?chē)粋€(gè)小區(qū)域的平均降水量，也就是通常所說(shuō)的某地下了多大的雨。雨量是分布得極不均勻的一種氣象要素，在大氣候環(huán)境相似的條件下，山區(qū)雨量多于平原，高地雨量多于河谷低地。而對(duì)于江河湖泊來(lái)說(shuō)，孤立的點(diǎn)狀的雨量數(shù)據(jù)并不能代表整個(gè)流域的降雨情況，流域范圍內(nèi)單位面積上的平均降水量才能客觀(guān)反映真實(shí)降雨情況[1]。

3 面雨量計(jì)算

用rM表示面雨量，用ri表示各點(diǎn)實(shí)測(cè)雨量，用Ai表示該點(diǎn)所代表的面積，面雨量的計(jì)算方法如公式（1）所示：

如果簡(jiǎn)單地求算術(shù)平均，是不能代表整個(gè)地區(qū)的真正面平均雨量的。于是，便發(fā)展了一些計(jì)算面雨量的方法。以不同的方法決定Ai，就得到不同的面平均雨量計(jì)算方法。面雨量常見(jiàn)計(jì)算方法有等值線(xiàn)法、數(shù)值法、算術(shù)平均法等。

3.1 等值線(xiàn)法

等值線(xiàn)法是以相鄰兩條等雨量線(xiàn)的平均值ri以及兩線(xiàn)之間的面積Ai，代入公式（1）即可，但由于雨量等值線(xiàn)往往是不規(guī)劃的曲線(xiàn)，用此方法計(jì)算兩條曲線(xiàn)所包含區(qū)間的面積Ai，其過(guò)程相當(dāng)復(fù)雜[2]。

3.2 數(shù)值法

數(shù)值法是指泰森多邊形法或三角形法。泰森多邊形法是將流域內(nèi)每3個(gè)最靠近的雨量站連成一個(gè)三角形，流域邊緣的站可利用流域外相鄰的站來(lái)連接三角形。對(duì)三角形各邊分別作垂直平分線(xiàn)，這些線(xiàn)組成若干個(gè)多邊形，使得每個(gè)多邊形內(nèi)有且只有一個(gè)雨量站。每個(gè)多邊形的面積為Ai，其雨量為ri，代入公式（1），便可求得rM。用此方法計(jì)算多邊形面積Ai過(guò)程更加復(fù)雜和煩瑣[3]。

三角形法如泰森多邊形法那樣，將相鄰的3個(gè)站連成三角形。任意一個(gè)三角形的面積為Ai，3個(gè)頂點(diǎn)是3個(gè)站，雨量分別為ri1、ri2、ri3，3個(gè)雨量的平均值是ri，代入公式（1）即可。可見(jiàn)，三角形法相對(duì)于泰森多邊形法來(lái)說(shuō)，算法上要方便一些，但是計(jì)算過(guò)程依然煩瑣。

3.3 算術(shù)平均法

算術(shù)平均法相對(duì)比較簡(jiǎn)單易行，將區(qū)域內(nèi)所有站點(diǎn)的雨量計(jì)算其平均值即可。用此方法計(jì)算面雨量結(jié)果的可用性要依賴(lài)于區(qū)域內(nèi)站點(diǎn)數(shù)量的多寡，而現(xiàn)實(shí)情況是區(qū)域自動(dòng)氣象站不可能特別密布，這就需要引入格點(diǎn)插值計(jì)算方法了。

3.4 格點(diǎn)插值

插值法是離散函數(shù)逼近的重要方法，即在一定的區(qū)域范圍內(nèi)按照一定的間隔，根據(jù)離散函數(shù)在有限個(gè)點(diǎn)處的取值狀況估算出函數(shù)在其他點(diǎn)處的近似值[4]。其實(shí)現(xiàn)原理如圖1所示。

但是，由于格點(diǎn)插值計(jì)算的區(qū)域是取江河湖泊流域邊界上下左右四個(gè)頂點(diǎn)所組成的區(qū)域（如圖2紅色外框所包圍的區(qū)域），是規(guī)則的四邊形，而實(shí)際江河湖泊流域邊界是不規(guī)則的區(qū)間（如圖2黑色線(xiàn)條所包圍的區(qū)域），位于格點(diǎn)區(qū)域內(nèi)，即被包含在格點(diǎn)區(qū)域內(nèi)。因此要對(duì)計(jì)算出的所有格點(diǎn)進(jìn)行判斷是否位于不規(guī)則區(qū)間內(nèi)，只有位于不規(guī)則區(qū)間內(nèi)的格點(diǎn)才是計(jì)算面雨量所需要的。

3.5 射線(xiàn)法判斷格點(diǎn)位置

通常采用射線(xiàn)法即用水平掃描線(xiàn)法或垂直線(xiàn)法即射線(xiàn)法來(lái)判斷一點(diǎn)是否在區(qū)域內(nèi)。在不考慮非歐空間的情況下，對(duì)于平面內(nèi)任意閉合曲線(xiàn)把平面分割成了多邊形區(qū)域內(nèi)、外兩部分。當(dāng)射線(xiàn)穿越多邊形區(qū)域邊界時(shí)，要么是進(jìn)入多邊形區(qū)域要么是穿出多邊形區(qū)域[5]。

如果某個(gè)點(diǎn)在多邊形區(qū)域內(nèi)部，射線(xiàn)第一次穿越邊界一定是穿出多邊形區(qū)域，當(dāng)再次穿越邊界一定是穿入多邊形區(qū)域；如果某個(gè)點(diǎn)在多邊形外部，射線(xiàn)第一次穿越邊界一定是穿入多邊形區(qū)域，當(dāng)再次穿越邊界一定是穿出多邊形區(qū)域。如圖3所示。

由于射線(xiàn)是由線(xiàn)段的一端無(wú)限延伸的直的線(xiàn)，而多邊形區(qū)域的邊界是固定的，因此射線(xiàn)最后一次穿越多邊形邊界，一定是穿出多邊形，到達(dá)外部。假定某點(diǎn)在多邊形外部，射線(xiàn)進(jìn)入和穿出多邊形至少需要零次或兩次即偶數(shù)次才能到達(dá)多邊形外部。假定某點(diǎn)在多邊形內(nèi)部，射線(xiàn)穿出或穿出、進(jìn)入、再穿出多邊形至少需要一次或三次即奇數(shù)次才能到達(dá)多邊形外部。因此，由射線(xiàn)和多邊形邊界相交的次數(shù)是奇數(shù)還是偶數(shù)就可以推斷出某點(diǎn)是在多邊形內(nèi)部還是在多邊形外部。

4 格點(diǎn)插值編程實(shí)現(xiàn)

Python作為大多數(shù)平臺(tái)上編寫(xiě)腳本和快速開(kāi)發(fā)應(yīng)用的編程語(yǔ)言，除自身提供豐富的標(biāo)準(zhǔn)庫(kù)外，其解釋器還易于擴(kuò)展。NumPy（Numerical Python）是Py?thon最為重要的一個(gè)擴(kuò)展模塊，特點(diǎn)針對(duì)數(shù)組運(yùn)算提供大量的數(shù)學(xué)函數(shù)庫(kù)，常被用來(lái)做數(shù)組、矩陣甚至多維運(yùn)算。其中的meshgrid函數(shù)功能是根據(jù)傳入的兩個(gè)一維數(shù)組參數(shù)生成兩個(gè)數(shù)組元素的列表，即將兩個(gè)一維數(shù)組生成一個(gè)二維矩陣，對(duì)應(yīng)兩個(gè)一維數(shù)組中所有的（x，y）對(duì)。scipy.interpolate是Python實(shí)現(xiàn)各種插值法的擴(kuò)展模塊，其griddata函數(shù)可以方便實(shí)現(xiàn)二維插值。插值方法有最鄰近插值nearest、階梯插值z(mì)ero、線(xiàn)性插值slinear或linear、樣條曲線(xiàn)插值quadratic或cu?bic，由于樣條曲線(xiàn)插值讓插值結(jié)果具有更好的平滑性，所以本文選用樣條曲線(xiàn)插值[6]。

4.1 計(jì)算區(qū)域雨量格點(diǎn)值

根據(jù)前文面雨量算法思路，先計(jì)算出四邊形格點(diǎn)，然后判斷所有格點(diǎn)是否在不規(guī)則區(qū)間內(nèi)，剔除在不規(guī)則區(qū)間外的格點(diǎn)后，即取得在不規(guī)則區(qū)間內(nèi)的格點(diǎn)，最后計(jì)算這些剩余格點(diǎn)雨量的平均值，即是最終的面雨量。具體實(shí)現(xiàn)過(guò)程如下。

4.1.1 計(jì)算生成經(jīng)緯度網(wǎng)格

先根據(jù)該區(qū)域的經(jīng)緯度范圍計(jì)算生成經(jīng)緯度網(wǎng)格，再根據(jù)經(jīng)緯度網(wǎng)格插值計(jì)算雨量的格點(diǎn)值。假定某區(qū)域內(nèi)站點(diǎn)雨量實(shí)況文件內(nèi)容如下（格式為經(jīng)度，緯度，雨量）：

117.33，32.39，0

......

# 生成經(jīng)緯度網(wǎng)格，其中minLon，maxLon，min?Lat，maxLat分別是某區(qū)域的最小經(jīng)度、最大經(jīng)度、最小緯度和最大緯度

LonLatArea=[minLon，maxLon，minLat，maxLat]

LonX=np. linspace（minLon， maxLon， int（（maxLonminLon）/經(jīng)度分辨率））

LatY=np.linspace（minLat，maxLat，int（（maxLat-min?Lat）/緯度分辨率））

LonGrid，LatGrid=np.meshgrid（LonX，LatY）

4.1.2 插值計(jì)算出雨量格點(diǎn)值

# 先將雨量文件按行讀取到列表PreTemp中，然后分別提取經(jīng)緯度和站點(diǎn)雨量值到Points和Pres列表中

for i in range（len（PreTemp））：

line=str（PreTemp[i][0]）

preinfo=line.split（'，'）

Points. append（[float（preinfo[0]）， float（prein?fo[1]）]）

Pres.append（float（preinfo[2]））

# 根據(jù)經(jīng)緯度網(wǎng)格插值計(jì)算出雨量格點(diǎn)值

PreGrid=griddata（Points， Pres， （LonGrid， LatGrid），method='cubic'，fill_value=0）

# 如果插值計(jì)算出的雨量小于0，則賦值為0

PreGrid[PreGrid<0]=0

4.2 判斷邊界內(nèi)格點(diǎn)

根據(jù)射線(xiàn)法基本原理，假若有一疑問(wèn)點(diǎn)P（x，y），要判斷它是否在多邊形內(nèi)，可從該疑問(wèn)點(diǎn)向左或向右引水平掃描線(xiàn)（即射線(xiàn)）并計(jì)算此線(xiàn)段與區(qū)域邊界的相交次數(shù)c。如果c為奇數(shù)，認(rèn)為疑問(wèn)點(diǎn)在多邊形內(nèi)；如果c為偶數(shù)，則疑問(wèn)點(diǎn)在多邊形外。

假定某區(qū)域邊界文件經(jīng)緯度文件內(nèi)容如下（格式為經(jīng)度，緯度）：

117.40178，31.69473

117.40212，31.69395

......

117.40178，31.69473

4.2.1 判斷坐標(biāo)點(diǎn)是否在外包矩形內(nèi)先要判斷某個(gè)格點(diǎn)坐標(biāo)是否在外包矩形內(nèi)，即是否在minLon，maxLon，minLat，maxLat范圍內(nèi)，以便快速剔除多余格點(diǎn)，然后利用射線(xiàn)法判斷是否在多邊形范圍內(nèi)。

# poi是格點(diǎn)坐標(biāo)，sbox是邊界經(jīng)緯度列表

def isPointInRect（poi，sbox，toler=0.0001）：

if poi[0]>sbox[0][0] and poi[0]<sbox[1][0] and poi[1]>sbox[0][1] and poi[1]<sbox[1][1]：

return True

if toler>0：

pass

return False

4.2.2 判斷射線(xiàn)與邊界是否有交點(diǎn)

# 判斷點(diǎn)，邊起點(diǎn)，邊終點(diǎn)，都是[lon，lat]格式列表def isRayIntersectsSegment（poi，s_poi，e_poi）：

# 排除與射線(xiàn)平行、重合，線(xiàn)段首尾端點(diǎn)重合的情況

if s_poi[1]==e_poi[1]：

return False

# 線(xiàn)段在射線(xiàn)上邊

if s_poi[1]>poi[1] and e_poi[1]>poi[1]：

return False

# 線(xiàn)段在射線(xiàn)下邊

if s_poi[1]<poi[1] and e_poi[1]<poi[1]：

return False

# 交點(diǎn)為下端點(diǎn)，對(duì)應(yīng)spoint

if s_poi[1]==poi[1] and e_poi[1]>poi[1]：

return False

# 交點(diǎn)為下端點(diǎn)，對(duì)應(yīng)epoint

if e_poi[1]==poi[1] and s_poi[1]>poi[1]：

return False

# 線(xiàn)段在射線(xiàn)左邊

if s_poi[0]<poi[0] and e_poi[1]<poi[1]：

return False

# 求交點(diǎn)

xseg=e_poi[0]-（e_poi[0]-s_poi[0]）*（e_poi[1]-poi

[1]）/（e_poi[1]-s_poi[1]）

# 交點(diǎn)在射線(xiàn)起點(diǎn)的左側(cè)

if xseg<poi[0]：

return False

return True

4.2.3 判斷格點(diǎn)坐標(biāo)是否在邊界內(nèi)

def isPointInPoly（poi，poly，toler=0.0001）：

# 先判斷格點(diǎn)坐標(biāo)是否在外包矩形內(nèi)，如果不在，函數(shù)返回False if not isPointInRect（poi，RectLatLon，toler）：

return False

# sinsc是交點(diǎn)個(gè)數(shù)

sinsc=0

# 對(duì)格點(diǎn)坐標(biāo)進(jìn)行判斷是否在多邊形范圍內(nèi)

for i in range（len（poly）-1）：

s_poi=poly[i]

e_poi=poly[i+1]

if isRayIntersectsSegment（poi，s_poi，e_poi）：

sinsc+=1

return True if sinsc % 2==1 else False

4.3 計(jì)算格點(diǎn)面雨量

將上述計(jì)算出的屬于邊界范圍內(nèi)的格點(diǎn)雨量進(jìn)行算術(shù)平均即可求得面雨量。

# 提取并計(jì)算邊界范圍內(nèi)的格點(diǎn)雨量

gribcnt=presum=0

for m in range（len（PreGrid））：

for n in range（len（PreGrid[m]））：

if isPointInPoly（[LonGrid[m][n]，LatGrid[m][n]]，

AreaPloy，toler=0.0001）：

presum+=PreGrid[m][n]

gribcnt+=1

PreAvg=round（presum/gribcnt，1）

5 結(jié)束語(yǔ)

由于不同的面雨量計(jì)算方法有不同的出發(fā)點(diǎn)和目標(biāo)，也有其適用條件。在一定條件下有某一種方法比其他方法準(zhǔn)確，因此不能籠統(tǒng)地說(shuō)某種方法最好。隨著氣象部門(mén)觀(guān)測(cè)站點(diǎn)安裝越來(lái)越密集，站點(diǎn)布局不再是以往零散分布，也已呈現(xiàn)網(wǎng)格化，這些觀(guān)測(cè)站點(diǎn)雨量數(shù)據(jù)再經(jīng)格點(diǎn)插值后計(jì)算的面雨量可以很好地作為氣象預(yù)報(bào)服務(wù)和水文水位預(yù)報(bào)的重要數(shù)據(jù)支撐。