El casco cóncavo

Creando un borde de clúster utilizando un enfoque de vecinos K más cercanos

Hace un par de meses, escribí aquí en Medium un artículo sobre el mapeo de los puntos calientes de accidentes de tráfico del Reino Unido. Me preocupaba sobre todo ilustrar el uso del algoritmo de agrupación en clústeres DBSCAN en datos geográficos. En el artículo, utilicé información geográfica publicada por el gobierno del Reino Unido sobre accidentes de tráfico reportados. Mi propósito era ejecutar un proceso de agrupación basado en la densidad para encontrar las áreas donde los accidentes de tráfico se reportan con mayor frecuencia. El resultado final fue la creación de un conjunto de geo-cercas que representan estos puntos calientes de accidentes.

Al recopilar todos los puntos de un clúster determinado, puede hacerse una idea de cómo se vería el clúster en un mapa, pero carecerá de una información importante: la forma exterior del clúster. En este caso, estamos hablando de un polígono cerrado que se puede representar en un mapa como una geo-valla. Se puede suponer que cualquier punto dentro de la geo-valla pertenece al clúster, lo que hace que esta forma sea una información interesante: se puede usar como función discriminadora. Se puede suponer que todos los puntos de muestreo nuevos que caen dentro del polígono pertenecen al clúster correspondiente. Como insinué en el artículo, podría usar dichos polígonos para afirmar su riesgo de conducción, usándolos para clasificar su propia ubicación GPS muestreada.

La pregunta ahora es cómo crear un polígono significativo a partir de una nube de puntos que conforman un clúster específico. Mi enfoque en el primer artículo era algo ingenuo y reflejaba una solución que ya había utilizado en código de producción. Esta solución implicaba colocar un círculo centrado en cada uno de los puntos del clúster y luego fusionar todos los círculos para formar un polígono en forma de nube. El resultado no es muy agradable, ni realista. Además, al usar círculos como forma de base para construir el polígono final, estos tendrán más puntos que una forma más aerodinámica, lo que aumentará los costos de almacenamiento y hará que los algoritmos de detección de inclusión sean más lentos de ejecutar.

Por otro lado, este enfoque tiene la ventaja de ser computacionalmente simple (al menos desde la perspectiva del desarrollador), ya que utiliza la función cascaded_union de Shapely para fusionar todos los círculos. Otra ventaja es que la forma del polígono se define implícitamente utilizando todos los puntos del clúster.

Para enfoques más sofisticados, necesitamos identificar de alguna manera los puntos fronterizos del clúster, los que parecen definir la forma de la nube de puntos. Curiosamente, con algunas implementaciones de DBSCAN, en realidad puede recuperar los puntos fronterizos como un subproducto del proceso de agrupación en clústeres. Desafortunadamente, esta información (aparentemente) no está disponible en la implementación de SciKit Learn, por lo que tenemos que conformarnos.

El primer enfoque que se me ocurrió fue calcular el casco convexo del conjunto de puntos. Este es un algoritmo bien entendido, pero sufre del problema de no manejar formas cóncavas, como esta:

Esta forma no captura correctamente la esencia de los puntos subyacentes. Si se usara como discriminador, algunos puntos se clasificarían incorrectamente como dentro del grupo cuando no lo están. Necesitamos otro enfoque.

La Alternativa de casco Cóncavo

Afortunadamente, hay alternativas a este estado de cosas: podemos calcular un casco cóncavo. Así es como se ve el casco cóncavo cuando se aplica al mismo conjunto de puntos que en la imagen anterior:

O tal vez este:

Como puede ver, y contrariamente al casco convexo, no hay una definición única de lo que es el casco cóncavo de un conjunto de puntos. Con el algoritmo que presento aquí, la elección de cuán cóncavo desea que sean sus cascos se realiza a través de un solo parámetro: k — el número de vecinos más cercanos considerados durante el cálculo del casco. Veamos cómo funciona esto.

El Algoritmo

El algoritmo que presento aquí fue descrito hace más de una década por Adriano Moreira y Maribel Yasmina Santos de la Universidad de Minho, Portugal . Del resumen:

Este artículo describe un algoritmo para calcular la envolvente de un conjunto de puntos en un plano, que genera cascos convexos en cascos no convexos que representan el área ocupada por los puntos dados. El algoritmo propuesto se basa en un enfoque de k-vecinos más cercanos, donde el valor de k, el único parámetro del algoritmo, se utiliza para controlar la “suavidad” de la solución final.

Debido a que aplicaré este algoritmo a la información geográfica, se tuvieron que hacer algunos cambios, especialmente al calcular ángulos y distancias . Pero estos no alteran de ninguna manera la esencia del algoritmo, que se puede describir ampliamente en los siguientes pasos:

1. Encuentre el punto con la coordenada y (latitud) más baja y conviértalo en el actual.
2. Encuentre los puntos k más cercanos al punto actual.
3. En los puntos k más cercanos, seleccione el que corresponda al giro a la derecha más grande desde el ángulo anterior. Aquí utilizaremos el concepto de rodamiento y comenzaremos con un ángulo de 270 grados (hacia el oeste).
4. Compruebe si al agregar el nuevo punto a la cadena de línea en crecimiento, no se cruza a sí misma. Si lo hace, seleccione otro punto del k-más cercano o reinicie con un valor mayor de k.
5. Haga que el nuevo punto sea el punto actual y elimínelo de la lista.
6. Después de k iteraciones, agregue el primer punto a la lista.
7. Bucle al número 2.

El algoritmo parece ser bastante simple, pero hay una serie de detalles que deben ser atendidos, especialmente porque estamos tratando con coordenadas geográficas. Las distancias y los ángulos se miden de una manera diferente.

El Código

Aquí estoy publicando es una versión adaptada del código del artículo anterior. Todavía encontrará el mismo código de agrupación en clústeres y el mismo generador de clústeres en forma de nube. La versión actualizada ahora contiene un paquete llamado geomath.hulls donde puede encontrar la clase ConcaveHull. Para crear sus cascos cóncavos, haga lo siguiente:

En el código anterior, points es una matriz de dimensiones (N, 2), donde las filas contienen los puntos observados y las columnas contienen las coordenadas geográficas (longitud, latitud). El array resultante tiene exactamente la misma estructura, pero contiene solo los puntos que pertenecen a la forma poligonal del clúster. Una especie de filtro.

Debido a que manejaremos matrices, es natural llevar a NumPy a la refriega. Todos los cálculos se vectorizaron debidamente, siempre que fue posible, y se hicieron esfuerzos para mejorar el rendimiento al agregar y eliminar elementos de matrices (spoiler: no se mueven en absoluto). Una de las mejoras que faltan es la paralelización de código. Pero eso puede esperar.

Organicé el código en torno al algoritmo expuesto en el documento, aunque se realizaron algunas optimizaciones durante la traducción. El algoritmo está construido alrededor de una serie de subrutinas que están claramente identificadas por el papel, así que quitémoslas del camino ahora mismo. Para su comodidad de lectura, usaré los mismos nombres que se usan en el periódico.

CleanList – La limpieza de la lista de puntos se realiza en el constructor de clases:

Como puede ver, la lista de puntos se implementa como una matriz NumPy por razones de rendimiento. La limpieza de la lista se realiza en la línea 10, donde solo se guardan los puntos únicos. El conjunto de datos está organizado con observaciones en filas y coordenadas geográficas en las dos columnas. Tenga en cuenta que también estoy creando una matriz booleana en la línea 13 que se utilizará para indexar en la matriz del conjunto de datos principal, facilitando la tarea de eliminar elementos y, de vez en cuando, agregarlos de nuevo. He visto esta técnica llamada “máscara” en la documentación de NumPy, y es muy poderosa. En cuanto a los números primos, los discutiré más tarde.

FindMinYPoint-Esto requiere una pequeña función:

Esta función se llama con la matriz de conjuntos de datos como argumento y devuelve el índice del punto con la latitud más baja. Tenga en cuenta que las filas están codificadas con longitud en la primera columna y latitud en la segunda.

RemovePoint
AddPoint – Estos son no cerebros, debido al uso de la matriz indices. Esta matriz se utiliza para almacenar los índices activos en la matriz del conjunto de datos principal, por lo que eliminar elementos del conjunto de datos es muy sencillo.

Aunque el algoritmo descrito en el documento requiere agregar un punto a la matriz que compone el casco, esto en realidad se implementa como:

Más tarde, la variable test_hull se asignará de nuevo a hull cuando se considere que la cadena de línea no se cruza. Pero me estoy adelantando al juego. Eliminar un punto de la matriz de conjuntos de datos es tan sencillo como:

self.indices = False

Volver a agregarlo es solo cuestión de voltear el valor de la matriz en el mismo índice a verdadero. Pero toda esta comodidad viene con el precio de tener que controlar los índices. Más sobre esto más adelante.

NearestPoints-Las cosas comienzan a ponerse interesantes aquí porque no estamos tratando con coordenadas planas, así que fuera con Pitágoras y dentro con Haversine:

Tenga en cuenta que el segundo y tercer parámetros son matrices en el formato de conjunto de datos: longitud en la primera columna y latitud en la segunda. Como puede ver, esta función devuelve una matriz de distancias en metros entre el punto del segundo argumento y los puntos del tercer argumento. Una vez que tengamos estos, podemos obtener los vecinos k más cercanos de la manera fácil. Pero hay una función especializada para eso y merece algunas explicaciones:

La función comienza creando una matriz con los índices base. Estos son los índices de los puntos que no se han eliminado de la matriz de conjunto de datos. Por ejemplo, si en un clúster de diez puntos comenzamos eliminando el primer punto, la matriz de índices base sería . A continuación, calculamos las distancias y ordenamos los índices de matriz resultantes. Las primeras k se extraen y luego se utilizan como máscara para recuperar los índices base. Es un poco retorcido, pero funciona. Como puede ver, la función no devuelve una matriz de coordenadas, sino una matriz de índices en la matriz del conjunto de datos.

SortByAngle-Hay más problemas aquí porque no estamos calculando ángulos simples, estamos calculando rodamientos. Estos se miden como cero grados hacia el Norte, con ángulos que aumentan en el sentido de las agujas del reloj. Aquí está el núcleo del código que calcula los rodamientos:

La función devuelve una matriz de rodamientos medidos desde el punto cuyo índice está en el primer argumento hasta los puntos cuyos índices están en el tercer argumento. La clasificación es sencilla:

En este punto, la matriz de candidatos contiene los índices de los puntos k más cercanos ordenados por orden descendente de rodamiento.

IntersectQ: En lugar de girar mis propias funciones de intersección de líneas, recurrí a Shapely para obtener ayuda. De hecho, mientras construimos el polígono, esencialmente estamos manejando una cadena de líneas, anexando segmentos que no se cruzan con los anteriores. Las pruebas para esto son simples: recogemos la matriz de casco en construcción, la convertimos en un objeto de cadena de línea bien formada y probamos si es simple (no auto-intersectable) o no.

En pocas palabras, una cadena de líneas bien formada se vuelve compleja si se cruza automáticamente, por lo que el predicado is_simple se convierte en falso. Sencillo.

PointInPolygon-Este resultó ser el más complicado de implementar. Permítanme explicarlo mirando el código que realiza la validación final del polígono de casco (comprueba si todos los puntos del clúster están incluidos en el polígono):

Las funciones de Shapely para probar la intersección y la inclusión deberían haber sido suficientes para verificar si el polígono de casco final se superpone a todos los puntos del clúster, pero no lo fueron. ¿Por qué? Shapely es independiente de las coordenadas, por lo que manejará las coordenadas geográficas expresadas en latitudes y longitudes exactamente de la misma manera que las coordenadas en un plano cartesiano. Pero el mundo se comporta de manera diferente cuando vives en una esfera, y los ángulos (o rodamientos) no son constantes a lo largo de una geodésica. El ejemplo de referencia de una línea geodésica que conecta Bagdad con Osaka ilustra perfectamente esto. Sucede que en algunas circunstancias, el algoritmo puede incluir un punto basado en el criterio de rodamiento, pero más tarde, utilizando los algoritmos planos de Shapely, se considera ligeramente fuera del polígono. Eso es lo que la corrección de distancia pequeña está haciendo allí.

Me llevó un tiempo descifrar esto. Mi ayuda de depuración fue QGIS, una gran pieza de software libre. En cada paso de los cálculos sospechosos, enviaba los datos en formato WKT a un archivo CSV para leerlos como una capa. Un verdadero salvavidas!

Finalmente, si el polígono no cubre todos los puntos del clúster, la única opción es aumentar k e intentarlo de nuevo. Aquí agregué un poco de mi propia intuición.

Prime k

El artículo sugiere que el valor de k se incremente en uno y que el algoritmo se ejecute de nuevo desde cero. Mis primeras pruebas con esta opción no fueron muy satisfactorias: los tiempos de ejecución serían lentos en clústeres problemáticos. Esto se debió al lento aumento de k, así que decidí usar otro horario de aumento: una tabla de números primos. El algoritmo ya comienza con k=3, por lo que fue una extensión fácil para hacerlo evolucionar en una lista de números primos. Esto es lo que ves que sucede en la llamada recursiva:

Me gustan los números primos, ya sabes

Blow Up

Los polígonos cóncavos del casco generados por este algoritmo todavía necesitan un procesamiento adicional, porque solo discriminarán puntos dentro del casco, pero no cerca de él. La solución es agregar un poco de relleno a estos racimos delgados. Aquí estoy usando exactamente la misma técnica que se usó antes, y esto es lo que parece:

Aquí, he utilizado bien formada del bufferfunción para hacer el truco.

La función acepta un polígono bien formado y devuelve una versión inflada de sí misma. El segundo parámetro es el radio en metros del relleno añadido.

Ejecutar el código

Comience por extraer el código del repositorio GitHub a su máquina local. El archivo que desea ejecutar es ShowHotSpots.py en el directorio principal. Tras la primera ejecución, el código se leerá en los datos de accidentes de tráfico del Reino Unido de 2013 a 2016 y se agrupará. Los resultados se almacenan en caché como un archivo CSV para las ejecuciones posteriores.

Se le presentarán dos mapas: el primero se genera utilizando los clústeres en forma de nube, mientras que el segundo utiliza el algoritmo de clústeres cóncavos que se discute aquí. Mientras se ejecuta el código de generación de polígonos, es posible que vea que se reportan algunos errores. Para ayudar a entender por qué el algoritmo no puede crear un casco cóncavo, el código escribe los clústeres en archivos CSV en el directorio data/out/failed/. Como de costumbre, puede usar QGIS para importar estos archivos como capas.

Esencialmente, este algoritmo falla cuando no encuentra suficientes puntos para” rodear ” la forma sin intersectarse automáticamente. Esto significa que debe estar listo para desechar estos grupos o para aplicarles un tratamiento diferente (casco convexo o burbujas fusionadas).

Concavidad

Es una envoltura. En este artículo presenté un método de posprocesamiento de clústeres geográficos generados por DBSCAN en formas cóncavas. Este método puede proporcionar un polígono exterior que se ajuste mejor a los clústeres en comparación con otras alternativas.

¡Gracias por leer y disfrutar retocando el código!

Kryszkiewicz M., Lasek P. (2010) TI-DBSCAN: Agrupación con DBSCAN por medio de la Desigualdad Triangular. In: Szczuka M., Kryszkiewicz M., Ramanna S., Jensen R., Hu Q. (eds) Rough Sets and Current Trends in Computing. RSCTC 2010. Lecture Notes in Computer Science, vol 6086. Springer, Berlin, Heidelberg

Scikit-learn: Aprendizaje automático en Python, Pedregosa et al., JMLR 12, pp. 2825-2830, 2011

Moreira, A. y Santos, M. Y., 2007, Casco Cóncavo: Un enfoque K-vecinos más cercanos para el cálculo de la región ocupada por un conjunto de puntos

Calcular la distancia, el rumbo y más entre los puntos de Latitud/Longitud

repositorio GitHub