Algorithm for clipping large volume of points to irregular geometries?

Question

I am currently using a home-brewed-hacked-together method for clipping a large number of points (~500k, could be more or less depending on the day) by a circle with a fixed radius. The output is a list containing only the points within the circle

This is my current method using Python & OGR bindings:

plot_points = []
circle = #< Some sort of circle geom>
for point in point_array:
     thing = ogr.Geometry(ogr.wkbPoint)
     thing.AddPoint(point[1], point[2])
     if thing.Within(circle):
         plot_points.append(point[0])

This method works, but is very slow.

I would like to avoid adding another dependency, as my list of dependencies is growing ever-longer.

Is there some sort of sorting algorithm I could implement that would speed up this process that uses ogr?

I have been looking into implementing quad-trees somehow, but have been unable to find any good examples of implementation in this context.

The points are initially stored in a 2D numpy array. It would be useful from an efficiency standpoint to not have to convert them into geometries.

I would prefer to clip to irregular polygons, not just circles.

Points are currently in UTM.

Answer 1

Here's a working example with the GDAL Python bindings using ogr2ogr (gdal.VectorTranslate) as library function and the Spatialite spatial index:

from osgeo import gdal

srcDS1 = gdal.OpenEx('/path/to/your/points.shp')
gdal.VectorTranslate('clip.sqlite', srcDS1, format = 'SQLite', layerName = 'points', datasetCreationOptions = ['SPATIALITE=YES'], callback=gdal.TermProgress)
srcDS2 = gdal.OpenEx('/path/to/your/clipper/polygon.shp')
gdal.VectorTranslate('clip.sqlite', srcDS2, format = 'SQLite', accessMode = 'update', layerName = 'clipper_polygon', datasetCreationOptions = ['SPATIALITE=YES'], callback=gdal.TermProgress)
gdal.VectorTranslate('clip.sqlite', gdal.OpenEx('clip.sqlite'), format = 'SQLite', accessMode = 'update', layerName = 'clipped_points', datasetCreationOptions = ['SPATIALITE=YES'], SQLStatement='SELECT g1.* FROM points g1, clipper_polygon g2 WHERE ST_Intersects(g1.geometry, g2.geometry) AND g1.rowid IN (SELECT rowid FROM SpatialIndex WHERE f_table_name = "points" AND search_frame = g2.geometry)', SQLDialect = 'SQLITE', callback=gdal.TermProgress)

Tested with a ~1.5M points dataset using an irregular polygon as clipper, it clips after ~35s (without spatial index) and ~4s (with spatial index) on my pc!

Answer 2

If the shapes are on a linear coordinate system, it's rather easy to check if a point is inside a polygon. The idea is to draw an infinite (or sufficiently long) line from the point towards right, and see how many times this line intersects segments in the polygon. If this happens odd times, then the point is inside the polygon.

Please see for example How to check if a given point lies inside or outside a polygon?

I will not copy the entire text into this answer, as it is a generic method and there should be multiple examles of it in the net.

Note that as a first optimisation, one could for each point to check exclude all polygon line segments that both start and end either above or below that point. If it is a closed polygon, it doesn't matter how many segments were left out - the polygon always "comes back". This however means that the 1/2 million points need to be checked against a smaller amount of polygon segments.

Distance does matter. If you work with local areas rather than countries, and allow some degree of approximation, and have sufficiently dense (and evenly distributed) polygon vertices (like a dozen or a few dozen), you could try to work directly with the UTM coordinates. However be aware that there will be errors, as the polygon line segments are not straight, out in nature. I take it this is a performance vs. quality issue. To get accuracy, one would have to convert all coordinates (both the 1/2 million points and the polygon vertices) into x,y-coordinates on a linear coordinate system.

As a game coder i might optimise. You didn't mention if the polygon(s?) is static, eg. changes only once a day. Neither did you mention at what rate the 1/2 million points drop in. If there is a performance demand, then pre-calculating when there is time is essential, even at the cost of larger data. If the points drop in along the day, one could already at that stage solve a coordinate trandform cache for each individual point.

A further optimisation would be a way to conclude, at a rough level, which points are most certainly inside a polygon, and only leave uncertain points for further investigation. I might simply draw a shape onto a 512x512 surface (texture) with a thick (eg. 32 pixel) red pen and a green fill. Then it would be easy to lookup any point. The point coordinate would snap onto some grid, but this doesn't matter as the painted shape has a fat borderline. Then, if a sampled pixel is green, it's most certainly inside the polygon. If it's red, it is somewhere near the bordeline and needs further work. Etc. Naturally one could grab the pixels into numeric data and do the lookups in an array. The convenient lazyness here is correctly drawing a closed polygon with a thick pen and getting all pixel colors right; this is often native functionality in computer languages.