[GIS] How to avoid self-intersections in geometry created by GDAL Polygonize in Python

gdalpolygonizepythonshapefile

I am creating geometry from a raster using GDAL Polygonize in Python. My issue is that doing this creates a lot of self intersection errors. This seems to occur due to the nature of the pixel grid.

I have noticed that this occurs mostly at edges, but can sometimes occur inside the geometry as well, as shown below.

This is part of an automated pipeline so I don't want to manually fix self-intersections in QGIS or the like; the solution needs to be able to be implemented in Python code.

Is there a way to either stop the self intersections being created in the first place, or a way to clean the generated geometry in Python?

I have tried to use shapely to clean the self intersections but the result is very poor.

Best Answer

Yes, this is because pixels always share corners with 8 others. There's no silver bullet for turning the invalid polygons into valid polygons or multipolygons. That said, in your examples above it looks like the self intersections would be solved by eroding (negative buffer, in other words) the invalid polygons by a small fraction, turning them into multipolygons.

Related Solutions

shapefile – Exploring Self-Intersecting Polygons in Shapefile Data Structures

Here are relevant excerpts from pages 8 and 9 of the Shapefile Technical Description (emphasis mine):

A polygon consists of one or more rings. A ring is a connected sequence of four or more points that form a closed, non-self-intersecting loop.

And:

Polygons stored in a shapefile must be clean. A clean polygon is one that

Has no self-intersections. This means that a segment belonging to one ring may not intersect a segment belonging to another ring. The rings of a polygon can touch each other at vertices but not along segments. Colinear segments are considered intersecting.

So, as far as the shapefile specification is concerned, the answer is no, polygons may not intersect themselves.

This is a conceptual rule, not a limitation of the file format. It is quite possible to construct shapefiles containing self-intersecting polygons. However, since self-intersecting geometry is invalid according to the specification, applications may differ in how they handle such cases: some may attempt to clean up the features, while others may [justifiably, if less helpfully] consider the shapefile malformed. Most problematic may be cases where the software assumes a shapefile conforms to the specification (ie, contains no self-intersections) but does not validate it before performing geometric operations, which may yield wrong results.

It's unclear exactly what you want to do, but if you're generating shapefiles, I'd advise against outputting polygons with self-intersections, as it may cause difficulty for your data users. Likewise, if you're writing some sort of algorithm to process shapefiles, you do not need to support self-intersecting polygons to comply with the specification - but if possible, I'd advise you check the input files are not malformed before proceeding.

GDAL and Python – Getting Coordinates for Specific Raster Values

In GDAL you can import the raster as a numpy array.

from osgeo import gdal
import numpy as np

r = gdal.Open("path/to/raster")
band = r.GetRasterBand(1) #bands start at one
a = band.ReadAsArray().astype(np.float)

Then using numpy it is a simple matter to get the indexes of an array matching a boolan expression:

(y_index, x_index) = np.nonzero(a > threshold)
#To demonstate this compare a.shape to band.XSize and band.YSize

From the raster geotransform we can get information such as the upper left x and y coordinates and the cell sizes.

(upper_left_x, x_size, x_rotation, upper_left_y, y_rotation, y_size) = r.GetGeoTransform()

The upper left cell corresponds to a[0, 0]. The Y size will always be negative, so using the x and y indices you can calculate the coordinates of each cell based on the indexes.

x_coords = x_index * x_size + upper_left_x + (x_size / 2) #add half the cell size
y_coords = y_index * y_size + upper_left_y + (y_size / 2) #to centre the point

From here it's a simple enough matter to create a shapefile using OGR. For some sample code see this question for how to generate a new dataset with point information.

Best Answer

Related Solutions

shapefile – Exploring Self-Intersecting Polygons in Shapefile Data Structures

GDAL and Python – Getting Coordinates for Specific Raster Values

Related Question