那些驚艷的 GIS 輪子

2019 年 10 月 29 日
筆記

一、前言

GIS 涉及測繪、幾何拓撲、人文社科等多方面的科學知識。在 .Net 平台下有着許多優秀的開源產品，比如：MapWindow、SharpMap、WorldWind等。而在這其中，CoordinateSharp與NetTopologySuite是兩款極其令人驚艷的中間開發組件產品。直到最近，我才遇到它們。

真的懊惱早沒有人告訴我這些優秀的作品的存在。此前都一直在調用 c/c++的接口，雖說其效率很高，但最終產品還是 .Net 桌面的產品，其間各種相互調用後誰也不能保證效率的優勢所在。並且出了問題還得反饋到底層開發組去重新修改編譯發佈一番。更別說調用 Python 的shapely、geopandas，或者 Java 的JTS Topology Suite、GeoTools等。正如聰明的讀者想到的那樣，可以將業務服務架構在這些優秀的產品之上。為此，有很長一段時間我都在研究wcf、asp.net core、Django、aspnet-microservices和docker等。的確也出了一些效果和性能均令人滿意的服務。但被告知後台業務將由 Java 組的人接手。於是，又開始了 Java 的研究，spring boot、spark、hbase等等。也寫了一些 Java 的服務端業務。但仍然避免不了高速實時數據處理，並且面向不同終端用戶要計算不同需求的問題。最終還是會有一些定製化的業務留在了桌面端。這就像有了雲計算後，還需要霧計算、邊緣計算作為有益的補充。不可避免的，還得使用 .Net 的實現。

以上都是我用過的各個平台上的優秀產品，沒有厚此薄彼的意思。這些也僅僅是為了具體的業務解決問題。下面特別地介紹一下CoordinateSharp與NetTopologySuite。二者皆是可以跨平台的 .net core 產品。

二、CoordinateSharp

CoordinateSharp 是一個簡單易用的進行地理坐標轉換、空間天體計算的產品庫。其強大與便捷之處我將以幾個簡單的小列子進行展示，僅拋磚引玉。

1.地理坐標轉換

# 北京天安門廣場的經緯度  CoordinateSharp.Coordinate.TryParse("N 39° 54' 27" E 116° 23' 17"",new DateTime(2019, 10, 1), out var c);  Console.WriteLine($"{c.Latitude.ToDouble()},{c.Longitude.ToDouble()}");//轉換結果：39.9075,116.38805555555555

這裡有一點令人疑惑的地方就是：為什麼會有時間信息。這正是它的獨到之處，不僅僅進行坐標轉換，還帶有計算日、月升落時間，位置等天體信息的能力：

Console.WriteLine($"{c.CelestialInfo.SunRise},{c.CelestialInfo.SunSet}");  # 10/1/2019 10:12:00 PM,10/1/2019 10:00:08 AM

由於時差原因，我們得加上 8 小時（東八區比格林尼治早 8 小時），於是結果變為10/2/2019 6:12:00 AM,10/1/2019 06:00:08 PM，日出時間變為第二天的早上了。日出減去 24 小時後為10/1/2019 6:12:00 AM。而日落仍然為10/1/2019 06:00:08 PM。查閱網上實時的信息：

日期	日出	日中	日落	晝長	天亮	天黑
2019 年 10 月 01 日周二	06:10:11	12:04:10	17:58:08	11:47:57	05:43:19	18:25:00

基本一致，誤差的原因由地球是橢球體、自轉公轉速率緩慢改變等引起。這個差別能夠讓人接受。

2.空間計算

計算球面距離

// 北京首都機場經緯度  var c1 = new CoordinateSharp.Coordinate(40.0760979329, 116.5953477768);  // 上海虹橋機場經緯度  var c2 = new CoordinateSharp.Coordinate(31.1982791377, 121.3356526703);  var d0 = new CoordinateSharp.Distance(c1, c2);  Console.WriteLine(d0.Kilometers);// 1074.0736250617888  Console.WriteLine(c1.ECEF);// 地球中心地球固定坐標（世界坐標系） -2188.311 km, 4369.881 km, 4084.559 km

而我們在網上查到的信息如下（注意，飛機真實飛行路線並非都是直線）：

上海到北京的空中距離為 1084 公里，上海虹橋機場到首都機場飛行距離為 1160 公里

這與計算結果 1074.0736250617888 接近（沒有辦法排除各各家數據差異對距離計算結果產生的影響，比如，谷歌、百度、騰訊高德的坐標都不同，上述兩個機場的位置坐標就是網上谷歌地球上取得的）。

此外，Coordinate 類有一個極其有用的方法void Move(double distance, double bearing, Shape shape)。其作用是計算往某個方向移動某個距離後的新坐標。

// 在橢球上向正北（方位角bearing正北為0）移動十公里,與c1原值比起來，經度基本沒有變化  c1.Move(10000,0,Shape.Ellipsoid);  Console.WriteLine($"{c1.Latitude.ToDouble()},{c1.Longitude.ToDouble()}");// 40.16739008225999,116.60039000000003

由此看來，地理空間能力基本上解決了距離與位置的相互轉換。更多功能歡迎探索https://github.com/Tronald/CoordinateSharp。去 star。

三、NetTopologySuite

NetTopologySuite 是一個快速可靠的 .Net GIS 解決方案。它提供了JTS Topology Suite所有功能的直接接口。JTS 是用於建模和平面幾何計算，並且遵循Open GIS的 SQL 簡單特性規範。而 NTS 基本上擁有這些能力，並且microsoft用其來為EF Core（This feature was added in EF Core 2.2.）提供空間計算能力。詳情可以參看Spatial Data。可以說 JTS 是幾何拓撲工具的 Java 實現，而 NTS 是 .NET 下的實現。

1.WKT 讀寫

var wkt = new WKTReader();  var gm = wkt.Read("POLYGON ((0 0,100 0, 100 100,0 100,0 0))");  //邊長為100的正方形  Console.WriteLine(new WKTWriter().Write(gm));// POLYGON ((0 0, 100 0, 100 100, 0 100, 0 0))

2.幾何圖形計算

// 幾何創建工廠（也可不使用工廠模式直接創建幾何圖形）  var gf = new GeometryFactory();  var pg3 = gf.CreatePolygon(new[]  {      new Coordinate(612, 612),      new Coordinate(144, 355),      new Coordinate(165, 188),      new Coordinate(277, 328),      new Coordinate(612, 612)  });  var pg4 = gf.CreatePolygon(new[]  {      new Coordinate(412, 612),      new Coordinate(555, 455),      new Coordinate(655, 188),      new Coordinate(577, 328),      new Coordinate(412, 612)  });

綠色是 pg3，金色是 pg4

求並集：var union = pg3.Union(pg4);

求交集：var intersection = pg3.Intersection(pg4);

求差集：var difference = pg3.Difference(pg4);
中間的連線是繪製時導致的，但是計算結果正確。我們查看其 WKT 描述為：MULTIPOLYGON (((478.68239179148486 538.78926215899912, 612 612, 499.14366946688551 516.32478247341942, 478.68239179148486 538.78926215899912)), ((478 498.4, 277 328, 165 188, 144 355, 460.37522887113056 528.73596970059953, 478 498.4)))的確是兩個多邊形。

其他能力，比如計算幾何間距、面積、凸包、判斷是否相交、是否覆蓋等可查看項目的介紹或者 test 例子。詳情訪問https://github.com/NetTopologySuite/NetTopologySuite。去 star。

四、CoordinateSharp 與 NetTopologySuite

這麼多優秀產品為何單獨介紹 CoordinateSharp 與 NetTopologySuite？

假設有需求為畫出某一城區的緩衝區，間距為 10km。這可怎麼辦？

在 NetTopologySuite 中直接提供緩衝區的計算函數public Geometry Buffer(double distance),效果也非常好：

紅色為原始直線，按彩虹色依次相距為 10 的緩衝區

但是我們卻忽略了一個重要的事情，NetTopologySuite 的計算的距離為平面坐標系下的歐氏距離。二維平面歐式距離的計算為Math.Sqrt(Math.Pow(x1-x2,2),Math.Pow(y1-y2,2)),直接用經緯度計算首都機場與虹橋機場的二維歐式距離為：10.06。而這個值顯然對應着球面距離1074.07千米。考慮到隨着維度的不同，兩者之間的比例也並非是定值。最簡單的例子就是，在赤道附近，一個經度差對應球面距離為 111.19 千米，而在 80° 緯度上則縮小到 19.31 千米，而在 90° 極點則幾乎為 0 千米。

這時，我們就可以利用 CoordinateSharp 中點的移動功能，計算出給定球面距離對應的經緯度，然後利用移動前後的經緯度再計算歐式距離，得出的結果才較為準確。

還是以虹橋機場為例，繪製其半徑為 100 千米的緩衝區：

// 對move進行簡單的封裝  public static double[] ConvertEarthDToPlaneXY(double lat, double lon, double distance, double bearing,      Shape shape = Shape.Ellipsoid)  {      var c0 = new CoordinateSharp.Coordinate(lat, lon);      c0.Move(distance, bearing, shape);      return new[] { c0.Latitude.ToDouble(), c0.Longitude.ToDouble() };  }    // 虹橋的位置  var hongqiao = new Coordinate(31.1982791377, 121.3356526703);  var hqPoint = gf.CreatePoint(hongqiao);  // 得出正西方向的100千米的位置hqMove  var hqMove = ConvertEarthDToPlaneXY(hongqiao.X, hongqiao.Y,100000, 270);  // 計算移動前後的歐式距離  var hqR = hongqiao.Distance(new Coordinate(hqMove[0], hqMove[1]);  // 計算buffer  var hqCircle = hqPoint.Buffer(hqR);  //驗證得出的buffer上的點與虹橋機場的位置  hqCircle      .Coordinates      .ToList()      .ForEach(t =>      {          var c1 = new CoordinateSharp.Coordinate(t.X, t.Y);          var c2 = new CoordinateSharp.Coordinate(31.1982791377, 121.3356526703);          var d0 = new CoordinateSharp.Distance(c1, c2, Shape.Sphere);          Console.WriteLine(d0.Kilometers);      });

驗證結果如下：

116.66885608202652  116.05364820797979  114.28764035412392  111.60513634687186  108.37944785253909  105.09000981559511  102.2637861494519  100.38650671049055  99.79581077763879  100.59276661692931  102.61595683755391  105.4929026708211  108.73914573002186  111.85882184494989  114.41831613054764  116.0891677559819  116.66885608202652  116.0891677559819  114.41831613054846  111.85882184495108  108.73914573002347  105.4929026708211  102.61595683755391  100.59276661692915  99.7958107776412  100.38650671049055  102.26378614945223  105.09000981559511  108.37944785253909  111.60513634687305  114.28764035412605  116.05364820797979  116.66885608202652

大部分結果都超過了 100 千米，並且誤差超過了 10 千米。比較妥善的方式就是計算多個方向的距離，然後取一個平均值。

public static double AverageDistance(double lat,double lon,double distance,int part)  {      var seg = 360 / part;      var rs = new List<double>();      var raw = new Coordinate(lat, lon);      for (var i = 0; i < 360; i += seg)      {          var move = ConvertEarthDToPlaneXY(lat, lon, distance, i);          var r = raw.Distance(new Coordinate(move[0], move[1]));          rs.Add(r);      }      return rs.Average();  }

上述緩衝區的距離計算改為：

// 取4個方向（0°、90°、180°和270°）對應的距離，然後求均值  var hqR = AverageDistance(hongqiao.X, hongqiao.Y, 100000,4);

得出的結果如下：

108.4796420377765  107.90879824241405  106.26991268291549  103.77978023046698  100.78401282732902  97.7268731801004  95.09732271673099  93.34698897092896  92.79099574532  93.52530908781738  95.401788528952  98.07519029469032  101.09498625442704  103.9991018274525  106.38288739908288  107.93950641121064  108.4796420377765  107.93950641121104  106.38288739908288  103.9991018274531  101.09498625442863  98.07519029469222  95.40178852895636  93.5253090878196  92.7909957453224  93.34698897093132  95.09732271673535  97.72687318010274  100.78401282733111  103.7797802304676  106.26991268291549  107.90879824241446  108.4796420377765

從結果中可以看出，誤差控制在 10 千米之內。

上述計算的是一個點的緩衝區，如果是 LineString 或者 Polygon 的緩衝區，則儘可能的取其上不同的點進行距離轉換後求取均值。

五、總結

總的來說，結合 CoordinateSharp 與 NetTopologySuite 可以進行許多有用的計算。但是誤差也不可避免，特別是緯度較高的時候。

如果各位有其他更好的解決方案，希望提交評論。