Rasterio: rasterize produces holes for lists of GeometryCollection

[astrojuanlu]

Expected behavior and actual behavior.

I expected rasterize(GeometryCollection([shp1, shp2])) to produce the same result than rasterize([GeometryCollection([shp1, shp2])]). Instead, the latter produces a hole in the overlap, regardless of the direction of the exterior of the two polygons.

Steps to reproduce the problem.

shp1 = Polygon.from_bounds(0, 0, 5, 5)
shp2 = Polygon.from_bounds(2, 2, 7, 7)

gc = GeometryCollection([shp1, shp2])

| Code | rasterize([shp1, shp2], (10, 10)) | rasterize(gc, (10, 10)) | rasterize([gc], (10, 10)) |
|--------|-----------------------------------------------------------------------------------------------------------------|-----------------------------------------------------------------------------------------------------------------|-----------------------------------------------------------------------------------------------------------------|
| Result | | | |

Operating system

Linux Mint 18.3

Rasterio version and provenance

• libgdal 1.11.2 installed with conda
• geos 3.6.2 installed with conda from conda-forge
• shapely 1.6.4.post1 installed with pip 9.0.1
• rasterio 1.0a12 manylinux1 wheel installed from PyPI using pip 9.0.1.

[sgillies]

Okay, rasterize([shp1, shp2]) and rasterize(gc) are the same because rasterize takes an interator over geometries and a shapely GeometryCollection is an iterator over its parts: list(gc) == [shp1, shp2]. I, too, would expect rasterize([gc]) to give the same results.

I think the only way this is not an upstream issue (in GDAL) is if there is a bug in Rasterio's OGR geometry factory specific to geometry collections. I'll take a look, but probably not until after 1.0, because we can document a work around for this in the short term.

@Juanlu001 would you be willing to see if updating GDAL to 2.2 fixes this?

Lastly, I'd like to point out that mixing conda packages, conda-forge packages, and the binary wheels I publish to PyPI is a perilous thing to do. My binary wheels are not compatible with conda. In general, I think C extension modules on PyPI can't be installed into a conda environment using pip.

[astrojuanlu]

I get the same behavior with GDAL 2.2. Just to make sure it's not a result of mixing binaries as you say, I installed rasterio from source to compile it against my libgdal:

$ conda install libgdal -c conda-forge
$ pip install --pre rasterio --no-binary :all:

And nothing changes.

Lastly, I'd like to point out that mixing conda packages, conda-forge packages, and the binary wheels I publish to PyPI is a perilous thing to do.

I know... "don't do this at home" :)

My binary wheels are not compatible with conda. In general, I think C extension modules on PyPI can't be installed into a conda environment using pip.

I do it all the time in Linux, and so far I have not found problems. Would you please point me to some past issue that was due to this?

[brendan-ward]

@Juanlu001 I checked in an x-failing test for this in the issue1253 branch.

I was able to isolate this from Shapely and reproduce the issue using GeoJSON dict objects.

It turns out that when we are iterating over a Shapely GeometryCollection object, it is equivalent to a list of the underlying features.

This explains why these are identical:

rasterize([shp1, shp2], (10, 10))
rasterize(gc, (10, 10))

I was able to reproduce the hole in the middle of the overlapping polygons in the geometry collection directly using gdal_rasterize, so my conclusion at this point is that it is either a bug or design choice in GDAL, and not a bug in rasterio.

To get around this, we could probably iterate over the contained geometries when we detect a geometry collection object (effectively same outcome as your rasterize(gc, (10, 10)) for a Shapely GeometryCollection). @sgillies any objections to that idea?

[astrojuanlu]

Thank you! 😊