臺灣各級行政區域(縣市/鄉鎮/村里)邊界座標檔， 自選解析度、 存成 csv、 畫成 svg

政府資料開放平臺提供全國 縣(市)行政區域界線、 鄉(鎮、市、區)行政區域界線、 村里界圖 三個資料集。 如果你想從中萃取 「宜蘭縣各鄉鎮邊界.json」 或是 「高雄市鼓山區各里邊界」， 該怎麼做呢? 要如何降低解析度、 減少資料量? 如何存成文字檔， 用 gnuplot 畫成 svg 圖檔， 或是匯入其他一般「非地理資訊系統」的軟體裡面? 也許 Quantum GIS 之類的專業地理資訊系統本來就能處理 (如果看到教學文， 請分享一下網址囉)； 不過我就是喜歡用命令列的文字檔工具來處理。 感覺學習門檻比較低； 每個中間步驟都看得見， 比較踏實； 若需要批次處理也比較簡單； 學這些通用小工具， 將來可能在更多地方可以派上用場。

一、 成果預覽及軟體安裝

請先到 我的 umap 地圖， 進入左側 「圖層」 圖示， 關掉三個 「眼睛」 圖示當中的兩個， 然後按 「更多」， 選 「將地圖內嵌並分享」， 分別把三個圖層以 geojson 的格式下載回來研究。 第二三四節的目標， 就是要從官方資料產生這幾個 geojson 檔。

(只有這一小段需要用 root 身份做事) 好的， 那就開始備料囉：

1. 請抓回 multi2maxpolygon.py 跟 geojson2csv.py， 放在 /usr/bin 底下。 (對， 我正在學 python。) 記得要 chmod a+x multi2maxpolygon.py geojson2csv.py。
2. 安裝 mapshaper： apt-get install npm ; (cd /usr/bin ; ln -s nodejs node ) ; npm install -g mapshaper。
3. 安裝 json 轉檔萬用瑞士刀 jq
4. 安裝地理資訊超強轉檔套件： apt-get install gdal-bin 這樣才能使用 ogr2ogr 這個指令。

回到普通用戶身份。 假設你把抓回來的 zip 檔放在家目錄， 並且建立一個 admin/ 工作目錄用來存放我們轉檔時所產生的中間步驟檔案。

二、 產生縣市界座標 geojson 檔

以下是縣市界座標轉檔分解動作：

1. 解壓縮。 其實我們只需要 zip 壓縮檔裡面的 rar 檔即可。 解出來檔名是亂碼。 不管它， 繼續解 rar 檔， 得到 dbp、 prj、 shp、 shx 等等多個檔案。 目前我們只對 shp 有興趣。
2. 用 ogr2ogr 把 shp 轉成 geojson。
3. 挑出最大塊。 轉出來的 geojson 檔裡面有很多 MultiPolygon 而不是 Polygon。 因為很多縣市都管轄了離島， 例如龜山島隸屬於宜蘭縣頭城鎮、 彭佳嶼隸屬於基隆市中正區。 但是很多軟體 (例如 umap) 無法處理 MultiPolygon。 此外， 平常製作簡要的示意圖時， 其實好像也不必把所有細節畫出來。 所以我寫了 multi2maxpolygon.py 用來把 geojson 檔裡面的每一個 MultiPolygon 都轉成它轄區內面積最大塊的單一 Polygon。 (不要控告我丟棄國土啊!) 轉檔時， 可以順便用 -n C_Name 指定要把每個縣市的 properties 欄位內的 C_Name 欄位拿來當做 name 欄位， 方便其他軟體顯示每一塊多邊形的名稱。
4. 降低解析度。 光是縣市界的原始 shp 檔就有 5MB； 鄉鎮村里界更大； 轉成適合人眼閱讀的 geojson 檔就更嚇死人了， 所以一定要降低解析度。 我趁著颱風天把程式寫好之後， 才發現人家早就寫好 mapshaper， 而且效果比我做的更好。 可以用 -simplify 選項指定要保留多少百分比的資料， 後面的數字可以是多少百分比， 或是零到一之間的一個小數。

最後產生的 geojson 檔就可以上傳到 umap 檢視囉!

以下是完整指令。

rm -f *.rar ; unzip ../縣\(市\)行政區域界線.zip '*1040808.rar' ; rar x *1040808.rar 
ogr2ogr -f GeoJSON 縣市界.json 縣市界\(經緯度\)1040808.shp 
multi2maxpolygon.py -n C_Name 縣市界.json | mapshaper - -simplify 5% -o 臺灣各縣市邊界.json

三、 產生某縣的鄉鎮界座標 geojson 檔

產生宜蘭縣各鄉鎮邊界的 geojson 檔， 大致步驟跟產生全國縣市邊界檔一樣。 為了只挑出宜蘭縣的資料， 還必須多一句： 用 jq 挑出隸屬於宜蘭縣的所有鄉鎮。 如果只用 map 抓出一堆 features， 產生出來的 geojson 檔， umap 是可以看得懂； 但 mapshaper 就看不懂了， 所以必須再把最外層的 .features 重新包回去。 jq 瑞士刀真的超級好用啊!

rm -f *.rar ; unzip ../鄉\(鎮、市、區\)行政區域界線.zip '*1040808.rar' ; rar x *1040808.rar
ogr2ogr -f GeoJSON 鄉鎮界.json 鄉鎮界\(經緯度\)1040808.shp
jq '.features = (.features | map(select(.properties.C_Name=="宜蘭縣")))' 鄉鎮界.json | multi2maxpolygon.py -n T_Name | mapshaper - -simplify 10% -o - > 宜蘭縣各鄉鎮邊界.json

四、 產生某鄉鎮的村里界座標 geojson 檔

產生高雄市鼓山區各里邊界 geojson 檔， 大致步驟跟產生某縣市各鄉鎮邊界檔一樣。 不過官網提供資料時， 在村里等級， 採用的座標單位是公尺而不是經緯度。 (我一開始沒注意， 轉出來的 geojson 檔裡面的 "經緯度" 座標都大得不像話、 繞著地球轉好幾圈!) 所以拿 ogr2ogr 把村里界的 shp 轉成 geojson 的同時， 應以 -t_srs 選項來指定 「要採用縣市界 (或鄉鎮界) 的投影及座標轉換」。

rm -f *.rar ; unzip ../村里界圖\(TWD97_121分帶\).zip '*1040804.rar' ; rar x *1040804.rar
ogr2ogr -f GeoJSON 村里界.json -t_srs 縣市界\(經緯度\)1040808.prj Village_NLSC_121_1040804.shp 
jq '.features = (.features | map(select(.properties.T_Name=="鼓山區")))' 村里界.json | multi2maxpolygon.py -n V_Name | mapshaper - -simplify 30% -o - > 高雄市鼓山區各里邊界.json

五、 轉成簡單文字檔、 用 gnuplot 畫出來

如果要餵給 gnuplot 或是其他 「非地理資訊系統」 的軟體吃， 最簡單的格式還是 csv 或原始的純文字格式。 搜尋不到現成的工具， 只好自己寫了 geojson2csv， 把 geojson 檔裡面所看到的所有座標， 管它三七二十一 (LineString 或 MultiLineString 或 Polygon 或 MultiPolygon) 全部以每列一對 x,y 座標的格式印出來。 每塊線段或多邊形跟相鄰的線段或多邊形之間會多印至少一個空列。 預設是以逗點分格 x y 座標 (所以才叫做 comma-separated values 嘛)； 但 gnuplot 喜歡吃空格分格的資料， 所以用 -d 選項指定分隔字元： geojson2csv.py -d ' ' 臺灣各縣市邊界.json > tw-counties.txt 。 然後把這段指令剪貼存檔， 叫做 tw-counties.gpt：

set size ratio 1/cos(24.0/180*pi)
set style data line
set term svg
set output 'tw-counties.svg'
plot 'tw-counties.txt'

再用 gnuplot tw-counties.gpt 產生 tw-counties.svg ， 最後用瀏覽器或 inkscape 打開 tw-counties.svg 就大功告成了! 其中 set size ratio 的用意是設定高矮胖瘦的比例。 高度必須拉長 1/cos(lat) 倍， 才會是我們習慣看的 麥卡托投影地圖。 (lat 是此圖大約的緯度)

如果你成功地用本文的方法把自製的地圖放進你的文章/海報/軟體/...等等作品， 請順便幫忙宣傳本文連結， 讓你的讀者/用戶也一起來學自製地圖吧! 也歡迎留言分享你的作品的連結哦!