Functional code suddenly breaking at rasterio.mask.mask(src, shapes)

[jonnyzwi]

Working code (as of two weeks ago) has stopped working. Neither the script itself nor the input has changed.

However, it appears that there was a new merge of rasterio into Main last week. Perhaps this is causing the sudden error?

After inspecting rasterio/_features.pyx in rasterio._features._bounds() that was mentioned in the error, my guess is that north_up is the empty sequence that ValueError: min() arg is an empty sequence is complaining about.

However, it is unclear as to how my call to rasterio.mask.mask(src, shapes) interacts with bounds() or _bounds(), or where north_up is supposed to be getting passed in via rasterio.mask.mask.

Any assistance would be appreciated.

ValueError                                Traceback (most recent call last)
[<ipython-input-18-a1f872635462>](https://localhost:8080/#) in <module>
     39 
     40 with rasterio.open(inputTIF) as src:
---> 41     out_image, out_transform = rasterio.mask.mask(src, shapes)
     42     out_meta = src.meta
     43 

4 frames
[/usr/local/lib/python3.8/site-packages/rasterio/features.py](https://localhost:8080/#) in bounds(geometry, north_up, transform)
    394         )
    395 
--> 396     return _bounds(geom, north_up=north_up, transform=transform)
    397 
    398 

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

ValueError: min() arg is an empty sequence

[snowman2]

@jonnyzwi are you able to provide code/files to reproduce this issue?

[jonnyzwi]

This portion of the code is really only a few lines long:

with fiona.open(parkingSHP, "r") as shapefile:
    shapes = [feature["geometry"] for feature in shapefile]
    shapes = [ i for i in shapes if i is not None]

with rasterio.open(parkingTIF) as src:
    out_image, out_transform = rasterio.mask.mask(src, shapes, crop=True)
    out_meta = src.meta

If you run the below in Google Colab, it will import all the necessary input files:

!pip install -q fiona
import fiona

!pip install -q rasterio
import rasterio
import rasterio.mask

import matplotlib.pyplot as plt
import numpy as np

### ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- 
!pip install --upgrade -q gdown
tncParking_zip = '1_67WOXuUoQYaF15-BbFE_0AXkpRx8rLa&confirm=t'
!gdown $tncParking_zip

!unzip TNC_PARKING.zip
!rm -r /content/TNC_PARKING.zip

parkingSHP = 'TNC_PARKING/yesParkingLots.shp'
notParkSHP = 'TNC_PARKING/notParkingLots.shp'

fullRaster = '17jmAmQD_UM6QtnFIazDVM0QGF1i-SbY3&confirm=t'
!gdown $fullRaster
parkingTIF = 'fullRasterMerge.TIF'

### ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- 

with fiona.open(parkingSHP, "r") as shapefile:
    shapes = [feature["geometry"] for feature in shapefile]
    shapes = [ i for i in shapes if i is not None]

with rasterio.open(parkingTIF) as src:
    out_image, out_transform = rasterio.mask.mask(src, shapes, crop=True)
    out_meta = src.meta

[jonnyzwi]

Hi @snowman2 — any ideas on this problem?

[snowman2]

Unfortunately, I haven't had a free moment to take a look. I haven't used Google Colab before and so it will likely take longer to dig in. If you can upload the files directly here, that may help save time.

[vincentsarago]

@jonnyzwi I just tested, and it worked for me 🤷‍♂️

print(rasterio.__version__)
print(rasterio.__gdal_version__)
print(rasterio.__geos_version__)
print(rasterio.__proj_version__)

Here are the files @snowman2 (note I created a COG from the tif because it was 600Mb)
TNC_PARKING.zip
fullRasterMerge_cog.tif.zip

[snowman2]

Thanks @vincentsarago 👍. @jonnyzwi mind providing the output of rasterio.show_versions()?

[jonnyzwi]

@snowman2 — indeed I had not uploaded because the files were too large.

Here is a link to a colab notebook containing the cnode snippet, which you can simply run via the browser

https://colab.research.google.com/drive/1tt-NRhPXvgsF9Y3nxWeco2ddCehZeN7R

This notebook includes your request to print

print(rasterio.__version__)
print(rasterio.__gdal_version__)
print(rasterio.__geos_version__)
print(rasterio.__proj_version__)

which result in

1.3.5.post1
3.5.3
3.11.1
9.0.1

@vincentsarago — did you successfully execute the code in a colab notebook, or elsewhere?

I should note that I also attempted to run my code with version 1.3.4 (which was the latest version when this code last worked) by calling !pip install rasterio==1.3.4 at the top of my script. That too had failed.

[groutr]

The bounds computation is happening because you have crop=True. Rasterio will compute the bound of the shapes and crop the raster to that window.

The line that is throwing the exception: https://github.com/rasterio/rasterio/blob/1.3.5.post1/rasterio/_features.pyx#L413

North_up has no relevance here. In fact, @snowman2, rasterio.features.geometry_window needs to be issuing deprecation warnings for a couple keyword arguments (per the docstring).

Inspecting the geointerface of your first shape, we see

{'coordinates': [[(-115.11462306899995, 36.184749054000065),
   (-115.11462315699998, 36.18473065500007),
   (-115.11462347799994, 36.18466354700007),
   (-115.11462245399997, 36.18459860100006),
   (-115.11462144599994, 36.18453040600008),
   (-115.11462043099999, 36.18446329400007),
   (-115.11461941099998, 36.18439726400004),
   (-115.11461973899998, 36.18432907500005),
   (-115.11498695699999, 36.18433239500007),
   (-115.11590297999999, 36.18434501200005),
   (-115.11590642199997, 36.18470146800007),
   (-115.11590630299997, 36.18470616600007),
   (-115.11590421599999, 36.18471701800007),
   (-115.11590003199996, 36.18472746400005),
   (-115.11589385599996, 36.18473723900007),
   (-115.11588584599997, 36.184746096000026),
   (-115.11587620299997, 36.184753811000064),
   (-115.11586517299997, 36.18476018800004),
   (-115.11585303499999, 36.18476506600007),
   (-115.11584009599994, 36.184768321000035),
   (-115.11582668299997, 36.18476987100007),
   (-115.11582160699999, 36.18476991500006),
   (-115.11581991099996, 36.18476989100003),
   (-115.11462306899995, 36.184749054000065)]],
 'type': 'Polygon',
 'geometries': []}

The presence of that geometries key is throwing _bounds down the wrong branch. Removing that key results in the correct bounds computation. I don't believe any of the code in this call stack has changed in a couple of years.

There are two possibilities I can think of:

1. Fiona (or one of its dependencies) was updated and somehow 'geometries' is now being added to fiona geometry object's geointerface.
2. The shapefile changed or is somehow now being interpreted by fiona differently.

Solutions:

1. perhaps the logic in _bounds is outdated? Should it check for empty containers?

elif 'geometries' in geometry and geometry['geometries']:
    ...

[snowman2]

to rasterio.features.geometry_window needs to be issuing deprecation warnings for a couple keyword arguments (per the docstring).

Good catch. A PR or issue to track is welcome.

[snowman2]

@groutr thanks for digging in. You are likely correct. Fiona 1.9 was recently released with new interfaces to Features/Geometry.

https://github.com/Toblerity/Fiona/blob/master/CHANGES.txt

[jonnyzwi]

@groutr thank you for catching that. The problem occurs even when crop is not set but, indeed, it seems geometries is causing the error.

I have filtered out the geometries item from each shape

with fiona.open(parkingSHP, "r") as shapefile:
    shapes = [ i['geometry'] for i in shapefile]
    shapes = [ i for i in shapes if i is not None]
    shapes = [ [i for i in shape.items() if 'geometries' not in i] for shape in shapes]

Unfortunately, a new issue arises:

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
[<ipython-input-100-4618ef0b5d41>](https://localhost:8080/#) in <module>
     13 # list(shape1.items())
     14 with rasterio.open(parkingTIF) as  src:
---> 15     out_image, out_transform = rasterio.mask.mask(src, shapes)
     16     out_meta = src.meta
     17 print(f'{type(shapes) = }')

4 frames
[/usr/local/lib/python3.8/dist-packages/rasterio/features.py](https://localhost:8080/#) in bounds(geometry, north_up, transform)
    385         return tuple(geometry['bbox'])
    386 
--> 387     geom = geometry.get('geometry') or geometry
    388 
    389     # geometry must be a geometry, GeometryCollection, or FeatureCollection

AttributeError: 'list' object has no attribute 'get'

[groutr]

@jonnyzwi yeah, that's my mistake, it computes the bounds all the time.
https://github.com/rasterio/rasterio/blob/1.3.5.post1/rasterio/features.py#L447

You need to remove the geometries key like this:

shapes = []
for i in shapefile:
    shape = i['geometry']
    if shape is not None:
        geom = dict(shape.__geo_interface__)
        del geom['geometries']
        shapes.append(geom)  # <-- shapes should be a list of dictionaries

[jonnyzwi]

@groutr Thank you! Things appear to be working again.

[sgillies]

This is a Fiona bug, introduced in 1.9.0: Toblerity/Fiona#1197.

[sgillies]

Wait, we can also fix this in rasterio for 1.3.6. I'll do that now.