SDSC23: Usando SQL para escalar tu analítica espacial

Prerequisitos

Una cuenta trial de CARTO. Todos los datos usados en este taller están disponibles en abierto dentro de la plataforma.

Vamos a utilizar CARTO en este taller y el motor de base de datos que lo respalda es BigQuery, así que el código SQL hace uso de las funciones y eventuales peculiaridades de dicho motor y de las extensiones de CARTO.

1. Algunos conceptos básicos

Vamos a comenzar explorando datos de puntos de interés (POIs) de OSM. Puedes copiar y pegar este código en la consola de CARTO

select
    *
from
    `cartobq.docs.sds_bootcamp_nyc_osm_points`

Ahora, vamos a filtrar esos POIs y nos vamos a quedar conlos que estén a una distancia específica de nuestro foco de estudio

WITH study_area AS (
    SELECT
    ST_BUFFER(
        ST_GEOGPOINT(-73.97421, 40.75976),
        5000
    ) AS geom
)
SELECT 
    osm.*
FROM
    `cartobq.docs.sds_bootcamp_nyc_osm_points` AS osm,
    study_area
WHERE 
    ST_INTERSECTS(osm.geom, study_area.geom)

Veamos cuales de ellos son restaurantes

WITH study_area AS (
    SELECT ST_BUFFER(ST_GEOGPOINT(-73.97421, 40.75976), 5000) AS geom
)
SELECT 
    osm.*
FROM
    `cartobq.docs.sds_bootcamp_nyc_osm_points` AS osm,
    study_area
WHERE 
    ST_INTERSECTS(osm.geom, study_area.geom)
    AND poi_type IN ('fast_food', 'restaurant')

Y si quisiéramos saber la densidad de POIs por cada "grupo de bloques censales"? Es interesante ver que el joinse hace basado en una condicion espacial

SELECT
    block.geo_id,
    ANY_VALUE(block.geom) AS geom,
    COUNT(osm.osm_id) AS num_osm_points,
    ANY_VALUE(ST_AREA(block.geom)) AS area,
    COUNT(osm.osm_id) / ANY_VALUE(ST_AREA(block.geom)) AS points_per_sq_m
FROM
    `cartobq.docs.sds_bootcamp_nyc_census_block_groups` AS block
JOIN
    `cartobq.docs.sds_bootcamp_nyc_osm_points` AS osm
ON
    ST_INTERSECTS(block.geom, osm.geom)
GROUP BY
    block.geo_id;

Cuando queremos estudiar la distribución de una métrica en el espacio continuo, puede ser muy interesante proyectar el dato a un teselado que sea independiente de decisiones administrativas o urbanísticas y que no varíe con el tiempo. Este tipo de teselado nos permiten almacenar, procesar y servir nuestros datos con gran agilidad y compatibilidad. Los más comunes actualmente

Probemos a generalizar los POIs anteriores a un teselado H3 de nivel 9 (0.1KM²) usando la función H3_FROMGEOGPOINT

WITH 
h3_points AS (
    SELECT
    `carto-un`.carto.H3_FROMGEOGPOINT(geom, 9) AS h3
    FROM
    `cartobq.docs.sds_bootcamp_nyc_osm_points`
)
SELECT
    h3,
    COUNT(h3) as num_points
FROM
    h3_points
GROUP BY
    h3;

2. Un caso de uso: analizando los datos de criminalidad de Chicago

Empecemos echando un vistazo a los datos de 2019

SELECT 
    primary_type, 
    COUNT(unique_key) AS ncount
FROM 
    `bigquery-public-data.chicago_crime.crime`
WHERE 
    year = 2019
GROUP BY 
    primary_type
ORDER BY 
    ncount DESC

Quedémonos con los datos de robos en domicilios, que por agilizar ya están recogidos en la tabla cartobq.docs.sdsb_chicago_burglary_crime_sample

SELECT 
    date, 
    description, 
    location_description, 
    ST_GEOGPOINT(longitude, latitude) AS geom
FROM `cartobq.docs.sdsb_chicago_burglary_crime_sample`

Añadiendo un widget categórico podemos ver la distribución por tipos de robo.

2.1. Distribución

Supongamos que la ciudad de Chicago quiere ubicar 5 nuevas comisarías de modo que estén lo más cerca posible a las zonas donde ser producen más robos en domicilio.

La aproximación más sencilla para encontrar las zonas requeridas, sería un análisis por k-means usando la función ST_CLUSTERKMEANS

WITH 
clustered_points AS
(
    SELECT 
        `carto-un`.carto.ST_CLUSTERKMEANS(
            ARRAY_AGG(ST_GEOGPOINT(longitude, latitude)), 
            5
        ) AS cluster_arr
    FROM `cartobq.docs.sdsb_chicago_burglary_crime_sample`
)
SELECT 
    cluster_element.cluster, 
    cluster_element.geom AS geom 
FROM 
    clustered_points,
    UNNEST(cluster_arr) AS cluster_element

Ahora, para estos clusters ya etiquetados, podemos calcular su centroide o la mediana de las ubicaciones con la función ST_CENTERMEDIAN, para definir las ubicaciones de esas potenciales comisarías.

WITH 
clustered_points AS
(
    SELECT 
        `carto-un`.carto.ST_CLUSTERKMEANS(
            ARRAY_AGG(ST_GEOGPOINT(longitude, latitude)),
            5
        ) AS cluster_arr
    FROM `cartobq.docs.sdsb_chicago_burglary_crime_sample`
)
SELECT 
    cluster_element.cluster,
    `carto-un`.carto.ST_CENTERMEDIAN(ST_UNION_AGG(cluster_element.geom)) AS geom
FROM
    clustered_points, 
    UNNEST(cluster_arr) AS cluster_element
GROUP BY 
    cluster_element.cluster

2.2. Concentración

Para otros casos de uso (inmobiliaria, seguros, administración pública, etc.) puede ser interesante encontrar las zonas de alta concentración de robos. Para ello, podemos hacer uso de un algoritmo de clustering basado en densidad como puede ser DBSCAN (ST_CLUSTERDBSCAN)

WITH 
crimes AS (
    SELECT 
        ST_GEOGPOINT(longitude, latitude) AS geom
    FROM `cartobq.docs.sdsb_chicago_burglary_crime_sample`
)
SELECT 
    geom, 
    ST_CLUSTERDBSCAN(geom, 100, 10) OVER () AS cluster_num
FROM crimes

2.3. Zonificación y ubicaciones atípicas

Otro posible análisis partiendo de estos datos es la zonificación de la ciudad en áreas similares en cuanto a número de robos y la identificación de ubicaciones con valores atípicos de robos (outliers) respecto a su contexto espacial. Empezaremos con el dato generalizado a celdas H3 de nivel 8 (H3_FROMLONGLAT)

SELECT 
    h3, 
    COUNT(h3) as n_crimes
FROM (
    SELECT 
        `carto-un`.carto.H3_FROMLONGLAT(longitude, latitude, 8) as h3
    FROM `cartobq.docs.sdsb_chicago_burglary_crime_sample`
)
GROUP BY h3

Vamos a calcular un índice local de asociación espacial (LISA), en concreto Local Moran's I, así que necesitamos rellenar el espacio con celdas H3 aunque su cuenta de robos haya sido cero. Para ello necesitamos identificar el área de estudio (en nuestro caso el límite municipal de Chicago) y muestrearla con celdas H3 del nivel deseado.

WITH 
chicago_polyfilled AS (
    SELECT 
        h3
    FROM UNNEST(
        `carto-un`.carto.H3_POLYFILL(
            (SELECT geom h3rray FROM `cartobq.docs.sdsb_chicago_boundaries`), 
            8
        )
    ) AS h3
),
h3_crime as(
    SELECT 
        `carto-un`.carto.H3_FROMLONGLAT(longitude, latitude, 8) as h3
    FROM `cartobq.docs.sdsb_chicago_burglary_crime_sample`
),
aggregated_crime AS (
    SELECT h3, COUNT(h3) AS n_crimes
    FROM 
        h3_crime
    GROUP BY h3
),
chicago_agg_crime AS (
    SELECT 
        cp.*, 
        IFNULL(aggt.n_crimes, 0) AS n_crimes
    FROM chicago_polyfilled cp
    LEFT JOIN aggregated_crime aggt
    USING(h3)
)
SELECT 
    i.* EXCEPT(index), 
    index AS h3,
    ['HH', 'LL', 'LH', 'HL'][ORDINAL(i.quad)] as clust
FROM UNNEST(
    (`carto-un`.carto.LOCAL_MORANS_I_H3(
        ARRAY(SELECT AS STRUCT h3, CAST(n_crimes AS FLOAT64) AS value FROM chicago_agg_crime),
        1, 
        'uniform', 
        100
    ))
) AS i

AbelVM/readme.md