手寫地理信息組件系列 第7篇
多邊形實(shí)體的面積計(jì)算
難度指數(shù)：★★★☆☆

“這個(gè)世界會(huì)好嗎？”

第一次聽到這一句，是在一個(gè)某已"行為不端"的民謠歌手的作品里聽到的。究起源頭，是出自于“中國(guó)最后一位大儒家”梁漱溟的父親，梁濟(jì)之口。

梁濟(jì)所處清末，朝野無(wú)能，江河日下。梁濟(jì)半生為國(guó)奔走呼號(hào)，奈何仍國(guó)之不國(guó)。一天，梁濟(jì)將自己25歲的兒子梁漱溟叫到身邊，問了兒子：

“這個(gè)世界會(huì)好嗎”？ “我相信世界是一天一天往好里去的”。
“能好就好啊”。

幾日，在梁濟(jì)距自己生辰三日之時(shí)，為警醒世人，在北京的凈業(yè)湖投湖自盡。
一句世紀(jì)之問，今日尤聽，仍振聾發(fā)聵。
你說這個(gè)世界會(huì)好嗎？

——讀《梁漱溟傳》

前情回顧

我們已經(jīng)掌握了Shapefile的線面實(shí)體讀取，了解線面shp是由文件頭+記錄頭+記錄內(nèi)容組成的。其中記錄內(nèi)容是由Parts首坐標(biāo)索引和Points點(diǎn)坐標(biāo)對(duì)構(gòu)成。依次的讀取這些記錄，就構(gòu)成了一個(gè)個(gè)的空間實(shí)體對(duì)象。

當(dāng)然我們不能止于文件讀取和繪圖，今天來(lái)走出GIS分析功能的第一步——多邊形實(shí)體的面積計(jì)算。

多邊形的構(gòu)成

相信關(guān)注過以往內(nèi)容的讀者會(huì)知道，我們當(dāng)初在定義點(diǎn)線面時(shí)，已經(jīng)為線面對(duì)象預(yù)先定義了Vertex屬性。因?yàn)樗菢?gòu)成多邊形對(duì)象最重要的屬性(環(huán)狀線段也視為多邊形)。

而Vertex被稱為節(jié)點(diǎn)，是構(gòu)成多邊形對(duì)象的骨架，每個(gè)Vertex帶有坐標(biāo)，而面積就是從中計(jì)算而來(lái)。

值得注意的是，我們平時(shí)在紙上畫多邊形并不強(qiáng)調(diào)點(diǎn)的順序，但這里Vertex的順序關(guān)系一定是必要的，因?yàn)椴煌墓?jié)點(diǎn)順序，可能構(gòu)成不同的多邊形，從而導(dǎo)致面積計(jì)算的不一致。說來(lái)可別不信，希望這張圖能給你留下印象。

圖中，AB兩點(diǎn)在多邊形節(jié)點(diǎn)中是兩個(gè)位置完全相同的點(diǎn)，而在多邊形環(huán)序中卻是兩個(gè)順序相反的點(diǎn)。所以說順序會(huì)決定形狀，而形狀可以決定面積。坐標(biāo)點(diǎn)一定要按序給出。

一般來(lái)說，對(duì)于凹凸多邊形，Shapefile都是按照順時(shí)針順序存儲(chǔ)多邊形節(jié)點(diǎn)。而對(duì)于內(nèi)部帶孔洞多邊形，節(jié)點(diǎn)存儲(chǔ)相對(duì)外層來(lái)說，順序是完全相反的。

向量與叉乘

回顧高中數(shù)學(xué)必修4當(dāng)中，第一次出現(xiàn)了向量。據(jù)此“史料”記載：向量是一個(gè)既有大小又有方向的量。

記得當(dāng)初，用這本數(shù)學(xué)書上學(xué)到的向量數(shù)量積知識(shí)，還解決了不少物理題上木塊移動(dòng)的做功問題。。

OK，還是回到我們當(dāng)前的問題。多邊形的面積計(jì)算用向量來(lái)解決，是一種非常科學(xué)高效的方法。但并非使用高中數(shù)學(xué)的數(shù)量積(點(diǎn)乘)方法，而是大學(xué)數(shù)學(xué)中與線性代數(shù)相關(guān)的叉乘。

叉乘，又稱叉積、向量積。在數(shù)值上等于以兩個(gè)向量的模長(zhǎng)為鄰邊，構(gòu)成的平行四邊形的面積。

向量點(diǎn)乘的結(jié)果是一個(gè)標(biāo)量，而叉乘的結(jié)果仍然是一個(gè)向量。在三維空間中，其位置處在兩不共線向量構(gòu)成平面的法線之上，并滿足右手定則。

結(jié)合圖中，簡(jiǎn)單的說右手定則：手掌四指指向OA方向，并將四指彎曲至OB方向，此時(shí)大拇指的指向，即為OA和OB叉乘結(jié)果的方向。

觀察圖中結(jié)果向量OC，可知OA和OB的叉乘為負(fù)數(shù)，而OB和OA的叉乘為正數(shù)。兩者的絕對(duì)值相等，都為OC的長(zhǎng)度。

多邊形面積計(jì)算

多邊形面積計(jì)算可由三角形計(jì)算引出：

依照向量叉乘定義，圖中平行四邊形的面積即為叉乘結(jié)果，因?yàn)橄蛄坎娉隧樞驗(yàn)?*a**×_b_，根據(jù)右手定則可知，結(jié)果為負(fù)數(shù)。三角形面積即為1/2平行四邊形面積。

而多邊形面積計(jì)算的思想，就是將多邊形分割為多個(gè)三角形，然后進(jìn)行加和：

值得注意的是，凸多邊形（左）和凹多邊形（右）不會(huì)因其凹凸性而影響面積的計(jì)算。右圖凹進(jìn)去的三角形（d、e為邊）看似會(huì)增大多邊形的面積，但實(shí)際上在加和中，會(huì)扣減這部分三角形的面積。你可以在上圖看出，左邊d×e為負(fù)數(shù)，右邊d×_e_為正數(shù)。

嗯是的，看起來(lái)還合理。但是，如果像這樣，遇到這種內(nèi)部帶有孔洞的多邊形怎么辦呢？

多邊形存在“洞”時(shí)，需要知道一點(diǎn)。孔洞的環(huán)序與外環(huán)順序是相反的，因此形成的面積會(huì)相減。所以帶孔多邊形的面積計(jì)算，只需將洞的環(huán)構(gòu)成的面積代入計(jì)算即可。這里也再次體現(xiàn)了節(jié)點(diǎn)順序的重要性。

以上，我們分析了利用叉積計(jì)算多邊形面積的幾種情況。可以發(fā)現(xiàn)，我們將多邊形劃分為多個(gè)三角形來(lái)計(jì)算面積的思想是可行的。但其中一個(gè)細(xì)節(jié)是，如何選擇第一個(gè)分割點(diǎn)開始分割呢？

我們本能的會(huì)在頭腦中想到多邊形左上角的那個(gè)節(jié)點(diǎn)，因?yàn)檫@來(lái)源于平時(shí)的繪圖習(xí)慣。而實(shí)際的答案是：

坐標(biāo)系內(nèi)的任意點(diǎn)

以任意點(diǎn)開始分割的推理，可以在文末點(diǎn)擊閱讀****原文。這里需要知道的是，以任意點(diǎn)作為起點(diǎn)分割多邊形，依此來(lái)計(jì)算面積是完全成立的。由于叉積的方向，在多邊形外面的部分會(huì)抵消，得到的仍是多邊形的面積。通常，我們將原點(diǎn)視為分割的第一點(diǎn)，這是出于計(jì)算方便的考慮。

到這里，我們利用這個(gè)原點(diǎn)，基本可以整理出利用叉積計(jì)算多邊形面積的等式：

其中：

既然公式有了，是不是馬上就可以開始用代碼算了呢？

且慢！這里存在一個(gè)問題！

這種計(jì)算方式，在數(shù)學(xué)上是完全成立的，但是理論需要結(jié)合實(shí)際。我們知道，計(jì)算機(jī)在計(jì)算小數(shù)的時(shí)候，是會(huì)有舍入誤差的。拿當(dāng)前這種以原點(diǎn)作為切割第一點(diǎn)的情況來(lái)說，如果圖形距離原點(diǎn)非常遠(yuǎn)，那這種誤差就會(huì)變得非常大。例如下面這種情況：

本初子午線和赤道的交點(diǎn)位于中非西部，它就是Web墨卡托投影坐標(biāo)系的原點(diǎn)(0,0)。如果以這個(gè)點(diǎn)來(lái)計(jì)算中國(guó)境內(nèi)的某個(gè)圖形面積，在數(shù)學(xué)上是可行的，但在計(jì)算機(jī)中就會(huì)出現(xiàn)誤差，而且還有可能發(fā)生溢出錯(cuò)誤，所以不推薦以原點(diǎn)進(jìn)行分割。而以圖形上的節(jié)點(diǎn)作為分割的第一點(diǎn)是很合適的。

這樣，向量的叉積就不會(huì)參照原點(diǎn)來(lái)計(jì)算，而是參照?qǐng)D形自身節(jié)點(diǎn)來(lái)計(jì)算，這樣，無(wú)論圖形距離原點(diǎn)多遠(yuǎn)，都不會(huì)造成誤差。

實(shí)現(xiàn)叉積算法

以上我們分析了多邊形面積計(jì)算的叉乘方法，現(xiàn)在結(jié)合代碼實(shí)現(xiàn)多邊形面積的計(jì)算。

在C#工程中單獨(dú)定義一個(gè)類，算法類Algorithm，用于提供統(tǒng)一的算法相關(guān)函數(shù)。在類中定義一個(gè)函數(shù)SignedArea(有符號(hào)面積計(jì)算)，以節(jié)點(diǎn)列表為參數(shù)。

public class Algorithm{public static double SignedArea(List<Vertex> vertexes){double area = 0.0d;//第一點(diǎn)坐標(biāo)（x0,y0）double x0 = vertexes[0].x;double y0 = vertexes[0].y;//從第二點(diǎn)開始遍歷節(jié)點(diǎn)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;//向量adouble va_x = ax - x0;double va_y = ay - y0;//向量bdouble vb_x = bx - x0;double vb_y = by - y0;//叉乘area += va_x * vb_y - vb_x * va_y;}return area / 2;} }

在面實(shí)體中新增面積屬性，面對(duì)象就可以計(jì)算自己的面積了。這里為了展示方便，在面繪制的方法中加入了面積數(shù)值的顯示：

//面實(shí)體public class Polygon : Geometry{private List<Vertex> vertexes =new List<Vertex>();public List<Vertex> Vertexes{get { return vertexes; }}public double Area{get{//返回面積的絕對(duì)值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);//在界面繪制面積數(shù)值graphics.DrawString(Area.ToString(),new Font("宋體", 15),new SolidBrush(Color.Navy),new PointF(points[0].X, points[0].Y));} }

面積計(jì)算

在Web墨卡托投影的底圖上，我選定了天津北部的一塊區(qū)域，手繪一個(gè)面shp，觀察河岸的曲線可以看出這是一個(gè)凹多邊形：

結(jié)合之前的Shpfile讀取代碼，現(xiàn)在面積可以計(jì)算并顯示出來(lái)了：

面積：3085928.87246157

再來(lái)對(duì)比ArcMap的計(jì)算結(jié)果：

可以認(rèn)為：我們的計(jì)算是準(zhǔn)確的

但是需要注意的是，涉及到面積計(jì)算，一定要將數(shù)據(jù)轉(zhuǎn)換為投影坐標(biāo)系下進(jìn)行。

我們的面積計(jì)算方法，本身是基于平面向量的，所以一定要將經(jīng)緯度這種球面坐標(biāo)轉(zhuǎn)換為平面坐標(biāo)后計(jì)算才有意義。同時(shí)，這也是ArcGIS在計(jì)算面積時(shí)，要求數(shù)據(jù)必須轉(zhuǎn)換為投影坐標(biāo)系的原因。

附：

如果你對(duì)文中，以“以原點(diǎn)為分割第一點(diǎn)會(huì)造成誤差”的分析有更多興趣，后臺(tái)回復(fù)【誤差】，查看對(duì)比結(jié)果。

看好關(guān)注，下期見。

上一篇：GIS底層｜線面shp讀取中，那些容易被遺落的重要細(xì)節(jié)

---

總結(jié)

以上是生活随笔為你收集整理的多边形面积计算的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。