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)
    srcband = src_ds.GetRasterBand(1)
except RuntimeError:
    print ('Band 1 not found')

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()
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
# > produces a file, but shows nothing

  • have you identified what causes the invalidity? – Leo Jul 16 '19 at 15:14
  • @Leo when I run it through QGIS GUI it says self-intersection – Perfalcon Jul 18 '19 at 8:58

Maybe you could try the polygonization with rasterio?


| improve this answer | |
  • 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

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)):
            if shape.is_valid != True:
                faultyshapefile = faultyshapefile.drop([i], axis=0)
                nonecontroller = 0

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.

| improve this answer | |

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.