基于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服务，直接调用数据库实时的数据，通过一整套的空间分析工具，结合脚本的计算语言，并在系统前端生成厦门市各行政区的泰森多边形及各区的区域平均雨量，同时可自主选择某一时段内的区域平均雨量，既可实时计算产生图文成果，又可统计时段内的平均雨量值，大大提高了平均雨量的效率和精度，该方法在防汛指挥、日常雨情、水文分析等工作中都具有很大的突破。