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@...> 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


wsdaphne.wang@...
 

Hi,

I am encountering a similar problem with getting polygons with holes from a raster image. Hoping someone can help.
The image data that I am working with is a multi-page label matrix in tiff format, I am working with the first page. I have attached it for reprex. My data has no crs.

from skimage import io
import numpy as np
from rasterio import features
import shapely
import geopandas as gpd
import matplotlib.pyplot as plt

# read the first page of the multipage tiff
img = io.imread(seg_map_path, img_num=0)
img_array = np.array(img).astype('int32')

# display the label matrix; value=0 is displayed as purple; (0,0) is in the top left corner
plt.imshow(img_array)
plt.show()


# attempt to get the multipolygon for value=0
(purple region)
mask = (img_array == 0)
polygons = []
for shape, value in features.shapes(img_array, mask):
    polygons.append(shapely.geometry.shape(shape))
polygons = shapely.geometry.MultiPolygon(polygons)
# polygons = polygons.apply(make_valid)
polygons = gpd.GeoSeries(polygons)
polygons.plot()

The resulting geoseries is from the purple geometry from above, reflected across the y axis (because (0,0) is now in the lower left corner), but with no holes.

I would really appreciate any help on this.


Fahim Hasan
 

Can you provide the data??--
Md Fahim Hasan 

Graduate Research Assistant
Geological Engineering
Missouri University of Science and Technology


wsdaphne.wang@...
 

The image was included in the original post! Thanks.