Skip to content
Coding

NetTopologySuite 處理 GIS (3)：座標系統轉換

By Anthea / 芫荽 3 min read
NetTopologySuite 處理 GIS (3)：座標系統轉換

前兩篇介紹了如何使用 NetTopologySuite 來讀寫 Shapefile 檔案，這篇要介紹的是如何將讀取出來的 Geometry 進行座標系統的轉換。

  • GIS 座標系統簡介
  • 使用 NetTopologySuite 搭配 ProjNet 進行座標系統轉換

GIS 座標系統

一組完整的座標系統會包含兩個資訊：大地基準（Datum）和座標格式（Format/Grid）。

大地基準

和地球形狀最相似的數學模型（參考橢球體），最常見的是由美國所制定的世界測量座標系統 WGS84。由於地球表面是凹凸不平的，如果在偏離參考橢球體太多的地方依舊使用 WGS84 進行參考，會有不準確的情況發生，因此大多數國家都會選擇使用比較符合自己地區的大地基準，以台灣為例，就會比較常使用 TWD97。

座標格式

通常會需要標記座標的地方會是在平面地圖上，將橢圓形的地球攤平到地圖上的動作就稱為投影。因此，座標格式也就分為兩種，一種是經緯度座標（Geographic Coordinate System），通常用在球面上；一種是投影座標（Projected Cooredinate System），通常用在地圖上。

以台灣來說，最常使用的就是二度分帶 TM2，完整的座標系統名稱就會是 TWD97 / TM2 zone 121。

利用 ProjNet 進行座標轉換

using NetTopologySuite.Geometries;
using NetTopologySuite.IO;
using ProjNet.CoordinateSystems;
using ProjNet.CoordinateSystems.Transformations;

// 座標系統的 WKT，來源：https://epsg.io/
const string Wgs84Wkt = "GEOGCS[\"WGS84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.01745329251994328,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4326\"]]";
const string Twd97Wkt = "PROJCS[\"TWD97 / TM2 zone 121\",GEOGCS[\"TWD97\",DATUM[\"Taiwan_Datum_1997\",SPHEROID[\"GRS 1980\",6378137,298.257222101],TOWGS84[0,0,0,0,0,0,0]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"3824\"]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"latitude_of_origin\",0],PARAMETER[\"central_meridian\",121],PARAMETER[\"scale_factor\",0.9999],PARAMETER[\"false_easting\",250000],PARAMETER[\"false_northing\",0],UNIT[\"metre\",1],AXIS[\"Easting\",EAST],AXIS[\"Northing\",NORTH],AUTHORITY[\"EPSG\",\"3826\"]]";

Console.WriteLine("請輸入經度：");
var lon = Console.ReadLine();

Console.WriteLine("請輸入緯度：");
var lat = Console.ReadLine();

var wktBeforeConvert = $"POINT({lon} {lat})";
Console.WriteLine($"轉換前的 WKT 為：{wktBeforeConvert}");

// 將 WKT 讀成 Geometry
var wktReader = new WKTReader();
var originPoint = wktReader.Read(wktBeforeConvert);

var coordinateSystemFactory = new CoordinateSystemFactory();
var coordinateTransformationFactory =  new CoordinateTransformationFactory();
    
var twd97 = coordinateSystemFactory.CreateFromWkt(Twd97Wkt);
var wgs84 = coordinateSystemFactory.CreateFromWkt(Wgs84Wkt);

// CreateFromCoordinateSystems(wgs84, twd97) = 從 WGS84 轉成 TWD97
var fromCoordinateSystems = coordinateTransformationFactory.CreateFromCoordinateSystems(wgs84, twd97);
var convertedPoint = Transform(originPoint, fromCoordinateSystems.MathTransform);

var wktAfterConvert = convertedPoint.AsText();
Console.WriteLine($"轉換後的 WKT 為：{wktAfterConvert}");


static Geometry Transform(Geometry geom, MathTransform transform)
{
    geom = geom.Copy();
    geom.Apply(new MTF(transform));
    return geom;
}

sealed class MTF : ICoordinateSequenceFilter
{
    private readonly MathTransform _mathTransform;

    public MTF(MathTransform mathTransform) => _mathTransform = mathTransform;

    public bool Done => false;
    public bool GeometryChanged => true;
    public void Filter(CoordinateSequence seq, int i)
    {
        // 如果有需要轉換 z 座標的話可以把註解打開
        
        double x = seq.GetX(i);
        double y = seq.GetY(i);
        // double z = seq.GetZ(i);
        
        _mathTransform.Transform(ref x, ref y);
        // _mathTransform.Transform(ref x, ref y, ref z);
        
        seq.SetX(i, x);
        seq.SetY(i, y);
        // seq.SetZ(i, z);
    }
}

輸出的結果如下：

從線上的轉換工具可以檢查我們轉換的結果：

或是利用 PostgreSQL 搭配 PostGIS 檢查：

參考文章