Topics

rasterio.features.shapes with holes in polygons


Paolo Corti
 

Hello users and devs

First: thanks a lot for this wonderful project, I am really enjoying using it for my geospatial needs.

For a specific process I am using rasterio.features.shapes and it is working greatly for my purpose. However I have a few cases where I am experiencing strange results.
For example, see this 1 band raster: 

Screen Shot 2020-12-15 at 5.33.13 PM.png
I would like to extract the polygons for the yellow area (value: 3).

When running the following code:

with rasterio.open(cover_path) as src:
cover_arr = src.read(1)
mask = (cover_arr == 3)
polygons = shapes(cover_arr, mask=mask, transform=src.transform)

I get these results (in blue):

Screen Shot 2020-12-15 at 5.35.37 PM.png

Apparently holes are not correctly parsed.

It would be great to know if I am doing something wrong here or if it is something which is not working properly.

For your reference, I am attaching a copy of the raster I am using here so you can reproduce this (sample_cover.tif)

Thanks in advance!
Paolo

--
Paolo Corti
Geospatial software developer
web: http://www.paolocorti.net
twitter: @capooti
skype: capooti
#drt3jc1


Sean Gillies
 

Hi Paolo,

Welcome! I ran the your file through rio-shapes and jq

rio shapes --bidx 1 ~/Downloads/sample_cover.tif | jq -c 'select(.properties.val == 3.0)' | fio collect

and then uploaded to a gist:


It looks like the holes aren't lost there. How did you make the second image?

On Tue, Dec 15, 2020 at 3:47 PM Paolo Corti <pcorti@...> wrote:
Hello users and devs

First: thanks a lot for this wonderful project, I am really enjoying using it for my geospatial needs.

For a specific process I am using rasterio.features.shapes and it is working greatly for my purpose. However I have a few cases where I am experiencing strange results.
For example, see this 1 band raster: 

Screen Shot 2020-12-15 at 5.33.13 PM.png
I would like to extract the polygons for the yellow area (value: 3).

When running the following code:

with rasterio.open(cover_path) as src:
cover_arr = src.read(1)
mask = (cover_arr == 3)
polygons = shapes(cover_arr, mask=mask, transform=src.transform)

I get these results (in blue):

Screen Shot 2020-12-15 at 5.35.37 PM.png

Apparently holes are not correctly parsed.

It would be great to know if I am doing something wrong here or if it is something which is not working properly.

For your reference, I am attaching a copy of the raster I am using here so you can reproduce this (sample_cover.tif)

Thanks in advance!
Paolo

--
Paolo Corti
Geospatial software developer
web: http://www.paolocorti.net
twitter: @capooti
skype: capooti
#drt3jc1



--
Sean Gillies


Paolo Corti
 

On Tue, Dec 15, 2020 at 6:37 PM Sean Gillies via groups.io
<sean=mapbox.com@groups.io> wrote:

Hi Paolo,

Welcome! I ran the your file through rio-shapes and jq

rio shapes --bidx 1 ~/Downloads/sample_cover.tif | jq -c 'select(.properties.val == 3.0)' | fio collect

and then uploaded to a gist:

https://gist.github.com/sgillies/96ea784a00540d174e02060dde9e1a5a

It looks like the holes aren't lost there. How did you make the second image?

Hi Sean!

Thanks a lot for your hint. My problem was caused by the fact that I
was merging the polygons returned by rasterio.features.shape using
unary_union. I changed the approach and I created a multipolygon from
the returned polygons using shapely.geometry.shape and this made the
trick.
However one thing which I have noticed is that the returned geometry
is invalid in some cases. So I had to buffer to avoid errors when writing to a
shapefile. Is this a good approach or you suggest something different?

MultiPolygon([shape(geom) for geom in polygons].buffer(0)

Thanks again for your help and for this wonderful project - and for
its siblings too :)


Paolo