1

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:

Original raster file

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):

vectorized result

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:

Ooops...export went wrong...

It seems correct, but...

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

differences

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):

invalid geometries

And a closeup:

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

invalid 0: dark blue

invalid 1: red

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:

Shapely's geometry fix method

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":

corner-paired polygons; one is kept, the other is lost

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

3
  • 1
    Can you ensure that all polygons are valid before calling dissolve? You can check with gdf.is_valid.all() and potentially fix with gdf.geometry = gdf.buffer(0). Jul 9, 2021 at 10:47
  • Updated main post. Thanks a lot for the hints! Tried gdf.geometry = gdf.buffer(0) with success (on my first case). Is it due to luck? Or is it entirely explainable? Especially regarding this "loopy" pattern I can notice on every single shape that is actually invalid...? But then, after dissolve, the problem persists: I miss those geometries that were invalid but were actually fixed... That's strange. Jul 9, 2021 at 17:17
  • Can you report it to github.com/geopandas/geopandas/issues? buffer(0) fixes the order of coordinates (difference between CW and CCW), so that bit is explainable. The issue with dissolve less so. Jul 12, 2021 at 8:20

1 Answer 1

0

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!

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.