NetTopologySuite 處理 GIS (1)：讀寫 Shapefile

上次使用 GDAL 進行 Shapefile 檔案的讀取之後，在實務上又碰到其他問題了。有個需求是要針對點位去找出他屬於哪一塊土地，原本想說把檔案內的資料一次讀出來，在透過程式去處理就好，但是 GDAL 有一個問題：沒有註冊、沒有開啟 Shapefile，貌似就不能進行 Geometry 的相關操作，如判斷交集、新增 Buffer 之類的，直接把我的美好幻想打碎 💔

為了解決這個問題，我使用了另外一個也相當有名的套件 — — NetTopologySuite，針對我從 Shapefile 裡取出的 WKT，去進行轉換和比對。但一個專案裡面使用兩個性質相同的套件，我個人是覺得很沒必要，所以決定統一使用 NetTopologySuite 作為地理資訊操作的主要套件。接下來就一起看看如何使用 NetTopologySuite 吧。

安裝 NetTopologySuite.IO.Esri

如果你使用 NetTopologySuite 和 Shapefile 這兩個關鍵字去 google ，第一個應該會看到的是 NetTopologySuite.IO.ShapeFile 這個套件，但可以在 README 的地方看到團隊沒有要繼續更新了，建議轉而使用 NetTopologySuite.IO.Esri。但經過搜尋之後，發現他還沒有被發佈到 Nuget 上，所以必須手動安裝。

安裝的方法很簡單，在 GitHub 頁面選擇下載成 zip 檔，或是 clone 下來之後，按照個人習慣（我都用 dotnet publish -c Release ）將他編譯成 dll 檔案，接著在你原本的專案內新增這些 dll 檔案的參考即可。

圖為編譯後的全部內容，選取起來的三個 dll 中，最主要是第三個 NetTopologySuite.IO.Esri.Shapefile，其他兩個也會使用到，可以從這裡加入參考，或是從 Nuget 中加入也行。

新增參考

如果你是使用 Visual Stadio 可以看看官方文件；如果你和我一樣是使用 Rider，對著你的 project 右鍵，點選 Add，就可以看到 Add Reference 的選項。點擊之後會打開一個視窗，在視窗最下面點選 Add From … 就能自行新增參考了。

Shapefile 操作

將 NetTopologySuite.IO.Esri 加入我們的專案之後，就可以開始對我們的 Shapefile 進行操作了。

讀取

using NetTopologySuite.IO.Esri;

foreach (var feature in Shapefile.ReadAllFeatures(@"C:\Users\emct\Downloads\taiwan\All\4326\COUNTY_MOI_1090820_4326.shp"))
{
    foreach (var attrName in feature.Attributes.GetNames())
    {
        Console.WriteLine($"{attrName,10}: {feature.Attributes[attrName]}");
    }

    // Console.WriteLine($"     SHAPE: {feature.Geometry.AsText()}");
    Console.WriteLine($"      SRID: {feature.Geometry.SRID}");
    Console.WriteLine();
}

根據 GitHub 上的寫法，只要簡單的這幾行，馬上就能讀取 Shapefile 了！執行後的結果如下：

可以看到之前讓我們在 GDAL 中痛苦不已的中文，在 NetTopologySuite 上居然很正常的顯示了！（痛哭流涕）

另外也可以看到，NetTopologySuite 在寫法上相對來說更為簡潔，直接讀取 Feature 就好，再透過欄位名稱從 Atrributes 取得 Field 的值，不像 GDAL 還需要從 Layer 先取得 Layer 的數量跑回圈，接著從 Layer 拿 layer.GetNextFeature()，看起來好看多了！

輸出

using System.Text;
using NetTopologySuite.Features;
using NetTopologySuite.IO;
using NetTopologySuite.IO.Esri;

var features = new List<Feature>();
var wktReader = new WKTReader();

var fcuGeometry = wktReader.Read(
    "POLYGON((214094.18733422138 2675204.397772933,214613.04806256315 2675203.096234374,214612.0331897791 2674795.5200602347,214093.15758024747 2674796.8214499434,214094.18733422138 2675204.397772933))");

var fcuAttributes = new AttributesTable();
fcuAttributes.Add("日期", DateTime.Today);
fcuAttributes.Add("地點", "逢甲大學");

var fcuFeature = new Feature(fcuGeometry, fcuAttributes);
features.Add(fcuFeature);

var nightMarketGeometry = wktReader.Read(
    "POLYGON((213945.3342767608 2674876.347737663,214097.93845940504 2674873.793433358,214099.4943339883 2674626.6362772766,214474.48457545802 2674625.694102451,214511.39058478415 2674562.7289472898,214223.11502519806 2674368.326804782,213983.77470059253 2674559.719947559,213945.3342767608 2674876.347737663))");
var nightMarketAttributes = new AttributesTable();
nightMarketAttributes.Add("日期", DateTime.Today);
nightMarketAttributes.Add("地點", "逢甲夜市");
var nightMarketFeature = new Feature(nightMarketGeometry, nightMarketAttributes);
features.Add(nightMarketFeature);

// 指定的 EPSG 字串，可以在 https://epsg.io/ 找到
var wkt = "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\"]]";

Shapefile.WriteAllFeatures(features, @"files\create.shp", Encoding.UTF8, wkt);

根據 GitHub 上的範例，使用 Shapefile.WriteAllFeatures() 來建立，執行之後遇到了一個錯誤：

要解決這個問題需要花一點篇幅進行，所以會另外開一篇來講解，包含為甚麼會有這個錯誤，以及如何修改原始碼讓我們可以成功建立一個 Shapefile。

參考文章

• NetTopologySuite.IO.Esri GitHub