基于ArcGIS的GP服務實現動態平均雨量的計算方法研究

肖惠珍

(廈門精圖信息技術有限公司，福建 廈門 361008)

0 引 言

通常雨量站所觀測到的降雨量，只能代表該雨量站周圍較小范圍的降雨情況，不能代表全流域或某一區域的平均降雨量，因此單獨某個雨量站的降雨數據不能作為洪澇預報與評估的依據。因此，需要采用全流域或某一區域的所有雨量站的降雨數據來計算區域平均降雨量。

目前，計算區域平均降雨量的方法很多，常用的有算術平均法、數值法、等值線法、泰森多邊形法等。在這些方法中，泰森多邊形法最適合流域或區域內雨量站或降雨量分布不均勻的情況，是我國計算區域平均雨量常用的方法，被水利、氣象、環境等部門廣泛應用。根據廈門的區域特點，雨量觀測站分散，降雨量時空分布不均，采用泰森多邊形法計算區域平均雨量非常適合。

1 研究方法

1.1 構建Delaunay三角網

基于ArcGIS平臺的泰森多邊形構建，是通過雨量監測點的雨量值計算區域平均雨量的方法，來實現離散點構建泰森多邊形，從而獲得區域的連續值。

構建Delaunay三角形網格是構建泰森多邊形的先決條件，就是由離散數據點構建三角網，即確定哪3個數據點構成一個三角形，也稱為自動聯接三角網。ArcGIS在繪制三角形時，要通過程序判斷得出最優三角形。從所有雨量站點任取一點作為第一個三角形的第一個頂點，找出離該點最近的一點作為第二個頂點， 然后再利用斜三角形的余弦定理，找出與第一、第二頂點形成夾角最大的一點作為第三個頂點，連接這3個頂點便得到最優三角形，最終得到優選三角網格。

1.2 泰森多邊形的裁剪

根據廈門市所有雨量監測站點構建出的泰森多邊形，以廈門市的最大經緯度和最小經緯度為邊界進行構建泰森多邊形，因此計算出來的邊界為長方形，需要對泰森多邊形進行裁剪。能否精準求出裁剪結果，是保證計算結果準確率的關鍵。利用ArcGIS的相交工具，將廈門市的行政區劃與泰森多邊形進行疊加，將區劃屬性附給泰森多邊形，從而進行行政區劃內泰森多邊形的面積和區域總面積的計算。

1.3 區域平均雨量計算方法

在泰森多邊形計算平均雨量時，非常重要的一步就是面積的計算，通過 ArcGIS 中計算幾何功能可快速得出各多邊形面積，面積可表示為：

在計算幾何工具中，坐標系采用廈門小92，單位選擇平方千米。之后用自定義公式計算出各多邊形與區域面積的比即雨量權重系數，通過實時雨情數據庫提出雨量再乘以權重系數， 最后按行政分區匯總得出所需的平均雨量。

設每個雨量觀測點的降雨量為xi，每個對應的泰森多邊形的面積為fx，則區域平均降雨量可按下式求得：

式中：xi為雨量觀測點的降雨量，mm；fx為泰森多邊形的面積，km2；n為區域內雨量觀測點或泰森多邊形的個數；F為區域的總面積，km2；Ai為雨量站權重系數。

2 基于GP服務的廈門市區域平均雨量

GP服務是ArcGIS一個自帶的特殊功能，通過Model Builder，通過拖動數據、工具，設置配置參數等操作，將整個運算過程建成一個完整的模型。同時可添加腳本指令語言，最后將模型發布成GP服務，并供系統調用，形成前后交互，實現實時、快速的運算模式。

2.1 構建模型

首先利用ArcGIS軟件的添加執行工具(圖1)，編寫計算雨量監測站點的平均降雨量的Python腳本，在模型“計算面積”工具中導入編寫好的腳本。運行該模型，ArcGIS自動執行一系列操作，包括獲得實時監測數據、建泰森多邊形、行政區劃裁剪、計算區域平均雨量、發布服務，最終的計算結果及泰森多邊形服務都在系統前端展示(圖2)，表明該模型構建成功。通過ArcGIS將該模型發布共享為地理處理服務，與系統進行交互，通過系統前端時間設置，ArcGIS執行模型獲取結果。

圖1 廈門市區域平均雨量GP模型Fig.1 GP Model of Regional Average Rainfall in Xiamen City

圖2 廈門市區域平均雨量成果圖Fig.2 Regional average precipitation map of Xiamen

2.2 計算站點平均降雨量Python代碼

Python對縮進要求嚴格，在sublime中縮進用tab，用空格會報錯，遇到錯誤可以去搜下有相關資料說明，主要用到arcpy這個類；

以下是代碼：

# -*- coding: utf-8 -*-

# 需加編碼，不然print中文字符會報錯

import arcpy

import time

import sys

reload(sys)

sys.setdefaultencoding( ″gbk″ )

#arcpy.env.workspace = r″D: est″

arcpy.env.overwriteOutput = True

addFieldName = ″AVG_RAIN″

#獲取參數方法

inFeature = arcpy.GetParameterAsText(0)

outFeature = arcpy.GetParameterAsText(1)

# 直接通過路徑去獲取輸入要素類和輸出要素類

#inFeature = ″calculateArea.shp″

#outFeature = ″copy.shp″

#復制輸入要素，后面進行編輯

arcpy.CopyFeatures_management(inFeature,outFeature)

print ″要素復制成功！″

print ″正在添加字段......″+ addFieldName

arcpy.AddField_management(outFeature, addFieldName, ″DOUBLE″,″″,″″,18);

print ″字段添加成功......″+ addFieldName

fields = [″RAIN″, ″F_AREA″, ″QHDM″]

HC_Area = 0

JM_Area = 0

TA_Area = 0

XA_Area = 0

HL_Area = 0

SM_Area = 0

#計算廈門各個行政區的面積

with arcpy.da.SearchCursor(outFeature, fields) as cursor:

for row in cursor:

if row[2] == ′350203′:

SM_Area = SM_Area + row[1]

elif row[2] == ′350206′:

HL_Area = HL_Area + row[1]

elif row[2] == ′350205′:

HC_Area = HC_Area + row[1]

elif row[2] == ′350211′:

JM_Area = JM_Area + row[1]

elif row[2] == ′350212′:

TA_Area = TA_Area + row[1]

elif row[2] == ′350213′:

XA_Area = XA_Area + row[1]

else:

print ″無符合條件″

print ″海滄區的面積為：{0},集美區的面積為：{1},同安區的面積為：{2},翔安區的面積為：{3},湖里區的面積為：{4},思明區的面積為：{5}″.format(HC_Area,JM_Area,TA_Area,XA_Area,HL_Area,SM_Area)

#把需要用到的字段加到下面數組當中，用游標去查詢更新

updateFields = [″RAIN″, ″F_AREA″,″AVG_RAIN″, ″QHDM″]

with arcpy.da.UpdateCursor(outFeature, updateFields) as cursor:

for row in cursor:

if row[3] == ′350203′:

row[2] = row[0]*(row[1]/SM_Area)

elif row[3] == ′350206′:

row[2] = row[0]*(row[1]/HL_Area)

elif row[3] == ′350205′:

row[2] = row[0]*(row[1]/JM_Area)

elif row[3] == ′350211′:

row[2] = row[0]*(row[1]/JM_Area)

elif row[3] == ′350212′:

row[2] = row[0]*(row[1]/TA_Area)

elif row[3] == ′350213′:

row[2] = row[0]*(row[1]/XA_Area)

else:

print ″無符合條件″

cursor.updateRow(row)

print ″站點平均雨量計算執行完畢！″。

3 結 語

本文中采用GP服務，直接調用數據庫實時的數據，通過一整套的空間分析工具，結合腳本的計算語言，并在系統前端生成廈門市各行政區的泰森多邊形及各區的區域平均雨量，同時可自主選擇某一時段內的區域平均雨量，既可實時計算產生圖文成果，又可統計時段內的平均雨量值，大大提高了平均雨量的效率和精度，該方法在防汛指揮、日常雨情、水文分析等工作中都具有很大的突破。