项目需要分析过滤osm上的数据,  在

http://download.geofabrik.de/

可以下载到按大洲,和国别 的 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__':
    #driver=ogr.GetDriverByName('OSM')
    #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到底啥关系?!