多边形面积计算

手写地理信息组件系列 第7篇
多边形实体的面积计算
难度指数：★★★☆☆

“这个世界会好吗？”

第一次听到这一句，是在一个某已"行为不端"的民谣歌手的作品里听到的。究起源头，是出自于“中国最后一位大儒家”梁漱溟的父亲，梁济之口。

梁济所处清末，朝野无能，江河日下。梁济半生为国奔走呼号，奈何仍国之不国。一天，梁济将自己25岁的儿子梁漱溟叫到身边，问了儿子：

“这个世界会好吗”？ “我相信世界是一天一天往好里去的”。
“能好就好啊”。

几日，在梁济距自己生辰三日之时，为警醒世人，在北京的净业湖投湖自尽。
一句世纪之问，今日尤听，仍振聋发聩。
你说这个世界会好吗？

——读《梁漱溟传》

前情回顾

我们已经掌握了Shapefile的线面实体读取，了解线面shp是由文件头+记录头+记录内容组成的。其中记录内容是由Parts首坐标索引和Points点坐标对构成。依次的读取这些记录，就构成了一个个的空间实体对象。

当然我们不能止于文件读取和绘图，今天来走出GIS分析功能的第一步——多边形实体的面积计算。

多边形的构成

相信关注过以往内容的读者会知道，我们当初在定义点线面时，已经为线面对象预先定义了Vertex属性。因为它是构成多边形对象最重要的属性(环状线段也视为多边形)。

而Vertex被称为节点，是构成多边形对象的骨架，每个Vertex带有坐标，而面积就是从中计算而来。

值得注意的是，我们平时在纸上画多边形并不强调点的顺序，但这里Vertex的顺序关系一定是必要的，因为不同的节点顺序，可能构成不同的多边形，从而导致面积计算的不一致。说来可别不信，希望这张图能给你留下印象。

图中，AB两点在多边形节点中是两个位置完全相同的点，而在多边形环序中却是两个顺序相反的点。所以说顺序会决定形状，而形状可以决定面积。坐标点一定要按序给出。

一般来说，对于凹凸多边形，Shapefile都是按照顺时针顺序存储多边形节点。而对于内部带孔洞多边形，节点存储相对外层来说，顺序是完全相反的。

向量与叉乘

回顾高中数学必修4当中，第一次出现了向量。据此“史料”记载：向量是一个既有大小又有方向的量。

记得当初，用这本数学书上学到的向量数量积知识，还解决了不少物理题上木块移动的做功问题。。

OK，还是回到我们当前的问题。多边形的面积计算用向量来解决，是一种非常科学高效的方法。但并非使用高中数学的数量积(点乘)方法，而是大学数学中与线性代数相关的叉乘。

叉乘，又称叉积、向量积。在数值上等于以两个向量的模长为邻边，构成的平行四边形的面积。

向量点乘的结果是一个标量，而叉乘的结果仍然是一个向量。在三维空间中，其位置处在两不共线向量构成平面的法线之上，并满足右手定则。

结合图中，简单的说右手定则：手掌四指指向OA方向，并将四指弯曲至OB方向，此时大拇指的指向，即为OA和OB叉乘结果的方向。

观察图中结果向量OC，可知OA和OB的叉乘为负数，而OB和OA的叉乘为正数。两者的绝对值相等，都为OC的长度。

多边形面积计算

多边形面积计算可由三角形计算引出：

依照向量叉乘定义，图中平行四边形的面积即为叉乘结果，因为向量叉乘顺序为**a**×_b_，根据右手定则可知，结果为负数。三角形面积即为1/2平行四边形面积。

而多边形面积计算的思想，就是将多边形分割为多个三角形，然后进行加和：

值得注意的是，凸多边形（左）和凹多边形（右）不会因其凹凸性而影响面积的计算。右图凹进去的三角形（d、e为边）看似会增大多边形的面积，但实际上在加和中，会扣减这部分三角形的面积。你可以在上图看出，左边d×e为负数，右边d×_e_为正数。

嗯是的，看起来还合理。但是，如果像这样，遇到这种内部带有孔洞的多边形怎么办呢？

多边形存在“洞”时，需要知道一点。孔洞的环序与外环顺序是相反的，因此形成的面积会相减。所以带孔多边形的面积计算，只需将洞的环构成的面积代入计算即可。这里也再次体现了节点顺序的重要性。

以上，我们分析了利用叉积计算多边形面积的几种情况。可以发现，我们将多边形划分为多个三角形来计算面积的思想是可行的。但其中一个细节是，如何选择第一个分割点开始分割呢？

我们本能的会在头脑中想到多边形左上角的那个节点，因为这来源于平时的绘图习惯。而实际的答案是：

坐标系内的任意点

以任意点开始分割的推理，可以在文末点击阅读****原文。这里需要知道的是，以任意点作为起点分割多边形，依此来计算面积是完全成立的。由于叉积的方向，在多边形外面的部分会抵消，得到的仍是多边形的面积。通常，我们将原点视为分割的第一点，这是出于计算方便的考虑。

到这里，我们利用这个原点，基本可以整理出利用叉积计算多边形面积的等式：

其中：

既然公式有了，是不是马上就可以开始用代码算了呢？

且慢！这里存在一个问题！

这种计算方式，在数学上是完全成立的，但是理论需要结合实际。我们知道，计算机在计算小数的时候，是会有舍入误差的。拿当前这种以原点作为切割第一点的情况来说，如果图形距离原点非常远，那这种误差就会变得非常大。例如下面这种情况：

本初子午线和赤道的交点位于中非西部，它就是Web墨卡托投影坐标系的原点(0,0)。如果以这个点来计算中国境内的某个图形面积，在数学上是可行的，但在计算机中就会出现误差，而且还有可能发生溢出错误，所以不推荐以原点进行分割。而以图形上的节点作为分割的第一点是很合适的。

这样，向量的叉积就不会参照原点来计算，而是参照图形自身节点来计算，这样，无论图形距离原点多远，都不会造成误差。

实现叉积算法

以上我们分析了多边形面积计算的叉乘方法，现在结合代码实现多边形面积的计算。

在C#工程中单独定义一个类，算法类Algorithm，用于提供统一的算法相关函数。在类中定义一个函数SignedArea(有符号面积计算)，以节点列表为参数。

 public class Algorithm
 {
   public static double SignedArea(List<Vertex> vertexes)
   {
     double area = 0.0d;
 
     //第一点坐标（x0,y0）
     double x0 = vertexes[0].x;
     double y0 = vertexes[0].y;

    //从第二点开始遍历节点
    for (int i = 1; i < vertexes.Count - 1; i++)
    {
      double ax = vertexes[i].x;
      double ay = vertexes[i].y;

      double bx = vertexes[i + 1].x;
      double by = vertexes[i + 1].y;

      //向量a
      double va_x = ax - x0;
      double va_y = ay - y0;

      //向量b
      double vb_x = bx - x0;
      double vb_y = by - y0;

      //叉乘
      area += va_x * vb_y - vb_x * va_y;
    }
    return area / 2;
  }
}

在面实体中新增面积属性，面对象就可以计算自己的面积了。这里为了展示方便，在面绘制的方法中加入了面积数值的显示：

//面实体
 public class Polygon : Geometry
 {
   private List<Vertex> vertexes =
     new List<Vertex>();
 
   public List<Vertex> Vertexes
   {
     get { return vertexes; }
  }

  public double Area
  {
    get
    {
      //返回面积的绝对值
      return Math.Abs(
        Algorithm.SignedArea(vertexes));
    }
  }

  public override void Draw(Graphics graphics, Map map)
  {
    System.Drawing.Point[] points =
      new System.Drawing.Point[vertexes.Count];

    for (int i = 0; i < points.Length; i++)
    {
      points[i] = map.ToScreenPoint(vertexes[i]);
    }
    graphics.FillPolygon(
      new SolidBrush(Color.LightSeaGreen), points);

    //在界面绘制面积数值
    graphics.DrawString(Area.ToString(),
    new Font("宋体", 15),
    new SolidBrush(Color.Navy),
    new PointF(points[0].X, points[0].Y));
  }
}

面积计算

在Web墨卡托投影的底图上，我选定了天津北部的一块区域，手绘一个面shp，观察河岸的曲线可以看出这是一个凹多边形：

结合之前的Shpfile读取代码，现在面积可以计算并显示出来了：

面积：3085928.87246157

再来对比ArcMap的计算结果：

可以认为：我们的计算是准确的

但是需要注意的是，涉及到面积计算，一定要将数据转换为投影坐标系下进行。

我们的面积计算方法，本身是基于平面向量的，所以一定要将经纬度这种球面坐标转换为平面坐标后计算才有意义。同时，这也是ArcGIS在计算面积时，要求数据必须转换为投影坐标系的原因。

附：

如果你对文中，以“以原点为分割第一点会造成误差”的分析有更多兴趣，后台回复【误差】，查看对比结果。

看好关注，下期见。

上一篇：GIS底层｜线面shp读取中，那些容易被遗落的重要细节

---