【ArcGIS】基于泰森多边形求流域面降水量

泰森多边形（Thiessen Polygon）法

泰森多边形又叫冯洛诺伊图（Voronoi diagram），得名于Georgy Voronoi，是一组由连接两邻点线段的垂直平分线组成的连续多边形。一个泰森多边形内的任一点到构成该多边形的控制点的距离小于到其他多边形控制点的距离。

1 原理

在开展流域水文分析时，往往流域范围仅有少量雨量站，为较为准确地将雨量站点实际观测的降雨量反应到区域或子流域上，需要进行数据的邻近分析，而泰森多边形方法被广泛应用与计算流域平均面雨量。
流域上各点的雨量用离该点最近雨量站的降雨量代表。用泰森多边形法计算流域的平均降雨量，是以各雨量站之间连线的垂直平分线，把流域划分为若干个多边形，然后以各个多边形的面积为权数，计算各站雨量的加权平均值，并把它作为流域的平均降雨量，一般来说结果比单纯算术平均法更为精确。

泰森多边形法特点：

• 每个泰森多边形内仅含有一个离散点数据
• 泰森多边形内的点到相应离散点的距离最近
• 位于泰森多边形边上的点到其两边的离散点的距离相等

面雨量计算方法：

• a.先计算每个泰森多形内的平均雨量，就是该多边形内的雨量站雨量乘以权重，权重即该多边形面积值除以流域面积。
• b.把所有多边形内的平均雨量相加，再除以多边形个数即是结果。

2 操作步骤

1.首先，在ArcMap中导入站点.shp（含降雨量数据）和流域边界.shp。需要注意的是，二者均为shp矢量数据，如下图所示。

2.接下来，利用ArcToolbox工具中的Create Thiessen Polygons工具（Toolboxes → Analysis Tools → Proximity → Create ThiessenPolygons）进行泰森多边形的创建。

点击Create Thiessen Polygons弹出窗口后所做的参数设置如下图所示，在Input Features中输入RainfallStation数据，在Output Features Class中自定义输出路径（一般选择之前已经定义好的默认路径），在Output Fields (optional)中选择All（即输出所有属性字段）。

【另】报错
ERROR：Input RainfallStation does not hace OIDs.
需要进行一步操作，先给降水数据表添加object_ID字段！！！

实现方法1： 添加新字段

操作步骤如下：在属性表中，如果没有 Object ID 字段，可以添加一个新字段。

• 右键单击内容列表中的表或图层，然后选择打开属性表。

• 单击表窗口中的表选项按钮 表选项。
  即使您未处于编辑会话中，也可以进行计算；但在这种情况下无法撤消计算结果。

• 单击添加字段。

• 输入字段的名称。

• 单击类型箭头，然后单击字段类型。

• 根据需要设置任何其他字段属性。

• 单击确定。

实现方法2： 将数据导出为.shp文件，会自动添加Object ID 字段。

3.设置Create Thiessen Polygons里的Environments

输入界面设置完成后进行环境变量设置，选择Create Thiessen Polygons窗口下面的【Environments】按钮，进入环境设置窗口，设置Output Coordinate System，选择Same as Input，也可以选择与untitled_poly保持一致的坐标系，不过选择与与untitled_poly保持一致的坐标系很有可能生不成泰森多边形，因此可以先选择与输入一致的坐标系，后面需要修改坐标系的时候再修改一下就可以了。
然后对【Processing Extent】进行设置，设置生成泰森多边形的四周边界，此处选择Same as Layer untitled_poly，其余保持默认。如图所示：

以上需要设置的地方都设置完成后，点击OK，在Create Thiessen Polygons窗口再点击OK，则生成的泰森多边形如图所示：

4.现在生成的泰森多边形是一个将untitled_poly流域包含在内的大四边形，不能直接用于untitled_poly流域的面雨量计算，因此需要按照untitled_poly流域的形状对新生成的泰森多边形数据进行裁剪。
采用ArcToolbox工具中的Clip工具（Toolboxes → Analysis Tools → Extract → Clip）进行裁剪，在弹出的窗口中如下图进行设置：单位选择Meters

点击OK后，裁剪后得到的图如下图所示：

5.接下来计算裁剪后生成的每一个多边形的面积，打开裁剪后的Rainfall_Station_CreateThies_Clip数据的属性表，并添加Area字段，并计算面积。
添加面积字段方法：单击Table Options → Add Field…，添加面积字段后，右击Area字段，点击Calculate Geometry…，然后按雨量站名称类别显示如下图所示：

各子区域面积如下：

【另】出错

直接设置图框为投影坐标系即可，此处设置为WGS_UTM_Zone_46N。

【另】出错长精度
可选择单位为【m】或【km】。

6.然后将属性表中所有数据全部选中，右击如图所示位置属性列表条件field——几何计算——导出dbf。

点击Copy Selected，粘贴到excel表格中进行面积权重的计算，最后根据各雨量站点所测的降雨量进行加权平均，就可以计算出untitled_poly流域的降雨量了。计算结果如图所示：

另：利用MATLAB计算面降水量

clc
close all
clear
load('P.mat')
load('areaRatio.mat')
% 求各站点年平均降水量
% ------------------------------------------------------------------
% 按年划分降水
PAnnual = cell(nStation, 1);                          % zeros(nyear,365);
PAnnually = zeros(nStation,nYear);
PAveAnnual = zeros(1,nStation);
% 按月划分降水
PMonth = cell(nStation,2);                 % 元胞数组（月降水数据）
PAveMonth = zeros(nStation,12);           % 元胞数组（月平均降水）
for in=1:nStation
     PAnnual{in,1} = year_data_simple( yearStart, yearEnd, P(:,in));
     PAnnually(in,:) = sum(PAnnual{in,1},2)';
     PAveAnnual(1,in) = mean( sum( PAnnual{in,1}, 2) );          % 多年平均降水
     
     PMonth{in,1} = month_data_simple( yearStart, yearEnd, P(:,in));
     PMonth{in,2} = sum( PMonth{in,1}, 2 );
     
     tempP = reshape( PMonth{in,2}, 12, nYear);
     PAveMonth(in,:) = mean(tempP,2)';
end
% 计算面降水量：各月降水占比和年降水量变化
PArea = sum (PAveMonth.*areaRatio, 1);
PAreaYear = sum (PAnnually.*areaRatio, 1);

成图如下所示：

参考

1.CSDN博客-泰森多边形的matlab实现
2.泰森多边形计算流域面雨量