1

I often use GDAL/OGR via their python bindings.

I know of the very useful ogr.Layer.SetAttributeFilter() command. However I am hoping to filter by feature ID's, rather than attributes. Also, I don't want to filter by just one value, but rather by a set of values. Ideally I could pass a list of feature ID's to the filter -- is there some method that can allow me to do this (either through SetAttributeFilter or otherwise)?

It occurred to me that I could copy the feature ID's to their own attribute field, but it would be nice to avoid having to do this.

CLARIFICATION: I am hoping to apply filters for read-only processes which would involve sorting, selecting & exporting/rasterizing the layer geometries. I primarily use ESRI Shapefile format.

  • It depends.. what is your feature source? If it's a shapefile this is a dangerous course to pursue, you'd be much better served to copy the features' ID to a new (static) field. – Michael Stimson Feb 8 '18 at 2:25
  • Yes, I am taking the features from a shapefile. Thanks for the input -- could you elaborate on why this would be a 'dangerous' move? – corvus Feb 8 '18 at 2:32
  • 1
    Feature IDs in a shapefile aren't static.. they can change at any time. – Michael Stimson Feb 8 '18 at 2:34
  • Interesting, I'm surprised I wasn't aware of that! Are they even liable to change if I am not doing any editing? I only need to run a handful of filters & then rasterize portions of my layer geometry to separate raster files. – corvus Feb 8 '18 at 2:43
  • IDs are always 0 based, unique and ascending, no they're not likely to change if you're not editing but the possibility exists that you (or someone else) could edit one attribute which would have a ripple effect on your query. The trouble is that $rownum isn't a real field - it's not in the DBF file (open it in Excel and you'll see) so it's not available for some functions. However if your feature source is file/personal geodatabase the OBJECTID is unchanging and a real field to base a query on. – Michael Stimson Feb 8 '18 at 2:49
1

For an ESRI Shapefile:

This can be accomplished using the ogr.layer.SetAttributeFilter() method. Feature ID's in a shapefile are stored in the FID field. This field cannot be accessed via OGR's GetField method (it is not a 'proper' attribute field in that it is not stored in the .dbf file), however it is apparently accessible through SQL queries. I was able to do this as follows:

## Open shapefile geometry
driver = ogr.GetDriverByName('ESRI Shapefile')
datasource = driver.Open(shapefile_path, 0)
layer = datasource.GetLayer()

## Store ID's of desired features in a list
fid_list = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]

## Pass ID's to a SQL query as a tuple, i.e. "(1, 2, 3, ...)"
layer.SetAttributeFilter("FID IN {}".format(tuple(fid_list)))

## The filter is now active & we can work with the layer as we like!

## For example, we can rasterize the layer & only the features
## that we specified in 'fid_list' will be included in the output

## Initialize output raster
dataset_out = gdal.GetDriverByName('GTiff').Create(
                   output_path, cols, rows, 1, gdal.GDT_Byte)
dataset_out.SetGeoTransform(gt)
dataset_out.SetProjection(srs)

## Rasterize filtered layer
gdal.RasterizeLayer(output_path, [1], layer,
    options=["ATTRIBUTE={}".format(field_name), "ALL_TOUCHED=TRUE"])

NOTE: the FID field of a shapefile can change under some circumstances. For example, if a feature is deleted, the feature ID's will be automatically renumbered so that the numbering is sequential and without gaps. So, this approach would probably not be advisable unless the file is on your local network & you do not intend to do any editing while the filter is in place (ideally with the shapefile opened in read-only mode).

Of course, you could just create a new static attribute field which could hold the FID values or some other unique identifier. The approach detailed above would still work in this case.

0

The ogr.Layer.GetFeature(i) function returns the feature with index (FID) i.

So you could get the feature thorugh a loop as follows:

FIDs = [1, 10, 20, 50, ...]  # FIDs example

for fid in FIDs:
    lyr.GetFeature(fid)
    # Keep this feature

Alternatively, when looping through a Layer object, the next avaiable feature (and thus the next FID) is accesed every time. You could then loop through the layer and control what features you are going to access as follows:

FIDs = [1, 10, 20, 50, ...]  # FIDs example

for i, feat in enumerate(lyr, start=1):
    if i in FIDs:
        # Keep this feature
  • Do you think you could elaborate on how I would go about passing the features on to some process, e.g. gdal.RasterizeLayer() or gdal.Rasterize()? – corvus Feb 9 '18 at 1:22
  • You could copy the features to a new shapefile and then use either of the Rasterize functions. It is probably not the most appropiate solution but it should work. – Marcelo Villa Feb 9 '18 at 1:45

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.