1

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

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

4 Answers 4

0

Maybe you could try the polygonization with rasterio?

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

1
  • 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
    Commented Jul 16, 2019 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.

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.

1
0

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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