那些驚豔的 GIS 輪子

時間 2020-04-18

標籤那些驚豔 gis 輪子简体版

原文原文鏈接

1、前言

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

真的懊惱早沒有人告訴我這些優秀的做品的存在。此前都一直在調用 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 的實現。html

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

2、CoordinateSharp

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

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

這裏有一點使人疑惑的地方就是：爲何會有時間信息。這正是它的獨到之處，不只僅進行座標轉換，還帶有計算日、月升落時間，位置等天體信息的能力：github

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

基本一致，偏差的緣由由地球是橢球體、自轉公轉速率緩慢改變等引發。這個差異可以讓人接受。spring

2.空間計算

計算球面距離docker

// 北京首都機場經緯度
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。

補充

雖然CoordinateSharp計算功能豐富，可是不少功能卻攪和在一塊兒。好比說計算距離，不少時候咱們僅僅須要距離只一個結果，而CoordinateSharp卻給了其餘豐富的信息，這是優勢。固然，這是以犧牲效率爲代價的。我在32核機子上併發計算幾百個點的距離居然花費了2秒左右的時間。這的確讓人難以忍受。若是讀者也僅僅須要距離計算，那麼不妨看看下面的算法：

public static double CalculateGeoDistance(double sLatitude, double sLongitude, double eLatitude, double eLongitude)
{
    var radiansOverDegrees = Math.PI / 180.0;

    var sLatitudeRadians = sLatitude * radiansOverDegrees;
    var sLongitudeRadians = sLongitude * radiansOverDegrees;
    var eLatitudeRadians = eLatitude * radiansOverDegrees;
    var eLongitudeRadians = eLongitude * radiansOverDegrees;

    var dLongitude = eLongitudeRadians - sLongitudeRadians;
    var dLatitude = eLatitudeRadians - sLatitudeRadians;

    var result1 = Math.Pow(Math.Sin(dLatitude / 2.0), 2.0) +
                  Math.Cos(sLatitudeRadians) * Math.Cos(eLatitudeRadians) *
                  Math.Pow(Math.Sin(dLongitude / 2.0), 2.0);

    // 地球半徑，單位：米
    var result2 = 6371000.0 * 2.0 * Math.Atan2(Math.Sqrt(result1), Math.Sqrt(1.0 - result1));

    return result2;
}

一樣的計算上海虹橋機場到北京首都機場的結果基本和CoordinateSharp一致，可是計算效率卻極高。單線程計算1000次CoordinateSharp耗費13秒左右，而該算法僅僅0.0003969秒左右。效率天然不言而喻，有着3萬多倍的差距。不過，最終看要解決什麼樣的問題，這不影響CoordinateSharp仍然成爲一個優秀的開源庫。

3、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。

4、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 的緩衝區，則儘量的取其上不一樣的點進行距離轉換後求取均值。