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

@Leo when I run it through QGIS GUI it says self-intersection — Perfalcon, Jul 18 '19 at 08:58

score 0 · Answer 1 · answered Jul 16 '19 at 12:46

0

Maybe you could try the polygonization with rasterio?

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

answered Jul 16 '19 at 12:46

Leo

952
6
23

Thanks for your comment. It is a good way to polygonize, but the Fix Geometries still do not work. Also, the link you posted is out-of-date, I think. Here is one that works with Python 3.7 – Perfalcon Jul 16 '19 at 14:00

score 0 · Answer 2 · answered Jul 18 '19 at 09:58

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.

score 0 · Answer 3 · edited Sep 17 '21 at 22:37

0

A little trick solved the problem to me: Openning 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.

edited Sep 17 '21 at 22:37

Kadir Şahbaz

76,800
56
247
389

answered Jul 27 '21 at 18:43

Rafael Matias

1

gdal version https://gis.stackexchange.com/a/263721/2856 – user2856 Sep 18 '21 at 03:29

score 0 · Answer 4 · answered Mar 01 '23 at 13:32

Му 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)

Invalid geometry after gdal.polygonize

4 Answers4