Introducción
En este trabajo intentaremos analizar el comportamiento geo-espacial de delitos denunciados en CABA durante 2022. Se trabajará con el dataset provisto por el portal de datos abiertos del GCBA (disponible aquí y un dataset de radios censales. El análisis es exploratorio e intentará detectar zonas con alto volumen de delitos y zonas con comportamientos particulares.
Importo librerías
# !pip install mapclassify folium pointpats contextily pygeos h3 pysalimport pandas as pd
import numpy as np
import geopandas as gpd
import os
import matplotlib.pyplot as plt
from matplotlib.patches import Polygon as PolygonM #noten como usamos otro alias
from matplotlib.patches import Patch
from shapely.geometry import Polygon
from splot.esda import plot_moran_bv_simulation, plot_moran_bv, moran_scatterplot
import seaborn as sns
import contextily as ctx
import libpysal as lps
import esda
from esda.moran import Moran_BV, Moran_Local_BVCargo data
data = pd.read_csv('./delitos/delitos_2022.csv', delimiter=",", encoding="Latin1")print(data.shape)## (112961, 16)print(data.dtypes)## Unnamed: 0 int64
## id.mapa int64
## anio int64
## mes object
## dia object
## fecha object
## franja int64
## tipo object
## subtipo object
## uso_arma object
## uso_moto object
## barrio object
## comuna float64
## latitud float64
## longitud float64
## cantidad int64
## dtype: objectdata.head()## Unnamed: 0 id.mapa anio mes ... comuna latitud longitud cantidad
## 0 1 1 2022 OCTUBRE ... 15.0 -34.584136 -58.454704 1
## 1 2 2 2022 OCTUBRE ... 4.0 -34.645043 -58.373194 1
## 2 3 3 2022 NOVIEMBRE ... 15.0 -34.589982 -58.446471 1
## 3 4 5 2022 NOVIEMBRE ... 2.0 -34.596748 -58.413609 1
## 4 5 6 2022 MAYO ... 9.0 -34.640978 -58.480254 1
##
## [5 rows x 16 columns]# Leemos los datos por radio censal
radios = gpd.read_file('./radios/radios.gpkg')print(radios.shape)## (3554, 17)print(radios.dtypes)## ID int64
## CO_FRAC_RA object
## COMUNA int32
## FRACCION int64
## RADIO int32
## TOTAL_POB int64
## T_VARON int64
## T_MUJER int64
## T_VIVIENDA float64
## V_PARTICUL int64
## V_COLECTIV float64
## T_HOGAR float64
## H_CON_NBI float64
## H_SIN_NBI float64
## univComP float64
## univComp_1 float64
## geometry geometry
## dtype: objectradios.head()## ID CO_FRAC_RA ... univComp_1 geometry
## 0 1 1_1_1 ... 7.0 MULTIPOLYGON (((6374084.199 6171829.069, 63740...
## 1 2 1_12_1 ... 27.0 MULTIPOLYGON (((6372815.703 6170430.589, 63728...
## 2 3 1_12_10 ... 29.0 MULTIPOLYGON (((6373471.888 6170345.81, 637359...
## 3 4 1_12_11 ... 32.0 MULTIPOLYGON (((6374523.865 6170322.842, 63745...
## 4 5 1_12_2 ... 39.0 MULTIPOLYGON (((6372943.061 6170441.384, 63730...
##
## [5 rows x 17 columns]Transformo df a geo
Para trabajar con estadística espacial, debemos pasar nuestro dataframe a geodataframe, es decir, debemos tener al menos una variable con datos georreferenciados. Vamos a construir esta variable a partir de las variables latitud y longitud y le daremos un sistema de coordenadas al dataframe.
def convertir_a_float(valor):
#elimino todos los puntos del string y luego aplico un punto después de los dos primeros números,luego convierto a float
valor = str(valor).replace('.', '') # elimino todos los puntos
if len(valor) > 3:
valor = valor[:3] + '.' + valor[2:] # inserto un punto después de los dos primeros números
try:
return float(valor)
except ValueError:
return None
# aplico la función a las columnas de latitud y longitud
data['latitud2'] = data['latitud'].apply(convertir_a_float)
data['longitud2'] = data['longitud'].apply(convertir_a_float)
# elimino filas con valores no válidos en latitud o longitud
data = data.dropna(subset=['latitud', 'longitud'])data_geo = gpd.GeoDataFrame(data, geometry = gpd.GeoSeries.from_xy(x = data.longitud, y = data.latitud, crs = 4326))print(data_geo.dtypes)## Unnamed: 0 int64
## id.mapa int64
## anio int64
## mes object
## dia object
## fecha object
## franja int64
## tipo object
## subtipo object
## uso_arma object
## uso_moto object
## barrio object
## comuna float64
## latitud float64
## longitud float64
## cantidad int64
## latitud2 float64
## longitud2 float64
## geometry geometry
## dtype: objectdata.head()## Unnamed: 0 id.mapa anio ... cantidad latitud2 longitud2
## 0 1 1 2022 ... 1 -34.458414 -58.845470
## 1 2 2 2022 ... 1 -34.464504 -58.837319
## 2 3 3 2022 ... 1 -34.458998 -58.844647
## 3 4 5 2022 ... 1 -34.459675 -58.841361
## 4 5 6 2022 ... 1 -34.464098 -58.848025
##
## [5 rows x 18 columns]Exploro
Veamos primero la densidad poblacional en CABA por radio censal
radios['densidad'] = (radios.TOTAL_POB / (radios.geometry.area / 1000000)).round().astype(int)radios.explore(column = 'densidad', scheme = 'boxplot', tiles = 'CartoDB positron', cmap = 'PuOr', legend = True)Ahora veamos la distribución espacial del volumen de delitos denunciados. Para esto vamos a crear un mapa de calor
f, ax = plt.subplots(1, figsize=(9, 9))
# Generamos un "heatmap"
sns.kdeplot(
x="longitud",
y="latitud",
data=data[1:50000], #Uso una porción de datos porque es algo muy costoso de computar
n_levels=20, #pueden ajustar este parametro
fill=True,
alpha=0.55,
cmap="coolwarm",
)
ctx.add_basemap(
ax, source=ctx.providers.CartoDB.Positron, crs="EPSG:4326"
)
ax.set_axis_off()
plt.show() 
En prinicpio podemos observar que hay cierta coincidencia entre la densidad porblacional y el volumen de delitos denunciados. Esta obervación tiene lógica a priori, ya que donde más personas se encuentran, mayor probabilidad de delitos puede haber (lo mismo se puede pensar para otro tipo de unidad de análisis, donde mayor concentración de personas haya, mayor cantidad de casos de una variable podremos encontrar). Sin embargo, en el mapa de calor, también obervamos algunos focos en zonas con menor densidad poblacional, como Liniers o Constitución. Indagaremos más sobre este comportamiento.
Procesamiento
Join espacial
Primero crearemos un join espacial entre ambos dataframes, a fin de ubicar los puntos que representan a los delitos dentro de cada polígono de los radios censales. Para esto primero debo analizar si el sistema de coordenadas de cada dataframe coincide
print(data_geo.crs)## EPSG:4326print(radios.crs)## PROJCS["POSGAR_2007_Argentina_6",GEOGCS["GCS_POSGAR 2007",DATUM["Posiciones_Geodesicas_Argentinas_2007",SPHEROID["GRS 1980",6378137,298.257222101],AUTHORITY["EPSG","1062"]],PRIMEM["Greenwich",0],UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",-90],PARAMETER["central_meridian",-57],PARAMETER["scale_factor",1],PARAMETER["false_easting",6500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH]]Normalizo el sistema de coordenadas y procedo con el join espacial
radios = radios.to_crs(epsg=4326)# realizo el join espacial
data_geo_radios = gpd.sjoin(data_geo, radios, how="left", predicate="within")
# verifico el resultado
print(data_geo_radios.head())## Unnamed: 0 id.mapa anio ... univComP univComp_1 densidad
## 0 1 1 2022 ... 7.0833 51.0 19632.0
## 1 2 2 2022 ... 9.1916 83.0 23777.0
## 2 3 3 2022 ... 12.4855 108.0 27084.0
## 3 4 5 2022 ... 15.8756 194.0 45067.0
## 4 5 6 2022 ... 5.5604 64.0 16813.0
##
## [5 rows x 37 columns]Ahora haremos un conteo de delitos por cada radio censal y sumaremos esa información como una variable más al dataframe de radios censales
# contar filas por CO_FRAC_RA
conteo_por_fraccion = data_geo_radios.groupby('CO_FRAC_RA')['CO_FRAC_RA'].count()
# agregar la columna de conteo al df radios
radios = radios.merge(conteo_por_fraccion.rename('conteo_delitos'), left_on='CO_FRAC_RA', right_index=True, how='left')
# verificaamos el resultado
print(radios.head())## ID CO_FRAC_RA ... densidad conteo_delitos
## 0 1 1_1_1 ... 187 696.0
## 1 2 1_12_1 ... 18360 72.0
## 2 3 1_12_10 ... 6667 62.0
## 3 4 1_12_11 ... 1441 263.0
## 4 5 1_12_2 ... 12465 36.0
##
## [5 rows x 19 columns]# removemos filas con valores nulos o NA para no entorpecer el posterior análisis
radios = radios.dropna(subset=['conteo_delitos'])
# verificaamos el resultado
print(radios.head())## ID CO_FRAC_RA ... densidad conteo_delitos
## 0 1 1_1_1 ... 187 696.0
## 1 2 1_12_1 ... 18360 72.0
## 2 3 1_12_10 ... 6667 62.0
## 3 4 1_12_11 ... 1441 263.0
## 4 5 1_12_2 ... 12465 36.0
##
## [5 rows x 19 columns]MORAN
Vamos a intentar analizar la relación entre los radios censales y los delitos denunciados en cada uno, aplicando una técnica de autocorrelación espacial y el indicador I de Moran. Primero construiremos una matriz de adyacencia de polígonos con un formato Queen (reina).
# Construimos la matriz
wq = lps.weights.Queen.from_dataframe(radios,ids ='CO_FRAC_RA')## E:\Programas\Anaconda\Lib\site-packages\libpysal\weights\contiguity.py:347: UserWarning: The weights matrix is not fully connected:
## There are 2 disconnected components.
## W.__init__(self, neighbors, ids=ids, **kw)wq.transform = 'r'Ahora crearemos el indicador y evaluaremos si existe una estructura espacial de nuestros datos.
#calcular estadistico global para saber si existe una estructura espacial
y = radios.conteo_delitos
np.random.seed(12345)
mi = esda.moran.Moran(y, wq)
print('Moran Global: ', mi.I, 'p-value: ',mi.p_sim)## Moran Global: 0.26091980771631124 p-value: 0.001sns.kdeplot(mi.sim, fill=True)
plt.vlines(mi.I, 0, 1, color='r')
plt.vlines(mi.EI, 0,1)
plt.xlabel("Moran's I")
plt.show();
Si bien el indicador quedó con un valor bajo, el p-value indica que existe correlación entre los polígonos por lo que seguiremos con el análisis planteado.
# Calculamos un Moran local
li = esda.moran.Moran_Local(y, wq)fig, ax = moran_scatterplot(li, p=0.05)
ax.set_xlabel('Delitos')
ax.set_ylabel('Lag espacial de Delitos')
plt.show();
rojo = '#d7191c'
celeste = '#abd9e9'
azul = '#2c7bb6'
naranja = '#fdae61'
legend_elements = [
Patch(facecolor=rojo, edgecolor='w',label='Alto - Bajo'),
Patch(facecolor=naranja, edgecolor='w',label='Hotspot'),
Patch(facecolor=celeste, edgecolor='w',label='Bajo - Alto'),
Patch(facecolor=azul, edgecolor='w',label='Coldspot')
]# Este análisis devuelve una clasificación de cada celda en los 4 clusters
sig = li.p_sim < 0.05
radios['cuadrante'] = sig * li.q
radios['color'] = radios['cuadrante'].replace({1:rojo, 2:celeste, 3:azul, 4:naranja})
radios.loc[radios.cuadrante == 0,'color'] = np.nan
radios['cuadrante'] = radios['cuadrante'].replace({1:'hotspot',2:'bajo_alto',3:'coldspot',4:'alto_bajo'})
radios.cuadrante.value_counts(normalize = True)## cuadrante
## 0 0.781047
## coldspot 0.104668
## hotspot 0.066761
## bajo_alto 0.039604
## alto_bajo 0.007921
## Name: proportion, dtype: float64Mapeamos los resultados del Moran local
import folium
# Crear un mapa base centrado en la ubicación deseada
m = folium.Map(location=[-34.60, -58.44], zoom_start=12) # Coordenadas de ejemplo para CABA
# Añadir capa base
folium.TileLayer('CartoDB.Positron').add_to(m)
# Agregar capas según cuadrantes, respetando el color asignado## <folium.raster_layers.TileLayer object at 0x00000000932AB910>for cuadrante, color in radios[['cuadrante', 'color']].dropna().drop_duplicates().values:
subset = radios[radios['cuadrante'] == cuadrante].dropna()
# Crear una capa para cada cuadrante con el color correspondiente y un label
folium.GeoJson(
subset,
name=f"{cuadrante}", # Nombre de la capa basado en el cuadrante
style_function=lambda feature, col=color: {'color': col, 'weight': 2, 'fillOpacity': 0.6},
tooltip=folium.features.GeoJsonTooltip(
fields=["cuadrante", "COMUNA", "TOTAL_POB", "densidad", "conteo_delitos"], # Ajusta los nombres de las variables que quieras mostrar en el label
aliases=["Cuadrante", "Comuna", "Habitantes", "Densidad Poblacional", "Total delitos denunciados"], # Ajusta los alias según cómo quieras que se vean en el tooltip
localize=True
)
).add_to(m)
# Añadir el control de capas al mapa después de añadir las capas## <folium.features.GeoJson object at 0x0000000093542D90>
## <folium.features.GeoJson object at 0x00000000939122D0>
## <folium.features.GeoJson object at 0x000000009EEDDCD0>
## <folium.features.GeoJson object at 0x00000000932A8D90>folium.LayerControl(collapsed=False).add_to(m) # `collapsed=False` para intentar expandir el control
# Mostrar el mapa## <folium.map.LayerControl object at 0x000000009EEDC650>mConclusión
Finalmente, si observamos los radios que quedaron categorizados con el valor alto_bajo (es decir, polígonos con alta cantidad de delitos y rodeado de polígonos con baja cantidad de delitos) podemos ver que en general se encuentran en zonas residenciales, de baja densidad poblacional pero mucho movimiento o circulación de personas, como ser el caso del polígono de Barrancas de Belgrano o la estación de Villa del Parque, ambos radios con altos niveles de circulación y tránsito en una zona residencial. Esta obervación podría sugerir un futuro análisis comparativo de radios censales teniendo en cuenta la densidad poblacional y la circulación o tránisto de personas, aplicando un análisis Moran bivariado.