肖惠珍

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

0 引 言

通常雨量站所觀測(cè)到的降雨量，只能代表該雨量站周圍較小范圍的降雨情況，不能代表全流域或某一區(qū)域的平均降雨量，因此單獨(dú)某個(gè)雨量站的降雨數(shù)據(jù)不能作為洪澇預(yù)報(bào)與評(píng)估的依據(jù)。因此，需要采用全流域或某一區(qū)域的所有雨量站的降雨數(shù)據(jù)來計(jì)算區(qū)域平均降雨量。

目前，計(jì)算區(qū)域平均降雨量的方法很多，常用的有算術(shù)平均法、數(shù)值法、等值線法、泰森多邊形法等。在這些方法中，泰森多邊形法最適合流域或區(qū)域內(nèi)雨量站或降雨量分布不均勻的情況，是我國(guó)計(jì)算區(qū)域平均雨量常用的方法，被水利、氣象、環(huán)境等部門廣泛應(yīng)用。根據(jù)廈門的區(qū)域特點(diǎn)，雨量觀測(cè)站分散，降雨量時(shí)空分布不均，采用泰森多邊形法計(jì)算區(qū)域平均雨量非常適合。

1 研究方法

1.1 構(gòu)建Delaunay三角網(wǎng)

基于ArcGIS平臺(tái)的泰森多邊形構(gòu)建，是通過雨量監(jiān)測(cè)點(diǎn)的雨量值計(jì)算區(qū)域平均雨量的方法，來實(shí)現(xiàn)離散點(diǎn)構(gòu)建泰森多邊形，從而獲得區(qū)域的連續(xù)值。

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

1.2 泰森多邊形的裁剪

根據(jù)廈門市所有雨量監(jiān)測(cè)站點(diǎn)構(gòu)建出的泰森多邊形，以廈門市的最大經(jīng)緯度和最小經(jīng)緯度為邊界進(jìn)行構(gòu)建泰森多邊形，因此計(jì)算出來的邊界為長(zhǎng)方形，需要對(duì)泰森多邊形進(jìn)行裁剪。能否精準(zhǔn)求出裁剪結(jié)果，是保證計(jì)算結(jié)果準(zhǔn)確率的關(guān)鍵。利用ArcGIS的相交工具，將廈門市的行政區(qū)劃與泰森多邊形進(jìn)行疊加，將區(qū)劃屬性附給泰森多邊形，從而進(jìn)行行政區(qū)劃內(nèi)泰森多邊形的面積和區(qū)域總面積的計(jì)算。

1.3 區(qū)域平均雨量計(jì)算方法

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

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

設(shè)每個(gè)雨量觀測(cè)點(diǎn)的降雨量為xi，每個(gè)對(duì)應(yīng)的泰森多邊形的面積為fx，則區(qū)域平均降雨量可按下式求得：

式中：xi為雨量觀測(cè)點(diǎn)的降雨量，mm；fx為泰森多邊形的面積，km2；n為區(qū)域內(nèi)雨量觀測(cè)點(diǎn)或泰森多邊形的個(gè)數(shù)；F為區(qū)域的總面積，km2；Ai為雨量站權(quán)重系數(shù)。

2 基于GP服務(wù)的廈門市區(qū)域平均雨量

GP服務(wù)是ArcGIS一個(gè)自帶的特殊功能，通過Model Builder，通過拖動(dòng)數(shù)據(jù)、工具，設(shè)置配置參數(shù)等操作，將整個(gè)運(yùn)算過程建成一個(gè)完整的模型。同時(shí)可添加腳本指令語言，最后將模型發(fā)布成GP服務(wù)，并供系統(tǒng)調(diào)用，形成前后交互，實(shí)現(xiàn)實(shí)時(shí)、快速的運(yùn)算模式。

2.1 構(gòu)建模型

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

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

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

2.2 計(jì)算站點(diǎn)平均降雨量Python代碼

Python對(duì)縮進(jìn)要求嚴(yán)格，在sublime中縮進(jìn)用tab，用空格會(huì)報(bào)錯(cuò)，遇到錯(cuò)誤可以去搜下有相關(guān)資料說明，主要用到arcpy這個(gè)類；

以下是代碼：

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

# 需加編碼，不然print中文字符會(huì)報(bào)錯(cuò)

import arcpy

import time

import sys

reload(sys)

sys.setdefaultencoding( ″gbk″ )

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

arcpy.env.overwriteOutput = True

addFieldName = ″AVG_RAIN″

#獲取參數(shù)方法

inFeature = arcpy.GetParameterAsText(0)

outFeature = arcpy.GetParameterAsText(1)

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

#inFeature = ″calculateArea.shp″

#outFeature = ″copy.shp″

#復(fù)制輸入要素，后面進(jìn)行編輯

arcpy.CopyFeatures_management(inFeature,outFeature)

print ″要素復(fù)制成功！″

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

#計(jì)算廈門各個(gè)行政區(qū)的面積

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 ″海滄區(qū)的面積為：{0},集美區(qū)的面積為：{1},同安區(qū)的面積為：{2},翔安區(qū)的面積為：{3},湖里區(qū)的面積為：{4},思明區(qū)的面積為：{5}″.format(HC_Area,JM_Area,TA_Area,XA_Area,HL_Area,SM_Area)

#把需要用到的字段加到下面數(shù)組當(dāng)中，用游標(biāo)去查詢更新

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 ″站點(diǎn)平均雨量計(jì)算執(zhí)行完畢！″。

3 結(jié) 語

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