python解析KML檔

Table of contents

背景

緣起

  • 網格模式在進行污染來源分析時，經常需要排放區域的分區定義，一個網格化的點陣圖(raster)：
    • 不論是CAMx-OSAT、CAMx-PSAT或者是
    • CMAQ-ISAM模組，都需要這樣的raster檔案。
    • 具體問題：因為缺乏大陸方面的地理分區定義資料，在反軌跡計算時無法歸因於省份、排放區域等定性標籤。
  • 然而一般網路公開的GIS檔案多是向量圖檔，是省界的polygon檔案，不正好是raster檔案。
    • 造成這種網路現象主要的理由，是因為raster檔案有其解析度的限制，提供者不知該提供到什麼樣的解析度才恰當，畢竟使用者定義的範圍解析度樣態真的太多樣了。
  • 本作業聚焦在省界shape(或KML)檔，轉成模式解析度(d1)的raster檔案。
    • KML(Keyhole markup language)是Google發展的地理資訊語言檔案格式，經常使用在google地圖、open street 地圖平台上。

目標

  • 建立範圍解析度(d1) raster的值(整數分區代碼)
  • 除此之外，此次作業的重點在建立合理的分區方式。
    • 這是因為大陸地區現行的省、市、自治區等，有57個之多，超過合理的處理範圍，需要進一步整併。

方案與檢討

  • 向量檔案轉成點陣圖檔的方式有很多種，除了在GIS內操作，也有獨立的軟體如Corel，也有業者提供網站服務。
    • 而在我們的應用中，2個檔案的座標系統也並非是完全一樣，不能直接套用。
  • 而一般的GIS軟體(gdal程式)也可以做到shape to raster轉換，然而不知道是大陸方面檔案格式(層次或資料庫)的差異、或者數據太多(約有170個polygon)，QGIS無法作動!
  • python可以直接讀取向量檔(shape及KML)，然而因為對KML格式較為熟悉，因此選擇以KML轉raster之進路。
    • shape檔案之讀取及處理詳參shape_to_raster應用在d4範圍、d5解析度案例。
  • 大陸行政區之整併方式，除了考量地理氣候方面的特徵，因為是空氣污染領域，除了文獻中常見的京津冀、長三角、珠三角 之外，也有人用華北、華中、華南來分大區。
    • 如果分區越細，軌跡機率就越少、排放量分配也較少。
    • 此外因為距離台灣很遠，d1範圍的解析度也有限(81K)，如果超過10區，每區網格數會太少。
    • 以美國而言，有所謂climate region，臺灣地區也有所謂空品區的概念。 建議還是以大陸官方空氣質量管理之分區方式為宜。

檔案來源與解壓縮

...
    <Placemark id="ID_00000">
      <name>Anhui Province</name>
...
      <styleUrl>#PolyStyle00</styleUrl>
      <MultiGeometry>
        <Polygon>
          <extrude>0</extrude><altitudeMode>clampToGround</altitudeMode>
          <outerBoundaryIs><LinearRing><coordinates> 
               116.935206535,34.08702338100009,0 
               116.930544829,34.08849598800008,0 
               116.918619955,34.09428245300003,0                
...

rd_kml.py

程式分段說明

  • 模組之調用
    • 此處使用到pykml的parser
1	from pykml import parser
2	from os import path
3	from pandas import *
4
  • 讀取doc.kml：來源見前述。
5	kml_file = path.join('doc.kml')
6	with open(kml_file) as f:
7	  doc = parser.parse(f).getroot()
  • 讀取所有的位置標籤(Placemark)，及其名稱(names)
    • 名稱names(會出現在KML檔案的Placemark(plms)的名稱標籤)，此名稱與網站提供的xls資料表名稱一致。共57筆，有重複。(line 8~9)
8	plms=doc.findall('.//{http://www.opengis.net/kml/2.2}Placemark')
9	names=[i.name for i in plms]
10
  • pykml並未內設一制式的從屬關係層次，而是用getparent 方法找出來物件之間的親子關係。
    • 在本案例中MultiGeometry與Polygon之間是有關係的，而MultiGeometry與Placemark的順序相同，可以引用相同的名稱，因此就可以由下向上，串連找到每一個Polygon的名稱。
      1. 先找到MultiGeometry物件，以及其個別的標籤(MultiGeometry 沒有名稱，只有標籤) (line 11~12)
      2. 找到Polygon物件，及其上層MultiGeometry 物件的標籤(line 14~15)
  • 讀取幾何形狀(mtgs)與其標籤(mtg_tag)
11	mtgs=doc.findall('.//{http://www.opengis.net/kml/2.2}MultiGeometry')
12	mtg_tag=[str(i.xpath).split()[-1][:-2] for i in mtgs]
13
  • 讀取多邊形(plgs)與其來源(plg_prt)
14	plgs=doc.findall('.//{http://www.opengis.net/kml/2.2}Polygon')
15	plg_prt=[str(i.getparent().values).split()[-1][:-2] for i in plgs]
16	

17	lon,lat,num,nam=[],[],[],[]
18	n=0
19	for plg in plgs:
20	  iplg=plgs.index(plg)
21	  imtg=mtg_tag.index(plg_prt[iplg])
22	  name=names[imtg]
23	  coord=plg.findall('.//{http://www.opengis.net/kml/2.2}coordinates')
24	  c=coord[0].pyval.split()
25	  for ln in c:
26	    if n%3==0:
27	      lon.append(ln.split(',')[0])
28	      lat.append(ln.split(',')[1])
29	      num.append('n='+str(n))
30	      nam.append(name)
31	    n+=1
  • 將結果寫進csv檔案(txt.write版本)
32	#with open('doc.csv','w') as f:
33	#  for i in range(len(num)):
34	#    f.write(str(lon[i])+','+str(lat[i])+','+str(num[i])+','+nam[i]+'\n')
35	#   if n%100==0:print(n)
  • DataFrame.to_csv版本
    • 功能主要在將每一個座標點輸出成csv檔，而將其名稱寫在點位的說明欄，可以應用csv2kml.py程式再次轉成kml檔案，檢視其多邊形及名稱解析是否正確。
36	df=DataFrame({'lon':lon,'lat':lat,'num':num,'nam':nam})
37	df.set_index('lon').to_csv('doc.csv')

程式下載

應用

  • 除了此處

    Reference

  • google,KML Tutoria,developers.google.com, Last updated 2021-09-07 UTC.
  • wiki, Keyhole標記語言, wiki, 2021年2月7日.
  • Tyler Erickson, pyKML v0.1.0 documentation,pythonhosted, 2011
  • 純淨天空, pykml.parser.fromstring函數,vimsky
  • contributors, Python fromstring Examples, python.hotexamples