Is there an efficient way to convert a set of node/relations in a PBF file, to the bounding boxes for each polygon, without going through a conversion to another GIS format?

Right now, my pipeline consists of:

  • Extract nodes/relations of interest (with osmium)
  • Convert to gpkg (costly)
  • Load with geopandas (costly for large files)
  • Extract bounding box for each polygon and save new file

Is there by any chance a program that could work directly with PBF?

Context: I am trying to supplement existing landusage data with building data, in order to then run quick estimates on land use, for any arbitrary region… But do not need the precision (and ensuing size/speed hit) of having exact building polygons, so would instead like to extract only bounding boxes and use that.

  • What is the data format after step 1 in your pipeline? If you have the raw coodinates there it would be possible to just extract the min and max values per feature and be done. – bugmenot123 Nov 22 '20 at 13:00
  • @bugmenot123 osmium generally filters out to another (smaller) PBF file. But even after exporting to gpkg, parsing the file to extract the bounds, while not extraordinarily difficult, is not something I'd want to do by hand… – Dave Nov 22 '20 at 13:50
  • I'd pipe the osmium_extract into osmium_export with boudary filter and import into PostgreSQL (you can get a GeoJSON, too); I know you explicitly want to use no intermediate storage, but its trivial to go from there. – geozelot Nov 22 '20 at 14:26
  • @geozelot intermediate storage isn't a big deal (though having to write/load/parse large files might be), but unless I'm missing something, the boundary filter on osmium_export will not result in a file that contains bounding boxes of the original polygons… – Dave Nov 22 '20 at 14:33
  • duplicate of : stackoverflow.com/questions/44975952/… ? – Spacedman Nov 23 '20 at 15:19

This way one command is enough but it may not be extremely fast because GDAL must in any case create temporary database structures from the whole .pbf file for resolving nodes-ways-relations.

ogrinfo -dialect sqlite -sql "select ST_Envelope(geometry) from multipolygons where building is not null limit 5" albania-latest.osm.pbf

Had to open data source read-only.
INFO: Open of `albania-latest.osm.pbf'
      using driver `OSM' successful.

Layer name: SELECT
Geometry: Unknown (any)
Feature Count: 5
Extent: (19.470962, 40.289919) - (20.153678, 41.329786)
Layer SRS WKT:
Data axis to CRS axis mapping: 2,1
Geometry Column = ST_Envelope(geometry)
  POLYGON ((19.4709618 40.7198672,19.4713055 40.7198672,19.4713055 40.7201239,19.4709618 40.7201239,19.4709618 40.7198672))

  POLYGON ((20.153198 40.297601,20.153678 40.297601,20.153678 40.2979293,20.153198 40.2979293,20.153198 40.297601))

  POLYGON ((19.8165231 41.3290525,19.8177414 41.3290525,19.8177414 41.3297856,19.8165231 41.3297856,19.8165231 41.3290525))

  POLYGON ((19.8073262 41.3292145,19.8078217 41.3292145,19.8078217 41.3296297,19.8073262 41.3296297,19.8073262 41.3292145))

  POLYGON ((20.0099954 40.2899189,20.0104853 40.2899189,20.0104853 40.2904093,20.0099954 40.2904093,20.0099954 40.2899189))

Edit the SQL if you want to include more attributes and use it with ogr2ogr when you are happy with it.

You did not tell the method that you used for parsing the bounding boxes from GeoPackage, but SQLite dialect and ST_Extent() can be used similarly with GeoPackage data.

  • As you suspected, the DB creation seems to add quite a cost on large files. But I've also come to realise that trying to work directly with the PBF file is bound to fail (no matter what, there needs to be a full pass to generate the individual polygon rows). Nonetheless, this helped put me toward the right direction! – Dave Nov 23 '20 at 8:36

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.