2012年11月6日 星期二

[Python] 讀取 Shapefile 格式(*.dbf, *.prj, *.shp, *.shx, *.sbn, *.sbx) 以臺北市公車站牌位置圖為例

前陣子參加 2012 Yahoo! Open Hack 時,終於有空把玩台北市政府公開資料,碰到了之前一直想解卻沒空解的問題,以 臺北市公車站牌位置圖 為例,其資料格式:


$ ls bus
busstop.dbf busstop.shx 公車站牌匯入.sbx
busstop.prj 公車站牌匯入.dbf 公車站牌匯入.shp
busstop.shp 公車站牌匯入.sbn 公車站牌匯入.shx


一堆很奇妙的格式,感覺就像給你一個加密過的資料,讓人十分不解,這大概就是門外漢的滋味 Orz 回想起來,這種感觸還真像當年自己去刻圖書館 MARC 格式的心情 XD


所幸參加 Open hack 才有時間把玩,當下不知哪根筋被打通,忽然間明瞭這種 *.shp 格式就叫 Shapefile ,如此一來有關鍵字就能解了,馬上找到 pyshp - Python Shapefile Library


廢話不多,用法:


#!/usr/bin/env python
# -*- coding: utf8 -*-
import shapefile
import sys
reload(sys)
sys.setdefaultencoding('utf-8')


#sf = shapefile.Reader("bus/公車站牌匯入") # big5 & TWD67 format
sf = shapefile.Reader("bus/busstop") # utf8 & WGS84 format
#fields = sf.fields
print sf.fields 


records = sf.records()
shapeRecs = sf.shapeRecords()
print "Total: ",len(shapeRecs)
bushash = {}
for i in range(len(shapeRecs)):
   #busstop = records[i][1].decode('big5', 'ignore') # decode from big5
   busstop = records[i][1].decode('utf8', 'ignore') # decode from utf8
   if busstop in bushash:
      continue
   else:
      bushash[busstop] = ''


   print busstop
   print shapeRecs[i].shape.points[0]
   #print = FromTWD67TM2ToWGS84( shapeRecs[i].shape.points[0][0], shapeRecs[i].shape.points[0][1] )


成果:


[('DeletionFlag', 'C', 1, 0), ['SATOPID', 'N', 20, 0], ['STOPNAME', 'C', 20, 0], ['CITYNAME', 'C', 5, 0]]
Total: 5540
101國際購物中
[121.56524, 25.03415]
101購物中心
[121.565389, 25.034183]
8號水門
[121.53166335449, 25.072300254772]
一女中
[121.51507, 25.03828]
二二八和平公
[121.51303, 25.04215]
二二八紀念館
[121.513913, 25.030761]
....


眼尖的大概可以看到,第一筆有漏字,資料原始為 "101國際購物中<E5><BF>",不知是不是 pyshp 讀取錯誤還是原始資料有誤。此外,同一個站名會有多個座標點,在此僅先過濾掉,有興趣的在自行把玩吧!


其他補充:


安裝 GDAL (MacPorts: $ sudo port install gdal)
$ ogr2ogr -f geojson /tmp/output.json busstop.shp
$ enca /tmp/output.json
Unrecognized encoding
$ cat /tmp/output.json


{
"type": "FeatureCollection",

"features": [
{ "type": "Feature", "properties": { "SATOPID": 2028.0, "STOPNAME": "101å<9c><8b>é<9a><9b>è³¼ç<89>©ä¸­å¿", "CITYNAME": "å<8f>°å<8c>" }, "geometry": { "type": "Point", "coordinates": [ 121.56524, 25.03415 ] } }
,
...
]


$ ogr2ogr -f geojson /tmp/output2.json 公車站牌匯入.shp
$ enca /tmp/output2.json
Universal transformation format 8 bits; UTF-8
$ cat /tmp/output2.json


{
"type": "FeatureCollection",

"features": [
{ "type": "Feature", "properties": { "SATOPID": 2028.0, "STOPNAME": "101國際購物中??", "CITYNAME": "台??", "stopNAMEB5": null }, "geometry": { "type": "Point", "coordinates": [ 306245.046938922838308, 2769874.23026815475896 ] } }
,
...
]


沒有留言:

張貼留言