library(tidyverse)
library(sf)

Vamos a trabajar con un dataset de inspecciones y fiscalizaciones de obras y comercios. El dataset original se descargó de BADATA

inspecciones <- read.csv("https://github.com/paulavidela/utdt_cienciadedatos/raw/main/data/inspecciones2021_caba.csv", encoding = "latin1")

Exploremos los datos:

head(inspecciones)
##         ticket                       area     direccion_entidad Seccion Manzana
## 1 DGFyC-643955  SGO Fiscalización Crítica          BACACAY 1709      57      20
## 2 DGFyC-645881 SGO Fiscalización Nocturna        GAONA AV. 1537      59     169
## 3 DGFyC-645877 SGO Fiscalización Nocturna   CORRIENTES AV. 4508      17      14
## 4 DGFyC-643004  SGO Fiscalización Crítica MONTES DE OCA AV. 616       8       2
## 5 DGFyC-642602 SGO Fiscalización Nocturna        NICARAGUA 6048      35      63
## 6 DGFyC-644029  SGO Fiscalización Crítica           ESTOMBA 851      49     117
##   Parcela fecha_inspeccion
## 1    024E       20/03/2021
## 2      31       20/03/2021
## 3    001E       20/03/2021
## 4    004j       20/03/2021
## 5       7       21/03/2021
## 6       7       20/03/2021

Como podemos ver, la variable direccion_entidad indica la calle y altura del lugar de la inspección. Dentro del dataset, no se encuentran otros datos para ubicar espacialmente la inspección.

Geocodificando con OSM

Open Street Map. es un proyecto abierto de mapas, información geográfica y datos geoespaciales. Esta plataforma cuenta con conocimiento local y es impulsado por la comunidad, que contribuye y mantiene la información de OSM. Además, cuenta con datos abiertos: pueden utilizarse libremente para cualquier propósito, y por eso muchos sitios web y aplicaciones consumen los datos de OSM.
Nominatim es una API que permite geocodificar direcciones a partir de OSM.


Vamos a utilizar las librerías de tmap y tmaptools

#install.packages("tmap")
library(tmap)
#install.packages("tmaptools")
library(tmaptools)

La función geocode_OSM() realiza una query (consulta) a la API de Nominatim. Vamos a utilizar esta función para localizar la Universidad Torcuato Di Tella. La universidad se encuentra en Figueroa Alcorta 7350, CABA

geocode_OSM("Figueroa Alcorta 7350, CABA")
## $query
## [1] "Figueroa Alcorta 7350, CABA"
## 
## $coords
##         x         y 
## -58.44667 -34.54778 
## 
## $bbox
##      xmin      ymin      xmax      ymax 
## -58.44751 -34.54861 -58.44586 -34.54716

La query nos devuelve 3 elementos:
- la búsqueda que realizamos
- las coordenadas de la ubicación obtenida
- la caja de coordenadas (bounding box) refiere a las coordenadas de un rectángulo que abarca la zona de interés

Como se puede ver, el output generado no es un dataframe. Es una lista con 3 elementos, que a su vez tienen otros elementos.


Ajustando parámetros dentro de la función, podemos obtener el dato en formato sf.

utdt_geo <- geocode_OSM("Figueroa Alcorta 7350, CABA", as.sf =  TRUE, geometry = "point")

¿Cuál es el formato del objeto obtenido?

str(utdt_geo)
## Classes 'sf' and 'data.frame':   1 obs. of  9 variables:
##  $ query  : chr "Figueroa Alcorta 7350, CABA"
##  $ lat    : num -34.5
##  $ lon    : num -58.4
##  $ lat_min: num -34.5
##  $ lat_max: num -34.5
##  $ lon_min: num -58.4
##  $ lon_max: num -58.4
##  $ bbox   :sfc_POLYGON of length 1; first list element: List of 1
##   ..$ : num [1:427, 1:2] -58.4 -58.4 -58.4 -58.4 -58.4 ...
##   ..- attr(*, "class")= chr [1:3] "XY" "POLYGON" "sfg"
##  $ point  :sfc_POINT of length 1; first list element:  'XY' num  -58.4 -34.5
##  - attr(*, "sf_column")= chr "point"
##  - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA NA NA NA NA
##   ..- attr(*, "names")= chr [1:8] "query" "lat" "lon" "lat_min" ...

Es un sf que cuenta con 9 variables. Contiene la información del punto y de la caja de coordenadas.


Podemos graficar el resultado en un mapa utilizando ggplot().

ggplot() + 
  geom_sf(data = utdt_geo)

Sin embargo, al no tener contexto, no logramos entender si se geolocalizó correctamente.

Mapas interactivos

Leaflet es una librería open-source escrita en JavaScript para realizar mapas interactivos.
Por ejemplo, el Instituto Geográfico Nacional de Argentina utiliza para sus mapas componentes de Leaflet. Podemos observar la referencia a Leaflet debajo a la derecha.

Uso de Leaflet Existe un paquete para R que nos permite utilizar las funcionalidades de leaflet.

#install.packages('leaflet')
library(leaflet)
## Warning: package 'leaflet' was built under R version 4.0.5

Vamos a generar un mapa interactivo. En primer lugar, utilizamos la función leaflet() para iniciar el contenido interactivo.

leaflet()

La función leaflet() nos permite iniciar el widget del mapa. Si bien ya tenemos ciertas funcionalidades (botones de zoom), aún no vemos nada. Para eso es necesario agregar imágenes de mapas.
Las imágenes de base se conocen como mosaicos (tiles), y las podemos agregar con la función addTiles()

leaflet() %>%
  addTiles()

¡Excelente! Ya generamos nuestro primer mapa interactivo. Ahora vamos a incorporar la geometría de la ubicación de UTDT. A diferencia de geom_sf que detecta automáticamente el tipo de geometría (puntos, líneas o polígonos), en el caso de leaflet debemos precisar que tipo de geometría es. Las geometrías más utilizadas son:
- puntos: addMarkers() o addCircleMarkers()
- líneas: addPolylines()
- polígonos: addPolygons()
Vamos a utilizar leaflet para visualizar la ubicación del punto de la universidad.

leaflet() %>%
  addTiles() %>% 
  addMarkers(data = utdt_geo)

Podemos ver en el mapa la ubicación de la facultad.

Probemos con otra función de leaflet. En este caso vamos a incorporar un pop-up que indique que se trata de la universidad.

leaflet() %>%
  addTiles() %>% 
  addCircleMarkers(data = utdt_geo, popup = "Universidad Torcuato Di Tella")

Inspecciones de obras y comercios en GCBA

Geolocalizar direcciones lleva tiempo. Por ese motivo, vamos a geolocalizar las direcciones de los 10 lugares donde más inspecciones se realizaron.

insp_agrupadas10 <- inspecciones %>% 
  group_by(direccion_entidad) %>%
  summarise(cant_insp = n()) %>%
  top_n(10)
## Selecting by cant_insp
head(insp_agrupadas10, 10)
## # A tibble: 10 × 2
##    direccion_entidad           cant_insp
##    <chr>                           <int>
##  1 AVELLANEDA AV. 1240                25
##  2 AVENIDA DEL LIBERTADOR 4103        94
##  3 AVENIDA DEL LIBERTADOR 7395        43
##  4 BRANDSEN 805                       26
##  5 CORRIENTES AV. 2322                27
##  6 CORRIENTES AV. 2330                28
##  7 GUARDIA VIEJA 3360                 24
##  8 HUMAHUACA 4310                     23
##  9 HUMBOLDT 450                       24
## 10 PARAGUAY 4979                      58

Para simplificar la búsqueda de dirección en OSM, incorporamos la ciudad y el país.

insp_agrupadas10 <- insp_agrupadas10 %>%
  mutate(direccion = paste0(direccion_entidad, ", CABA, Argentina"))

Vamos a utilizar la función geocode_OSM() para georrefenciar las 10 direcciones

direcciones_localizadas <- geocode_OSM(insp_agrupadas10$direccion, 
                           as.data.frame = TRUE)

Veamos los resultados

head(direcciones_localizadas)
##                                          query       lat       lon   lat_min
## 1         AVELLANEDA AV. 1240, CABA, Argentina -34.61864 -58.44942 -34.61869
## 2 AVENIDA DEL LIBERTADOR 4103, CABA, Argentina -34.56911 -58.42648 -34.56916
## 3 AVENIDA DEL LIBERTADOR 7395, CABA, Argentina -34.54475 -58.45927 -34.54480
## 4                BRANDSEN 805, CABA, Argentina -34.63568 -58.36486 -34.63659
## 5         CORRIENTES AV. 2322, CABA, Argentina -34.60412 -58.40534 -34.60912
## 6         CORRIENTES AV. 2330, CABA, Argentina -34.60412 -58.40534 -34.60912
##     lat_max   lon_min   lon_max
## 1 -34.61859 -58.44947 -58.44937
## 2 -34.56906 -58.42653 -58.42643
## 3 -34.54470 -58.45932 -58.45922
## 4 -34.63464 -58.36584 -58.36378
## 5 -34.59912 -58.41034 -58.40034
## 6 -34.59912 -58.41034 -58.40034

Ahora tenemos las 10 direcciones con sus respectivas coordenadas. Vamos a unir la tabla de direcciones a la tabla de inspecciones agrupadas ‘insp_agrupadas’. Ambas tablas tienen la variable ‘direccion’ en comun. Podemos utilizar distintos tipos de uniones (join). En este caso, vamos a utilizar left_join().

insp_agrupadas10_geo <- insp_agrupadas10 %>%
  left_join(direcciones_localizadas, by = c("direccion" = "query"))

Veamos el resultado

insp_agrupadas10_geo
## # A tibble: 10 × 9
##    direccion_entidad cant_…¹ direc…²   lat   lon lat_min lat_max lon_min lon_max
##    <chr>               <int> <chr>   <dbl> <dbl>   <dbl>   <dbl>   <dbl>   <dbl>
##  1 AVELLANEDA AV. 1…      25 AVELLA… -34.6 -58.4   -34.6   -34.6   -58.4   -58.4
##  2 AVENIDA DEL LIBE…      94 AVENID… -34.6 -58.4   -34.6   -34.6   -58.4   -58.4
##  3 AVENIDA DEL LIBE…      43 AVENID… -34.5 -58.5   -34.5   -34.5   -58.5   -58.5
##  4 BRANDSEN 805           26 BRANDS… -34.6 -58.4   -34.6   -34.6   -58.4   -58.4
##  5 CORRIENTES AV. 2…      27 CORRIE… -34.6 -58.4   -34.6   -34.6   -58.4   -58.4
##  6 CORRIENTES AV. 2…      28 CORRIE… -34.6 -58.4   -34.6   -34.6   -58.4   -58.4
##  7 GUARDIA VIEJA 33…      24 GUARDI… -34.6 -58.4   -34.6   -34.6   -58.4   -58.4
##  8 HUMAHUACA 4310         23 HUMAHU… -34.6 -58.4   -34.6   -34.6   -58.4   -58.4
##  9 HUMBOLDT 450           24 HUMBOL… -34.6 -58.4   -34.6   -34.6   -58.4   -58.4
## 10 PARAGUAY 4979          58 PARAGU… -34.6 -58.4   -34.6   -34.6   -58.4   -58.4
## # … with abbreviated variable names ¹​cant_insp, ²​direccion

Grafiquemos los resultamos en un mapa interactivo

leaflet() %>%
  addTiles() %>% 
  addCircleMarkers(data = insp_agrupadas10_geo, lng = ~lon, lat = ~lat )

Si bien obtenemos el resultado deseado, la función lleva su tiempo y no sabemos que está ocurriendo. Nos gustaría imprimir las direcciones que vamos geolocalizando en la consola.
Vamos a generar una nueva función para eso con un output simplificado - con geocode_OSM(as.data.frame = TRUE) obtenemos un output con las siguientes variables: query, coordenadas y bbox. Al agregar el parámetro as.data.frame = TRUE, el output es un dataframe y no un lista con subelementos.
- extraemos la latitud y la longitud
- imprimos la direccion geolocalizada

geocode_verbose <- function(x){
  
  geocode_direcciones <- geocode_OSM(x, as.data.frame = TRUE) 
  geocode_direcciones$direccion <- geocode_direcciones$query
  print(x)
  df <- select(geocode_direcciones, direccion, lat, lon)

  
}

Esta funcion recibe como input una dirección x y no un vector (una variable entera de un dataframe). Para repetir la operación sobre un vector, se pueden utilizar múltiples métodos. En este caso, vamos a utilizar una función de mapeo como map y realizar una iteración con un for loop.


Iterar una función mediante mapeo

La función map() toma como input una lista y una función (se repite la función sobre cada elemento de la lista).
Primero, se arma la lista con las direcciones.

direcciones <- list(insp_agrupadas10$direccion)

El segundo argumento que toma map es una función. Vamos a utilizar la función que ya generamos. Luego, realizamos la operación con map(). Esta función devuelve una lista. Para transformar la lista en un dataframe reducimos la lista y vamos uniendo las filas con rbind() (row bind)

direcciones_geo <- map(direcciones, geocode_verbose) %>% 
  reduce(rbind) 
##  [1] "AVELLANEDA AV. 1240, CABA, Argentina"        
##  [2] "AVENIDA DEL LIBERTADOR 4103, CABA, Argentina"
##  [3] "AVENIDA DEL LIBERTADOR 7395, CABA, Argentina"
##  [4] "BRANDSEN 805, CABA, Argentina"               
##  [5] "CORRIENTES AV. 2322, CABA, Argentina"        
##  [6] "CORRIENTES AV. 2330, CABA, Argentina"        
##  [7] "GUARDIA VIEJA 3360, CABA, Argentina"         
##  [8] "HUMAHUACA 4310, CABA, Argentina"             
##  [9] "HUMBOLDT 450, CABA, Argentina"               
## [10] "PARAGUAY 4979, CABA, Argentina"

Veamos el resultado.

head(direcciones_geo, 10)
##                                       direccion       lat       lon
## 1          AVELLANEDA AV. 1240, CABA, Argentina -34.61864 -58.44942
## 2  AVENIDA DEL LIBERTADOR 4103, CABA, Argentina -34.56911 -58.42648
## 3  AVENIDA DEL LIBERTADOR 7395, CABA, Argentina -34.54475 -58.45927
## 4                 BRANDSEN 805, CABA, Argentina -34.63568 -58.36486
## 5          CORRIENTES AV. 2322, CABA, Argentina -34.60412 -58.40534
## 6          CORRIENTES AV. 2330, CABA, Argentina -34.60412 -58.40534
## 7           GUARDIA VIEJA 3360, CABA, Argentina -34.60275 -58.41224
## 8               HUMAHUACA 4310, CABA, Argentina -34.60155 -58.42573
## 9                 HUMBOLDT 450, CABA, Argentina -34.59435 -58.44810
## 10               PARAGUAY 4979, CABA, Argentina -34.58032 -58.42743

Ahora tenemos las 10 direcciones con sus respectivas coordenadas. Vamos a unir la tabla de direcciones a la tabla de inspecciones agrupadas ‘insp_agrupadas’. Ambas tablas tienen la variable ‘direccion’ en comun. Podemos utilizar distintos tipos de uniones (join). En este caso, vamos a utilizar left_join().

insp_agrupadas10_geo <- left_join(insp_agrupadas10, direcciones_geo)
## Joining, by = "direccion"

Se realiza la unión por la variable ‘direccion’.

head(insp_agrupadas10_geo)
## # A tibble: 6 × 5
##   direccion_entidad           cant_insp direccion                      lat   lon
##   <chr>                           <int> <chr>                        <dbl> <dbl>
## 1 AVELLANEDA AV. 1240                25 AVELLANEDA AV. 1240, CABA, … -34.6 -58.4
## 2 AVENIDA DEL LIBERTADOR 4103        94 AVENIDA DEL LIBERTADOR 4103… -34.6 -58.4
## 3 AVENIDA DEL LIBERTADOR 7395        43 AVENIDA DEL LIBERTADOR 7395… -34.5 -58.5
## 4 BRANDSEN 805                       26 BRANDSEN 805, CABA, Argenti… -34.6 -58.4
## 5 CORRIENTES AV. 2322                27 CORRIENTES AV. 2322, CABA, … -34.6 -58.4
## 6 CORRIENTES AV. 2330                28 CORRIENTES AV. 2330, CABA, … -34.6 -58.4

Vamos a transformar el dataset, en un dataset espacial.

insp_agrupadas10_geo <- st_as_sf(insp_agrupadas10_geo, coords = c("lon", "lat"), crs = 4326)

Ahora vamos a mapear los lugares de inspecciones. Podemos ajustar el parámetro radius (tamaño de radio) según la cantidad de inspecciones. Como el parámetro depende de una variable, es necesario utilizar ~ para indicar que se trata de una variable.

leaflet() %>%
  addTiles() %>% 
  addCircleMarkers(data = insp_agrupadas10_geo, radius = ~cant_insp)

Iterar una función mediante un loop

Como mencionamos anteriormente, nuestra función geocode_OSM() toma como input una variable de tipo caracter. Para iterar la función sobre un dataset, podemos realizar un loop a fin de repetir la función para cada valor del vector.
Un for loop nos permite definir una iteración. La estructura del for loop es la siguiente:
En primer lugar se establece el intervalo en el cual realizar el loop. En este caso vamos a definir que el loop debe realizarse desde la primer observación hasta la última, es decir desde la observación 1 hasta la n. La ‘n’ observación puede calcularse con la función nrow() que permite calcular la cantidad de filas del dataset.

1:nrow(insp_agrupadas10)
##  [1]  1  2  3  4  5  6  7  8  9 10

Una vez establecido el rango, es necesario definir que realizar dentro del loop. En este caso, vamos a generar una nueva función que tome nuestro dataset original como input y genere un dataset como resultado.

geocode_df <- function(x, direccion = y){
  ## x es un dataframe, y es el nombre de la variable
  vector_direcciones <- x$direccion
  geocode_direcciones <- geocode_OSM(vector_direcciones, as.data.frame = TRUE) 
  x$lat <- geocode_direcciones$lat
  x$lon <- geocode_direcciones$lon
  x

}

Verifiquemos que nuestra función se comporta como esperamos con el primer valor del dataset

  direccion_test <- insp_agrupadas10[1,]
  direccion_test_geo <- geocode_df(direccion_test, direccion)
  direccion_test_geo
## # A tibble: 1 × 5
##   direccion_entidad   cant_insp direccion                              lat   lon
##   <chr>                   <int> <chr>                                <dbl> <dbl>
## 1 AVELLANEDA AV. 1240        25 AVELLANEDA AV. 1240, CABA, Argentina -34.6 -58.4

Finalmente, dentro del loop, vamos a unir las filas con rbind(). Se realiza dentro del loop para poder guardar cada output de cada iteración. Se incorpora la función print() para ir viendo cuales son las direcciones que se está geolocalizando. Para eso generamos un nuevo valor ‘estado’ que indica la posición que estamos geolocalizando.

direcciones_geo2  <- NULL
total <- nrow(insp_agrupadas10)
for(i in 1:total) {
  direccion_a_geolocalizar <- insp_agrupadas10[i,]
  direccion_unica_geo <- geocode_df(direccion_a_geolocalizar, direccion)
  estado <- paste0(i,"/", total, ": ", direccion_unica_geo$direccion )
  print(estado)
  direcciones_geo2 <- rbind(direcciones_geo2, direccion_unica_geo)
}
## [1] "1/10: AVELLANEDA AV. 1240, CABA, Argentina"
## [1] "2/10: AVENIDA DEL LIBERTADOR 4103, CABA, Argentina"
## [1] "3/10: AVENIDA DEL LIBERTADOR 7395, CABA, Argentina"
## [1] "4/10: BRANDSEN 805, CABA, Argentina"
## [1] "5/10: CORRIENTES AV. 2322, CABA, Argentina"
## [1] "6/10: CORRIENTES AV. 2330, CABA, Argentina"
## [1] "7/10: GUARDIA VIEJA 3360, CABA, Argentina"
## [1] "8/10: HUMAHUACA 4310, CABA, Argentina"
## [1] "9/10: HUMBOLDT 450, CABA, Argentina"
## [1] "10/10: PARAGUAY 4979, CABA, Argentina"

Veamos el resultado.

head(direcciones_geo2, 10)
## # A tibble: 10 × 5
##    direccion_entidad           cant_insp direccion                     lat   lon
##    <chr>                           <int> <chr>                       <dbl> <dbl>
##  1 AVELLANEDA AV. 1240                25 AVELLANEDA AV. 1240, CABA,… -34.6 -58.4
##  2 AVENIDA DEL LIBERTADOR 4103        94 AVENIDA DEL LIBERTADOR 410… -34.6 -58.4
##  3 AVENIDA DEL LIBERTADOR 7395        43 AVENIDA DEL LIBERTADOR 739… -34.5 -58.5
##  4 BRANDSEN 805                       26 BRANDSEN 805, CABA, Argent… -34.6 -58.4
##  5 CORRIENTES AV. 2322                27 CORRIENTES AV. 2322, CABA,… -34.6 -58.4
##  6 CORRIENTES AV. 2330                28 CORRIENTES AV. 2330, CABA,… -34.6 -58.4
##  7 GUARDIA VIEJA 3360                 24 GUARDIA VIEJA 3360, CABA, … -34.6 -58.4
##  8 HUMAHUACA 4310                     23 HUMAHUACA 4310, CABA, Arge… -34.6 -58.4
##  9 HUMBOLDT 450                       24 HUMBOLDT 450, CABA, Argent… -34.6 -58.4
## 10 PARAGUAY 4979                      58 PARAGUAY 4979, CABA, Argen… -34.6 -58.4

De ambas maneras se obtienen los mismos resultados.




Localizando todas las inspecciones

Primero vamos a incluir la ciudad y el país a nuestras direcciones.

inspecciones <- inspecciones %>% 
  mutate(direccion = paste0(direccion_entidad, ", CABA, Argentina"))

¿Cuántas inspecciones tenemos que geolocalizar?

dim(inspecciones)
## [1] 47445     8

Vemos que nuestro dataset tiene 47455 inspecciones. Vamos a generar un dataset que contenga solamente las calles, y vamos a quedarnos con registros únicos, para evitar geolocalizar más de una vez la misma dirección.

direcciones_insp <- inspecciones %>%
  select(direccion) %>%
  unique()
dim(direcciones_insp)
## [1] 28405     1

La cantidad de direcciones únicas es mucho menor. Esto va a ahorrar tiempo de cómputo.

Para geolocalizar todas las direcciones de las inspecciones, se utilizó la siguiente función.
Se incorporaron los parámetros de if / else para evitar que la función falle en caso de que no se encuentre la dirección correspondiente.

geolocalizador <- function(x){
  geocode_direcciones <- geocode_OSM(x, as.data.frame = TRUE) 
  if (is_null(geocode_direcciones)) {
    df <- data.frame(direccion = x, 
                     lat = NA, 
                     lon = NA)
  }
  else {
    geocode_direcciones$direccion <- geocode_direcciones$query
     df <- select(geocode_direcciones, direccion, lat, lon)
  }
}

Se genera el dataset vacío con las direcciones geolocalizadas

direcciones_geolocalizadas  <- NULL
#este código no corre en el Rmarkdown dado que tardaría mucho tiempo en geolocalizar las direcciones
x <- direcciones_insp
total_direcciones <- nrow(x)

Para el Rmarkdown, se corre a modo de ejemplo una muestra de 20 direcciones.

set.seed(10)
x <- sample_n(direcciones_insp, 20)
total_direcciones <- nrow(x)
for(i in 1:nrow(x)) {
  direccion_unica <- x$direccion[i]
  direccion_unica_geo <- geolocalizador(direccion_unica)

  
  print(paste0(i, "/", total_direcciones, ": ", direccion_unica))
  direcciones_geolocalizadas <- rbind(direcciones_geolocalizadas, direccion_unica_geo)
}
## [1] "1/20: URUGUAY 250, CABA, Argentina"
## [1] "2/20: CARACAS 4778, CABA, Argentina"
## [1] "3/20: OLAZABAL 3199, CABA, Argentina"
## [1] "4/20: VALLE 602, CABA, Argentina"
## [1] "5/20: VARELA AV. 912, CABA, Argentina"
## [1] "6/20: MOSCONI GENERAL AV. 3075, CABA, Argentina"
## [1] "7/20: AVENIDA LA PLATA 239, CABA, Argentina"
## [1] "8/20: TRES ARROYOS 564, CABA, Argentina"
## [1] "9/20: AZCUENAGA 1150, CABA, Argentina"
## No results found for "LINCOLN AV. 3899, CABA, Argentina".
## [1] "10/20: LINCOLN AV. 3899, CABA, Argentina"
## [1] "11/20: CASTILLO 1366, CABA, Argentina"
## [1] "12/20: CORDOBA AV. 5285, CABA, Argentina"
## [1] "13/20: BILLINGHURST 1689, CABA, Argentina"
## [1] "14/20: DESAGUADERO 3441, CABA, Argentina"
## [1] "15/20: PASO 770, CABA, Argentina"
## [1] "16/20: OLAZABAL AV. 4466, CABA, Argentina"
## [1] "17/20: DIAZ VELEZ AV. 4770, CABA, Argentina"
## [1] "18/20: COCHABAMBA 2934, CABA, Argentina"
## [1] "19/20: LOPE DE VEGA AV. 3238, CABA, Argentina"
## [1] "20/20: 20 de Septiembre 411, CABA, Argentina"
head(direcciones_geolocalizadas, 20)
##                                    direccion       lat       lon
## 1               URUGUAY 250, CABA, Argentina -34.60590 -58.38650
## 2              CARACAS 4778, CABA, Argentina -34.58144 -58.49595
## 3             OLAZABAL 3199, CABA, Argentina -34.56458 -58.46593
## 4                 VALLE 602, CABA, Argentina -34.62470 -58.43637
## 5            VARELA AV. 912, CABA, Argentina -34.63963 -58.45828
## 6  MOSCONI GENERAL AV. 3075, CABA, Argentina -34.58763 -58.50451
## 7      AVENIDA LA PLATA 239, CABA, Argentina -34.61765 -58.42881
## 8          TRES ARROYOS 564, CABA, Argentina -34.60811 -58.45893
## 9            AZCUENAGA 1150, CABA, Argentina -34.59462 -58.39975
## 10         LINCOLN AV. 3899, CABA, Argentina        NA        NA
## 11            CASTILLO 1366, CABA, Argentina -34.58858 -58.44097
## 12         CORDOBA AV. 5285, CABA, Argentina -34.59920 -58.40394
## 13        BILLINGHURST 1689, CABA, Argentina -34.59108 -58.41068
## 14         DESAGUADERO 3441, CABA, Argentina -34.61024 -58.51750
## 15                 PASO 770, CABA, Argentina -34.59954 -58.40327
## 16        OLAZABAL AV. 4466, CABA, Argentina -34.55933 -58.45893
## 17      DIAZ VELEZ AV. 4770, CABA, Argentina -34.60841 -58.41684
## 18          COCHABAMBA 2934, CABA, Argentina -34.62590 -58.40604
## 19    LOPE DE VEGA AV. 3238, CABA, Argentina -34.62678 -58.50868
## 20     20 de Septiembre 411, CABA, Argentina -34.63054 -58.36366

Podemos cargar las direcciones ya geolocalizadas

direcciones_geolocalizadas <- read.csv("https://raw.githubusercontent.com/paulavidela/utdt_cienciadedatos/main/data/direcciones_geolocalizadas.csv", encoding = "UTF-8")

Vamos a unir las direcciones a las inspecciones y vamos a transformar el dataset a espacial.

inspecciones_geo <- inspecciones %>%
  left_join(direcciones_geolocalizadas, by = "direccion" ) %>%
  filter(!is.na(lat) & !is.na(lon)) %>%
  st_as_sf(coords = c("lon", "lat"), crs = 4326)

Grafiquemos los resultados

ggplot() + 
  geom_sf(data = inspecciones_geo )

Ajustes en leaflet

¿Cómo son las inspecciones nocturas?

table(inspecciones$area)
## 
##                                   Antenas 
##                                        35 
##         Auditorias Integrales Programadas 
##                                      1998 
##                                       AVO 
##                                      4314 
##                                Eléctricas 
##                                       406 
##                                Elevadores 
##                                      1354 
## Establecimientos y Productos Alimenticios 
##                                      8260 
##                                  Incendio 
##                                      1535 
##      Mercadería en Tránsito y Vía Pública 
##                                         1 
##                                     Obras 
##                                      7438 
##                            Sanitarias AVO 
##                                         1 
##                       SGO Est. Deportivos 
##                                       956 
##                       SGO Est. Educativos 
##                                       752 
##              SGO Est. Hoteles y Pensiones 
##                                      1082 
##              SGO Est. Salud y Geriátricos 
##                                       582 
##                    SGO Eventos Deportivos 
##                                       591 
##                       SGO Eventos Masivos 
##                                       329 
##                 SGO Fiscalización Crítica 
##                                      2687 
##                SGO Fiscalización Nocturna 
##                                      4000 
##                   SGO Grandes Superficies 
##                                      1430 
##                  SGO Inspecciones simples 
##                                      5943 
##               SGO Suspensiones y Rechazos 
##                                      1115 
##                        SGO Verificaciones 
##                                      1577 
##                                  Térmicas 
##                                      1059

Queremos calcular la cantidad de días de inspecciones, y la cantidad de inspecciones totales. Luego unimos al dataset resultante la geolocalización.

insp_nocturnas_total <- inspecciones %>%
  filter(area == "SGO Fiscalización Nocturna") %>%
  group_by(direccion, fecha_inspeccion) %>%
  summarise(cantidad = n()) %>%
  ungroup() %>%
  group_by(direccion) %>%
  summarise(dias_inspeccion = n(), 
            total_inspecciones = sum(cantidad, na.rm = TRUE)) %>%
  ungroup() %>%
  left_join(direcciones_geolocalizadas, by = "direccion") %>%
  filter(!is.na(lat) & !is.na(lon)) %>%
  st_as_sf(coords = c("lon", "lat"), crs = 4326)
## `summarise()` has grouped output by 'direccion'. You can override using the
## `.groups` argument.

El dataset de insp_nocturas_total tiene 2024 observaciones. Para facilitar el trabajo con leaflet vamos a seleccionar aquellos lugares que cuentan con más de 3 inspecciones

insp_nocturnas <- insp_nocturnas_total %>%
  filter(total_inspecciones > 3)

Veamos los resultados

head(insp_nocturnas)
## Simple feature collection with 6 features and 3 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: -58.48319 ymin: -34.61341 xmax: -58.41092 ymax: -34.59929
## geographic CRS: WGS 84
## # A tibble: 6 × 4
##   direccion                            dias_…¹ total…²              geometry
##   <chr>                                  <int>   <int>           <POINT [°]>
## 1 24 DE NOVIEMBRE 225, CABA, Argentina       5       5 (-58.41152 -34.61341)
## 2 AGUERO 445, CABA, Argentina               12      12  (-58.4118 -34.60465)
## 3 AGUERO 449, CABA, Argentina                6       6 (-58.41179 -34.60462)
## 4 AGUERO 806, CABA, Argentina               10      10 (-58.41092 -34.60037)
## 5 ALVAREZ JONTE 2294, CABA, Argentina        5       6 (-58.46822 -34.59929)
## 6 ALVAREZ JONTE AV. 3402, CABA, Argen…       5       5 (-58.48319 -34.60904)
## # … with abbreviated variable names ¹​dias_inspeccion, ²​total_inspecciones

Hagamos un primer mapa interactivo

leaflet() %>%
  addTiles() %>% 
  addCircleMarkers(data = insp_nocturnas)

Ajustes básicos

En el mapa anterior es muy difícil distinguir concentraciones, podemos ajustar el radio de los círculos y la transparencia.

leaflet() %>%
  addTiles() %>% 
  addCircleMarkers(data = insp_nocturnas, 
                   radius = 1.5, 
                   fillOpacity = 0.4, 
                   opacity = 0.5)

Si bien se visualiza un poco mejor, aún no se logran distinguir concentraciones.

leaflet() %>%
  addTiles() %>% 
  addCircleMarkers(data = insp_nocturnas, 
                   radius = 1, 
                   fillOpacity = 0.5, 
                   opacity = 0.6, 
                   clusterOptions = markerClusterOptions())

Como tenemos las inspecciones nocturnas por cantidad y dirección, queremos visualizar en un mapa como varía la cantidad. Probemos con el radio.

leaflet() %>%
  addTiles() %>% 
  addCircleMarkers(data = insp_nocturnas, 
                   radius = ~total_inspecciones, 
                   color = "red")

¿Qué pasa si ajustamos el color?

leaflet() %>%
  addTiles() %>% 
  addCircleMarkers(data = insp_nocturnas, 
                   color = ~total_inspecciones, 
                   radius = 4)

El color no se ajusta…

Ajustes de colores

Para utilizar escalas de color en leaflet primero tenemos que configurarlas. Leaflet cuenta con funciones según el tipo de escala, por ejemplo colorNumeric(), colorFactor(), colorQuartile(), etc. En nuestro caso se trata de una escala numérica, entonces vamos a utilizar colorNumeric().

paleta <- colorNumeric(palette = "viridis", domain =  insp_nocturnas$total_inspecciones)

El objeto que creamos es una función. Vamos a aplicar esta función a nuestra variable ‘cantidad’ para que varíe el color.

leaflet() %>%
  addTiles() %>% 
  addCircleMarkers(data =  insp_nocturnas, 
                   color = ~paleta(total_inspecciones), 
                   radius = 4)

Los colores más claros se ven muy lavados. Vamos a seleccionar nosotros nuestra escala.

paleta <- colorNumeric(palette = c("#c7e9b4", "#1d91c0", "#081d58" ), 
                       domain =  insp_nocturnas$total_inspecciones)
leaflet() %>%
  addTiles() %>% 
  addCircleMarkers(data =  insp_nocturnas,  
                   color = ~paleta(total_inspecciones), 
                   radius = 2, 
                   fillOpacity = 0.5, 
                   opacity = 0.7)

El fondo no colabora. Vamos a cambiar el fondo por un fondo en blanco y negro. Para eso vamos a cambiar la función addTiles() por addProviderTiles(). Se puede elegir el proveedor desde aquí

leaflet() %>%
  addProviderTiles(provider = "CartoDB.Positron" ) %>%
  addCircleMarkers(data = insp_nocturnas, 
                   color = ~paleta(total_inspecciones), 
                   radius = 2)

Ahora si, mucho mejor!

head(insp_nocturnas)
## Simple feature collection with 6 features and 3 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: -58.48319 ymin: -34.61341 xmax: -58.41092 ymax: -34.59929
## geographic CRS: WGS 84
## # A tibble: 6 × 4
##   direccion                            dias_…¹ total…²              geometry
##   <chr>                                  <int>   <int>           <POINT [°]>
## 1 24 DE NOVIEMBRE 225, CABA, Argentina       5       5 (-58.41152 -34.61341)
## 2 AGUERO 445, CABA, Argentina               12      12  (-58.4118 -34.60465)
## 3 AGUERO 449, CABA, Argentina                6       6 (-58.41179 -34.60462)
## 4 AGUERO 806, CABA, Argentina               10      10 (-58.41092 -34.60037)
## 5 ALVAREZ JONTE 2294, CABA, Argentina        5       6 (-58.46822 -34.59929)
## 6 ALVAREZ JONTE AV. 3402, CABA, Argen…       5       5 (-58.48319 -34.60904)
## # … with abbreviated variable names ¹​dias_inspeccion, ²​total_inspecciones

Pop-ups y etiquetas

Vamos a incorporar pop-ups. Podemos configurar los pop-ups a partir de una variable. Nos gustaría que el texto del pop-up sea:
- El establecimiento fue inspeccionado dias_inspeccion días distintos
- En total se realizaron total_inspecciones inspecciones

Para lograr este pop-up necesitamos incluir un ‘enter’ entre ambas frases. Para eso incorporamos código en HTML.
</br> permite ir a la línea
<b> pasa la letra a negrita. Con </b> se termina la negrita
<i> pasa la letra a itálica. Con </i> se termina la itálica

insp_nocturnas <- insp_nocturnas %>%
  mutate(leyenda = paste0( "El establecimiento fue inspeccionado ", dias_inspeccion, " días distintos", "</br>", 
                           "En total se realizaron " , "<b>",total_inspecciones, " inspecciones"  , "</b>"
    
  ))
leaflet() %>%
  addProviderTiles(provider = "CartoDB.Positron" ) %>%
  addCircleMarkers(data = insp_nocturnas, 
                   color = ~paleta(total_inspecciones), 
                   popup = ~leyenda,
                   radius = 4)

También podemos incorporar información cuando el cursor se encuentra sobre el establecimiento.

leaflet() %>%
  addProviderTiles(provider = "CartoDB.Positron" ) %>%
  addCircleMarkers(data = insp_nocturnas, 
                   color = ~paleta(total_inspecciones), 
                   popup = ~leyenda,
                   label = ~direccion,
                   radius = 4)

Elementos adicionales

Podemos incorporar una referencia a nuestra escala de color. Para eso, utilizaremos la función addLegend

paleta <- colorNumeric(palette = c("#c7e9b4", "#1d91c0", "#081d58" ), 
                       domain = insp_nocturnas$total_inspecciones)

leaflet() %>%
  addProviderTiles(provider = "CartoDB.Positron" ) %>%
  addCircleMarkers(data = insp_nocturnas, 
                   color = ~paleta(total_inspecciones), 
                   popup = ~leyenda,
                   label = ~direccion,
                   radius = 4) %>%
    addLegend(data = insp_nocturnas,
              "bottomright", 
              pal = paleta, 
              values = ~total_inspecciones,
              title = "Inspecciones",
              opacity = 1)

Otras geometrías

En leaflet, es posible agregar otras geometrías. Vamos a incluir los barrios de la ciudad de Buenos Aires.

barrios <- st_read("https://raw.githubusercontent.com/paulavidela/utdt_cienciadedatos/main/data/barrios.geojson")
## Reading layer `barrios2' from data source `https://raw.githubusercontent.com/paulavidela/utdt_cienciadedatos/main/data/barrios.geojson' using driver `GeoJSON'
## Simple feature collection with 48 features and 6 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -58.53152 ymin: -34.70529 xmax: -58.33515 ymax: -34.52649
## geographic CRS: WGS 84
leaflet() %>%
  addPolygons(data = barrios, fillColor = NA, color = "black") %>%
  addProviderTiles(provider = "CartoDB.Positron" ) %>%
  addCircleMarkers(data = insp_nocturnas, 
                   color = ~paleta(total_inspecciones), 
                   popup = ~leyenda,
                   label = ~direccion,
                   radius = 4) %>%
 
    addLegend(data = insp_nocturnas,
              "bottomright", 
              pal = paleta, 
              values = ~total_inspecciones,
              title = "Inspecciones",
              opacity = 1)

Bonus! Geocodificando con Here

Here es una plataforma orientada a soluciones de localización y mobilidad. A diferencia de OSM, no es una plataforma abierta, sin embargo cuenta con un nivel gratuito que permite consultar realizar consultas limitadas a la API.
Al igual que OSM, cuenta con una API para geocodificar direcciones. Además, se desarrolló en R el paquete hereR que permite utilizar los servicios desde R.

#install.packages("hereR")
library(hereR)

Para utilizar HERE es necesario contar con una llave (KEY). Para eso, primero hay que crearse una cuenta en la plataforma de desarrollador. No es necesario agregar datos de tarjeta. Dentro de la plataforma, hay que crear generar una nueva app. Allí se obtiene la API KEY.

# esta key es ficticia
set_key("llavejemplo-aqui-va-la-llave")

Al igual que con Nominatim, la librería hereR cuenta con una función para geolocalizar: geocode(). Probemos localizar nuevamente la Universidad Torcuato Di Tella. La universidad se encuentra en Figueroa Alcorta 7350, CABA.

utdt_here <- geocode("Figueroa Alcorta 7350, CABA")

La función geocode() también existe en el paquete ggmap. En caso de tener ambos paquetes activos, es recomendable declarar de que paquete es la función.

Veamos el output que generamos con la función de Here.

head(utdt_here)
## Simple feature collection with 1 feature and 12 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: -58.44668 ymin: -34.54719 xmax: -58.44668 ymax: -34.54719
## geographic CRS: WGS 84
##   id
## 1  1
##                                                                            address
## 1 Avenida Presidente Figueroa Alcorta 7350, 1428 Ciudad de Buenos Aires, Argentina
##          type                              street house_number postal_code
## 1 houseNumber Avenida Presidente Figueroa Alcorta         7350        1428
##   district                   city                 county
## 1 Belgrano Ciudad de Buenos Aires Ciudad de Buenos Aires
##                             state   country                     access
## 1 Ciudad Autónoma de Buenos Aires Argentina POINT (-58.44665 -34.5471)
##                      geometry
## 1 POINT (-58.44668 -34.54719)

A diferencia de Nomitatim, Here devuelve por defecto un sf. Incluye más información como el código postal.

leaflet() %>%
  addTiles() %>% 
  addCircleMarkers(data = utdt_here)

Existen otras formas de geocodificar, por ejemplo utilizando la API de Google. Para utilizarla, es necesario utilizar una tarjeta de crédito. Aquí hay un tutorial de como utilizar Google. Otra opción, es utilizar OpenCage, un servicio similar a Here. Aquí hay un tutorial.