logo

Objetivo

Realizar mapas con la mayor granularidad posible, a escala intraurbana, de vulnerabilidad social ante la pandemia de COVID-19.

Con ello se busca contribuir en la toma de decisiones por parte de entidades municipales para la planificación y distribución de atención a nivel intraurbano (víveres con bancos de comida, atención médica in situ, etc.).

En esta instancia se utilizarán solo fuentes abiertas con cobertura global, de modo que el procedimiento sea reproducible en cualquier ciudad de Latinoamérica sin requerir fuentes locales específicas. Por supuesto, los datos locales de gran valor y podrán ser incorporados a éste procedimiento para refinar o extender los resultados.

Esta prueba de concepto toma a Quito, Ecuador, como ubicación a analizar.


Mapa online con los resultados: https://vulnerabilidad.netlify.com

Fuentes de datos

Paquetes a utilizar

library(tidyverse) # funciones útiles en general para manipulación y visualización
library(sf) # para datos georeferenciados en formato vectorial
library(osmdata) # para acceso a datos de OpenStreetMap
library(osrm) # para calcular rutas a pie utilizando la grilla de calles local
library(leaflet) # para visualización de datos geográficos
library(nngeo) # para identificar elementos próximos en el espacio  
               # No disponible en CRAN, usar devtools::install_github("michaeldorman/nngeo")

Descarga de datasets

Estimados demográficos

Población general

Descargamos los datos de estimados demográficos para Ecuador producidos por Facebook

url <- "https://data.humdata.org/dataset/58c3ac3f-febd-4222-8969-59c0fe0e7a0d/resource/c05a3c81-a78c-4e6c-ac05-de1316d4ba12/download/population_ecu_2018-10-01.csv.zip"

zipfile <- tempfile()
download.file(url, zipfile)
filename <- unzip(zipfile, list = TRUE)$Name

poblacion <- read_csv(unz(zipfile, filename))

unlink(temp)

El dataset contiene puntos definidos por pares de coordenadas en proyección Mercator, y estimados de población en torno a esa posición en 2015 y 2020, para todo Ecuador:

head(poblacion)

Cada par latitud/longitud representa un área que corresponde a 1 segundo de arco de resolución (aproximadamente 30m x 30m)

Población mayor a 60 años

De la misma fuente descargamos la cantidad estimada de personas de más de 60 años

url <- "https://data.humdata.org/dataset/58c3ac3f-febd-4222-8969-59c0fe0e7a0d/resource/904d8988-d18d-41a5-a7f7-22668204cefe/download/ecu_elderly_60_plus_2019-06-01_csv.zip"

zipfile <- tempfile()
download.file(url, zipfile)
filename <- unzip(zipfile, list = TRUE)$Name

personas_mayores <- read_csv(unz(zipfile, filename))

unlink(temp)

En este caso el dataset contiene el estimado de personas mayores a 60 años viviendo en el área que corresponde a cada punto:

head(personas_mayores)

Límites de la ciudad

El polígono de las fronteras administrativas de la ciudad puede encontrarse en Nominatim: https://nominatim.openstreetmap.org/

Debido a la infinidad de marcos territoriales y legales que existen en el mundo para definir las fronteras de las ciudades, en general existen varias formas de interpretar cuales son los límites de una ciudad.

Nominatim es una base de datos global de nombres propios de lugares. Al realizarse una búsqueda, por ejemplo por “Quito”, se encuentran distintas entidades geográficas con ese nombre. La interfaz de Nominatim permite comparar opciones y verificar los límites en el mapa, identificando la entidad buscada.

En este caso, buscando por “Quito, Ecuador”, el resultado deseado es el quinto resultado: Quito, Pichincha, Ecuador (County), que coincide con los límites visibles en el SIG de la secretaría de Territorio de Quito

Identificado el nombre, descargamos los límites.

ciudad <- "Quito, Ecuador"
# La posición en que aparece entre los resultados de Nominatim 
posicion <- 5

url <- URLencode(paste0("https://nominatim.openstreetmap.org/search.php?q=",
                        ciudad, "&polygon_geojson=1&format=geojson"))

destfile = paste0(tempdir(), "/limits.json")

download.file(url, destfile)

limites <- st_read(destfile)[posicion,]
## Reading layer `limitesquito' from data source `/home/havb/Downloads/limitesquito.geojson' using driver `GeoJSON'
## Simple feature collection with 1 feature and 9 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: -78.94768 ymin: -0.5956779 xmax: -78.1649 ymax: 0.2562735
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
leaflet(limites) %>%
    addProviderTiles(providers$CartoDB.Positron) %>% 
    addPolygons() 

Puntos de Interés

Para obtener datos georeferenciados locales realizaremos consultas a Overpass (http://overpass-api.de/), una interfaz que permite extraer información de la base de datos global de OpenStreetMap.

Overpass requiere que se especifique una “bounding box”, es decir las coordenadas de un rectángulo que abarque la zona de interés (en la práctica, los valores máximos y mínimos de latitud y longitud).

La generamos:

bbox <- matrix(st_bbox(limites), nrow = 2, dimnames = list(c("x", "y"), c("min", "max")))

bbox
##           min         max
## x -78.9476819 -78.1649035
## y  -0.5956779   0.2562735

Para descargar entidades georeferenciadas se requiere conocer las palabras clave con las que se identifican los registros en la base de OSM. Existe gran detalle para el tipo de datos georeferenciados disponibles: áreas de parques públicos, posición de oficinas de correo o cajeros automáticos, vías de ferrocarril, etc. La nomenclatura se puede consultar en https://wiki.openstreetmap.org/wiki/Map_Features

En este caso vamos a solicitar todas las vías de circulación (calles, avenidas, autopistas, etc) de la ciudad. En la base de datos de OSM ciertas entidades agrupadas bajo la categoría “Healthcare” resultan de inmediato interés:

key value description
amenity clinic A medium-sized medical facility or health centre.
amenity hospital A hospital providing in-patient medical treatment
amenity pharmacy Pharmacy: a shop where a pharmacist sells medications

Para esta prueba de concepto descargaremos clínicas y hospitales.

clinicas <-  opq(bbox) %>% 
  add_osm_feature(key = "amenity", value = "clinic") %>% 
  osmdata_sf() 

hospitales <-  opq(bbox) %>% 
  add_osm_feature(key = "amenity", value = "hospital") %>% 
  osmdata_sf() 

Dependiendo del nivel de detalle con el que hayan sido registrados, la mayoría de las clínicas y hospitales están representados por puntos, pero en algunos casos por los polígonos de su planta. Para nuestros fines sólo necesitamos los puntos. Convertiremos los polígonos en puntos, tomando su centroide.

todo_a_puntos <- function(osm_query_sf) {
  

  # extraemos los elementos de la lista que refieren a entidades OSM. Son dataframes con nombre
  # "osm_points", "osm_lines", "osm_polygons", "osm_multilines", "osm_multipolygons"
  puntos <- osm_query_sf[startsWith(names(osm_query_sf), "osm_")] %>% 
    # retiramos los que estan vacíos
    compact() %>% 
    # nos quedamos solo con el campo "name" 
    map(~select(., name)) %>% 
    # combinamos en un solo dataframe
    reduce(rbind) %>% 
    #tomamos los centroides
    st_centroid() 
  
  
  # Solo retenemos los puntos con nombre, ya que los que tienen ese campo vacío
  # parecen ser ubicaciones repetidas
  filter(puntos, !is.na(name))
}
  
hospitales <- cbind(todo_a_puntos(hospitales), clase = "hospital")
clinicas <- cbind(todo_a_puntos(clinicas), clase = "clinica")

hospitales_y_clinicas <- rbind(clinicas, hospitales)
## Reading layer `hospitales_y_clinicas_quito' from data source `/home/havb/Downloads/hospitales_y_clinicas_quito.geojson' using driver `GeoJSON'
## Simple feature collection with 457 features and 2 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: -78.89925 ymin: -0.5168509 xmax: -78.20187 ymax: 0.2484693
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs

El resultado:

hospitales_y_clinicas
leaflet(hospitales_y_clinicas) %>%
    addProviderTiles(providers$OpenStreetMap) %>% 
    addMarkers(popup = ~paste(clase, name)) %>% 
    addPolygons(data = limites, fill = NA)

Procesamiento

División de la ciudad en zonas de análisis

Aquí es donde cada ciudad puede decidir su unidad geográfica de análisis. Lo ideal serían unidades estadísticas censales, con la mayor resolución (lo más pequeñas) que se disponga.

EL análisis a nivel de agregación similar al barrio (o distritos, comunas, etc) no es deseable debido a su baja resolución espacial. En ausencia de cartografía censal de alta granularidad, puede realizarse la partición de la superficie de la ciudad en celdas arbitrarias, lo suficientemente pequeñas.

Para este ejemplo particionaremos la superficie de la ciudad en celdas de unos 500m de radio (aproximando una superficie de 0.785 km2)

# Usamos una re-proyección equiareal para mayor precisión al calcular superficies

n_celdas <- limites %>% 
  st_transform(crs = "+proj=laea +ellps=WGS84 +units=m +no_defs") %>%
  st_area() %>% 
  {. / (pi * 500^2)} %>% 
  as.numeric() %>% 
  round()

n_celdas
## [1] 5364

Definiremos entonces 5364 celdas.

celdas <- limites %>%
  st_sample(size = 100000) %>%  # sampleamos al azar un número grande, a gusto
  st_coordinates() %>% # extraemos sendas columnas con long y lat
  kmeans(centers = n_celdas) %>%  # nuestros N grupos de puntos con separación máxima entre si
  .$centers %>% # sólo queremos los centroides
  as_tibble() %>% # convertimos en tibble que es lo que le gusta a sf
  st_as_sf(coords = c("X", "Y"), crs = 4326) %>% # pasamos los centroide a objeto espacial
  st_union() %>% # los combinamos en un sólo objeto multipunto
  st_voronoi() %>% # particionamos en voronoi el espacio donde estan
  st_collection_extract("POLYGON") %>% # del resultado extraemos los polígonos
  st_intersection(limites) %>% # recortamos el cuadrado de los voronoi de acuerdo a los límites
  st_sf %>% # Convertimos en dataframe espacial
  mutate(id = rev(row_number())) # agregamos columna con id
## Reading layer `celdas' from data source `/home/havb/Downloads/celdas.geojson' using driver `GeoJSON'
## Simple feature collection with 5364 features and 1 field
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -78.94768 ymin: -0.5956779 xmax: -78.1649 ymax: 0.2562735
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
ggplot(celdas) +
  geom_sf(aes(fill = NULL)) +
  theme_void()

Estimados demogŕaficos por área

Para el ejemplo usaremos los estimados de población mayor a 60 años, como grupo de riesgo. También podría considerase tanto la población general como los mayores, combinando ambas variables un un índice agregado.

En aras de la performance, primero extraemos del dataset nacional los datos que caen dentro de la bounding box de la ciudad (esta operación es muy rápida). Habiendo limitado así los puntos a clasificar, toma mucho menos tiempo verificar a que celda de la ciudad corresponde cada punto del dataset.

personas_mayores_quito <- personas_mayores %>% 
  filter(between(latitude, bbox["y", "min"], bbox["y", "max"]), 
         between(longitude, bbox["x", "min"], bbox["x", "max"])) %>% 
  st_as_sf(coords = c("longitude", "latitude"), crs = 4326) %>%
  st_join(celdas, left = FALSE)
ggplot(personas_mayores_quito) +
  geom_sf(aes(color = population), size = .0001, alpha = .3) +
  scale_color_viridis_c() +
  theme_minimal() + 
  labs(title = "Población de mayores de 60 años por punto", color = "personas")

Agregamos la cantidad de personas por celda, en base a los punto que caen en cada una.

celdas <- personas_mayores_quito %>% 
  st_set_geometry(NULL) %>% 
  group_by(id) %>% 
  summarise(population = sum(population, na.rm = TRUE)) %>% 
  {left_join(celdas, .)}
## Reading layer `celdas_con_data' from data source `/home/havb/Downloads/celdas_con_data.geojson' using driver `GeoJSON'
## Simple feature collection with 5364 features and 2 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -78.94768 ymin: -0.5956779 xmax: -78.1649 ymax: 0.2562735
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
ggplot(celdas) +
  geom_sf(aes(fill = population), color = NA) +
  scale_fill_viridis_c() +
  theme_minimal() + 
  labs(title = "Población de mayores de 60 años por área", fill = "personas")

Efector de salud más cercano a cada área

Para estimar distancias entre cada celda y su efector de salud más cercano, crearemos una matriz Origen-Destino. Tomamos el centroide de cada celda, y determinamos el efector de salud (hospital o clínica) más cercano a su posición.

Como paso previo, descartamos celdas con población menor a un umbral dado:

umbral <- 10 # Solo celdas con al menos 10 personas

celdas <- celdas %>% 
  filter(population >= umbral)

Estimaremos distancia a pie a través de la grilla de calles local usando Open Source Routing Machine (OSRM), un motor de ruteo de alta performance de uso libre, que permite además adquirir grillas de calles de forma muy simple usando datos de OpenStreetMap.

La instalación del software, adquisición de la grilla de calles y puesta en marcha del servicio de ruteo puede resumirse en 4 pasos:

  1. Instalar Docker

  2. Descargar un extracto de los datos de OpenStreetMap para el país o región donde se estimarán los ruteos

  3. Preparar los datos de ruteo (preprocesamiento)

  4. Iniciar el servidor de ruteo

Detalles para la instalación y puesta en marcha de OSRM

Estas instrucciones fueron testeadas en un sistema corriendo Ubuntu Linux. Deberían funcionar sin modificaciones con cualquier otra versión de Linux.

En MacOS, las instrucciones podrían funcionar sin cambios o con ajustes mínimos.

En Windows, los pasos son los mismos pero los comandos podrían ser distintos.

Instalar Docker

Descargar el software de https://www.docker.com/get-started y seguir las instrucciones

Descargar extracto de datos de OpenStreetMap

Crear una carpeta de trabajo

Por ejemplo,

mkdir -p ~/data/osrm/

cd ~/data/osrm/

Es importante no usar espacios en blanco en los nombres de los directorios, ya que causa problemas con docker. Es decir, evitar directorios tipo “~/data/ruteo con osrm/”

Luego descargams allí la grilla de calles. Podemos usar la información geográfica contenida en la base de OpenStreetMap para obtener la grilla de calles de nuestra área de interés… o de todo el mundo.

Por ejemplo, para el contenido completo de todo Ecuador podemos usar el servicio OpenStreetMap Data Extracts provisto por GeoFabrik.

Descargamos el extracto más reciente para Ecuador con:

wget https://download.geofabrik.de/south-america/ecuador-latest.osm.pbf

Preprocesar los datos de ruteo

OSRM requiere un procesamiento previo de los datos extraidos de OpenStreetMap, construyendo un modelo optimizado de la grilla de calles. Para la optimización se debe elegir un modo de transporte. Las opciones disponibles por default son: car - foot - bicycle

Hacemos la extracción para el modo a pie (foot) con:

docker run -t -v $(pwd):/data osrm/osrm-backend osrm-extract -p /opt/foot.lua /data/ecuador-latest.osm.pbf

Esto inicia una instancia de Docker, corre el comando osrm-extract y luego desactiva la instacia.

La opción -v $(pwd):/data crea un directorio /data dentro del contenedor de Docker, y muestra allí el contenido del directorio desde donde corrimos el comando; es decir, el directorio donde dejamos los datos de OSM. Así es como queda accesible para la instancia en Docker como /data/ecuador-latest.osm.pbf.

Cuando corremos el comando por primera vez tomará unos cuantos minutos, ya que docker necesita descargar antes el container con OSRM. En lo sucesivo no será necesario, ya que tras la descaga guarda una copia local del container.

Quedan dos operaciones pendientes que OSRM necesita antes de poder realizar ruteos.

Dividir el grafo en celdas, ejecutando:

docker run -t -v $(pwd):/data osrm/osrm-backend osrm-partition /data/ecuador-latest.osrm

Obsérvese que ahora usamos ecuador-latest._osrm_ (ya no es “.pbf”), un archivo que quedó en el directorio de trabajo tras ejecutar el paso anterior.

… y por último asignar “peso” a cada celda del grafo con

docker run -t -v $(pwd):/data osrm/osrm-backend osrm-customize /data/ecuador-latest.osrm

Iniciar el servidor de ruteo

Con los datos ya preparados, sólo resta iniciar el servicio de ruteo:

docker run -t -i -p 5000:5000 -v $(pwd):/data osrm/osrm-backend osrm-routed --algorithm mld /data/ecuador-latest.osrm

Con ello nuestro sistema esta listo para ejecutar los siguientes pasos.

A partir de ahora, pasamos a R. El resto del documento muestra código ejecutado en el entorno R.


Disponiendo del servidor de ruteo, identificamos el elemento de la lista Y (sitios de salud) más cercano a cada elemento de la lista X (centroides de celdas), para luego poder calcular la ruta a pie.

# Calculamos el centroide
centroides <- st_centroid(celdas)

# Identificamos el índice (la posición en el dataframe) de hospital/clínica más cercano a cada uno
id_cercanos <- unlist(st_nn(centroides, hospitales_y_clinicas))

Luego creamos una tabla de trayectos, con las coordenadas del centroide de cada celda (los orígenes) junto las coordenadas de su hospital más cercano (los destinos):

# Creamos un dataframe de trayectos
trayectos <- tibble(
  origen_id = celdas$id,
  origen_X = st_coordinates(st_centroid(celdas))[, 1],
  origen_Y = st_coordinates(st_centroid(celdas))[, 2],
  destino_id = hospitales_y_clinicas[id_cercanos, "name"]$name,
  destino_X = st_coordinates(hospitales_y_clinicas[id_cercanos, ])[, 1],
  destino_Y = st_coordinates(hospitales_y_clinicas[id_cercanos, ])[, 2]
)
head(trayectos)

Probamos con el primer trayecto la obtención de la ruta a pie:

# Utilizamos nuestra instancia local de OSRM, en lugar de la opción por defecto que es el servidor público del proyecto OSRM

options(osrm.server = "http://127.0.0.1:5000/")

ruta <- osrmRoute(c(trayectos[1,"origen_id"], trayectos[1,"origen_X"], trayectos[1,"origen_Y"]),
                  c(trayectos[1,"destino_id"], trayectos[1,"destino_X"], trayectos[1,"destino_Y"]),
                  overview = "full", returnclass = "sf")
## Reading layer `ruta_osrm_quito' from data source `/home/havb/Downloads/ruta_osrm_quito.geojson' using driver `GeoJSON'
## Simple feature collection with 1 feature and 4 fields
## geometry type:  LINESTRING
## dimension:      XY
## bbox:           xmin: -78.29649 ymin: -0.10832 xmax: -78.28615 ymax: -0.09761
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs

El resultado:

ruta$duration
## [1] 25.33333
ruta$distance
## [1] 2.111
leaflet(ruta) %>% 
  addProviderTiles(provider = providers$OpenStreetMap) %>% 
  addPolylines(color = "red") %>% 
  addMarkers(data = trayectos[1,2:3], lng = ~origen_X, lat = ~origen_Y,
             label = "origen", labelOptions = labelOptions(noHide = T, direction = "bottom")) %>% 
  addMarkers(data = trayectos[1,4:6], lng = ~destino_X, lat = ~destino_Y, 
             label = ~destino_id, labelOptions = labelOptions(noHide = T, direction = "bottom"))

Ahora estimamos tiempo y distancia del viaje a pie para todos los trayectos, y agregamos a cada celda el tiempo de viaje desde su centro hasta el efector de salud más cercano:

# Función que toma una fila con origen y destino, y devuelve sus identificadores 
# junto a tiempo y distancia del viaje
obtener_ruta <- possibly(function(origen_id, origen_X, origen_Y, 
                                  destino_id, destino_X, destino_Y) {
  
                            ruta <- osrmRoute(src = c(origen_id, origen_X, origen_Y),
                                              dst = c(destino_id, destino_X, destino_Y),
                                              overview = FALSE)
                            
                            tibble(origen_id, destino_id, 
                                   duracion = ruta["duration"], distancia = ruta["distance"]) }, 
                        otherwise = NULL)


# Aplicamos la función a todas las filas del dataframe de trayectos

distancias <- pmap_df(trayectos, obtener_ruta) 


# Agregamos a cada celda su tiempo de viaje y distancia hasta efector de salud más cercano


celdas <- celdas %>% 
  left_join(distancias, by = c("id" = "origen_id"))
## Reading layer `celdas_con_data_distancia' from data source `/home/havb/Downloads/celdas_con_data_distancia.geojson' using driver `GeoJSON'
## Simple feature collection with 1161 features and 5 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -78.93785 ymin: -0.4184543 xmax: -78.27289 ymax: 0.2433535
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
celdas <- celdas %>% 
  mutate(rango = cut(duracion, breaks = c(0, 15, 30, 45, 60, 90, 120, max(duracion)), 
                   labels = c("menos de 15 minutos", "15 a 30 minutos", "30 a 45 minutos", 
                              "45 a 60 minutos", "60 a 90 minutos", "90 a 120 minutos",
                              "más dos horas")))

ggplot(celdas) +
  geom_sf(aes(fill = rango), color = NA) +
  scale_fill_viridis_d(option = "magma", direction = -1) +
  theme_minimal() + 
  labs(title = "Tiempo de viaje a pie hasta el efector de salud más cercano", fill = NULL)

Identificación de zonas de riesgo

Combinando los datos demográficos con los de cercanía a salud, resaltamos las zonas donde se concentra población mayor, y a la vez se verifican las mayores distancias.

ggplot(celdas) +
  geom_point(aes(x = population, y = duracion)) +
  labs(x = "personas mayores a 60 años") +
  theme_minimal()

En base al gráfico de dispersión, consideraremos casos crítico a las celdas con más de 500 residentes mayores que distan más de 1.5 km del centro de salud más cercano

celdas <- celdas %>% 
  mutate(sitio_critico = population > 500 & duracion > 30)
ggplot(celdas) +
  geom_point(aes(x = population, y = duracion)) +
  geom_point(data = filter(celdas, sitio_critico), 
             aes(x = population, y = duracion), color = "red") +
  labs(x = "personas mayores a 60 años") +
  theme_minimal()

En el mapa:

celdas %>% 
  filter(sitio_critico) %>% 
  leaflet() %>% 
  addProviderTiles(providers$OpenStreetMap) %>% 
  addPolygons(color = "purple", popup = ~paste("población > 60 años:", round(population),"</br>",
                                            "tiempo de viaje:", duracion, "minutos"))