Topics

Rasterio.features.geometry_window error when geojson has GeometryCollection


filthyhands@...
 



I'm attempting to isolate areas in a rasterio raster object, using Fiona to open a geojson and then window the parts I want to look at. 

Where a geojson entry is a geometry collection rasterio.features.geometry_window throws a key error. It's expecting each entry in the geojson to hold a 'coordinates' key but in this geojson the 'coordinates' are inside a list a level below 'geometryCollection'

A workaround could be to dismantle the geometry collections into separate objects. It would work but it's not ideal as the collections are there to group various styles of object together.

Does anyone have experience with this issue, or can suggest another solution? My code etc below.

Linestring: 
[{'type': 'LineString',
  'coordinates': [(144.937247997, -26.3759745069999, 0.0),

Geometry Collection: 
{'type': 'GeometryCollection',
  'geometries': [{'type': 'LineString',
    'coordinates': [(145.619783254, -38.5831492689999, 0.0),


My code:
with fiona.open('./Methane gas detection docs/NationalOnshoreGasPipelines.geojson') as fi: shapes = [each['geometry']for each in fi] # isolate a pipeline area pipes = rasterio.features.geometry_window(win, shapes, pad_x = 0.001, pad_y=0.001)



Error:
---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
<ipython-input-183-b9dff7a7db24> in <module>
     10 
     11     # isolate a pipeline area
---> 12     pipes = rasterio.features.geometry_window(win, shapes, pad_x = 0.001, pad_y=0.001)
     13 
     14     # print(shapes)

~\Anaconda3\envs\sat\lib\site-packages\rasterio\features.py in geometry_window(dataset, shapes, pad_x, pad_y, north_up, rotated, pixel_precision)
    419 
    420     if not rotated:
--> 421         all_bounds = [bounds(shape, north_up=north_up) for shape in shapes]
    422         lefts, bottoms, rights, tops = zip(*all_bounds)
    423 

~\Anaconda3\envs\sat\lib\site-packages\rasterio\features.py in <listcomp>(.0)
    419 
    420     if not rotated:
--> 421         all_bounds = [bounds(shape, north_up=north_up) for shape in shapes]
    422         lefts, bottoms, rights, tops = zip(*all_bounds)
    423 

~\Anaconda3\envs\sat\lib\site-packages\rasterio\features.py in bounds(geometry, north_up, transform)
    374 
    375     geom = geometry.get('geometry') or geometry
--> 376     return _bounds(geom, north_up=north_up, transform=transform)
    377 
    378 

rasterio\_features.pyx in rasterio._features._bounds()

KeyError: 'coordinates'


Brendan Ward
 

This looks like an oversight when we implemented bounds().

Rasterio currently supports FeatureCollection or individual geometry objects in this operation.

Logged here: https://github.com/mapbox/rasterio/issues/1914

Working on a quick fix now...


filthyhands@...
 

That's great Brendan, much appreciated. 


filthyhands@...
 

How can I patch my install of Rasterio to incorporate your changes? 


Brendan Ward
 

The target release date that will include this fix is 5/1/2020 (https://github.com/mapbox/rasterio/issues/1911)

If you need changes sooner, you have a few options:
1) clone the repo and checkout the maint-1.1 branch, and rebuild the Cython files.  This can be a bit involved depending on your OS / env; you'll need a compiler, GDAL, Cython, etc.  (https://github.com/mapbox/rasterio/blob/master/CONTRIBUTING.rst#development-environment)

2) coerce each GeometryCollection to FeatureCollection

3) use shapely to union each GeometryCollection into a single geometry object before passing into bounds()

4) calculate bounds on each individual geometry within a GeometryCollection, then union the bounds.  Here's how it was done in the Cython file (https://github.com/mapbox/rasterio/pull/1915/files#diff-83bc1acea1f8d435fd6c4798a58fc072R411-R426); you can do the same thing on your end.


filthyhands@...
 

Thanks again Brendan, very helpful.