3

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?

7

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"

  • 1
    The layer instance doesn't "know" that there's only a single feature and doesn't read all of its features into Python objects, so it makes sense to require a "next" call. If you wanted to get your single result by using a 0 index, do it like: feature = list(layer)[0]. – sgillies May 21 '16 at 16:12

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.