Invalid geometry after gdal.polygonize

Question

After converting a raster file to a vector file using gdal.polygonize, I cannot carry out further commands (e.g. Clip) because the vector is produced with invalid geometries (as determined by the QGIS GUI). Within GIS itself, I can correct this by using Fix Geometries, but cannot using Python code outside of the GUI. I can't work out why this happens.

Is there a way to check/make valid geometries in Python/GDAL framework? I have not been able to find anything that works.

src_ds = gdal.Open('input.tif')
if src_ds is None:
    print ('Unable to open %s' % src_filename)
    sys.exit(1)
try:
    srcband = src_ds.GetRasterBand(1)
except RuntimeError:
    print ('Band 1 not found')
    sys.exit(1)

dst_layername = "output/path/filename"
drv = ogr.GetDriverByName("ESRI Shapefile")
if os.path.exists(dst_layername + ".shp"):
    drv.DeleteDataSource(dst_layername + ".shp")

# create spatial reference for the new vector based on the raster
srs = osr.SpatialReference()
srs.ImportFromWkt(src_ds.GetProjection())
dst_layer = dst_ds.CreateLayer(dst_layername, srs = srs)

# create polygon layer
dst_ds = drv.CreateDataSource(dst_layername + ".shp")
gdal.Polygonize(srcband, srcband, dst_layer, -1, [], callback=None)

### fix geometries
processing.run("native:fixgeometries",
        {'INPUT':"output/path/filename.shp",
        'OUTPUT':"output_fixed.shp"})
# > produces a file, but shows nothing

Answer 1

Maybe you could try the polygonization with rasterio?

https://gist.github.com/sgillies/8465883

Answer 2

you will have to look for the geometry that is not valid and then delete it or change it so that it becomes valid. try it with a code like this:

faultyshapefile = gp.read_file('your shapefile.shp')
    nonecontroller = 0
    while nonecontroller == 0:
        nonecontroller = 1
        for i in range(len(faultyshapefile)):
            shape=faultyshapefile.iloc[i]['geometry']
            if shape.is_valid != True:
                faultyshapefile = faultyshapefile.drop([i], axis=0)
                nonecontroller = 0
                break

If you do not want to loose the shape, instead of "drop" you can do any other operation on the geometry that fixes it. for example buffering with a very small radius.

Answer 3

Му solution is also make buffer

gdal.Polygonize(rasterband, rasterband, mem_layer, 0)
for mem_feat in mem_layer:
    if mem_geom.GetGeometryName() == 'MULTIPOLYGON':
        for i in range(0, mem_geom.GetGeometryCount()):
            g = mem_geom.GetGeometryRef(i)
            if not g.IsValid(): 
                g = g.Buffer(0)
    elif mem_geom.GetGeometryName() == 'POLYGON':
        if not mem_geom.IsValid(): 
            mem_geom = mem_geom.Buffer(0)

Answer 4

A little trick solved the problem to me: Opening the gdal exported shapefile with Geopandas and fix geometries.

shape1['geometry'] = shape1.buffer(0), where shape1 is a GeoDataFrame.

Although it works with Geopandas I couldn't make it work directly with gdal.Polygonize.