rasterio.features.shapes with holes in polygons
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:
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):
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!
Geospatial software developer
toggle quoted messageShow quoted text
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:
On Tue, Dec 15, 2020 at 6:37 PM Sean Gillies via groups.io
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
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 :)
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 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
# 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 = shapely.geometry.MultiPolygon(polygons)
# polygons = polygons.apply(make_valid)
polygons = gpd.GeoSeries(polygons)
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.
Can you provide the data??--
Md Fahim Hasan
Graduate Research Assistant
Missouri University of Science and Technology
The image was included in the original post! Thanks.