DBSCAN: Machine Learning para detectar centros de actividad urbana

 In Blogs

Este post fue publicado originalmente en Bits & Bricks por avazquez@gmail.com (H. Antonio Vazquez Brust)

El algoritmo DBSCAN

DBSCAN es un algoritmo de machine learning diseñado para detectar en forma automática “clusters”, es decir elementos próximos entre si de acuerdo a sus atributos en varias dimensiones.

A diferencia de otros algoritmos de clustering como KMeans, DBSCAN resulta muy adecuado para buscar patrones de agrupación en el espacio físico. Por ejemplo, en la distribución espacial de actividades humanas.

Como se ilustra debajo, entre varias alternativas DBSCAN es la única cuyos resultados aproximan los de un analista humano que estuviera clasificando puntos aglomerados en un mapa (fuente).

Pongámoslo a trabajar entonces, para identificar los centros comerciales (y de otros tipos de actividad) en la ciudad de Mendoza y sus alrededores. Un ejercicio como este, si se repitiera con cierta frecuencia, permitiría detectar la aparición, desplazamiento o desaparición de distintos centros de actividad especializada en la ciudad: zonas con oferta concentrada de bares, de servicios profesionales, de tiendas de bicicletas o lo que fuere.

Paquetes a utilizar

Los herramientas generales de Tidyverse, siempre útiles para el proceso general de análisis de datos

library(tidyverse)

y cinco paquetes especializados:

  • sf para leer y procesar datos geoespaciales
  • geosphere para calcular distancias sobre la superficie terrestre
  • dbscan para emplear el algoritmo homónimo de machine learning que permite encontrar clusters de elementos próximos
  • ggmap para obtener y visualizar mapas detallados de la ciudad que analizaremos
  • ggrepel para etiquetar en forma prolija nuestros mapas con gran cantidad de puntos
library(sf)
library(geosphere)
library(dbscan)
library(ggmap)
library(ggrepel)

Los datos

Trabajaremos con un dataset de “puntos de interés”, o PoI en la jerga de aplicaciones espaciales (por “Points of Interest”).

mdz_poi <- read_csv('https://bitsandbricks.github.io/data/mendoza_poi.csv')

head(mdz_poi)
## # A tibble: 6 x 5
##   nombre                                        tipo  categoria   lat   lng
##   <chr>                                         <chr> <chr>     <dbl> <dbl>
## 1 Banco de la Nación Argentina                  atm   banca     -32.9 -68.8
## 2 Hospital Central BANCO DE LA NACION ARGENTIN… atm   banca     -32.9 -68.8
## 3 Banelco                                       atm   banca     -32.9 -68.8
## 4 Banelco                                       atm   banca     -32.9 -68.8
## 5 Cajero Automatico Banelco                     atm   banca     -32.9 -68.8
## 6 Citibank                                      bank  banca     -32.9 -68.8

Se trata de lugares en el área urbbana de Mendoza, de 93 tipos distintos, agrupados en 12 categorías generales:

count(mdz_poi, categoria, tipo)
## # A tibble: 93 x 3
##    categoria tipo                 n
##    <chr>     <chr>            <int>
##  1 banca     atm                 93
##  2 banca     bank               112
##  3 culto     church             248
##  4 culto     mosque               2
##  5 culto     place_of_worship     4
##  6 culto     synagogue            3
##  7 cultura   art_gallery         14
##  8 cultura   library             18
##  9 cultura   museum              55
## 10 educacion school             890
## # … with 83 more rows

Los datos fueron compildos a mediados del 2017, a partir de consultas a la base de Google Places. Google Places es la base de datos donde Google registra información sobre puntos de interés (comercios, prestadores de servicios, terminales de transporte y un largo etcétera) en cualquier parte del mundo. La forma habitual de acceder a partes de esta información es mediante la Google Maps, desde un navegador o una app de smartphone, pero también es posible consultar la base de forma automatizada, la forma en la cual se relevaron las datos que estamos usando.

Implementando DBSCAN

Para practicar, empecemor por detectar clusters de bares.

I. Realizar una matriz de distancias

El primer paso para aplicar el algorimto es obtener una matriz con la distancia que media entre todos los puntos a analiar. Dado que trabajamos con lugares situados sobre la superficie terrestre, debemos calcular distancias geográficas; es decir, tomar en cuenta la curvatura de la superficie. Para ellos podemos usar la función distm() del paquete geosphere.

Veamos como lucen nuestros bares en el mapa:

mdz_bares <- mdz_poi %>% 
    filter(tipo == "bar")

# Definir la "caja" de coordenadas donde entran los datos, para luego pedir el mapa
bbox <- c(min(mdz_poi$lng),
          min(mdz_poi$lat),
          max(mdz_poi$lng),
          max(mdz_poi$lat))

mendoza <- get_stamenmap(bbox = bbox, maptype = "toner-background", zoom = 13)


ggmap(mendoza) +
    geom_point(data = mdz_bares, aes(x = lng, y = lat), color = "orange", alpha = .5)

Ahora creamos la matriz con la distancia entre cada par de puntos (es decir, medimos la distancia de cada bar contra todos los demás), el insumo que requiere DBSCAN.

distancias <- mdz_bares %>%
    select(lng, lat) %>% # extraemos las columnas de longitud y latitud
    distm(fun = distGeo) %>% # Calculamos las distancias de acuerdo a la curvatura de la Tierra
    as.dist() # convertimos en una matriz de distancias (el tipo de objeto que DBSCAN espera)

II. Definir parámetros

Con la matriz de distancias a mano, es hora de decisiones. Tenemos que fijar dos parámetros, con el fin de indentificar zonas densas, medidas por la cantidad de objetos cercanos a cualquier punto dado:

  • epsilon es el radio de “vecindad” en torno a un punto.
  • minPts es la cantidad mínima de vecinos dentro del radio epsilon que alcanzan para definir a ese punto como miembro de un cluster

((fuente)[http://www.sthda.com/english/wiki/print.php?id=246#concepts-of-density-based-clustering])

En (a) vemos a x, un punto “núcleo”(core). Son puntos “núcleo” aquellos que dentro de su radio epsilon tienen la cantidad mínima de vecinos definida como minPts. En (b) vemos tambiéb a y y a z. y es un punto “borde” porque su cantiadad de vecinos es menor a minPts, pero de todas formas se encuentra en la vecindad de un punto núcleo. Los puntos núcleo y brode que comparten vecindad son parte del mismo clúster. z es considerado “ruido” ya que no es núcleo ni borde, y por tanto no pertenece a un cluster.

Ahora bien, ¿cómo decidimos que cantidad es buena para minPts, y que distancia para epsilon? No hay receta infalible. Depende de cada caso, y se pueden encontrar por prueba y error. Para nuestro problema, decidamoslo así: vamos a considerar como cluster aquellas zonas donde, en un radio de una cuadra a la redonda de un bar dado (1 cuadra = 128 m, le), se encuentren al menos otros 5 bares. Tendremos entonces epsilon = 128 y minPts = 5.

III. Aplicar el algoritmo

Alla vamos:

clusters <- dbscan(distancias, eps = 128, minPts = 5)

clusters
## DBSCAN clustering for 321 objects.
## Parameters: eps = 128, minPts = 5
## The clustering contains 9 cluster(s) and 207 noise points.
## 
##   0   1   2   3   4   5   6   7   8   9 
## 207   8  43   6   6  22   6  12   6   5 
## 
## Available fields: cluster, eps, minPts

Con nuestras reglas, encontramos 9 clusteres en la ciudad. Los bares cuyo valor asignado es 0 son “ruido”, los que no pertenecen a ningún cluster.

IV. Visualizar los resultados

Agreguemos a cada bar el cluster que le corresponde:

mdz_bares <- mdz_bares%>% 
    cbind(cluster = clusters[["cluster"]])

Y con eso a podemos proyectar los resultados en un mapa:

# Visualizamos en un capa los bares sin cluster, y luego los que estan agrupados en otra
ggmap(mendoza) +
    geom_point(data = filter(mdz_bares, cluster == 0), 
               aes(x = lng, y = lat), 
               alpha = .5) +
    geom_point(data = filter(mdz_bares, cluster > 0), 
               aes(x = lng, y = lat, color = factor(cluster)), 
               alpha = .5) +
    labs(title = "Mendoza: clusters de bares",
         color = "cluster")

Automatizando la detección

Vamos a crear una función que realice todos los pasos que vimos antes: realizar la matriz de distancia, fijar los parámetros de DBSCAN, aplicar el algoritmo, y agregar los identificadores de cluster a la data original.

Ya que vamos a lidiar con categorías generales (por ejemplo “servicios” incluye de todo, desde dentistas a funerarias) vamos a ponernos exigentes con los parámetros. Consideraremos que un punto pertenece a un cluster de actividad cuando tiene al menos 25 vecinos en una cuadra a la redonda. Con algunas excepciones: para banca, cultura, culto, educacion, entretenimiento, espacios_verdes, gobierno_y_serv_publicos, y “salud” dejaremos en 5 al umbral, ya que se trata de sitios de mayor “peso”. Seis bares en estrecha proximidad no representan una agrupación notable, pero con seis sedes de la administración pública, o seis museos, ya hablaríamos un centro de actividad.

Si no elegimos un valor par minPts, nuestra función se encargará de verificar el tipo de categoría a la que corresponden los puntos, y fijar el parámetro en forma acorde.

assign_clusters <- function(poi_df, minPts = NA) {
    
    # Para ciertas categorias, como educacion o salud, consideramos que a partir de n = 5 ya se 
    # forma clustering.
    # Para las demas (retail y servicios) el umbral se eleva a 25
    
    if(is.na(minPts)) {
        if(poi_df[1, "categoria"] %in% c("banca", "cultura", "culto", "educacion", "entretenimiento", 
                                         "espacios_verdes", "gobierno_y_serv_publicos", "salud")) {
            minPts <- 5
        } else minPts <- 25
    }
    
    
    eps <- 128 #  metros de distancia máxima entre miembros de un cluster
    
    poi_df[c("lng", "lat")] %>% 
        distm(fun = distHaversine) %>%
        as.dist() %>% 
        dbscan(eps = eps, minPts = minPts) %>% 
        .[["cluster"]] %>% 
        cbind(poi_df, cluster = .)
    
}

Ahora aplicamos la función a nuestros datos.

ATENCION! Con menos de 16 GB de RAM en el sistema, es posible que los recursos no alcancen para procesar el próximo bloque de código. El consumo de memoria de DBSCAN aumenta en forma exponencial (cuadrática) a medida que se incrementan la cantidad de puntos a analizar.

mdz_poi <- mdz_poi %>%
    # separamos el dataset en una lista con un dataframe por cada categoría
    split(mdz_poi$categoria) %>%  
    # Le aplicamos al dataframe de cada categoría la función de clustering
    map_df(assign_clusters) 

Listo!

Producto final: un gráfico facetado por categoría con todos los clusters hallados

Antes de llevar los resultados al mapa, extraigamos para cada cluster su convex hull.

La “envolvente convexa” de un cluster no es otra cosa que un polígono que encierra todos sus puntos. La idea us usarlas para mostrar en el mapa la extensión de las áreas que hemos identificado como centros de actividad.

get_hull<- function(df) {
    
    cbind(df$lng, df$lat) %>% 
        as.matrix() %>%
        st_multipoint() %>% 
        st_convex_hull() %>% 
        st_sfc(crs = 4326) %>% 
        {st_sf(categoria = df$categoria[1], cluster = df$cluster[1], geom = .)}
}

hulls <- function(df) {
    
    df %>%
        split(.$cluster) %>% 
        map(get_hull) %>% 
        reduce(rbind)
    
}

mdz_cluster_hulls <- mdz_poi %>%
    filter(cluster != 0) %>%
    select(lng, lat, categoria, cluster) %>% 
    split(.$categoria) %>% 
    map(hulls) %>% 
    reduce(rbind) 
#Para reducir la cantidad de etiquetas
hulls_por_cat <- mdz_cluster_hulls %>% 
    group_by(categoria) %>% 
    summarise()

# A graficar

ggmap(mendoza) + 
  coord_sf(crs = st_crs(3857)) + # Necesario cuando se combinan ggmap y geom_sf
    geom_sf(data = hulls_por_cat, fill = NA, size = 1.5,
            aes(color = factor(categoria)), inherit.aes = FALSE) +
    geom_point(data = filter(mdz_poi, cluster == 0), size = .1, alpha =.1,
               aes(x = lng, y = lat)) +
    geom_point(data = filter(mdz_poi, cluster != 0), size = .2, alpha =.2,
               aes(x = lng, y = lat, color = factor(categoria))) + 
    ggrepel::geom_label_repel(data = hulls_por_cat,
                              aes(label = categoria, color = factor(categoria), geometry = geom),
                              stat = "sf_coordinates",
                              inherit.aes = FALSE) +
    labs(y = "", x = "",
         title="Mendoza: clusters de actividad identificados",
         color = "categoría") +
    # Eliminamos las etiquetas de latitud y longitud de los ejes 
    theme(axis.text.x = element_blank(),
          axis.text.y = element_blank(),
          axis.ticks = element_blank()) +
    guides(colour = FALSE)

Recommended Posts

Start typing and press Enter to search