Fiona: Random access reads for Shapefile

[mccarthyryanc]

Is there way within Fiona to get records based on an index? I understand this is driver dependent, but does Fiona take advantage of the indexing that comes with Shapefiles (.shx)? Perhaps something like:

with fiona.open('data.shp', 'r') as src_dataset:
    record = src_dataset.get_record(42)

[geowurster]

@mccarthyryanc Kind of. Fiona supports Python's __getitem__() protocol to get features by their feature ID, which is often the index but not guaranteed to be.

with fiona.open('data.shp', 'r') as src_dataset:
    record = src_dataset[42]

However, slices operate on feature indexes so this is a safer (and uglier) method for grabbing a single feature by its index, although we should add a method for doing just this:

with fiona.open('data.shp', 'r') as src_dataset:
    record = src_dataset[42:43][0]

It looks like the ID is always the index for shapefiles. Fiona doesn't support setting a feature's ID directly like OGR, but the code below illustrates how to do it with OGR's Python bindings. SQLite supports setting the feature ID directly, so the code below prints 2 but if you change the driver to ESRI Shapefile it prints 0.

from osgeo import ogr

# Create a layer with a single feature
# Set the feature's ID to 2
# The feature does not have any properties or geometry
driver = ogr.GetDriverByName('SQLite')
ds = driver.CreateDataSource('data.sqlite')
layer = ds.CreateLayer('data', geom_type=ogr.wkbPoint)
feature_defn = layer.GetLayerDefn()
feature = ogr.Feature(feature_defn)
feature.SetFID(2)
layer.CreateFeature(feature)
layer = None
ds = None
driver = None

# Read the first (only) feature and print its ID
ds = ogr.Open('data.sqlite')
layer = ds.GetLayer()
feature = layer.GetNextFeature()
print(feature.GetFID())
layer = None
ds = None

[mccarthyryanc]

@geowurster Thanks for the clarification. I'll stick with the uglier and safer slicing operation for now.

[snorfalorpagus]

@geowurster I find it very surprising/confusing that the behaviour of slices is different to that of a single index. Is this something that needs to be improved? Looking at the documentation, OGR_L_SetNextByIndex isn't particularily efficient anyway. If possible we should provide methods to access features both by index or by FID. I'm not sure if __getitem__ should work on index or FID, but it should be consistent.

[sgillies]

@snorfalorpagus @geowurster @mccarthyryanc I'm not as certain about OGR FIDs as I should be, but my understanding of the OGR data model is that they are not guaranteed to be anything more than unique (within a layer) integer identifiers for features. Details may vary from format to format. For example, in the Shapefile format reference (http://www.gdal.org/drv_shapefile.html)

Also, .dbf files are required to have at least one field. If none are created by the application an "FID" field will be automatically created and populated with the record number.

The OGR shapefile driver supports rewriting existing shapes in a shapefile as well as deleting shapes. Deleted shapes are marked for deletion in the .dbf file, and then ignored by OGR. To actually remove them permanently (resulting in renumbering of FIDs) invoke the SQL 'REPACK ' via the datasource ExecuteSQL() method.

Starting with GDAL 2.0, REPACK will also result in .shp being rewritten if a feature geometry has been modified with SetFeature() and resulted in a change of the size the binary encoding of the geometry in the .shp file.)

Deleted features may remain in a shapefile, holding on to their FID, which means that iterating over a shapefile and printing the FID of features can result in a sequence with gaps.

Making Collection not directly sliceable could improve item access consistency, yes? We'd keep the existing dict-like semantics of item access, which is that you give an integer OGR FID like colxn[42] and you get the feature with OGR FID 42, not the 43rd record from an abstract list of records.

The results of list(colxn) and list(colxn.filter()) would remain sliceable, of course, and it would remain possible to import islice from itertools and do islice(colxn, start, stop).

[geowurster]

@sgillies @snorfalorpagus @mccarthyryanc I should have mentioned that I discovered this behavior when I started looking at the code to make sure I understood how Fiona handles __getitem__(), so I was surprised too.

Another potential issue for our FID handling is that we explicitly cast them to string when building a feature. GDAL produces a GIntBig. @sgillies is this GeoJSON related? I don't see a type requirement for id in the current spec on geojson.org but haven't checked the draft for the new spec.

I didn't know deleted features could remain in a shapefile so modifying __getitem__()'s behavior as defined above seems appropriate, especially since I'm sure more users will be bringing a list of needed FID's rather than a range of features.

[sgillies]

For 1.8.0 I plan to add a get() method to a collection, so that we can do, for example

with fiona.open('tests/data/coutwildrnp.shp') as src:
    print(src.get('0'))

In 2.0 we could change the behavior of __getitem__() to be consistent with this.

[snorfalorpagus]

I'm unclear on what the proposal is here. I think making Collection not directly sliceable is a good idea, as the current behavior is inconsistent. @sgillies What would src.get('0') return, if different from src[0]? Do we have a long term plan for int vs str FIDs?

[sgillies]

@snorfalorpagus src.get(fid) should return the feature with that integer fid using OGR_L_GetFeature (http://www.gdal.org/ogr__api_8h.html#abdad5bebd6b71d136ce21e70ef4c63c7). src[fid] should return the same thing in 2.0.

I think we (I, in particular) need to give up on string FIDs. I don't see any sign of OGR ever supporting them across the board.

[snorfalorpagus]

Another example of where this is an issue is GPKG, where FIDs start at 1 not 0.

I've also noticed that accessing an invalid index, e.g. collection[999999] returns None rather than raising an KeyError.

We already have src.session.get_feature(1) so src.get could be an alias?