read format .osm.pbf from geofabrik

#--------------------------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上的数据, 在

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__':

    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到底啥关系?!