#--------------------------20200227更新----------------------------------------
德國人第二天郵件就反饋了。確實嚴謹，也學到了。
1 osm 全部坐標信息只在points里， 后面的way 和relation 只有點id的索引。
他給出的讀取方法是：
l0 = osm.GetLayer(0) f0 = None fcount = 0 while l0.GetNextFeature() is not None: fcount += 1
必須順次讀取，先讀完全部的points再讀后面的
2 pbftoosm 不維護了。
Pbftoosm is not maintained any more. Please use Osmconvert, its direct
replacement, or Osmium instead.
osmium cat -o nk.osm north-korea-latest.osm.pbf
3 不兼容報錯，主要是因為刪除了上傳者個人信息。 他們驗證過的工具可以處理這個問題。
Our files follow the standard with the exception that personal metadata
fields are set to 0 or empty strings to protect the privacy of
OpenStreetMap contributors. Our tests pointed out that this is not an
issue for established tools such as Osmium, Osmosis and Osmconvert.
4 處理是基於導入pg的
We do not plan to support pbftoosm in future because its development
stopped nine years ago. We do not use GDAL to read PBF files ourselves
because it is not performant. Instead, we import them into a PostgreSQL
database and run
畢竟咱是外行，還是學到很多啊
#----------------------------------------------------------------------------
項目需要分析過濾osm上的數據, 在
可以下載到按大洲,和國別 的 osm數據, 所有國家提供 .osm.pbf格式 和 raw osm壓縮 .osm.bz2 格式, 數據少的國家 提供 . shp 格式的.
對shp的處理當然沒問題. 但在數據大的國家動輒幾G, pbf體積最小, 就下載了.但是讀取可坑死我了. 前后折騰了N天.終於明白了.
一.症狀
直接說結論:這個網站的pbf格式, 不是標准的pbf! 表現症狀有2:
1.1 不被 gdal ogr 庫支持
在py下
import gdal, ogr gdal.SetConfigOption('OGR_INTERLEAVED_READING', 'YES') #gdal.SetConfigOption('OSM_MAX_TMPFILE_SIZE', '2000') if __name__ == '__main__': fname = './data/osm_pbf/north-korea-latest.osm.pbf' osm = ogr.Open(fname) print(osm.GetLayerCount()) layer = osm.GetLayer(0) #layer = data.GetLayer('lines') print(layer.GetFeatureCount()) feature = layer.GetNextFeature() o1 =feature.ExportToJson(as_object=True) print(o1)
1 所有的對 layer 0 也就是points 是可以運行的,但是 lines multipolygons 一個features都讀不出來, 只返回None
layer.GetFeatureCount()
對points 在win下, 使用老外編譯好的gdal https://www.lfd.uci.edu/~gohlke/pythonlibs/#gdal
可以讀出數量
在ubuntu下, 自己辛苦源碼編譯的gdal下, 全是-1
但是, gdal 測試里自帶的pbf GDAL/autotest/ogr/data/test.pbf
是沒任何問題的 !
1.2 不被openstreetmap官網的pbf轉換工具支持
https://wiki.openstreetmap.org/wiki/Pbftoosm#Download
官網提供了pbf2osm 直接一個可執行文件 比如x64 linux 直接下載一個pbftoosm64 放到bin 之類的路徑下
對gdal自帶的pbf是可以正常轉型的 轉成osm的 注意尖括號不能省略
pbftoosm64 <test.pbf > test.osm
但是轉他家下載的.osm.pbf就報錯
➜ osm_pbf git:(master) ✗ pbftoosm64 <north-korea-latest.osm.pbf > nk.osm
pbftoosm Warning: header block element type unknown: 0x80 0x02.
pbftoosm Warning: header block element type unknown: 0xA6 0xC8.
pbftoosm: Format error: 0xA6.
pbftoosm: Number of bytes read: 24
也就是說,這TM根本不是 osm官網的pbf格式啊 https://wiki.openstreetmap.org/wiki/PBF_Format
這幾天,為了編譯gdal和他的依賴庫, 兜了個大圈子,坑死我了.
2 解決
但是,這畢竟是個天天更新全球數據的網站, 德國人開的網站,總不至於發布完全無法讀取的數據呀.
又仔細看見下載頁面http://download.geofabrik.de/asia.html, 發現一行小字
Commonly Used Formats
- asia-latest.osm.pbf, suitable for Osmium, Osmosis, imposm, osm2pgsql, mkgmap, and others. This file was last modified 6 hours ago and contains all OSM data up to 2020-02-25T21:59:02Z. File size: 7.7 GB; MD5 sum:
好了, 直接搜第一個, 發現有py版本, 按https://pypi.org/project/osmium/ 安裝
sudo apt-get install build-essential cmake libboost-dev \ libexpat1-dev zlib1g-dev libbz2-dev pip3.8 install osmium
然后根據https://docs.osmcode.org/pyosmium/latest/intro.html
寫個新的py腳本read_osm_pbf.py
'''https://download.geofabrik.de/asia.html asia-latest.osm.pbf, suitable for Osmium, Osmosis, imposm, osm2pgsql, mkgmap, and others. ''' import osmium class CounterHandler(osmium.SimpleHandler): def __init__(self): osmium.SimpleHandler.__init__(self) self.num_nodes = 0 def node(self, n): self.num_nodes += 1 if __name__ == '__main__': fname = './data/osm_pbf/north-korea-latest.osm.pbf' h = CounterHandler() h.apply_file(fname) print("Number of nodes: %d" % h.num_nodes)
運行一下:
python3.8 read_osm_pbf.py
果然能讀, 沒報錯
下面繼續研究具體怎么讀坐標,怎么關聯shapely
有時間真想問問德佬, 你們這pbf和osm官網的pbf到底啥關系?!