0

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

bumped to the homepage by Community yesterday

This question has answers that may be good or bad; the system has marked it active so that they can be reviewed.

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

Maybe you could try the polygonization with rasterio?

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

  • 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 at 14:00
0

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.

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.