Geocodificar con OSM y Google. Mapear con Leaflet - COVID-19 Mapa de la Red Solidaria en José C. Paz para personas Mayores

Hola Rers!

Durante el 2020 y en el contexto de la Pandemia por COVID-19, tuve la oportunidad de participar en el diseño de una estrategia de acompañamiento a adultos/as mayores y personas de alto riesgo en la localidad de José C. Paz, Provincia de Buenos Aires, Argentina. Consistió en “mapear” (o ubicar en un mapa) tanto a estas personas, como a: * voluntarios/as * distintas organizaciones sociales * Comedores populares * Huertas urbanas * Negocios solidarios (exp de alimentos)

y que en términos generales, se encargarían de hacerles llegar alimentos e insumos, especialmente en los momentos más duros del confinamiento o cuarentena obligatoria.

A los efectos de esta publicación, cuya finalidad es mostrar el proceso de geocodificación y cartografiado, omitiré todo dato personal y de contacto.

Geocodificar es…

…el proceso de asignar coordenadas (ergo, localización espacial) a uno o varios atributos alfanuméricos que constituyen una dirección. Se requiere de 4 argumentos fundamentales para ello: Calle, Altura o número, Localidad y Estado o Provincia. Otros atributos mejoran la calidad y la precisión: Barrio, País, Código Postal, Números de teléfono, etc.

Para geocodificar hacemos uso de APIs y servicios geográficos, en este caso, veremos cómo usar el de Google y el de OSM (Nominatim)

Las librerías:

library(ggmap)
library(tmaptools)
library(RCurl)
library(jsonlite)
library(tidyverse)
library(leaflet)
library(DT)
library(rgeos)
library(maps)
library(sf)
library(sp)
library(ggvoronoi)

Servicio de Google

Deberás crear un Token en Google Cloud Platform para usar la API de geocodificación. Ten en cuenta que al registrarte, te regalan 300 USD por un año como crédito, ya que el servicio de Google no es gratuito (aunque hace unos años lo era). Los costos dependerán de la cantidad de consultas que hagas. (aprox 5 USD cada 1.000 direcciones)

Para acceder a las funcionalidades de la API usaremos la libería ggmap

register_google(api_key, write = TRUE) # registramos una variable api_key que contendrá tu token: api_key <- AizIwahEsja82_288401sjjauwqa (esto es un ejemplo)

Preparando los datos

#Organizaciones sociales
oss <- read.csv("ORGANIZACIONES_SOCIALES.csv", sep= ";", encoding = "UTF-8")
DT::datatable(oss)

# Transformar la tabla y construir un solo campo de dirección
oss_dir <- oss %>% 
  mutate(Dir_completa = paste(Calle, Altura, ",", Localidad, ",", Provincia, ", Argentina")) #Completo con ARGENTINA como país para mejorar la precisión
head(oss_dir[11])

Dir_completa
ALSINA 5141 , JOSE C. PAZ , BUENOS AIRES , Argentina
ARREGUI 338 , JOSE C. PAZ , BUENOS AIRES , Argentina
SANTIAGO DE LINIERS 4158 , JOSE C. PAZ , BUENOS AIRES , Argentina
SUECIA 4400 , JOSE C. PAZ , BUENOS AIRES , Argentina
TRES DE FEBRERO 3500 , JOSE C. PAZ , BUENOS AIRES , Argentina
OTAWA 4599 , JOSE C. PAZ , BUENOS AIRES , Argentina

Asignar coordenadas

oss_dir <- oss_dir %>% 
  mutate_geocode(location = Dir_completa) #agregamos una nueva columna donde estarán las coordenadas

Resultados:

summary(oss_dir$lon) # no hay datos vacíos en LON, lo que indica que se encontraron las 31 direcciones

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  -58.81  -58.80  -58.79  -58.78  -58.76  -58.73

DT::datatable(head(oss_dir[10:13]))

# lon: X o longitud, mientras que lat: Y o latitud

De esta manera habremos obtenido las coordenadas a través de la API de Google.

---

El servicio de OpenStreetMap(Nominatim)

El servicio de OSM no tiene costo, y a diferencia de Google, no necesitarás una clave/token. OSM es una plataforma colaborativa, y Nominatim es el servicio de geolocalización. Deberás registrarte con un mail unicamente y usar el servicio “a consciencia”, esto quiere decir que si te pasas con la cantidad de consultas, te restrigen (normalmente por 24 horas). Mirá la Documentación

Para acceder a las funcionalidades de la API usaremos la libería tmaptools , que tiene una amplia funcionalidad, entre ellas, conectarse a los servicios de OSM. Para más info: Cran/tmaptools

Preparando los datos

# En esta oportunidad, geocodifcaremos las Huertas Urbanas

huertas <- read.csv("Huertas.csv", sep= ";", encoding = "UTF-8")
DT::datatable(head(huertas))

# Transformar la tabla y construir un solo campo de dirección
huertas_dir <- huertas %>% 
  mutate(Dir_completa = paste(Calle, Altura, ",", Localidad, ",", Provincia, ", Argentina"))
head(huertas_dir[12])

Dir_completa
SAN LORENZO 3814 , JOSE C. PAZ , BUENOS AIRES , Argentina
ANTARTIDA ARGENTINA 4990 , JOSE C. PAZ , BUENOS AIRES , Argentina
PEDRO AGUSTONI 1230 , PILAR , BUENOS AIRES , Argentina
PEDRO DE MENDOZA 3736 , JOSE C. PAZ , BUENOS AIRES , Argentina
MATHEO BOOTZ 1609 , JOSE C. PAZ , BUENOS AIRES , Argentina
ANTÁRTIDA ARGENTINA 4235 , JOSE C. PAZ , BUENOS AIRES , Argentina

Asignar coordendas

# El código contempla varios argumentos a saber:

huertas_dir <- huertas_dir %>% 
  mutate(geo=geocode_OSM(  #me devuelve las coordenadas (geocode_OSM) en un campo adicional a los originales
    huertas_dir$Dir_completa,  
    return.first.only = TRUE,  # me devuelve solo el primer resultado encontrado, FALSE si quieres lo contrario
    keep.unfound = TRUE,  # que mantenga aquellas direcciones que No pudo geolocalizar, con un null en el campo "geo"
    details = FALSE, # que no devuelva detalles, id de OSM u objetos cercanos, TRUE si quieres lo contrario
    as.data.frame = NA,
    as.sf = FALSE, # que lo deje como una tabla y no lo convierta en un spatial frame
    geometry = "point", # que la geometría sea de tipo PUNTO (es decir, un par de coordenadas)
    server = "https://nominatim.openstreetmap.org" # el servicio que deseo usar: Nominatim
  ))

Resultados:

summary(huertas_dir$geo) # de los 147 registros, se encontraron 107 direcciones, NA's = 40

##     query                lat              lon            lat_min      
##  Length:147         Min.   :-40.80   Min.   :-62.97   Min.   :-40.80  
##  Class :character   1st Qu.:-34.54   1st Qu.:-58.79   1st Qu.:-34.54  
##  Mode  :character   Median :-34.52   Median :-58.76   Median :-34.52  
##                     Mean   :-34.58   Mean   :-58.79   Mean   :-34.59  
##                     3rd Qu.:-34.50   3rd Qu.:-58.74   3rd Qu.:-34.50  
##                     Max.   :-34.45   Max.   :-58.41   Max.   :-34.45  
##                     NA's   :40       NA's   :40       NA's   :40      
##     lat_max          lon_min          lon_max      
##  Min.   :-40.80   Min.   :-62.97   Min.   :-62.97  
##  1st Qu.:-34.54   1st Qu.:-58.79   1st Qu.:-58.78  
##  Median :-34.52   Median :-58.76   Median :-58.76  
##  Mean   :-34.58   Mean   :-58.79   Mean   :-58.79  
##  3rd Qu.:-34.50   3rd Qu.:-58.74   3rd Qu.:-58.73  
##  Max.   :-34.45   Max.   :-58.41   Max.   :-58.41  
##  NA's   :40       NA's   :40       NA's   :40

head(huertas_dir$geo) 

query	lat	lon	lat_min	lat_max	lon_min	lon_max
SAN LORENZO 3814 , JOSE C. PAZ , BUENOS AIRES , Argentina	-34.51970	-58.73151	-34.52356	-34.51577	-58.73685	-58.72648
ANTARTIDA ARGENTINA 4990 , JOSE C. PAZ , BUENOS AIRES , Argentina	-34.52823	-58.76535	-34.53107	-34.52598	-58.76678	-58.76309
PEDRO AGUSTONI 1230 , PILAR , BUENOS AIRES , Argentina	-34.44894	-58.91098	-34.44895	-34.44881	-58.91114	-58.91083
PEDRO DE MENDOZA 3736 , JOSE C. PAZ , BUENOS AIRES , Argentina	-34.51506	-58.72724	-34.51605	-34.51396	-58.72906	-58.72514
MATHEO BOOTZ 1609 , JOSE C. PAZ , BUENOS AIRES , Argentina	NA	NA	NA	NA	NA	NA
ANTÁRTIDA ARGENTINA 4235 , JOSE C. PAZ , BUENOS AIRES , Argentina	-34.52823	-58.76535	-34.53107	-34.52598	-58.76678	-58.76309

# Base de huertas final, con coordenadas
huertas_dir <- huertas_dir %>% 
  mutate(lon=geo$lon, lat = geo$lat) %>% # el campo geo se arma como lista, deberás seleccionar una lon y una lat
  select(X.U.FEFF.ID, Exp_mantenimiento, Caracteristicas_Generales, Dir_completa, lon, lat)

DT::datatable(head(huertas_dir))

## Resultado en el mapa

Para visualizar los dos resultados (huertas y organizaciones sociales), utilizaré la libreria leaflet

# Exploremos las huertas primero 

mapa <- leaflet(huertas_dir) %>% 
    addTiles()  %>%  
    addProviderTiles("CartoDB") %>% # usar el mapa base de CARTO
  addMarkers(lng = ~lon,    # los marcadores serán las huertas 
               lat = ~lat,
               popup = ~Caracteristicas_Generales) %>%  # Esta información será la que aparezca en el popup de cada marcador
  setView(lng = -58.765954,
            lat = -34.515433, # ubicar un punto en el área de estudio: José C. Paz
            zoom = 12) %>% # de entrada, con zoom a nivel de partidos / localidades
    addMiniMap( toggleDisplay = TRUE)# display de un mapa chico para navegar

## Warning in validateCoords(lng, lat, funcName): Data contains 40 rows with either
## missing or invalid lat/lon values and will be ignored

mapa

# Notar que las 40 no geocodifcadas se omiten en la visualización. 
# Mejorar el número de geocodificación normalizando datos de calle, altura y localidad.

De la misma forma de las dos tablas anteriores, geocodificamos los merenderos, los negocios solidarios y a las personas (Voultarixs y adultxs mayores); información que verán a continuación, previamente anonimizada.

comercios_csv <- "https://docs.google.com/spreadsheets/d/e/2PACX-1vQOsaKBvAKvRIuszu77w7bq6o6xSwh6kBLWa-mII_vjr1ClQDpMsMsA7y0SUrwhWbQUm7hnNJF0hRLo/pub?output=csv"

merenderos_csv <- "https://docs.google.com/spreadsheets/d/e/2PACX-1vQz7A_ekC0DtAbvmzi4nlE9GZk3tDQUUNArDD_gDIkEkPwVYeaBhGt75-VcTJ6pq5AQrGpH1nYMdCTF/pub?output=csv"

voluntarios <- read.csv2("voluntarios_josecpaz.csv", dec = ".")

Adultos <- read.csv2("adultos_josecpaz.csv", dec = ".")

Mapa de la red solidaria: Comercios, Huertas, Merenderos y Organizaciones sociales

# con estos Iconos identificaré los merenderos y comercios, para más información visitar: https://www.glyphicons.com/

icon_merendero <- makeAwesomeIcon(icon = "cutlery",markerColor="lightblue", iconColor = "white", library="glyphicon" )
icon_comercio <- makeAwesomeIcon(icon = "shopping-cart", markerColor="red", iconColor = "white", library="fa", spin = F, markerOptions(fontSize=9))



# importando datos de comercios y merenderos

comercios <- read.csv(comercios_csv)
merenderos <- read.csv(merenderos_csv)

# Mapa con las organizaciones, comederos y huertas

mapa3 <- leaflet() %>% 
  addTiles()  %>%  
  addProviderTiles("CartoDB") %>% 
  addCircleMarkers(data= huertas_dir,
             lng = ~lon,    
             lat = ~lat,
             radius = 3,
             color = "#2ca25f",
             popup = ~Caracteristicas_Generales,
             group = "Huertas Urbanas") %>%
  addAwesomeMarkers(data= comercios,
             lng= ~LON,    
             lat = ~LAT,
             icon= icon_comercio,
             popup = ~Nombre.del.Comercio,
             group = "Comercios Solidarios") %>% 
  addAwesomeMarkers(data= merenderos,
             lng= ~longitude,    
             lat = ~latitude,
             icon= icon_merendero,
             popup = ~nombreinst,
             clusterOptions = markerClusterOptions(),
             group = "Merenderos y Comedores") %>% # como son muchos merenderos, decidí mostrarlos en clusters
  addMarkers(data = oss_dir,
             lng = ~lon,    
             lat = ~lat,
             popup = ~Nombre_de_Organizacion,
             group= "Organizaciones Sociales") %>% 
  setView(lng = -58.765954,
            lat = -34.515433, 
            zoom = 12) %>% 
  addMiniMap( toggleDisplay = TRUE) %>% 
  addLegend(values = 1, group = "Huertas Urbanas", position = "bottomleft", labels = "Huertas Urbanas", colors = "#2ca25f") %>% 
  addLegend(values = 2, group = "Merenderos y Comedores", position = "bottomleft", labels = "Merenderos/Comedores", colors = "lightblue") %>% 
  addLegend(values = 3, group = "Organizaciones Sociales", position = "bottomleft", labels = "Org. Sociales", colors = "#3182bd") %>% 
  addLegend(values = 4, group = "Comercios Solidarios", position = "bottomleft", labels = "Comercios Solidarios", colors = "red") %>%
  addLayersControl(overlayGroups = c("Huertas Urbanas", "Merenderos y Comedores", "Organizaciones Sociales", "Comercios Solidarios"),
                   options = layersControlOptions(collapsed = FALSE))

## Warning in validateCoords(lng, lat, funcName): Data contains 40 rows with either
## missing or invalid lat/lon values and will be ignored

mapa3

### Incorporando a las Personas que harán realidad el proyecto

mapa4 <- leaflet() %>% 
  addTiles()  %>%  
  addProviderTiles("CartoDB") %>% 
  addCircleMarkers(data= Adultos,
             lng = ~LON,    
             lat = ~LAT,
             radius = 1,
             color = "orange",
             group = "Adultos") %>%
  addCircleMarkers(data= voluntarios,
             lng = ~LON,    
             lat = ~LAT,
             radius = 1,
             color = "black",
             group = "Voluntarios") %>% 
  addCircles(data= voluntarios, # área de influencia de voluntarios a aprox 10 cuadras a la redonda. 
             lng = ~LON, 
             lat = ~LAT, 
             weight = 1,
             radius = ~1200, 
             fillColor = "gray",
             group = "buffer") %>% 
  setView(lng = -58.765954,
            lat = -34.515433, 
            zoom = 12) %>% 
  addMiniMap( toggleDisplay = TRUE) %>% 
  addLegend(values = 1, group = "Adultos", position = "bottomleft", labels = "Adultxs Mayores", colors = "orange") %>% 
  addLegend(values = 2, group = "Voluntarios", position = "bottomleft", labels = "Voluntarixs", colors = "darkblue") %>%
  addLegend(values = 3, group = "buffer", position = "bottomleft", labels = "Buffer", colors = "gray") %>% 
  addLayersControl(overlayGroups = c("Adultos", "Voluntarios", "buffer"),
                   options = layersControlOptions(collapsed = FALSE))



mapa4

Con un buffer (o área de influencia) establecido a una distancia promedio de 10 cuadras a la redonda, se puede facilmente la zona de cobertura, así como también, las personas que quedaron lejos de la red de voluntarios/redes de apoyo y comercios; Esto activaría el relevamiento de más voluntarios en zonas ya identificadas mediante inteligencia espacial, por ejemplo: Muñiz, Mayor del Pino, etc.

Patrón de distribución con Polígonos de Thiessen

También conocidos como Polígonos de Voronoi, se trata de un método de interpolación que permite crear áreas de influencia a través de un grupo de puntos y la distancia entre ellos (vecinos), donde el límite entre un área y otra será la intersección resultante por mediatriz entre dos puntos. El resultado, es una red de polígonos equidistantes a los puntos vecinos, de mucha utilidad para entender mejor las áreas de influencia especialmente en estudios cualitativos. Para crearlos, usé la libreria ggvoronoi

# Generamos los polígonos basados en Voluntarios
?voronoi_polygon()

## starting httpd help server ... done

voronoi_voluntarios <- voronoi_polygon(voluntarios, x="LON", y="LAT")

# Mapa ...

mapa5 <- addPolygons(mapa4, 
                     data= voronoi_voluntarios,
                     weight = 1,
                     fillColor = "lightgray",
                     group = "voronoi",
                     highlightOptions = highlightOptions(color = "red", weight = 1,
                                                  bringToFront = TRUE)) %>%
  addLabelOnlyMarkers(data= voluntarios,
                      label = ~ï..ID,
                      lng = ~LON, 
                      lat = ~LAT) %>% 
  addLegend(values = 4, group = "voronoi", position = "bottomleft", labels = "voronoi", colors = "lightgray") %>% 
  addLayersControl(overlayGroups = c("Adultos", "Voluntarios", "buffer", "voronoi"),
                   options = layersControlOptions(collapsed = FALSE))
mapa5

Sugerencia: apagar la capa de Buffers para observar mejor la distribución en los polígonos, o: la densidad de “Adultos/as” en cada área de influencia-voronoi de los voluntarios y voluntarias.

Este tipo de operaciones son de gran utilidad para hacer inteligencia espcacial, pudiendo entender gráfica (y cartográficamente) los elementos en el territorio. Cada voluntarix puede disponer de una red de Personas Mayores a quien apoyar en el momento de mayor dificultad, teniendo también a disposición, la red de Negocios Solidarios, Huertas, Establecimientos y organizaciones de apoyo.

Mapa Web Final

mapa_final <- leaflet() %>% 
  addTiles()  %>%  
  addProviderTiles("CartoDB") %>% 
  addCircleMarkers(data= huertas_dir,
             lng = ~lon,    
             lat = ~lat,
             radius = 3,
             color = "#2ca25f",
             popup = ~Caracteristicas_Generales,
             group = "Huertas Urbanas") %>%
  addAwesomeMarkers(data= comercios,
             lng= ~LON,    
             lat = ~LAT,
             icon= icon_comercio,
             popup = ~Nombre.del.Comercio,
             group = "Comercios Solidarios") %>% 
  addAwesomeMarkers(data= merenderos,
             lng= ~longitude,    
             lat = ~latitude,
             icon= icon_merendero,
             popup = ~nombreinst,
             clusterOptions = markerClusterOptions(),
             group = "Merenderos y Comedores") %>%
  addMarkers(data = oss_dir,
             lng = ~lon,    
             lat = ~lat,
             popup = ~Nombre_de_Organizacion,
             group= "Organizaciones Sociales") %>%
  addCircleMarkers(data= Adultos,
             lng = ~LON,    
             lat = ~LAT,
             radius = 1,
             color = "orange",
             group = "Adultos") %>%
  addCircleMarkers(data= voluntarios,
             lng = ~LON,    
             lat = ~LAT,
             radius = 1,
             color = "black",
             group = "Voluntarios") %>% 
  addCircles(data= voluntarios, # área de influencia de voluntarios a aprox 10 cuadras a la redonda. 
             lng = ~LON, 
             lat = ~LAT, 
             weight = 1,
             radius = ~1200, 
             fillColor = "gray",
             group = "buffer") %>% 
  addPolygons(data= voronoi_voluntarios,
              weight = 1,
              fillColor = "lightgray",
              group = "voronoi",
              highlightOptions = highlightOptions(color = "red", weight = 1,
                                                  bringToFront = TRUE)) %>%
  setView(lng = -58.765954,
            lat = -34.515433, 
            zoom = 12) %>% 
  addMiniMap( toggleDisplay = TRUE) %>% 
  addLegend(values = 1, group = "Huertas Urbanas", position = "bottomleft", labels = "Huertas Urbanas", colors = "#2ca25f") %>% 
  addLegend(values = 2, group = "Merenderos y Comedores", position = "bottomleft", labels = "Merenderos/Comedores", colors = "lightblue") %>% 
  addLegend(values = 3, group = "Organizaciones Sociales", position = "bottomleft", labels = "Org. Sociales", colors = "#3182bd") %>% 
  addLegend(values = 4, group = "Comercios Solidarios", position = "bottomleft", labels = "Comercios Solidarios", colors = "red") %>%
  addLegend(values = 5, group = "Adultos", position = "bottomleft", labels = "Adultxs Mayores", colors = "orange") %>% 
  addLegend(values = 6, group = "Voluntarios", position = "bottomleft", labels = "Voluntarixs", colors = "darkblue") %>%
  addLegend(values = 7, group = "buffer", position = "bottomleft", labels = "Buffer", colors = "gray") %>% 
  addLegend(values = 8, group = "voronoi", position = "bottomleft", labels = "voronoi", colors = "lightgray") %>% 
  addLayersControl(overlayGroups = c("Huertas Urbanas", "Merenderos y Comedores", "Organizaciones Sociales", "Comercios Solidarios", "Adultos", "Voluntarios", "buffer", "voronoi"),
                   options = layersControlOptions(collapsed = FALSE))

## Warning in validateCoords(lng, lat, funcName): Data contains 40 rows with either
## missing or invalid lat/lon values and will be ignored

mapa_final

A los efectos de esta publicación, vuelco en el mapa todas las capas generadas, sin embargo, considera que para un mapa final las capas analíticas no hace falta mostrarlas, ya que su finalidad no es de visualización en el producto final, sino que está inmersa en las decisiones tomadas o estrategias elaboradas para el objetivo del proyecto.

Próximamente…

Publicaré un modelo analítico de redes: conectar puntos entre sí mediante una red urbana de calles 😉

---

Encuéntrame en LinkedIn Daniela Rattia - Geógrafa., Especialista en Datos y GisAnalytics