arcpy泰森多边形法计算面雨量工具

   在水利部门或气象部门中面平均降水量是降雨中很重要的指标，传统计算多用使用算术平均法，泰森多边形法和等值线法，后两种计算方法在传统的计算中很难计算，但使用用GIS十分方便计算，可以是任意区域的任意一场降雨。

一、泰森多边形介绍：

荷兰气候学家A·H·Thiessen提出了一种根据离散分布的气象站的降雨量来计算平均降雨量的方法，即将所有相邻气象站连成三角形，作这些三角形各边的垂直平分线，于是每个气象站周围的若干垂直平分线便围成一个多边形。用这个多边形内所包含的一个唯一气象站的降雨强度来表示这个多边形区域内的降雨强度，并称这个多边形为泰森多边形。

特点：

1、每个泰森多边形内仅含有一个离散点数据；

2、泰森多边形内的点到相应离散点的距离最近；

3、位于泰森多边形边上的点到其两边的离散点的距离相等。

二、在ArcGIS中利用泰森多边形法分析流域降雨量步骤（长江上游为例）没有使用工具：

1、首先在ArcMap中加长江流域分区和雨量站点的.shp格式的数据（必须是.shp格式的）。若雨量站点信息为.xls或.txt格式的，则应该将其加载后先转成.shp格式，再进行以下的操作。加载数据结果如图：

2、在ArcToolbox工具中选择Analysis Tools—Proximity—Create Thiessen Polygons工具，打开Create Thiessen Polygons窗口，在Input Features中输入站点数据:雨量站点，在Output Features Class中设置输出路径，在Output Fields（optional）中选择All（即输出所有属性字段）。如图所示：

  然后设置其环境，即选择Create Thiessen Polygons窗口下面的Environments…按钮，进入环境设置窗口

选择General Settings进行设置：主要设置包括两项，第一项对Output Coordinate System设置，选取流域面矢量数据以和其保持一致的坐标系，此处选择Same as Layer “流域分区图”；第二项对Extent进行设置，设置生成泰森多边形的四周边缘，此处选择Same as Layer 流域分区图，其余保持默认。如图所示：

选择OK后，Create Thiessen Polygons窗口再次选择OK，则生成的泰森多边形如图所示：

3、因为需要对江流域范围内的每一块多边形的面积进行统计，所以对新生成的泰森多边形数据进行裁剪，由于是矢量数据，所以采用ArcToolbox工具中的Analysis Tools—Extract—Clip工具进行裁剪，设置如图：

点击OK，生成的Clip后的图如下图所示：

为了统计每一多边形的面积，打开Clip后的数据的属性表，并添加属性数据：Area计算面积，在excel中计算出面平均降水量。

三、以上是熟悉GIS工具的专业人员可以很快计算出区域泰森多边形面平均降水量，量对于水利或气象专业的人员不一定熟悉ARCGIS，以下是本人通过arcpy编写的小程序，只要有降水量站点及降水值和区域面积图就可以直接计算出面平均降水量。工具如下，选择降水量站点图，带有降水量数据字段。

1、添加降水量图层，选择年降水量字段。

2、添加区域图层，如果流域中有多个行政区选择行政区名称字段

3、点确定之后就计算出区域平均降水量如下表，流域内每个县年平均降水量计算完成

4、arcpy泰森多边形法计算面雨量工具界面如下

 

 

5、arcpy代码如下

# -*- coding: utf-8 -*-
import arcpy

Rian_stion = arcpy.GetParameterAsText(0)

Xzfq = arcpy.GetParameterAsText(1)

Outinter_Thiessen = arcpy.GetParameterAsText(2)

filedname = arcpy.GetParameterAsText(3)

outtable = arcpy.GetParameterAsText(4)

summation = arcpy.GetParameterAsText(5)

outshp = ""
Outclp_Thiessen = ""
Thiessen = Outinter_Thiessen
outshp = Outinter_Thiessen
xzq_areas = Thiessen

# Process: 创建泰森多边形
arcpy.CreateThiessenPolygons_analysis(Rian_stion, outshp, "ALL")

# Process: 裁剪
arcpy.Clip_analysis(outshp, Xzfq, Outclp_Thiessen, "1 Meters")

# Process: 相交
arcpy.Intersect_analysis("# #;# #", Outinter_Thiessen, "ALL", "1 Meters", "INPUT")

# Process: 添加字段
arcpy.AddField_management(Outinter_Thiessen, "Thiessen_Area", "FLOAT", "", "", "", "泰森多边形面积", "NULLABLE", "NON_REQUIRED", "")

# Process: 计算字段
arcpy.CalculateField_management(Thiessen, "Thiessen_Area", "!Shape.area!/1000000", "PYTHON", "")

# Process: 添加字段 (2)
arcpy.AddField_management(Outinter_Thiessen, "Pa", "FLOAT", "", "", "", "降水量", "NULLABLE", "NON_REQUIRED", "")

# Process: 计算字段 (2)
arcpy.CalculateField_management(outshp2, "Pa", "float(!字段名!)*!Thiessen_Area!", "PYTHON", "")

# Process: 汇总统计数据
arcpy.Statistics_analysis(outshp3, outtable, "Pa SUM", summation)