LA Crime: Ahora en 3D (no se requieren gafas) |  de Aleksei Rozanov |  junio de 2024

Visualizando geodatos del crimen en Python

Imagen por autor.

En En mi opinión, una de las mejores características que poseen los datos geográficos es su capacidad de visualización 3D. Sin embargo, debido a la cantidad de recursos computacionales necesarios para tales cálculos, rara vez se realiza en Python (a menudo javascript y sus bibliotecas utilizadas como alternativas). En uno de mis artículos anteriores, compartí seis paquetes de Python que permiten la creación de hermosos mapas estáticos e interactivos, pero solo en un espacio 2D.

Hoy quiero cerrar esta brecha e investigar con ustedes un marco realmente elegante y eficiente para visualizaciones basadas en web de alto rendimiento. cubierta.jlque también tiene una biblioteca de Python PyDeck.

Para explorar adecuadamente sus capacidades en Python necesitamos un gran conjunto de datos geoespaciales. Un candidato perfecto para eso es Datos sobre criminalidad de Los Ángeles 2010-2020 conjunto de datos de Kaggle. Por suerte, tiene un licencia abiertapor lo que podríamos utilizarlo libremente para nuestros fines.

Los autores distribuyen dos csv archivos, que vamos a fusionar en uno solo filtrando todas las columnas excepto longitud y latitud (coordenadas de los lugares donde se cometieron los crímenes).

🐍El código Python completo: GitHub.

import numpy as np
import pandas as pd
import xarray as xr
import geopandas as gpd
import random
import math

import matplotlib.pyplot as plt
from shapely import Point

import cartopy
import cartopy.crs as ccrs
import cartopy.feature as cfeature

import warnings

warnings.filterwarnings('ignore')

recent = pd.read_csv('./LA Crimes/Crime_Data_from_2020_to_Present.csv')[['LAT', 'LON']]
old = pd.read_csv('./LA Crimes/Crime_Data_from_2010_to_2019.csv')[['LAT', 'LON']]
df = pd.concat([old, recent])
df = df[(df.LON!=0) & (df.LAT!=0)] #zeros are Nans according to meta info

Después de cargarlo en pandas, quiero realizar una visualización 2D estática usando cartopía, sólo para tener una referencia confiable. Si simplemente trazamos los datos, obtendremos un montón de puntos de datos, que no nos sirven de nada.

Imagen por autor.

En su lugar, realicemos una interpolación espacial utilizando el método NN (puedes leer más sobre esto en mi otro artículo).

Básicamente, significa que transferimos observaciones escasas a una cuadrícula geográfica (PyDeck hará lo mismo, en este caso se puede llamar agregación de datos).

def coords(x,y, base=0.01):
x, y = round(base * math.ceil(abs(x)/base),2), round(base * math.ceil(y/base),2)
return (y,x)

def NN(data, LAT, LON):
array = np.zeros((LAT.shape[0], LON.shape[0]),dtype=int)
onGrid = data.apply(lambda row: coords(row.LAT, row.LON, 0.01), axis = 1).value_counts()
for coor in onGrid.index:
lon_idx, lat_idx = np.where(LON==coor[0]), np.where(LAT==coor[1])
array[lat_idx,lon_idx] = int(onGrid[coor])
return array

Una vez finalizado el algoritmo (deberá esperar un tiempo, ya que tenemos más de 2 mil filas para procesar), podemos empaquetar sus resultados en un matriz x conjunto de datos y mapearlo.

LAT, LON = np.arange(round(df.LAT.min()), round(df.LAT.max()), 0.01).astype(np.float32), np.arange(round(df.LON.min()), round(df.LON.max()), 0.01).astype(np.float32)
crimes = NN(df, LAT, LON)
ds = xr.Dataset(
{'Crimes': (['lat', 'lon'], crimes)},
coords={'lat': LAT, 'lon': LON})

fig, ax = plt.subplots(subplot_kw=dict(projection=ccrs.PlateCarree()), figsize=(16, 9))

ds.Crimes.plot(ax=ax, cmap='Reds')
ax.set_extent([-118.9, -118.1, 33.6, 34.5 ], crs=ccrs.PlateCarree())
ax.gridlines(draw_labels=True,linewidth=2, color='black', alpha=0.5, linestyle='--')
ax.add_feature(cfeature.BORDERS, edgecolor='black', linewidth=1)

ax.add_feature(cfeature.COASTLINE, edgecolor='black', linewidth=1)
ax.add_feature(cartopy.feature.RIVERS, edgecolor='blue', linewidth=0.5)
states_provinces = cfeature.NaturalEarthFeature(
category='cultural', name='admin_1_states_provinces',
scale='10m', facecolor='none')

plt.show()

Imagen por autor.

Desde mi punto de vista, parece agradable e informativo, pero si necesitas vender este proyecto a alguien, lo más probable es que fracases con ese mapa xD. ¡Así que instalemos PyDeck y veamos qué puede hacer por nosotros!

El primer tipo de capa de mapa que encontré hermosa fue la capa Hexágono. Debe especificarlo al crear la L.ayer variable en PyDeck. Hay varios otros parámetros cruciales para nosotros:

  • radio (el radio del hexágono en metros, es decir, resolución espacial en m);
  • elevación_scale (factor de escala para contenedores, cuanto mayor es el valor, mayores son los hexágonos);
  • rango_elevación (altura mínima y máxima);
  • seleccionable (mostrando los valores de forma interactiva);
  • extruido (elevación de celda).
layer = pdk.Layer(
'HexagonLayer',
df,
get_position=['LON', 'LAT'],
radius=500, #bin radius
auto_highlight=True,
elevation_scale=50, #scale factor for bins (the greater - the higher)
elevation_range=[0, 3000],
pickable=True,
extruded=True,#cell elevation
)

La segunda variable que necesitamos crear es el estado de Vista. Necesitamos alimentarlo con:

  • longitud y latitud;
  • zoom (zoom inicial);
  • min_zoom y max_zoom;
  • rumbo (ángulo de visión izquierdo/derecho);
  • tono (ángulo de visión arriba/abajo).
view_state = pdk.ViewState(
longitude=-118.3,
latitude=34.4,
zoom=8,
min_zoom=6,
max_zoom=15,
bearing=-20,#left/right angle
pitch=20, #up/down angle
)

Dado que tenemos un conjunto de datos grande, Google Colab no muestra todo el conjunto de datos, por lo que tiene dos opciones:

  • Muestra N filas del conjunto de datos;
  • Guarde el mapa en html y ábralo en su navegador.

Si eliges el segundo, obtendrás esto:

Imagen por autor.

Para ser honesto, me encanta cómo se ven los hexágonos, pero nunca los he visto en ningún artículo, presentación o conferencia científica, por lo que recomendaría usarlos conscientemente.

Ahora intentemos crear una visualización similar, pero usando columnas. Pero en este caso necesitaremos pasar a la función el conjunto de datos xarray que creamos anteriormente y especificar el color y la variable a visualizar:

layer = pdk.Layer(
'ColumnLayer',
ds.to_dataframe().reset_index(),
get_position=['lon', 'lat'],
get_elevation='Crimes',
elevation_scale=10,
radius=200,
get_fill_color=['Crimes', 220],
pickable=True,
extruded=True,
)
Imagen por autor.

En esencia, el diagrama de dispersión es una nube de puntos, pero los autores de PyDeck desarrollaron cilindros que parecen bastante únicos:

Imagen por autor.
layer = pdk.Layer(
'ColumnLayer',
df[:15000],
get_position=['LON', 'LAT'],
auto_highlight=True,
get_radius=200, # Radius is given in meters
get_fill_color=[180, 0, 200, 140], # Set an RGBA value for fill
pickable=True)
Imagen por autor.

Lo realmente interesante de PyDeck es que, al igual que plotly, geemap, folium y otras herramientas de mapeo interactivo, permite a los usuarios cambiar el mapa base, lo que significa que puedes diseñar un mapa de acuerdo con el contexto de tu proyecto:

r = pdk.Deck(layers=[layer],
initial_view_state=view_state,
map_style=pdk.map_styles.LIGHT, # ‘light’, ‘dark’, ‘road’, ‘satellite’, ‘dark_no_labels’, and ‘light_no_labels
)
Imagen por autor.

La siguiente característica, que encuentro muy útil, es cambiar la descripción de los datos interactivos. Al colocar un cursor en una columna/hexágono/punto, puede obtener metainformación, pero a menudo parece un poco absurdo. Pero en PyDeck puedes superarlo fácilmente:

Imagen por autor.
r = pdk.Deck(layers=[layer],
initial_view_state=view_state,
"html": "<b>Number of crimes:</b> {elevationValue}",
"style": {
"backgroundColor": "yellow",
"color": "black"
}
},
)
Imagen por autor.

Finalmente, la característica más sorprendente de esta biblioteca es su capacidad de cambiar el ángulo de visión en el mapa simplemente haciendo clic con el botón derecho del mouse:

from ipywidgets import HTML

text = HTML(value='Move the viewpoint')

def filter_by_bbox(row, west_lng, east_lng, north_lat, south_lat):
return west_lng < row['lng'] < east_lng and south_lat < row['lat'] < north_lat

def filter_by_viewport(widget_instance, payload):
west_lng, north_lat = payload['data']['nw']
east_lng, south_lat = payload['data']['se']
filtered_df = df[df.apply(lambda row: filter_by_bbox(row, west_lng, east_lng, north_lat, south_lat), axis=1)]

r.deck_widget.on_click(filter_by_viewport)

Imagen por autor.

Definitivamente disfruto PyDeck y planeamos profundizar en el cubierta.jl estructura. Con una sintaxis súper sencilla e intuitiva, permite a los usuarios crear visualizaciones impresionantes con eficiencia energética. Python limita mucho las capacidades de este paquete, así que consulte su galeríaes alucinante, especialmente su función experimental GlobalView…

¡Esperamos que este artículo haya sido informativo y revelador para usted!

=============================================

Todas mis publicaciones en Medium son gratuitas y de acceso abierto, ¡por eso te agradecería mucho que me siguieras aquí!

PD: Soy un apasionado de la ciencia (geo)datos, el aprendizaje automático/inteligencia artificial y el cambio climático. Entonces, si quieres trabajar juntos en algún proyecto, por favor contáctame en LinkedIn.

🛰️Sigue para más🛰️