格點插值計算面雨量方法的Python編程實現

李勝

關鍵詞：面雨量；格點；插值；射線法

中圖分類號：TP311.1 文獻標識碼：A

文章編號：1009-3044（2022）36-0039-04

1 前言

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

面雨量從定義上講是指計算統計區域內單位面積上降水量的平均值，能比較客觀真實地反映某個區域的整體降水情況。面雨量是水利部門開展洪水預報非常重要的基礎數據，開展面雨量計算能更好地為各級政府組織防汛抗洪以及水庫涵閘調度、攔蓄泄洪等重大決策提供依據，在氣象與水文水利部門的合作中發揮著十分重要的作用，是服務地方經濟建設和防災減災工作的一個重要手段。

2 面雨量簡介

氣象學中常用的降雨量，是指一個氣象觀測站點上測得的可代表觀測站點周圍一個小區域的平均降水量，也就是通常所說的某地下了多大的雨。雨量是分布得極不均勻的一種氣象要素，在大氣候環境相似的條件下，山區雨量多于平原，高地雨量多于河谷低地。而對于江河湖泊來說，孤立的點狀的雨量數據并不能代表整個流域的降雨情況，流域范圍內單位面積上的平均降水量才能客觀反映真實降雨情況[1]。

3 面雨量計算

用rM表示面雨量，用ri表示各點實測雨量，用Ai表示該點所代表的面積，面雨量的計算方法如公式（1）所示：

如果簡單地求算術平均，是不能代表整個地區的真正面平均雨量的。于是，便發展了一些計算面雨量的方法。以不同的方法決定Ai，就得到不同的面平均雨量計算方法。面雨量常見計算方法有等值線法、數值法、算術平均法等。

3.1 等值線法

等值線法是以相鄰兩條等雨量線的平均值ri以及兩線之間的面積Ai，代入公式（1）即可，但由于雨量等值線往往是不規劃的曲線，用此方法計算兩條曲線所包含區間的面積Ai，其過程相當復雜[2]。

3.2 數值法

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

三角形法如泰森多邊形法那樣，將相鄰的3個站連成三角形。任意一個三角形的面積為Ai，3個頂點是3個站，雨量分別為ri1、ri2、ri3，3個雨量的平均值是ri，代入公式（1）即可。可見，三角形法相對于泰森多邊形法來說，算法上要方便一些，但是計算過程依然煩瑣。

3.3 算術平均法

算術平均法相對比較簡單易行，將區域內所有站點的雨量計算其平均值即可。用此方法計算面雨量結果的可用性要依賴于區域內站點數量的多寡，而現實情況是區域自動氣象站不可能特別密布，這就需要引入格點插值計算方法了。

3.4 格點插值

插值法是離散函數逼近的重要方法，即在一定的區域范圍內按照一定的間隔，根據離散函數在有限個點處的取值狀況估算出函數在其他點處的近似值[4]。其實現原理如圖1所示。

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

3.5 射線法判斷格點位置

通常采用射線法即用水平掃描線法或垂直線法即射線法來判斷一點是否在區域內。在不考慮非歐空間的情況下，對于平面內任意閉合曲線把平面分割成了多邊形區域內、外兩部分。當射線穿越多邊形區域邊界時，要么是進入多邊形區域要么是穿出多邊形區域[5]。

如果某個點在多邊形區域內部，射線第一次穿越邊界一定是穿出多邊形區域，當再次穿越邊界一定是穿入多邊形區域；如果某個點在多邊形外部，射線第一次穿越邊界一定是穿入多邊形區域，當再次穿越邊界一定是穿出多邊形區域。如圖3所示。

由于射線是由線段的一端無限延伸的直的線，而多邊形區域的邊界是固定的，因此射線最后一次穿越多邊形邊界，一定是穿出多邊形，到達外部。假定某點在多邊形外部，射線進入和穿出多邊形至少需要零次或兩次即偶數次才能到達多邊形外部。假定某點在多邊形內部，射線穿出或穿出、進入、再穿出多邊形至少需要一次或三次即奇數次才能到達多邊形外部。因此，由射線和多邊形邊界相交的次數是奇數還是偶數就可以推斷出某點是在多邊形內部還是在多邊形外部。

4 格點插值編程實現

Python作為大多數平臺上編寫腳本和快速開發應用的編程語言，除自身提供豐富的標準庫外，其解釋器還易于擴展。NumPy（Numerical Python）是Py?thon最為重要的一個擴展模塊，特點針對數組運算提供大量的數學函數庫，常被用來做數組、矩陣甚至多維運算。其中的meshgrid函數功能是根據傳入的兩個一維數組參數生成兩個數組元素的列表，即將兩個一維數組生成一個二維矩陣，對應兩個一維數組中所有的（x，y）對。scipy.interpolate是Python實現各種插值法的擴展模塊，其griddata函數可以方便實現二維插值。插值方法有最鄰近插值nearest、階梯插值zero、線性插值slinear或linear、樣條曲線插值quadratic或cu?bic，由于樣條曲線插值讓插值結果具有更好的平滑性，所以本文選用樣條曲線插值[6]。

4.1 計算區域雨量格點值

根據前文面雨量算法思路，先計算出四邊形格點，然后判斷所有格點是否在不規則區間內，剔除在不規則區間外的格點后，即取得在不規則區間內的格點，最后計算這些剩余格點雨量的平均值，即是最終的面雨量。具體實現過程如下。

4.1.1 計算生成經緯度網格

先根據該區域的經緯度范圍計算生成經緯度網格，再根據經緯度網格插值計算雨量的格點值。假定某區域內站點雨量實況文件內容如下（格式為經度，緯度，雨量）：

117.33，32.39，0

......

# 生成經緯度網格，其中minLon，maxLon，min?Lat，maxLat分別是某區域的最小經度、最大經度、最小緯度和最大緯度

LonLatArea=[minLon，maxLon，minLat，maxLat]

LonX=np. linspace（minLon， maxLon， int（（maxLonminLon）/經度分辨率））

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

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

4.1.2 插值計算出雨量格點值

# 先將雨量文件按行讀取到列表PreTemp中，然后分別提取經緯度和站點雨量值到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]））

# 根據經緯度網格插值計算出雨量格點值

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

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

PreGrid[PreGrid<0]=0

4.2 判斷邊界內格點

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

假定某區域邊界文件經緯度文件內容如下（格式為經度，緯度）：

117.40178，31.69473

117.40212，31.69395

......

117.40178，31.69473

4.2.1 判斷坐標點是否在外包矩形內先要判斷某個格點坐標是否在外包矩形內，即是否在minLon，maxLon，minLat，maxLat范圍內，以便快速剔除多余格點，然后利用射線法判斷是否在多邊形范圍內。

# poi是格點坐標，sbox是邊界經緯度列表

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 判斷射線與邊界是否有交點

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

# 排除與射線平行、重合，線段首尾端點重合的情況

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

return False

# 線段在射線上邊

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

return False

# 線段在射線下邊

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

return False

# 交點為下端點，對應spoint

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

return False

# 交點為下端點，對應epoint

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

return False

# 線段在射線左邊

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

return False

# 求交點

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

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

# 交點在射線起點的左側

if xseg<poi[0]：

return False

return True

4.2.3 判斷格點坐標是否在邊界內

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

# 先判斷格點坐標是否在外包矩形內，如果不在，函數返回False if not isPointInRect（poi，RectLatLon，toler）：

return False

# sinsc是交點個數

sinsc=0

# 對格點坐標進行判斷是否在多邊形范圍內

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 計算格點面雨量

將上述計算出的屬于邊界范圍內的格點雨量進行算術平均即可求得面雨量。

# 提取并計算邊界范圍內的格點雨量

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 結束語

由于不同的面雨量計算方法有不同的出發點和目標，也有其適用條件。在一定條件下有某一種方法比其他方法準確，因此不能籠統地說某種方法最好。隨著氣象部門觀測站點安裝越來越密集，站點布局不再是以往零散分布，也已呈現網格化，這些觀測站點雨量數據再經格點插值后計算的面雨量可以很好地作為氣象預報服務和水文水位預報的重要數據支撐。