python解析GML檔

Table of contents

背景
GML檔案之讀取
polygons.csv的應用
程式下載
Reference

背景

GML(Geography Markup Language是臺灣地區的鄉鎮區界檔案格式，也是Open GIS社群通用的格式。
此處目標設定將GML檔案轉寫成raster檔。raster檔案在空品模擬系統之功用詳見python解析KML檔之背景說明。
目標raster解析度範圍為d4，GML雖然跟KML有很大的差別，但處理程序有很高的相似性。

目標

建立範圍解析度(d4) raster的值(整數分區代碼)

樣式範例

檔案來源：鄉鎮市區界線(TWD97經緯度))

<行政區域界線 xmlns="http://standards.moi.gov.tw/schema/pub" xmlns:gml="http://www.opengis.net/gml" xmlns:xlink="http://www.w3.org/1999/xlink" xmlns:ngis_primitive="http://standards.moi.gov.tw/schema/ngis_primitive" xmlns:gmd="http://www.isotc211.org/2005/gmd" xmlns:gco="http://www.isotc211.org/2005/gco" xmlns:utility="http://standards.moi.gov.tw/schema/utility" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xsi:schemaLocation="http://standards.moi.gov.tw/schema/pub pub.xsd">
  <gml:metaDataProperty>
    <ngis_primitive:NGIS_Primitive>
      <ngis_primitive:資料描述>鄉鎮市區界線(TWD97經緯度)</ngis_primitive:資料描述>
      <ngis_primitive:坐標參考系統識別碼>
        <gmd:RS_Identifier>
          <gmd:code>
            <gco:CharacterString>EPSG:3824</gco:CharacterString>
          </gmd:code>
        </gmd:RS_Identifier>
      </ngis_primitive:坐標參考系統識別碼>
      <ngis_primitive:坐標參考系統定義 xlink:herf="http://standards.moi.gov.tw/schema/epsg/3824.xml" />
      <ngis_primitive:資料內容對應時間>
        <gml:TimeInstant>
          <gml:timePosition indeterminatePosition="after">109-07-21</gml:timePosition>
        </gml:TimeInstant>
      </ngis_primitive:資料內容對應時間>
    </ngis_primitive:NGIS_Primitive>
  </gml:metaDataProperty>
  <gml:featureMember>
    <PUB_行政區域>
      <名稱>臺東縣成功鎮</名稱>
      <涵蓋範圍>
        <gml:MultiPolygon>
          <gml:polygonMember>
            <gml:Polygon>
              <gml:outerBoundaryIs>
                <gml:LinearRing>
                  <gml:coordinates>121.40981573700003,23.213692785000092 121.40984267700003,23.213661019000085 ...</gml:coordinates>
                </gml:LinearRing>
              </gml:outerBoundaryIs>
            </gml:Polygon>
          </gml:polygonMember>
        </gml:MultiPolygon>
      </涵蓋範圍>
      <行政區域代碼>10014020</行政區域代碼>
      <比例尺分母>5000</比例尺分母>
      <行政區域設置時間 />
    </PUB_行政區域>
  </gml:featureMember>
  <gml:featureMember>
    <PUB_行政區域>
      <名稱>屏東縣佳冬鄉</名稱>
      <涵蓋範圍>
        <gml:MultiPolygon>
          <gml:polygonMember>
...
 120.99327470700021,24.780576327000048</gml:coordinates>
                </gml:LinearRing>
              </gml:innerBoundaryIs>
            </gml:Polygon>
          </gml:polygonMember>
        </gml:MultiPolygon>
      </涵蓋範圍>
      <行政區域代碼>10004100</行政區域代碼>
      <比例尺分母>5000</比例尺分母>
      <行政區域設置時間 />
    </PUB_行政區域>
  </gml:featureMember>
</行政區域界線>

方案與檢討

python並沒有針對GML發展特殊的parser，如果用eTree來解讀，還算方便。但因為多邊形多達1036筆(群島)，再加上d4高解析度(由1公里加總成3公里)，因此須解決迴圈太多的問題。

GML檔案之讀取

GML是open GIS的共通語言，但相關軟體和討論似乎並不是很普遍。
以下為應用etree解讀的範例

GML檔案來源

檔案TOWN_MOI_1090727.gml：由內政地理資訊圖資雲整合服務平台TGOS網站下載。
或由內政部國土測繪圖資商城開放資料區下載(如下圖)
中文問題：此處為方便並沒有直接使用中文來搜尋，而是用嘗試錯誤法找出下列標籤。
- xzqy：「行政區域」(為字串序44以後)
- xzqudm：行政區域代碼(為新的8碼)，為序列中第2項。

rd_gml.py程式說明

GML雖然是open GIS的共通語言，但相關軟體和討論似乎並不是很普遍。還好etree還可以讀取。
- 獲取檔案的rootElement

import xml.etree.cElementTree as ET
from pandas import *
     3
rootElement = ET.parse("TOWN_MOI_1090727.gml").getroot()

次元件標籤中含有「行政區域」(xzqy)4個中文字，如果字串中有此4個中文字，則讀取「行政區域代碼」6個中文字(xzqydm)

s=set()
for subelement in rootElement.getiterator():
  for subsub in subelement:
    s=s|set([subsub.tag])
s=list(s)
s.sort()
xzqy=s[5][44:]
xzqydm=[i for i in s if xzqy in i][1]

讀取所有鄉鎮區的代碼all_town

all_town=set()
for subelement in rootElement.getiterator():
  for subsub in subelement:
    if xzqydm == subsub.tag:
      all_town=all_town|set([subsub.text])
all_town=list(all_town)
all_town.sort()

20~31：將多邊形座標讀出，存到wkt序列，共1036筆。

twnid={}
isq=0
wkt=[]
for subelement in rootElement.getiterator():
  for subsub in subelement:
    if xzqydm == subsub.tag:tid=subsub.text
    if subsub.tag == "{http://www.opengis.net/gml}coordinates":
      x = subsub.text
      point_for_pol =[(float(i.split(',')[0]),float(i.split(',')[1])) for i in x.split()]
      wkt.append(point_for_pol)
      twnid.update({isq:tid})
      isq+=1
    32

鄉鎮區名稱TOWN_MOI_1090727E.csv：由開放平台下載，去掉中文字以簡化過程，用代碼勾串。

    33  df_twn=read_csv('TOWN_MOI_1090727E.csv')
    34  ii=[int(twnid[i]) for i in range(len(twnid))]

儲存：存成csv檔，序列存在csv中會變成很長的字串，要調用會有一點麻煩(詳下)，但因為概念架構上比較方便，也就將就了。

df=DataFrame({'twnid':ii,'lonlats':wkt})
TOWNENG={i:j for i,j in zip(df_twn.TOWNCODE, df_twn.TOWNENG)}
df['TOWNENG']=[TOWNENG[i] for i in df.twnid]
df.set_index('twnid').to_csv('polygons.csv')

str2lst

DataFrame儲存格如果是一個字串，要改成序列，可以參考下列片段：

def str2lst(A):
    return [float(i) for i in A[1:-1].split(',')]

polygons.csv的應用

中央氣象局測站所在鄉鎮區之判別(town_cwb.py)
鄉鎮行政區多邊形轉1公里解析度tiff檔(多邊形csv檔)
附近行政區之定位中研判鄉鎮區是否相鄰

程式下載

github

Reference

wiki, 地理标记语言, wiki, 2021年3月1日