GDAL/OGR losing OSM closed way features that OGR2OGR imports

Question

I have managed to import an OSM PBF file into QGIS using ogr2ogr -> SQLite and then importing layers using the GUI- (essentially as per the first step in this answer :

ogr2ogr -f "SQLite" -dsco SPATIALITE=YES texas.db texas-latest.osm.pbf

My version of GDAL and OGR2OGR is 2.4.1 (as reported by both - I checked)

That appears to work fine. I'm trying to duplicate this in a Python plugin without the SQLite intermediate step. Points, lines, and multilines appear to be readable fine (exact same numbers of features). However, the total number of multipolygon features is significantly less. Further investigation shows the missing features are things like buildings which have an osm_way_id set. I.e. these have been encoded in the PBF file as closed ways.

Here is an example imported with ogr2ogr:

The GDAL library code does not see any features with osm_way_id set! (although it recognizes this is a valid field name). Here is my code:

    gdal.SetConfigOption("OGR_INTERLEAVED_READING", "YES")    

    input_pbfname = "C:/Data/QGIS/Texas/texas-latest.osm.pbf"
    ds = gdal.OpenEx(input_pbfname, gdal.OF_VECTOR)
    if ds is None:
        QMessageBox.information(None, "open failed", "Open failed")

    nly = ds.GetLayerCount()

    #version_num = int(gdal.VersionInfo('VERSION_NUM'))
    #QMessageBox.information(None, "Version",  str(version_num))

    nPoints = 0
    nLines = 0
    nMultiLines=0    
    nPolygons = 0

    # GetNextFeature handles interleaved reading correctly 
    for ily in range (0,nly):
        layer = ds.GetLayerByIndex(ily)

        layer.ResetReading()
        feat = layer.GetNextFeature()
        while feat is not None:
            if layer.GetName() == "points":
                nPoints = nPoints+1
            elif layer.GetName() == "lines":
                nLines = nLines+1
            elif layer.GetName() == "multilinestrings":
                nMultiLines = nMultiLines+1
            elif layer.GetName() == "multipolygons":
                nPolygons = nPolygons+1
                idx = feat.GetFieldIndex("osm_way_id")
                if idx>=0:
                    wayid = feat.GetFieldAsString(idx)
                    if (wayid is not None and wayid!="" and wayid!="0" ):
                        nWayId=nWayId+1

            idx = feat.GetFieldIndex("osm_way_id")
            if idx>=0:
                # tried both GetFieldAsInteger and GetFieldAsString
                # also tried counting all multipolygon features with osm_way_id set
                wayid = feat.GetFieldAsInteger(idx)
                if wayid == 646651935:
                    QMessageBox.information(None,"646651935 building found","Layer:"+layer.GetName() )

            feat = layer.GetNextFeature()

    QMessageBox.information(None, "Count", "Points: %d Lines:%d MultiLine:%d  Poly:%d" % (nPoints,nLines,nMultiLines, nPolygons) )


    # close the temporary layer
    ds = None

So why are these "closed ways" getting imported correctly by ogr2ogr but not the GDAL api? Is there a flag or setting that I'm missing? Or another subtlety with the interleaving that I have missed? (I'm aware of the reasons for this with OSM data - in the past I've coded an OSM PBF reader in C#...)

ogr2ogr is importing 3385834 multipolygons, but the gdal API is only importing 10896. The other layers (point, line, multiline all match)

Answer 1

Found the issue. Specifically this code:

# GetNextFeature handles interleaved reading correctly 
for ily in range (0,nly):
    layer = ds.GetLayerByIndex(ily)

I initially replaced it with:

# GetNextFeature handles interleaved reading correctly 
for ily in (0,1,2,4,3):
    layer = ds.GetLayerByIndex(ily)

...which appeared to work. It ensures Layer 4 ("other relations") is read before layer 3 (multi-polygons). It was clear the points and lines/linestrings had to be read first, but I thought the missing polygons (closed ways) weren't defined with relations.

However there can be other repeats of layer feature types - ie. there should be a third loop surrounding the layer loop. This should keep repeating the original range() layer loop until no more features are read.

---- edit ----

Starting with GDAL 2.2 there is a much simpler way to do this. That is to use the dataset's GetNextFeature() call. This can be used in a single level loop, eg:

        feat,layer = ds.GetNextFeature()

        while feat is not None:
            num_features = num_features+1
            geom = feat.GetGeometryRef()
            feat,layer = ds.GetNextFeature()

Note that features can still arrive out of order - there is no guarantee all points will be together,etc. Also, dataset.GetNextFeature() can also return a "percentage complete" value - see the docs for details.

---- end_of_edit ----

Although it wasn't the eventual cause, during my investigations, I found the osmconf.ini file wasn't properly defined. Or it was difficult to tell which was being used. This can be set explicitly with:

    gdal.SetConfigOption("OSM_CONFIG_FILE", "c:/my/path/to/osmconf.ini")

In my case, I point it to the installed plugin directory - ie. local to my plugin code. This can then be packaged with the installer. I have control over it, and the end user is free to modify it if they so wish.