How do I get the feature from a single-feature OGR layer in Python without looping?

Question

A while ago I asked why the 0 index in a OGR layer was not a feature and it turned out that OGR interprets the index as Feature ID and looks that up.

Now I filtered a layer to have only one feature in it. I would like to use that feature straight-away but cannot figure out how.

The setup:

from osgeo import ogr

driver = ogr.GetDriverByName("WFS")
wfs = driver.Open("WFS:https://geodienste.hamburg.de/HH_WFS_Statistische_Gebiete?REQUEST=GetCapabilities&SERVICE=WFS")  # This is a public WFS with open data
layer = wfs.GetLayerByName("Stadtteile")  # Some polygons, possibly multipolygons

# A random point inside one of them
point_in_hamburg = ogr.CreateGeometryFromWkt("POINT (566795 5935774)")

# Checking which features on layer intersect with our point
layer.SetSpatialFilter(point_in_hamburg)

# Due to the nature of Stadtteile (=districts) we
# will get only one intersecting feature:
print("One single feature: {}".format(len(layer) == 1))  # Prints "True"

Now I want to use that single feature. The only way I managed to is by using a loop over that length 1 layer:

for feature in layer:
    print(feature.GetFID())  # Prints "94"

I find that incredibly ugly and illogical. How can I directly access the feature instead without knowing anything about it?

layer[0] is not it. Strangely enough I can use layer[94] before the Spatial Filter and I get the feature (even though the WFS server returned a "Generic WFS service error") but if I try to use it afterwards, I get an IndexError from ogr. So even if the FID is 94 and is still 94 after filtering, layer[94] will fail at that point.

Is there some alternative?

Answer 1

Use layer.GetNextFeature()

from osgeo import ogr

driver = ogr.GetDriverByName("WFS")
wfs = driver.Open("WFS:https://geodienste.hamburg.de/HH_WFS_Statistische_Gebiete?REQUEST=GetCapabilities&SERVICE=WFS")  # This is a public WFS with open data

layer = wfs.GetLayerByName("Stadtteile")  # Some polygons, possibly multipolygons

# A random point inside one of them
point_in_hamburg = ogr.CreateGeometryFromWkt("POINT (566795 5935774)")

# Checking which features on layer intersect with our point
layer.SetSpatialFilter(point_in_hamburg)

# Due to the nature of Stadtteile (=districts) we
# will get only one intersecting feature:
print("One single feature: {}".format(len(layer) == 1))  # Prints "True"

feature = layer.GetNextFeature()
print(feature.GetFID())  # Prints "94"

I find that incredibly ugly and illogical.

I agree, but gdal/ogr is not a python library, it is a C++ library with python bindings. It is quite "unpythonic" and there are many "gotchas"