Algoritmo para recortar un gran volumen de puntos en geometrías irregulares?

Question

Actualmente estoy utilizando un método casero para recortar un gran número de puntos (~500k, podría ser más o menos dependiendo del día) por un círculo con un radio fijo. La salida es una lista que contiene sólo los puntos dentro del círculo

Este es mi método actual usando Python y 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])

Este método funciona, pero es muy lento.

Me gustaría evitar añadir otra dependencia, ya que mi lista de dependencias es cada vez más larga.

¿Existe algún tipo de algoritmo de ordenación que pueda implementar para acelerar este proceso y que utilice ogr?

He estado buscando alguna forma de implementar los quad-trees, pero no he podido encontrar ningún buen ejemplo de implementación en este contexto.

Los puntos se almacenan inicialmente en un array 2D de numpy. Sería útil desde el punto de vista de la eficiencia no tener que convertirlos en geometrías.

Preferiría recortar a los polígonos irregulares, no sólo a los círculos.

Los puntos están actualmente en UTM.

Answer 1

Aquí hay un ejemplo de trabajo con los enlaces de GDAL Python usando ogr2ogr ( gdal.VectorTranslate ) como función de biblioteca y el Índice espacial Spatialite :

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)

Probado con un conjunto de datos de ~1,5 millones de puntos utilizando un polígono irregular como cortador, ¡se corta después de ~35s (sin índice espacial) y ~4s (con índice espacial) en mi PC!

Answer 2

Si las formas están en un sistema de coordenadas lineal, es bastante fácil comprobar si un punto está dentro de un polígono. La idea es dibujar una línea infinita (o suficientemente larga) desde el punto hacia la derecha, y ver cuántas veces esta línea interseca segmentos en el polígono. Si esto sucede impar veces, entonces el punto está dentro del polígono.

Véase, por ejemplo ¿Cómo comprobar si un punto dado está dentro o fuera de un polígono?

No voy a copiar todo el texto en esta respuesta, ya que es un método genérico y debe haber múltiples exámenes del mismo en la red.

Tenga en cuenta que, como primera optimización, se podría, para cada punto a comprobar, excluir todos los segmentos de línea del polígono que empiezan y terminan por encima o por debajo de ese punto. Si se trata de un polígono cerrado, no importa cuántos segmentos se hayan dejado fuera: el polígono siempre "vuelve". Sin embargo, esto significa que el medio millón de puntos tiene que ser comprobado con una cantidad menor de segmentos de polígono.

La distancia sí importa. Si trabaja con áreas locales en lugar de países, y permite cierto grado de aproximación, y tiene vértices de polígonos suficientemente densos (y distribuidos uniformemente) (como una docena o unas pocas docenas), podría intentar trabajar directamente con las coordenadas UTM. Sin embargo, tenga en cuenta que habrá errores, ya que los segmentos de las líneas de los polígonos no son rectos, en la naturaleza. Supongo que se trata de una cuestión de rendimiento frente a la calidad. Para conseguir precisión, habría que convertir todas las coordenadas (tanto el medio millón de puntos como los vértices del polígono) en coordenadas x,y en un sistema de coordenadas lineal.

Como codificador de juegos podría optimizar. No mencionaste si los polígonos son estáticos, por ejemplo, si cambian sólo una vez al día. Tampoco mencionaste a qué ritmo caen los 1/2 millones de puntos. Si hay una demanda de rendimiento, entonces es esencial precalcular cuando hay tiempo, incluso a costa de datos más grandes. Si los puntos caen a lo largo del día, ya se podría resolver en ese momento un caché de trandform de coordenadas para cada punto individual.

Una optimización adicional sería una forma de concluir, a nivel aproximado, qué puntos están con toda seguridad dentro de un polígono, y sólo dejar los puntos inciertos para una investigación posterior. Podría simplemente dibujar una forma en una superficie (textura) de 512x512 con un lápiz rojo grueso (por ejemplo, 32 píxeles) y un relleno verde. Entonces sería fácil buscar cualquier punto. La coordenada del punto se ajustaría a alguna cuadrícula, pero esto no importa ya que la forma pintada tiene un borde gordo. Entonces, si un píxel muestreado es verde, lo más seguro es que esté dentro del polígono. Si es rojo, está en algún lugar cerca de la línea de borde y necesita más trabajo. etc. Naturalmente, uno podría tomar los píxeles en datos numéricos y hacer las búsquedas en una matriz. La pereza conveniente aquí es dibujar correctamente un polígono cerrado con un lápiz grueso y obtener todos los colores de los píxeles correctamente; esto es a menudo una funcionalidad nativa en los lenguajes de ordenador.