Geoinformática - Práctica 7¶
Patrones de Puntos¶
Los Puntos son entidades espaciales que pueden estudiarse de dos formas fundamentamentales diferentes. Por un lado, los puntos pueden ser vistos como objetos fijos en el espacio, esto es, su ubicación se encuentra determinada; desde esta perspectiva, el análisis de los puntos es muy similar al de otros tipos de datos espaciales, como son polígonos y líneas. Por otra parte, los puntos pueden ser vistos como la ocurrrencia de un evento que en teoría pudiese ocurrir en donde sea, pero únicamente se manifiesta en ciertos lugares; ésta es la aproximación que se utilizará para el resto de la práctica.
Dado que los puntos son vistos como eventos que pueden ocurrir en múltiples ubicaciones, pero únicamente suceden en unas cuantas, se tiene que una colección de dichos eventos se llama un Patrón de Puntos; en este caso, la ubicación de los puntos es un aspecto clave del análisis. Un buen ejemplo de un Patrón de Puntos son los crímenes en una ciudad; en teoría, podrían ocurrir en cualquier lugar dentro de la misma, pero normalmente se observa que únicamente suceden en algunos cuantos. Los Patrones de Puntos se consideran Marcados si más atributos además de la ubicación son proporcionados, o No Marcados si únicamente se conocen las coordenadas del evento. Continuando con el ejemplo de los crímenes, un Patrón No Marcado resultaría de únicamente utilizar el lugar en donde ocurren éstos para el análisis, mientras que se tendría un Patrón Marcado si otros atributos, como el tipo de crimen, la extensión del daño, entre otros, también se tuviesen a la mano.
Como tal, el Análisis de Patrones de Puntos se enfoca en la descripción, caracterización estadística y modelación de Patrones de Puntos, centrándose específicamente en metodologías que expliquen el comportamiento de los datos, a través de preguntas como ¿Cuál es la naturaleza de la distribución de los puntos?, ¿Hay alguna característica en la forma en que las ubicaciones se acomodan en el espacio que pueda discernirse estadísticamente?, ¿Por qué esos eventos ocurren en esos lugares y no en otros?
Esta práctica pretende ser un primer acercamiento a trabajar con Patrones de Puntos en Python; como tal, se estudiará cómo importar, procesar y transformar datos puntuales, así como múltiples formas en que pueden visualizarse los Patrones de Puntos.
# Importación de Librerías
%matplotlib inline
import numpy as np
import pandas as pd
import geopandas as gpd
import pysal as ps
import seaborn as sns
import matplotlib.pyplot as plt
from sklearn.cluster import dbscan # Instalada bajo el nombre 'scikit-learn'
from ipywidgets import interact, fixed
Datos¶
Para iniciar con el Análisis de Patrones Puntuales, se trabajará con todos los robos reportados a la Procuraduría General de Justicia de la CDMX durante el 2018, así como las AGEB's de la misma, con las que ya se ha trabajado en prácticas anteriores. Como es usual, es necesario definir la ubicación de estos archivos:
f = 'data/'
Dado que los datos se encuentran almacenados como ShapeFile, su método de importación es el ya practicado con GeoPandas:
%%time
# Importar los Datos
robos = gpd.read_file(f + 'robos_cdmx_2018.shp')
# Analizar rápidamente los datos
robos.info()
Puede notarse que se ha añadido el comando %%time en la parte superior de la celda; esto es pues, una vez ejecutados los comandos, permite ver una medida exacta del tiempo que le tomó a la computadora procesar el código. Se ha añadido esto pues, como puede observarse en la descripción de las columnas, se trata de una tabla bastante extensa, con 95,000 registros.
Dependidendo del tiempo de ejecución, es recomendable no utilizar los datos en su totalidad, sino únicamente una muestra aleatoria de éstos (la cual retiene las mismas propiedades). Si le tomó a la computadora más de 20s el importar el archivo (indicado por la palabra total al final del resultado), entonces es recomendable obtener la muestra aleatoria, de modo que el análisis y el resto de la práctica sean más sencillas de trabajar.
Muestra Aleatoria de Datos¶
Una vez importados los datos, obtener una muestra aleatoria es relativamente simple. Primero se ejecutará la operación completa, para analizar a detalle cada uno de los pasos:
robos = gpd.read_file(f + 'robos_cdmx_2018.shp')
# Establecer un "Estado de Semilla" de modo que siempre se obtengan los mismos números aleatorios
np.random.seed(5432)
# Crear una secuencia con el mismo número de filas que las de la tabla original
ri = np.arange(len(robos))
# Reorganizar aleatoriamente los valores
np.random.shuffle(ri)
# Redefinir los índices de la tabla utilizando únicamente los primeros 10,000 valores
# de la secuencia, la cual se encuentra aleatoriamente reorganizada.
robos = robos.iloc[ri[:10000], :]
# Analizar el resultado
robos.info()
Analizando cada uno de los pasos a detalle:
- Primero, se crea por separado una secuencia de números, iniciando desde el cero (Python siempre inicia sus cuentas desde cero, no uno) tan largo como el número de filas en la tabla de la cual se pretende tomar la muestra. Como tal, la lista inicia en 0 y continua como 1, 2, 3, 4, 5, ..., $N-1$ (siendo $N$ la longitud de la tabla, esto es, 95,000).
- Después, en la línea 10, la lista es reacomodada al azar. Como tal, la longitud aún se mantiene de 95,000, pero el orden de los números ha sido cambiado a uno completamente aleatorio.
- Finalmente, es posible obtener la muestra, lo cual se hace en la Línea 14. Este comando se compone de dos elementos: el primero (
ri[:10000]), en el que se conservan únicamente los primeros 10,000 elementos de la lista guardada aleatoriamente (si se quiere cambiar el tamaño de la muestra, éste es el número que debe de ser cambiado); segundo (robos.iloc), que es una de las funciones de búsqueda utilizadas hasta ahora de la libreríaPandas.
El truco reside en que, al consultar la tabla con el subconjunto de 10,000 números obtenidos de la reorganización aleatoria, únicamente se conservan las filas que poseen un índice entre esos números. Esto permite dos cosas: l primero, que únicamente se obtengan 10,000 observaciones en lugar de las 95,000 originales; segundo, que la muestra conservada sea completamente aleatoria, debido a que la serie utilizada para su obtención también lo es.
Visualización de Patrones de Puntos¶
El resto de la práctica se enfocará en aprender diferentes formas de visualizar Patrones de Puntos. En particular, se considerarán dos estrategias principales: una depende de la agregación de puntos en polígonos, mientras que la otra se basa en crear superficies contínuas utilizando Estimaciones de Densidad de Kernel (KDE).
De Puntos a Polígonos¶
Dado que en prácticas anteriores se ha aprendido sobre la visualización de datos poligonales, la forma más sencilla de visualizar Patrones de Puntos es "convirtiéndolos" en polígonos y aplicándoles técnicas conocidas, como el Mapeo de Coropletas (Choropleths), para visualizar su distribución espacial. Para hacer eso, se colocará una capa de Polígonos encima de un Patrón Puntual, se unirán estos dos al asignar a cada punto el polígono al que corresponden, y se crearán Coropletas a partir de las cuentas por polígono.
Esta aproximación, aunque bastante intuitiva, también genera la pregunta ¿Qué polígonos deben de utilizarse para agregar los puntos? Idealmente, se buscarían unos límites que empaten de la forma más exacta posible con el proceso que generó a los puntos, y divida al espacio en áreas con una intensidad similar de puntos; sin embargo, esto normalmente no es posible, principalmente dado el hecho de que una de las razones por la cual se busca visualizar Patrones Puntuales es aprender a cómo generar dicho proceso, de modo que no es posible conocer a priori los polígonos que le corresopnden.
Si no se tiene el conjunto de polígonos ideal desde un inicio, pueden adoptarse dos aproximaciones más realistas: Utilizar áreas irregulares preexistentes, o crear un conjunto artificial de polígonos regulares.
Polígonos Irregulares Preexistentes¶
Para ejemplificar esta aproximación, se utilizarán las AGEB's de la Ciudad de México utilizadas en prácticas anteriores. Así que, antes que cualquier otro paso, deben de importarse éstos:
# Importar y establecer índices en un sólo paso.
agebs = gpd.read_file(f + 'agebs_cdmx.shp').set_index('cvegeo')
El siguiente paso requiere asignar cada robo a la AGEB en donde ocurrió. Esto puede ser logrado a través de una operación de SIG común llamada Point-In-Polygon; a través de QGIS, esto puede ser conseguido fácilmente colocando ambas capas en un mismo proyecto y utilizando la herramienta correspondiente en el menú de Vector (Vectos -> Data Management Tools -> Join Attributes by Location). Otra forma, más rápida y eficiente, es a través de la función Spatial Join (.sjoin()) de la librería GeoPandas; la forma de ejecutarlo es la siguiente:
# Realizar el Spatial Join
robos = gpd.sjoin(robos, agebs, op = 'intersects')
# Revisar el resultado
robos.head()
Podrá observarse que el GeoDataFrame contenido en robos ahora posee una nueva columna llamada index_right, que corresponde a la Clave Geográfica del AGEB en donde se ocurrió el robo. La lógica detrás de esta operación reside en que .sjoin() realiza la unión de atributos a través de sus propiedades espaciales, y que se encuentra dictada por el argumento op; en este caso, debido a que se seleccionó intersects, a todos los robos que intersecten con un AGEB dado se les asignará la Clave Geográfica del mismo. Es importante notar que el orden en que las tablas son colocadas dentro de la función es importante, pues éste dictará qué tabla será a la que se le asignarán los atributos espaciales.
Como tal, lo único que hace falta es renombrar la nueva columna para saber en futuros procesos qué es la información contenida en ella, lo cual puede lograrse a través de la función .rename():
# Renombrar la columna de interés
robos = robos.rename(columns = {'index_right':'ageb_urbana_cvegeo'})
# Visualizar el resultado
robos.head()
Una vez registrada la Clave Geográfica del AGEB a la que pertenece cada robo, obtener la cuenta del número de robos por polígono es más sencillo. Nuevamente se depende de la función .groupby(), que toma todos los robos de la tabla y los agrupa por la Clave Geográfica; una vez agrupados, se aplica el método .size(), que cuenta los elementos de cada grupo y regresa una columna indexada por las Claves Geográficas con todas las cuentas como sus valores. Para facilitar la generación de mapas, se asignan las cuentas a una nueva columna en la tabla original de agebs:
# Crear las cuentas
robos_agebs = robos.groupby('ageb_urbana_cvegeo').size()
# Asignar las cuentas como una nueva columna de la tabla de AGEB's
agebs['no_robos'] = robos_agebs
El código de arriba ha creado una nueva columna en la tabla de agebs que contiene el número de robos que han ocurrido dentro de cada polígono. Revisando rápidamente la tabla:
agebs.info()
Podrá notarse un pequeño problema. Como la palabra Index indica, la tabla contiene un total de 2,432 registros; sin embargo, la columna no_robos únicamente posee 1,977 registros en los que existe un valor, esto es, existen 455 elementos sin un valor en esta columna o, en otras palabras, en las cuales no ocurrió ni un solo robo, y por ello la función .groupby() no fue capaz de asignarles un valor.
Para resolver este problema, primero deben de identificarse las AGEB's en las que no ocurrió robo alguno, esto es, sin un valor en la columna no_robos, utilizando el operador .isnull():
sin_robos = agebs.loc[agebs['no_robos'].isnull() , :]
sin_robos.head()
También es posible verificar que en la tabla de robos que no existe ningún registro ocurriendo en alguna de estas AGEB's:
# El método '.index' obtiene todas las Claves Geográficas de 'sin_robos', y '[0]' selecciona la primera
robos.loc[robos['ageb_urbana_cvegeo'] == sin_robos.index[0] , :]
Puede observarse que, la búsqueda de todos los robos en la que su columna ageb_urbana_cvegeo sea igual a uno de los índices de la tabla sin_robos regresa una tabla vacía.
En el contexto de esta práctica, un área sin alguna ocurrencia del evento podría derivarse directamente en un cero para la cuenta total; sin embargo, en otros casos, esta no es necesariamente la respuesta correcta, razón por la cual pandas decide dejar un Valor Nulo (NaN) cuando no es capaz de obtener una cuenta. Para corregir este detalle, existe una función llamada .fillna(), que buscará todos los Valores Nulos de una tabla y los reemplazará con el valor que se le indique que, en este caso, es 0:
agebs = agebs.fillna(0)
Como tal, si se revisa nuevamente la tabla de agebs, ya no se encontrarán valores nulos:
agebs.info()
Dado lo anterior, es posible generar un mapa con las cuentas. Desde el punto de vista técnico, se trata de un Mapa de Coropletas, tal y como se ha realizado en prácticas anteriores:
# Preparar la figura y sus filas
fig, filas = plt.subplots(1 , figsize = (10,10))
# Graficar una Coropleta dividia con Fisher-Jenks, con su leyenda
agebs.plot(column = 'no_robos' , scheme = 'fisher_jenks' , legend = True , ax = filas , cmap = 'BuPu' , edgecolor = 'grey', linewidth = 0.2)
# Eliminar los ejes
filas.set_axis_off()
# Asignar un título
filas.set_title('Número de Robos por AGEB en la CDMX')
# Mostrar el Resultado
plt.show()
Ejercicio Opcional¶
Crea un Mapa de Coropletas similar al anterior, pero esta vez utilizando la Clasificación por Cuantiles (quantiles) e Intervalos Equivalentes (equalinterval), y responde las preguntas: ¿Cuáles son las diferencias principales? ¿A qué se deben estas diferencias? ¿Cómo se relacionan con la Distribución de las cuentas en los polígonos?
El mapa anterior muestra una concentración de robos principalmente en la Alcaldía Cuauhtémoc de la Ciudad de México; sin embargo, también muestra una alta concentración en las AGEB's de gran extensión, detalle que lleva a una importante reflexión. Como se mencionó anteriormente, los robos son un evento que, en teoría, pueden ocurrir en cualquier parte de la ciudad y, debido a esta naturaleza, el tamaño de una AGEB puede influir en el hecho que dentro de ésta ocurra el crimen o no; mientras más grande sea un AGEB, más problable es que un robo ocurra en su interior. Como tal, si no se consideran factores como éste dentro del mapa, se puede caer en el error de simplemente estar mapeando la incidencia de robos sin realizar un análisis detallado.
Para obtener un mapeado más exacto, lo que se desearía observar es la intensidad de los robos, no las cuentas brutas de éstos. Para conseguir esto, pueden recurrirse a indicadores como el Número de Robos por cada 1,000 habitantes o, como el que se utilizará en esta práctica, el Número de Robos por km² dentro del AGEB.
Para hacer esto, puede recurrirse a la propiedad .area de la tabla agebs, la cual es posible obtener gracias a la librería de GeoPandas:
agebs.area.head()
Es importante notar que el área calculada por GeoPandas depende directamente de la proyección en la que se encuentran los datos; si se observa la propiedad .crs de la tabla:
agebs.crs
Puede observarse que ésta se encuentra en la proyección EPSG: 6362, que equivale a la Cónica Conforme de Lambert (LCC) para México; esta proyección utiliza unidades en metros y, como tal, el área calculada es medida en metros cuadrados (m²). Para convertir esta área en km², únicamente es necesario realizar la conversión correspondiente, para entonces poder añadir este valor como una nueva columna en la tabla original:
# Cálculo de Área y creación de nueva columna
agebs['area'] = agebs.area / 1000000
# Visualización de resultados
agebs.head()
A partir de esto, es sencillo obtener la tasa de robos por unidad de área en cada AGEB, esto es, un indicador de intensidad:
# Cálculo de Robos por km2
agebs['robos_int'] = agebs['no_robos'] / agebs['area']
# Visualización de Resultados
agebs.head()
Dada la intensidad, es posible crear un nuevo Mapa de Coropletas, esta vez utilizando la columna robos_int como medida para la clasificación:
# Preparar la figura y sus filas
fig, filas = plt.subplots(1 , figsize = (10,10))
# Graficar una Coropleta dividia con Fisher-Jenks, con su leyenda
agebs.plot(column = 'robos_int' , scheme = 'fisher_jenks' , legend = True , ax = filas , cmap = 'BuPu' , edgecolor = 'grey', linewidth = 0.2)
# Eliminar los ejes
filas.set_axis_off()
# Asignar un título
filas.set_title('Número de Robos por km² en cada AGEB de la CDMX')
# Mostrar el Resultado
plt.show()
Dado este ajuste, es posible observar que las AGEB's más grandes ya no necesariamente son las de mayor intensidad, pues el indicador permite deducir que su tamaño es tan grande que el hecho de que un robo haya ocurrido en su interior puede derivarse meramente de este detalle. Por otra parte, las AGEB's en la Alcaldía Cuahtémoc continúan con una coloración intensa, mostrando que, aún siendo pequeñas, el número de robos que han ocurrido en su interior puede derivarse de factores más allá de su tamaño.
Ejercicio Opcional¶
Nuevamente, crea un Mapa de Coropletas similar al anterior, pero esta vez utilizando la Clasificación por Cuantiles (quantiles) e Intervalos Equivalentes (equalinterval), y responde las preguntas: ¿Cuáles son las diferencias principales? ¿A qué se deben estas diferencias? ¿Cómo se relacionan con la Distribución de las cuentas en los polígonos?
Polígonos Regulares - Hex-Binning¶
En algunos casos, se tiene que los polígonos irregulares disponibles no son lo suficientemente adecuados para agregar los puntos en ellos, o simplemente no se cuenta con ningún tipo de área de agregación; bajo estas circunstancias, una alternativa válida es crear una topología artificial de polígonos que pueda ser utilizada para agregar los puntos. Existen múltiples formas de hacer esto, pero la más común es crear una Red de Hexágonos; ésta proporciona una topología regular, pues cada polígono es de la misma forma y tamaño, lo cual trae ventajas sobre trabajar con círculos, pues cubre todo el espacio sin sobreponer figuras, y posee más vértices que los cuadrados, resolviendo problemas relacionados a éstos.
A través de Python existe una forma simplificada de crear esta capa de hexágonos y agregar los puntos en ellos a través de una sola línea, gracias a la función .hexbin() de la librería matplotlib. Éstos son los comandos para crear la Capa Hexagonal:
# Preparación de la Figura y sus Filas
fig , filas = plt.subplots(1 , figsize = (10,10))
# Añadir capa de hexágonos que muestra la cuenta de puntos en cada polígono
hb = filas.hexbin(robos['geometry'].x , robos['geometry'].y, gridsize = 50, alpha = 0.8, cmap = 'BuPu')
# Añadir Barra de Color de referencia (Opcional)
plt.colorbar(hb)
Puede notarse que lo único que se necesita es llamar directamente a la función .hexbin() sobre las filas (filas) de la figura que se está creando, añadir las coordenadas de los puntos (llamadas con robos['geometry'].x y robos['geometry'].y); los argumentos adicionales manejan el número de hexágonos por cada eje (gridsize, teniendo 50 para una capa de 50x50 polígonos), la transparencia deseada (alpha) y los colores a utilizar en el mapa (cmap). Adicionalmente, se incluye una Barra de Colores para entender el significado de los colores; debe notarse que se realiza colocando el objeto en el que se ha almacenado el resultado de hexbin, siendo en este caso la variable hb.
Teniendo la base, es posible seguir la rutina de Generación de Mapas utilizada hasta ahora para generar un mapa completo del número de robos en la Ciudad de México:
# Preparación de la figura y sus filas
fig , filas = plt.subplots(1, figsize = (10,10))
# Añadir capa de hexágonos que muestra la cuenta de puntos en cada polígono
hb = filas.hexbin(robos['geometry'].x , robos['geometry'].y, gridsize = 50, alpha = 0.8, cmap = 'BuPu')
# Añadir Barra de Color de referencia (Opcional)
plt.colorbar(hb)
# Añadir una Capa Base con las AGEB's
agebs.plot(ax = filas, facecolor = '', linewidth = 0.2, edgecolor = 'grey', alpha = 0.5)
# Remover los ejes
filas.set_axis_off()
# Añadir título al mapa
filas.set_title('Mapa tipo Hex-Binning de los Robos en la CDMX')
# Mostrar el resultado
plt.show()
Estimación de Densidad de Kernel¶
¡Importante! - Es recomendable que, para esta sección, se utilice la muestra aleatoria de robos en lugar de todo el conjunto de 95,000 registros.
Utilizar la técnica de 'Hexagonal-Binning' puede resultar una solución rápida cuando no se tiene una buena capa de polígonos para agregar los puntos, y alguna de sus propiedades, como la equivalencia en tamaños de cada polígono, puede ayudar a aliviar algunos de los problemas asociados a una "mala" toplogía irregular (esto es, una que no se ajusta perfectamente al proceso que generó el Patrón de Puntos). Sin embargo, esto aún no resuelve del todo el dilema derivado del Problema de Unidad de Área Modificable (MUAP, por sus siglas en inglés, estudiado en prácticas posteriores); después de todo, aún se están imponiendo límites arbitrarios y agregando según éstos, de modo que la posibilidad de no acercarse del todo a la distribución correcta del Patrón de Puntos es bastante probable.
Una forma de resolver este problema es evitando la agregación completamente. En su lugar, puede intentarse estimar una variable contínua, como lo es la probabilidad de observar el evento de interés. El método más común utilizado para realizar esto es la Estimación de la Densidad de Kernel (KDE, por sus siglas en inglés), método que busca contar el número de puntos de forma contínua; en lugar de utilizar una cuenta discreta, en la que se incluye un punto en la cuenta si éste se encuentra dentro de los límites de un polígono, el método KDE utiliza funciones (kernel) que incluye a los puntos, pero les asigna diferentes pesos dependiendo de su distancia al punto desde el cual se está realizando la cuenta.
Aunque el algoritmo para realizar esta estimación no es nada trivial, dentro de Python se trata de una tarea muy simplificada gracias a la existencia de la función .kdeplot() de la librería seaborn. El comando para crear un mapa de Densidad de Kernel de forma muy sencilla es:
sns.kdeplot(robos['geometry'].x, robos['geometry'].y, n_levels = 50, shade = True, cmap = 'BuPu')
La librería seaborn sintetiza enormemente todo el proceso en una sola línea. La función .kdeplot() (que también puede utilizarse para crear una Estimación de Densidad de Kernel de una sola variable) toma las coordenadas X (robos['geometry'].x) y Y (robos['geometry'].y) de los puntos como sus únicos argumentos obligatorios; adicionalmente, se especifican el número de niveles en el que se busca que se coloree el gradiente (n_levels), si se busca colorear el espacio entre cada uno de los niveles (shade) y los colores de preferencia (cmap).
Una vez entendido el funcionamiento básico de la herramienta, puede ser colocada dentro de un mapa más completo, siendo la diferencia más importante del proceso que debe de indicársele a .kdeplot() dónde debe de ser colocado (en este caso, las filas de la figura):
# Preparación de figuras y filas
fig, filas = plt.subplots(1, figsize = (10,10))
# Generar el KDE
sns.kdeplot(robos['geometry'].x, robos['geometry'].y, n_levels = 50, shade = True, cmap = 'BuPu', ax = filas)
# Añadir una Capa Base con las AGEB's
agebs.plot(ax = filas, facecolor = '', linewidth = 0.2, edgecolor = 'grey', alpha = 0.75)
# Remover los ejes
filas.set_axis_off()
# Añadir título al mapa
filas.set_title('Mapa tipo KDE de los Robos en la CDMX')
# Mostrar el resultado
plt.show()
Agregación de Puntos¶
En esta última sección, se aprenderá el método para identificar Agregaciones (Clusters) de puntos, basándose en su densidad a lo largo del espacio. Para esto, se utilizará el conocido algoritmo DBSCAN; en éste, se considera como Cluster a una concentración de por lo menos $m$ puntos, cada uno con una distancia de $r$ a por lo menos uno de los puntos de la misma agregacion. A partir de esto, los puntos del conjunto son dividos en tres categorías:
- Ruido (Noise) - Aquellos puntos fuera del Cluster.
- Núcleos (Cores) - Puntos dentro de una agregación que cumpla con tener por lo menos $m$ puntos dentro de una distancia de $r$.
- Fronteras (Borders) - Puntos dentro de una agregación con menos de $m$ puntos dentro del mismo en una distancia de $r$.
Tanto $m$ y $r$ necesitan ser especificados con anterioridad por el usuario antes de ejecutar el algoritmo DBSCAN; esto es un factor clave, pues el valor de estos parámetros puede influir significativamente en el resultado final. Antes de continuar explorando a detalle, resulta pertinente realizar una primera ejecución del algoritmo.
Conceptos Básicos¶
Para facilitar la ejecución del algoritmo, es de utilidad adicionar al GeoDataFrame de robos una nueva columna que almacene las Coordenadas X y Y de cada uno de los registros; esto puede hacerse fácilmente gracias a la capacidad de GeoPandas de obtener las propiedades .x y .y de las geometrías:
# Obtener Coordenadas X
robos['X'] = robos['geometry'].x
# Obtener Coordenadas Y
robos['Y'] = robos['geometry'].y
# Visualizar resultados
robos.head()
Con lo anterior, ejecutar el algoritmo DBSCAN es muy sencillo gracias a la función dbscan(), presente en la librería scikit-learn, al punto que es necesaria una única línea:
# Ejecutar DBSCAN
nucleos , etiquetas = dbscan(robos[['X' , 'Y']])
La función requiere de dos objetos. El primero, llamado nucleos, contiene los índices ordenados, iniciando desde el cero, de cada uno de los puntos que han sido clasificados como Núcleo (Core). Es posible observar la apariencia de esta variable:
# Observar los primeros cinco elementos de 'nucleos'
nucleos[:5]
Lo anterior indica que el elemento Número 169 (que posee el índice 168, debido a que Python numera desde cero) del conjunto de datos es un Núcleo, como lo son el No. 176, 177, 178 y 180. El objeto nucleos siempre tendrá una dimensión variable, pues dependerá del número de núcleos que el algoritmo encuentre.
Ahora, si se observa el objeto etiquetas:
# Observar los primeros cinco elementos de 'etiquetas'
etiquetas[:5]
El objeto etiquetas siempre tendrá la misma dimensión que el número de puntos utilizado para ejecutar el alboritmo DBSCAN; cada valor de la lista representa el índice del clúster al que pertenece un punto. Si un punto es clasificado como Ruido (Noise), entonces recibirá un valor de -1; como tal, puede observarse que los primeros cinco puntos analizados por el algoritmo no pertenecen a algún agregado. Para facilitar el trabajo a futuro, se convertirá etiquetas en un objeto del tipo Series cuyos índices sean los mismos que los puntos introducidos al algoritmo:
etiquetas = pd.Series(etiquetas , index = robos.index)
etiquetas.head()
Una vez obtenidos los Clusters, es posible visualizarlos. Existen múltiples formas para conseguir esto, siendo la más sencilla el colorear todos los puntos (Núcleos) pertenecientes a un Clúster de rojo, y el Ruido de gris:
# Preparación de Figura y Filas
fig , filas = plt.subplots(1, figsize = (10,10))
# Subconjunto de Puntos que no son parte de algún Cluster (Ruido)
ruido = robos.loc[etiquetas == -1 , ['X' , 'Y']]
# Graficar el Ruido de color gris
filas.scatter(ruido['X'] , ruido['Y'], c = 'grey' , s = 5 , linewidth = 0)
# Graficar todos los puntos que no sean Ruido de color rojo
# Cabe destacar el uso de '.index.difference' en el que se obtienen el índice de todos los puntos en 'robos'
# y se les resta el de aquellos que pertenezcan a 'ruido'
filas.scatter(robos.loc[robos.index.difference(ruido.index) , 'X'] , robos.loc[robos.index.difference(ruido.index) , 'Y'] , c= 'red' , linewidth = 0)
# Mostrar el resultado
plt.show()
Aunque el resultado es informativo, no se trata de uno completamente satisfactorio, pues no puede asegurarse a simple vista que el tamaño de las Agregaciones cubra la totalidad de una zona en particular, además de que se muestran Clusters en zonas en las que aparentemente no deberían de existir algunos; esto se debe a que el algoritmo DBSCAN ha sido ejcutado utilizando los parámetros por defecto del mismo.
Si se escribe dbscan? en alguna de las celdas, se obtendrá información auxiliar de la función, donde puede verificarse la afirmación anterior: por defecto, se tiene un radio (argumento llamado eps) de 0.5 y un mínimo de 5 puntos por agregación (min_samples); debido a que la proyección de los datos utiliza metros, el radio de 0.5m únicamente arrojará Clusters en extremo locales, lo cual puede ser de utilidad para algunos proyectos, pero puede producir resultados extraños en otros.
Como tal, deben de modificarse estos parámetros para poder obtener patrones más generales. Por ejemplo, si se desea establecer que el Clúster contenga, por lo menos, al 1% del total de puntos en los datos:
# Obtener un mínimo de puntos igual al 1% del total de la muestra
minp = np.round(len(robos) * 0.01)
minp
Asimismo, puede aumentarse el Radio Máximo a 2.5km, siendo éste un valor arbitrario. Como tal, puede volverse a ejecutar algoritmo y graficar el resultado, todo dentro de la misma celda:
# Ejecutar DBSCAN
nucleos , etiquetas = dbscan(robos[['X' , 'Y']] , eps = 2500 , min_samples = minp)
# Convertir la lista 'etiquetas' en una Serie
etiquetas = pd.Series(etiquetas , index = robos.index)
# Preparación de Figura y Filas
fig , filas = plt.subplots(1, figsize = (10,10))
# Subconjunto de Puntos que no son parte de algún Cluster (Ruido)
ruido = robos.loc[etiquetas == -1 , ['X' , 'Y']]
# Graficar el Ruido de color gris
filas.scatter(ruido['X'] , ruido['Y'], c = 'grey' , s = 5 , linewidth = 0)
# Graficar todos los puntos que no sean Ruido de color rojo
# Cabe destacar el uso de '.index.difference' en el que se obtienen el índice de todos los puntos en 'robos'
# y se les resta el de aquellos que pertenezcan a 'ruido'
filas.scatter(robos.loc[robos.index.difference(ruido.index) , 'X'] , robos.loc[robos.index.difference(ruido.index) , 'Y'] , c= 'red' , linewidth = 0)
# Mostrar el resultado
plt.show()
Ahora, el resultado es completamente diferente, teniendo lo que parece múltiples Clústers en prácticamente toda la extensión de la Ciudad de México. Esto ejemplifica cómo los parámetros pueden modificar sustancialmente el resultado, aún teniendo los mismos datos y ejecutando el mismo algoritmo.
Ejercicio Opcional¶
Intenta crear mapas similares a los anteriores, pero esta vez también representando aquellos puntos que son clasificados como Fronteras (Borders) por el algoritmo.
Para Practicar...¶
Reproduce el Análisis de Patrones de Puntos realizado en la práctica, esta vez con un conjunto de datos de tu elección. Es necesario:
- Obtener los datos a analizar.
- Importar los datos.
- Si encuentras Polígonos Irregulares aptos para agregar los puntos:
- Utilizar la función
.sjoin()para realizar el Spatial Join entre Puntos y Polígonos. - Obtener las cuentas del número de puntos por polígono.
- Crear un Mapa de Coropletas con las cuentas en bruto.
- Generar una medida de intensidad y un nuevo Mapa de Coropletas.
- Utilizar la función
- Crear un Mapa de Hex-Binning.
- Calcular y mostrar la Estimación de Densidad de Kernel (KDE) de la distribución de los puntos.
- Obtener Clusters a través de DBSCAN
Como recomendación, puedes utilizar los siguientes datos:
- Tipo de Uso de Suelo por Vivienda en la Ciudad de México (Liga)
- Carpetas de Investigación de delitos en la Ciudad de México (Liga)
En ambos casos, el portal ofrece diversos filtros para descargar el subconjunto deseado, de modo que no sea necesario trabajar desde un inicio con todos los datos existentes. Es recomendable utilizar estos filtros, pues ambas bases de datos son bastante extensas.