Wrong features when exporting vectorized raster data to a file using GeoPandas after applying a filter and a dissolve operation

Question

I have a raster file with black areas having a value of 0 and white areas having a value of 1, resulting from a visibility analysis:

I vectorized this raster after loading it with rasterio:

import rasterio # version 1.2.6
from rasterio.features import shapes
import geopandas as gpd # version 0.9.0

with rasterio.Env():
    with rasterio.open('/path/to/geotiff_file.tif') as src:
        dataset = src.read(1)
        results = ({'properties': {'raster_value': v}, 'geometry': s}
                  for i, (s,v) in enumerate(
                          shapes(dataset,
                                 mask=None,
                                 transform=src.transform)
                          )
                  )

geometries = list(results)
gdf_polygons = gpd.GeoDataFrame.from_features(geometries)

Here's a sample of this GeoDataFrame:

In [4]: gdf_polygons.head()
Out[4]: 
                                            geometry  raster_value
0  POLYGON ((2679580.000 1342480.000, 2678980.000...           1.0
1  POLYGON ((2678840.000 1345340.000, 2677140.000...           1.0
2  POLYGON ((2681020.000 1344100.000, 2681220.000...           0.0
3  POLYGON ((2673260.000 1347900.000, 2678360.000...           1.0
4  POLYGON ((2682040.000 1344140.000, 2681240.000...           0.0
...

I export this result to a GeoPackage, which I then load in QGIS and where I set up the symbology so that it fits the different values of the raster_value attribute (which can be: 0, 1 or -99999):

Here the green color matches the black regions (0) of the original raster and the purple areas the white ones (1). Perfect.

But then, I only want to export the geometries where raster_value == 1 (white/purple):

dis = gdf_polygons[gdf_polygons['raster_value']==1].dissolve()
dis.to_file('/path/to/output.gpkg'), driver='GPKG')
# I also tried using sp_ops.unary_union() from shapely 1.7.1 but the same problem occurs

And here's the resulting GeoPackage, exported polygons are depicted in light purple here:

It seems correct, but...

...when overlaying this result with the original raster using transparency on the former:

The values of the "dark" purple regions are 1 on this new vector result but they were actual 0 on the original raster.

Why the differences and how to fix this problem so that, when exporting features by applying a filter where raster_value==1 + a dissolve operation, the result perfectly matches the original raster?

Update (2021-07-09) based on martinfleis' comment:

Apparently, not all geometries were valid:

In [32]: gdf_polygons.is_valid.all()
Out[33]: False

In [34]: gdf_polygons[~gdf_polygons.is_valid]
Out[34]: 
                                              geometry  raster_val
148  POLYGON ((2680060.000 1110880.000, 2680060.000...         0.0
241  POLYGON ((2680040.000 1110220.000, 2680740.000...         0.0
448  POLYGON ((2673000.000 1110060.000, 2670300.000...         1.0
462  POLYGON ((2679040.000 1110040.000, 2670340.000...         0.0
504  POLYGON ((2682080.000 1110540.000, 2680580.000...         0.0
625  POLYGON ((2678020.000 1100040.000, 2670920.000...         1.0
...

And when plotting these polygons I always notice the same "loop" pattern (circled in red):

gdf_polygons[~gdf_polygons.is_valid].plot()

Here a zoom on some of the 0 green polygons (:=black regions on the raster):

And a closeup:

Using gdf_polygons.geometry = gdf_polygons.buffer(0) apparently "did something with" the invalid geometries as now gdf_polygons.is_valid.all() evaluates to True, but when exporting the dissolved GeoPackage to a file and when importing it in QGIS, I have the same wrong shapes than in the figures above. Strangely.

And second strange thing, is that I'm actually loosing polygons which were valid as shown by the following figures where I added the invalid geometries:

• with a simple dark blue outline (no fill color) for shapes with a raster_value of 0 (black on the original raster) and
• with a red outline for those with a raster_value of 1 (white of the original raster).

In the background, there is still the black and white original raster with some transparency.

The polygons which are actually lost when dissolving the GeoDataFrame are the small ones with a light gray shading, not the darkest ones. One can easily see that they are not in the set of invalid geometries!

And as shown by how the original Shapely's geometries are fixed:

Verified facts:

The problem appears when using either the .dissolve() method alone on the whole GeoDataFrame, or the .buffer(0) method applied on all invalid features (rows) alone.

If I do nothing, every shape is present in the GeoPackage, but if I use one of these methods, the exact same small shapes actually disappear.

It seems that my largest feature (which is quite big and complex; encapsulating many small black regions) is "closing"* some of those encapsulated black regions... I'm not 100% sure why. The only thing that seems logical, is that it is closing one of those "corner-paired" encapsulated black shapes, keeping (i.e. not "closing") the other of the pair "alive":

*Note: what I mean by closing is a kind of morphological closing operation.

Answer 1

An upgrade from the Python package rasterio version 1.2.6 to its newest version, 1.2.8, has fixed the issue. All my geometries are now valid. This is maybe because this version has been build upon GDAL 3.3.2 which now also fixes the same issue when using gdal_polygonize.py, and which seems to be the default on Debian 11, release on 14 August 2021.

As it has been for at least >1 year I was having troubles with vectorization of complicated binary raster resulting from a viewshed analysis, I'd like to be particularly thankful to the GDAL and rasterio developers (+all the others of course) for their work! Really!