用QGIS + OpenStreetMap + Python 處理線形地理資料

本文內容是敘述如何用QGIS，結合開放街圖（OpenStreetMap）的圖資，以及Python資料處理，在QGIS上處理資料。

本次使用圖資：福爾摩沙高速公路（即我國國道三號，下稱福高）的路線圖。

取得圖資（XML）：

去OpenStreetMap尋找福高路線圖的XML檔。搜尋「福爾摩沙高速公路(北上)」，取得該路線圖編號（3876887）。並於畫面下方點選「下載XML」，進入XML顯示頁面。操作過程如下三張圖所示

然而，上述的XML不含完整的地理資訊。因此須在如上圖XML頁面的網址列加入字串，變成如：

https://www.openstreetmap.org/api/0.6/relation/3876887/full

方可使該XML顯示頁面完整顯示經緯度等資訊。正確結果如下圖：

在該頁面空白處點選右鍵 -> 「另存新檔」，即可下載該高速公路線形的XML檔。

匯入Python與處理

然而，XML並不好利用。因此必須將其轉換成CSV等更好處理的格式。較方便的做法，便是利用Python的語法與套件處理。

利用語法（參考並修改自此網頁）如下：

from xml.etree import ElementTree
import pandas as pd

tree = ElementTree.parse("full.xml") #請使用您下載的檔案的位置與檔名
root = tree.getroot()

rows = [] 
 
for i in root:        
   rows.append(i.attrib) 
cols = list(rows[1])
df = pd.DataFrame(rows, columns=cols) 
df.to_csv('Formosa_Freeway.csv') 

完成圖如下：

將CSV檔匯入QGIS

最後，進入QGIS將CSV檔匯入（Layer -> Add Layer -> Add Delimited Text Layer），過程如下圖，並開啟您後製完成的CSV檔案。並在Geometry Definition處，指定檔案的經緯度欄位，與座標參考系統（Geometry CRS, 即地圖投影方式；本例為WGS84）後完成匯入。

匯入結果如下圖所示：

從結果也可看出，相當符合福高的實際線形。

不過，用XML轉CSV匯入時是以點為單位，若您覺得一個點一個點看起來很醜，是可以透過一些工具/外掛程式去把點轉換成線或面。但我還要找找。但對於初階的地理資訊應用，應該尚稱足夠。

進階利用1 - 套疊地形圖

如果我想知道做出來的結果是不是對的，可以套疊其他加工良好的圖資，例如內政部公佈的「內政部20公尺網格數值地形模型資料」（具體作法可用Layer -> Add Layer -> Add Raster Layer 加入您下載回來的圖資）。完成結果如下圖所示。亦可進行其他加值應用。

進階利用2：用Python自動加工OpenStreetMap的XML

如果嫌上述手動操作Python和網頁的模式麻煩，其實可透過Python的語法與套件，自動擷取與存檔OpenStreetMap的XML。具體語法（參考並修改自此）如下：

import traceback
import urllib3
import xmltodict
import pandas as pd

def getxml():
   url = "https://www.openstreetmap.org/api/0.6/relation/3876887/full"

   http = urllib3.PoolManager()

   response = http.request('GET', url)
   try:
       data = xmltodict.parse(response.data)
   except:
       print("Failed to parse xml from response (%s)" % traceback.format_exc())
   return data

data=getxml()

my_list = data["osm"]["node"]
cols = list(my_list[1])
df = pd.DataFrame(my_list, columns=cols) 
df.to_csv('Formosa_Freeway_Northbound.csv') 

匯入QGIS的方式則大致與前一段相同。如此便可省卻前處理XML和存成CSV檔案的龐雜步驟與心力。然而運用該語法時，仍須留意著作權條款以及該網站是否有禁止自動化與批次下載等規定，以免觸法。