I am trying to come up with a pipeline for turning very large OSM PBF files (eg buildings layer for a whole country: ~100M polygons) into a more manageable density estimate (with arbitrary precision).

So far my approach has been to convert the PBF file to a more usable GPKG format with ogr2ogr:

ogr2ogr -select building -overwrite -progress -gt 5000000 -f GPKG buildings.gpkg buildings.osm.pbf multipolygons

followed by a geopandas script that loads chunks (in parallel) and sums up polygon areas for each tile (using polygon boundary centre as a quick tile assignment estimate):

precision_scale = 10

def parse_slice(s):
    gdf = gpd.read_file(gpkg_path, rows=s, driver='GPKG', mode='r')
    snapped = gdf.bounds.apply(lambda b: (
        int((b.minx+b.maxx)/2 * precision_scale),
        int((b.miny+b.maxy)/2 * precision_scale)
    ), axis=1)
    df = pd.DataFrame({'area': gdf.to_crs(CRS.from_epsg('6933')).area,
                       'snapped': snapped}).groupby('snapped').sum().reset_index()
    return df

merged_df = pd.DataFrame()
with multiprocess.Pool(processes = 4) as pool:
    for df in pool.imap(parse_slice, slices):
        merged_df = merged_df.append(df)

df = merged_df.groupby('snapped').sum().reset_index()

grid = gpd.GeoDataFrame(df.area, geometry=df.snapped.apply(
    lambda pt: Point(pt[0]/precision_scale, pt[1]/precision_scale)
), crs=gdf_crs)

This works fairly well for medium-sized data (PBF file < 2GB), but ogr2ogr seems to choke on anything bigger (writes up to ~9GB of GPKG file, then hangs forever without any error message).

Is there any method/tool/algorithm that is more cpu and/or memory efficient that I could use to achieve a similar result?

  • Decoding OSM protobufs and serializing to OSM objects is well supported in the Python eco-system, and very fast, as are geometric operations on the resulting geometries; also considering your other question (creating bboxes), I would suggest to add some code to parse the PBFs yourself if you can, for maximum performance and memory management; note that the OSM model will likely require you to hold the nodes in memory at least temporarily, to being able to refer them from ways or relations...that's the bane of the data model, and all tools have to do it. – geozelot Nov 23 '20 at 19:52
  • Given the somewhat "exotic" nature of PBF (as opposed to plain one-polygon-by-line formats), it feels like writing my own parser would be error-prone (and overkill). Dealing with a stream of polygons would be very feasible (though not sure what tool would give me this), but I've noticed much better performances when taking advantage of geopandas' vectorising to load and treat larger chunks, rather than individual polygons… – Dave Nov 23 '20 at 20:03
  • @geozelot I am currently testing whether I can get ogr2ogr to do the whole computation in one single sqlite query (likely), but I am afraid the conversion to a sqlite DB will be too costly (not sure). It doesn't seem doable with the plain SQL dialect option. – Dave Nov 23 '20 at 20:29
  • OSM data and tools are somewhat oriented towards maintaing an up-to-date database primarily for mapping purposes; flexibility in the direction of common GIS concepts is costly. However, filtering and decoding, while memory/CPU intense, is fast. The parsing is mostly done by the decoder implementation, you would be able to work with lists and dicts. At least extract buildings only with e.g. osmium and then run ogr2ogr. – geozelot Nov 23 '20 at 20:53
  • @geozelot I omitted it from above for brevity, but I actually do filter out the building layer first using osmium. Even then, the remaining PBF is large. I think I have a possible solution based on ogr2ogr + sqlite. Will post it as a self-answer, but would love feedback or improvements… – Dave Nov 23 '20 at 21:58

Thanks to some feedback and further experimentation, I was able to come up with a single ogr2ogr command that seems to do the job:

ogr2ogr -dialect sqlite -sql "select \
MakePoint(ROUND(X(Centroid(Envelope(GEOMETRY))), 1),ROUND(Y(Centroid(Envelope(GEOMETRY))), 1)) as pt,\
 SUM(ST_Area(GEOMETRY, 6933)) AS area \
from multipolygons \
where building is not null \
GROUP BY pt" \
-progress -f GPKG ./output.gpkg ./input.osm.pbf 
  • The ROUND function (with decimal precision argument) is an extremely coarse way to set the tiling (in particular, will give off very different tile sizes for different lats). But not sure it is worth projecting back/forth.

  • Oddly enough, using Envelope(GEOMETRY) before computing Area did not improve performance (even decreased it slightly). Guessing that has to do with a lot of the geometries being fairly simple to begin with.

  • Not sure whether there is a simple way to avoid computing the Centroid twice, but obviously the cost is near-zero when using the envelope anyway.

  • Not sure there is any way to provide some progress feedback (command can take a few minutes for large files).

  • This seems an ideal candidate for parallelised processing. Not sure if worth the gain.

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.