Mapas con R

Visualización de datos geoespaciales con R.

En este taller veremos como visualizar datos geoespaciales (vectoriales) con R, sin embargo, los paquetes de R permiten hacer mucho más que sólo visualizar. Para aprender más sobre esto ver Geocomputation with R (Lovelace, Nowosad, y Muenchow 2019).

💻 Preparación

Instala los paquetes y descarga el material del taller.

install.packages(c("usethis", "tidyverse", "sf", "leaflet"))
usethis::use_course("https://github.com/tereom/taller-mapas-mv/blob/master/material.zip?raw=true")

🗺 Estáticos

Usaremos los paquetes tidyverse (Wickham et al. 2019) y sf(Pebesma 2018).

library(tidyverse) # arreglar datos (dplyr) y graficar (ggplot2)
library(sf) # simple features para leer y manipular datos espaciales

🐶 Nuestro primer mapa

• Objetivo: Un mapa de población por departamento de Uruguay

• ¿Qué necesito?

  1. Datos que determinan la forma de los departamentos (shapefile),
     
  2. Datos de población por departamento (csv)

Los datos geoespaciales de los departamentos los obtenemos de mapas vectoriales del INE y usaremos la función st_read() para leerlos:

depart_uy <- st_read("datos/ine_depto")

Reading layer `ine_depto' from data source `/Users/tereom/Documents/mapas-taller-mv/datos/ine_depto' using driver `ESRI Shapefile'
Simple feature collection with 20 features and 5 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 366582.2 ymin: 6127919 xmax: 858252.1 ymax: 6671738
Projected CRS: WGS 84 / UTM zone 21S

La función st_read() recibe como argumento la ruta a la carpeta (puede ser al archivo con terminación .shp) y la asigna a depart_uy.

Los podemos graficar directo con plot(depart_uy) o con ggplot (❤️) en el segundo caso usamos geom_sf() para agregar la capa de coordenadas:

ggplot(depart_uy) +
  geom_sf()

Si inspeccionamos el objeto depart_uy descubrimos que es un poco familiar (y un poco nuevo), lo familiar es que es una tabla de datos de tipo data.frame, con una columna adicional de geometría que almacena la información de la geometría espacial.

class(depart_uy)

[1] "sf"         "data.frame"

glimpse(depart_uy)

Rows: 20
Columns: 6
$ AREA_KM2_  <int> 530, 11928, 4536, 6106, 11643, 10417, 10016, 1392…
$ PERIMETER  <dbl> 202136.96, 722385.94, 404534.80, 431056.43, 82132…
$ DEPTO      <int> 1, 2, 3, 5, 6, 8, 9, 11, 12, 13, 14, 15, 16, 17, …
$ NOMBRE     <chr> "MONTEVIDEO", "ARTIGAS", "CANELONES", "COLONIA", …
$ CDEPTO_ISO <chr> "UYMO", "UYAR", "UYCA", "UYCO", "UYDU", "UYFD", "…
$ geometry   <MULTIPOLYGON [m]> MULTIPOLYGON (((579197 6136..., MULT…

En la gráfica podemos colorear (como en cuaqluier ggplot) los polígonos utilizando con el valor de alguna de las columnas del data.frame.

ggplot(depart_uy) +
  geom_sf(aes(fill = AREA_KM2_))

También podemos unirle una tabla de datos externa, por ejemplo de población (datos de juntas departamentales) como lo haríamos de normal:

basicos <- read_csv("datos/datosbasicosjds.csv")
glimpse(basicos %>% select(1:10, 38))

Rows: 19
Columns: 11
$ `Junta Departamental` <chr> "Artigas", "Canelones", "Cerro Largo",…
$ `Tiene WEB`           <chr> "NO", "Si", "Si", "Si", "Si", "Si", "S…
$ `Link de la web`      <chr> NA, "http://www.juntadecanelones.gub.u…
$ Dirección             <chr> "Cnel. Carlos Lecueder 510", "Luis Alb…
$ `Teléfono 1`          <chr> "4772 4833", "4332 2420", "4642 2283",…
$ `Telefono 2`          <chr> NA, NA, "4642 3471", NA, "4362 4778", …
$ Email                 <chr> "juntadeartigas@gmail.com", "juntacan@…
$ Plenarios             <chr> "1er - 2do y 3er jueves de cada mes", …
$ `Hora Plenario`       <chr> "20.30", "18.00", "19.00", "20.00", "1…
$ `Local propio`        <chr> "Si", "Si", "Si", "Si", "Si", "No", "S…
$ Población             <chr> "73.377", "520.173", "84.698", "123203…

# un poco de limpieza para poder unir las tablas
depart_pob <- basicos %>% 
  select(nombre = `Junta Departamental`, pob = `Población`) %>% 
  mutate(nombre = janitor::make_clean_names(nombre), 
         pob = as.numeric(str_remove_all(pob, "\\.")))

depart_uy_pob <- depart_uy %>% 
  mutate(nombre = janitor::make_clean_names(NOMBRE)) %>% 
  left_join(depart_pob)

glimpse(depart_uy_pob)

Rows: 20
Columns: 8
$ AREA_KM2_  <int> 530, 11928, 4536, 6106, 11643, 10417, 10016, 1392…
$ PERIMETER  <dbl> 202136.96, 722385.94, 404534.80, 431056.43, 82132…
$ DEPTO      <int> 1, 2, 3, 5, 6, 8, 9, 11, 12, 13, 14, 15, 16, 17, …
$ NOMBRE     <chr> "MONTEVIDEO", "ARTIGAS", "CANELONES", "COLONIA", …
$ CDEPTO_ISO <chr> "UYMO", "UYAR", "UYCA", "UYCO", "UYDU", "UYFD", "…
$ nombre     <chr> "montevideo", "artigas", "canelones", "colonia", …
$ pob        <dbl> 1318755, 73377, 520173, 123203, 57084, 67047, 588…
$ geometry   <MULTIPOLYGON [m]> MULTIPOLYGON (((579197 6136..., MULT…

Utilizamos el paquete janitor para arreglar los nombres y que coincidan entre las dos tablas que vamos a unir. Unimos con la función left_join() como haríamos en una manipulación de datos ususal de data.frames.

Y graficamos nuestro mapa coloreado por población.

ggplot(depart_uy_pob) +
  geom_sf(aes(fill = pob))

😕 Lo no tan familiar

Shapefiles

Ahora inspeccionamos el lado nuevo de nuestro objeto, cuando leímos los datos apuntamos a una carpeta:

depart_uy <- st_read("datos/ine_depto")

esta carpeta de donde leímos los datos tiene muchos archivos.

fs::dir_tree("datos/ine_depto/")

datos/ine_depto/
├── ine_depto.dbf
├── ine_depto.prj
├── ine_depto.sbn
├── ine_depto.sbx
├── ine_depto.shp
├── ine_depto.shx
└── ine_depto_metadatos.pdf

Esto es porque, en este ejemplo, nuestros datos vectoriales están almacenados en un shapefile (un formato muy común para datos geoespaciales).

Un shapefile es un grupo de archivos que contienen geometrías e información de acuerdo a un estándar especificado por el Insituto de Investigación en Sistemas de Ecosistemas (ESRI). En nuestro ejemplo tenemos los siguientes:

• nombre.shp: es el archivo principal y contiene la geometría.

• nombre.dbf: es la base de datos asoiciada al objeto y almacena la información de los atributos de los objetos, en el caso de ine_depto hay un renglón por polígono, es decir, cada línea contiene las características de un departamento.

ine_depto_dbf <- foreign::read.dbf("datos/ine_depto/ine_depto.dbf")
head(ine_depto_dbf)

  AREA_KM2_ PERIMETER DEPTO     NOMBRE CDEPTO_ISO
1       530  202137.0     1 MONTEVIDEO       UYMO
2     11928  722385.9     2    ARTIGAS       UYAR
3      4536  404534.8     3  CANELONES       UYCA
4      6106  431056.4     5    COLONIA       UYCO
5     11643  821325.4     6    DURAZNO       UYDU
6     10417  669409.9     8    FLORIDA       UYFD

• nombre.prj (opcional): sistema de referencia de coordenadas asociado a los datos.

crs_depart <- st_crs(depart_uy)
crs_depart$input

[1] "WGS 84 / UTM zone 21S"

• nombre.shx: almacena el índice de las entidades geométricas.

• nombre.sbn y .sbx (opcional): almacena índice espacial.

Como vimos en el ejemplo, read_sf() recibe la ruta a la carpeta o a alguno de los archivos de la capa de interés (en caso de tener varias capas en una carpeta podemos indicar la ruta hasta alguno de los archivos): read_sf("datos/ine_depto/ine_depto.shp") o utilizar el argumento layer: read_sf("datos/ine_depto", layer = "ine_depto")). Aún cuando un shapefile este almacenado en un conjunto de archivos, ésta instrucción es lo único que necesitamos para cargar los datos en R.

Geometrías

El segundo aspecto nuevo, es la columna geometry en los datos,

depart_uy$geometry

Geometry set for 20 features 
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 366582.2 ymin: 6127919 xmax: 858252.1 ymax: 6671738
Projected CRS: WGS 84 / UTM zone 21S
First 5 geometries:

ésta puede almacenar distintos tipos de datos vectoriales.

Los tipos van aumentando en complejidad, el más simple serían puntos, cuya unión forma líneas, con estos dos formas polígonos. Algunas de las geometrías más comunes están representadas en la siguiente imagen:

Figure 1: Lovelace, Nowosad, y Muenchow (2019) (CC BY-NC-ND 4.0).

Tu turno

🌎 Grafica los polígonos de barrios de Montevideo (los archivos se encuentran en datos/ine_barrios).

🌐 Lee el archivo en datos/uptu_cultura y agrega las ubicaciones contenidas a la misma gráfica que los polígonos.

🌍 Describe la diferencia entre la geometría de los barrios de Montevideo y del archivo de cultura.

🤸 Un ejemplo más

Notemos que podemos utilizar las operaciones usuales para transformar tablas de datos, por ejemplo, si queremos crear una columna de tipo de establecimiento usamos la función mutate de dplyr.

glimpse(cultura_mv)

Rows: 398
Columns: 3
$ NOMBRE    <chr> "TEATRO ALFREDO MORENO", "TEATRO ALIANZA URUGUAY-E…
$ DIRECCION <chr> NA, "PARAGUAY 1217", "SAN JOSE 1426", "COLONIA 187…
$ geometry  <POINT [m]> POINT (582125 6138279), POINT (573785.3 6136…

cultura_mv <- cultura_mv %>% 
  mutate(tipo_establecimiento = case_when(
    str_detect(NOMBRE, "TEATRO") ~ "teatro", 
    str_detect(NOMBRE, "MUSEO") ~ "museo",
    str_detect(NOMBRE, "BIBLIOTECA") ~ "biblioteca",
    str_detect(NOMBRE, "CINE") ~ "cine",
    TRUE ~ "otro"
  ))

Usamos la función case_when() para categorizar los puntos de interés culturales de acuerdo al tipo, queremos esta variable para mostrar la información en un mapa.

glimpse(cultura_mv)

Rows: 398
Columns: 4
$ NOMBRE               <chr> "TEATRO ALFREDO MORENO", "TEATRO ALIANZ…
$ DIRECCION            <chr> NA, "PARAGUAY 1217", "SAN JOSE 1426", "…
$ geometry             <POINT [m]> POINT (582125 6138279), POINT (57…
$ tipo_establecimiento <chr> "teatro", "teatro", "teatro", "teatro",…

Veamos que tenemos disponibles las opciones típicas de ggplot como facet_wrap().

ggplot() +
  geom_sf(data = barrios_mv) +
  geom_sf(data = cultura_mv, color = "red", alpha = 0.8, size = 0.5) + 
  labs(subtitle = "Establecimientos culturales en Montevideo", 
       color = "tipo") +
  theme_void() + 
  facet_wrap(~tipo_establecimiento)

📄 Crear uno objeto sf de una tabla de datos

Imaginemos que tenemos una tabla de datos georreferenciados, pueden ser datos colectados con gps en una encuesta o ubicaciones obtenidas de google maps. Lo podemos almacenar en un data.frame que podemos a su vez convertir en un objeto de tipo sf, lo que debemos hacer es indicar cuáles son las columnas que incluyen información geográfica.

• 🖍 Utiliza google maps para ubicar un lugar especial para ti y agrégalo a este documento.

library(tidyverse)
library(googlesheets4)
url <- ("https://docs.google.com/spreadsheets/d/1h7n3dpjqGofWYtKf5ADCDAr8nT5U5UalMFq9nS3G_ew/edit#gid=0")
gs4_deauth()
colab <- read_sheet(url)
glimpse(colab)

Rows: 17
Columns: 3
$ nombre <chr> "Teresa", "Pablo", "Flo", "Mica", "Antoine", "Lucía",…
$ long   <dbl> 18.95584, -34.93513, -32.14170, -34.79537, -34.93536,…
$ lat    <dbl> -99.22654, -56.16043, -53.73057, -55.86009, -56.16036…

Por ahora colab es simplemente un data.frame,

class(colab)

[1] "tbl_df"     "tbl"        "data.frame"

Lo podemos convertir a un objeto de clase sf con la función st_as_sf, indicando cuales son las columnas con información geográfica. De manera opcional podemos indicar también la información de la proyección (Google Maps utiliza la 4326 para reportar).

colab_longlat <- colab %>% 
  st_as_sf(coords = c("lat", "long")) %>% 
  st_sf(crs = 4326)

class(colab_longlat)

[1] "sf"         "tbl_df"     "tbl"        "data.frame"

Y ahora lo podemos graficar con ggplot y geom_sf().

ggplot(colab_longlat) +
  geom_sf(aes(color = nombre))

📏 CRS y Proyecciones

Como vimos con nuestros primeros mapas, los objetos espaciales suelen tener asociado un valor en sistema de referencia de coordenadas (CRS), en esta sección explicamos a qué se refiere y porque es importante.

La Tierra tiene forma elipsoidal y para ubicar puntos en la Tierra debemos especificar un modelo de elipse, un punto de origen y dirección de los ejes, y una estrategia para proyectar a un plano.

1. Elipse: Describe la forma de la Tierra, es una aproximación que no ajusta a la perfección. Existen distintos elipsoides, algunos están diseñados para ajustar toda la Tierra (WGS84, GRS80) y algunos ajustan únicamente por regiones (NAD27). Los locales son más precisos para el área para la que fueron diseñados pero no funcionan en otras partes del mundo.

rgdal::projInfo(type = "ellps") %>% 
  rmarkdown::paged_table()

2. Datum: Define un punto de origen para el sistema de coordenadas de la Tierra y la dirección de los ejes. El datum define el elipsoide que se usará pero el sentido contrario no es cierto.

1. Las distancias deben tener unidades (metros en el ej.).
2. Se necesita un punto de referencia (datum).

Los conceptos anteriores son la base de un sistema de coordenadas, es importante recalcar que no todos los datos geoespaciales usan el mismo sistema de coordenadas. Y para ser representados en un mismo mapa debemos cuadrarlos.

3. Proyección: Nos permiten desplegar el elipsoide en un espacio de dos dimensiones, es decir pasar de longitud y latitud en el elipsoide a coordenadas de un mapa, esto facilita ciertas operaciones (por ejemplo. las longitudes y los ángulos son constantes en dos dimensiones, en contraste con una esfera donde las líneas de longitud convergen en los polos), sin embargo, conlleva necesariamente una distorsión. La estrategia usual para proyectar es utilizar superficies de desarrollo y después aplanar la superficie.

Figure 2: Comparación entre proyecciones tangente y secante cilíndirca, cónica y de azimut con licencia CC-BY 4.0

Es posible trabajar con datos no proyectados (longitud/latitud), a estos se les determina coordenadas geográficas.

Las distorsiones resultan en pérdidas de información que pueden ser en área, forma, distancia o dirección. Por ejemplo Mercator preserva dirección y es útil para navegación:

Show code

states <- map_data("state")
usmap <- ggplot(states, aes(x=long, y=lat, group=group)) +
  geom_polygon(fill="white", colour="black")
usmap + coord_map("mercator")

Por su parte Azimuthal Equal Area preserva area pero no dirección:

Show code

usmap + coord_map("azequalarea") 

En particular en investigación es común usar Universal Transverse Mercator (UTM) porque preserva ángulos y dirección (por lo que es común en sistemas de navegación). Para minimizar la distorsión se divide la Tierra en 60 regiones y utiliza una proyección (de secantes) en cada una. Esta es la proyección más común para Uruguay.

Figure 3: Husos y Zonas UTM con licencia CC-BY 3.0

📂 Sistemas de referencia de coordenadas (CRS)

Como vimos arriba, los conjuntos de datos geoespaciales tienen asociado un sistema de referencia de coordenadas (y proyección si es el caso), Hay varias maneras de hacer referencia al sistema asociado a un conjunto de datos, una notación estándar que se utiliza para describir un CRS es proj4string y se ven como:

st_crs(barrios_mv)$proj4string # cartesianas

[1] "+proj=utm +zone=21 +south +datum=WGS84 +units=m +no_defs"

st_crs(colab_longlat)$proj4string # geográficas

[1] "+proj=longlat +datum=WGS84 +no_defs"

También se puede usar el código ESPG asociado, los códigos EPSG. Los códigos EPSG (European Petroleum Survey Group) son códigos de 4 a 5 dígitos que representan definiciones de CRS:

st_crs(32721)$proj4string

[1] "+proj=utm +zone=21 +south +datum=WGS84 +units=m +no_defs"

Hay más maneras de documentar el sistema de coordenadas, sin embargo las dos mencionadas son muy comunes.

Tu turno

🗾 En un editor de texto abre el archivo .prj de alguno de los conjuntos de datos geoespaciales que hemos utilizado.

🗺 ¿Cuáles hemos usado hasta ahora?

👁 Observaciones

• Como usarios que consumimos datos geoespaciales producidos por otros, usualmente podemos trabajar con la proyección en la que vienen los datos, es decir, muchas veces no tenemos que preocuparnos por el detalle de la proyección. Como vimos, pudimos hacer los primeros mapas sin saber de esto.

• Hay algunas situaciones en que si es necesario transformar nuestros datos a otro sistema de coordenadas, por ejemplo, si queremos trabajar con conjuntos de datos con distinto CRS, o para utilizar ciertas herramientas como leaflet.

📍 Interactivos

Para hacer mapas interactivos usaremos el paquete leaflet (Cheng, Karambelkar, y Xie 2021):

library(leaflet)

Leaflet es una librería de JavaScript para realizar mapas interactivos, el paquete de R nos permite utilizarla desde R.

🖼 Puntos culturales con Leaflet

Leaflet recibe datos geoespaciales con CRS 4326 (al igual que las coordenadas que copiamos de GoogleMaps). Por ellos si queremos graficar los puntos culturales debemos comenzar reproyectando para que leaflet funcione como debe.

La función para reproyectar recibe un objeto de clase sf, y la referencia a la proyección que queremos obtener, en este caso vamos a reproyectar cultura_mv a la proyección asociada al código ESPG 4326:

cultura_mv_longlat <- st_transform(cultura_mv, 4326)

Ahora graficamos con leaflet, la sintaxis de leaflet es similar a la de ggplot, pero usamos %>% en lugar de +:

leaflet(data = cultura_mv_longlat) %>% 
  addTiles() %>% 
  addCircles()

La función addTiles() agrega las tejas del fondo. En este ejemplo agregagmos marcadores que muestran la ubicación de los sitios culturales de Montevideo, podemos añadir información de variables en los datos, por ejemplo, el nombre de la atracción. Y con addCircles() agregamos la geometría de puntos representados con círculos.

leaflet(data = cultura_mv_longlat) %>% 
  addTiles() %>% 
  addMarkers(popup = ~NOMBRE)

Con leaflet hay mucha flexibilidad, podemos meter html a las etiquetas o definir nuestros propios íconos, en el sitio se explica claramente con ejemplos.

rLadiesIcon <- makeIcon(
  iconUrl = "https://raw.githubusercontent.com/rladies/starter-kit/master/logo/R-LadiesGlobal_CMYK_offline_LogoOnly.svg",
  iconWidth = 38, iconHeight = 38
)

leaflet(data = colab_longlat) %>% 
  addTiles() %>% 
  addMarkers(popup = ~nombre, icon = rLadiesIcon) 

También podemos agregar objetos espaciales de otros tipos, para polígonos se utiliza addPolygons().

barrios_longlat <- barrios_mv %>% 
  st_transform(crs = 4326)
leaflet(data = barrios_longlat) %>% 
  addTiles() %>% 
  addPolygons() 

Tu turno

🌍 ¿Qué ocurre si utilizas leaflet con los datos no transformados, es decir si utilizas colab en lugar de colab_longlat?

🌍 Agrega al siguiente código la ubicación de los puntos de interés de Montevideo.

leaflet(data = colab_longlat) %>% 
  addTiles() %>% 
  addMarkers(popup = ~nombre, icon = rLadiesIcon) 

🌎 Experimenta con otros fondos, utiliza la función addProviderTiles() y elige un fondo de aquí.

leaflet(data = colab) %>% 
  addProviderTiles(providers$Stamen.Watercolor) %>% 
  addMarkers(popup = ~nombre, icon = rLadiesIcon) 

➕ Extra: Operaciones

Las operaciones con datos geoespaciales van más allá del alcance de este taller, pero aquí hay algunos ejemplos de lo que se puede hacer.

• Área de polígonos (si no incluyera en los datos).

barrios_mv %>% 
  mutate(area = st_area(geometry))

Simple feature collection with 62 features and 4 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 551982.7 ymin: 6133499 xmax: 589227.2 ymax: 6159810
Projected CRS: WGS 84 / UTM zone 21S
# A tibble: 62 x 5
   AREA_KM NOMBBARR   NROBARRIO                      geometry     area
 *   <dbl> <chr>          <int>            <MULTIPOLYGON [m]>    [m^2]
 1   2.11  Ciudad Vi…         1 (((573146.3 6137145, 573148.… 2297414…
 2   1.30  Centro             2 (((574497.2 6137037, 574506.… 1294508…
 3   0.694 Barrio Sur         3 (((574564 6136237, 574570.5 …  726727…
 4   2.28  Cordon             4 (((574397.4 6137451, 574382.… 2278114…
 5   0.798 Palermo            5 (((575149.1 6136881, 575144 …  832445…
 6   0.766 Parque Ro…         6 (((576530.7 6135354, 576532.…  778529…
 7   2.73  Punta Car…         7 (((576200.9 6135810, 576200.… 2818381…
 8   4.16  Buceo              9 (((580161.9 6137234, 580167.… 4164074…
 9   3.38  Pque. Bat…        10 (((578338.6 6137972, 578352.… 3371317…
10   3.50  Malvin            11 (((580642.4 6138757, 580590.… 3502102…
# … with 52 more rows

• En que barrio están los museos?

barrios_cultura_mv <- barrios_mv %>%
    st_join(cultura_mv, join = st_intersects) # usualmente st_within es más apropiado 

barrios_cultura_mv %>% 
  st_drop_geometry() %>% 
  group_by(NOMBBARR) %>% 
  summarise(n_lugares = sum(!is.na(NOMBRE))) %>% 
  arrange(-n_lugares) 

# A tibble: 62 x 2
   NOMBBARR                 n_lugares
   <chr>                        <int>
 1 Centro                          81
 2 Cordon                          57
 3 Ciudad Vieja                    50
 4 Pque. Batlle, V. Dolores        24
 5 Buceo                           14
 6 Palermo                         13
 7 Punta Carretas                  13
 8 Parque Rodo                     10
 9 Pocitos                         10
10 Aguada                           8
# … with 52 more rows

• Distancia entre puntos.

colab_longlat %>% 
  st_distance()

Units: [m]
             [,1]         [,2]       [,3]       [,4]         [,5]
 [1,]       0.000 7.502629e+06 7448775.05 7510874.54 7.502652e+06
 [2,] 7502629.366 0.000000e+00  383282.07   31536.82 2.602528e+01
 [3,] 7448775.047 3.832821e+05       0.00  354676.44 3.832982e+05
 [4,] 7510874.536 3.153682e+04  354676.44       0.00 3.154332e+04
 [5,] 7502652.006 2.602528e+01  383298.22   31543.32 0.000000e+00
 [6,] 7257670.754 3.373620e+05  239764.86  323808.85 3.373873e+05
 [7,] 7497463.328 5.432652e+03  380862.59   31727.21 5.458155e+03
 [8,]    1764.185 7.503542e+06 7449783.57 7511795.22 7.503565e+06
 [9,] 7500443.235 2.585246e+03  381589.35   30873.05 2.611242e+03
[10,] 7530147.568 1.193379e+05  276168.20   87826.71 1.193452e+05
[11,] 7498177.981 4.507945e+03  383331.81   33635.88 4.528303e+03
[12,] 3152959.565 6.068736e+06 5838781.74 6061168.81 6.068762e+06
[13,] 7505126.464 9.214092e+03  374916.81   22323.46 9.220354e+03
[14,] 7497823.883 4.831816e+03  383075.98   33578.93 4.853018e+03
[15,] 7610373.070 1.842366e+05  282417.27  155071.53 1.842339e+05
[16,] 7499714.033 2.987428e+03  382202.12   31806.81 3.012334e+03
[17,] 7436051.589 3.391470e+05   48394.74  311388.81 3.391643e+05
           [,6]        [,7]        [,8]        [,9]      [,10]
 [1,] 7257670.8 7497463.328    1764.185 7500443.235 7530147.57
 [2,]  337362.0    5432.652 7503542.253    2585.246  119337.87
 [3,]  239764.9  380862.590 7449783.571  381589.350  276168.20
 [4,]  323808.9   31727.212 7511795.225   30873.046   87826.71
 [5,]  337387.3    5458.155 7503564.890    2611.242  119345.23
 [6,]       0.0  332445.133 7258642.965  334822.361  296539.00
 [7,]  332445.1       0.000 7498376.647    2995.084  118996.53
 [8,] 7258643.0 7498376.647       0.000 7501356.476 7531090.06
 [9,]  334822.4    2995.084 7501356.476       0.000  118515.28
[10,]  296539.0  118996.534 7531090.064  118515.284       0.00
[11,]  334610.3    2496.943 7499090.685    3083.189  121113.16
[12,] 5741659.9 6063404.752 3154722.243 6066155.172 6036651.41
[13,]  333238.7   10505.787 7506041.627    8832.392  110133.98
[14,]  334210.4    2213.971 7498736.641    3226.223  121005.11
[15,]  362697.1  185941.907 7611323.898  184401.518   86769.67
[16,]  334788.2    2474.609 7500627.088    1030.070  119401.56
[17,]  203695.9  336466.294 7437048.540  337349.621  236505.79
             [,11]   [,12]       [,13]        [,14]      [,15]
 [1,] 7498177.9805 3152960 7505126.464 7497823.8827 7610373.07
 [2,]    4507.9447 6068736    9214.092    4831.8162  184236.59
 [3,]  383331.8147 5838782  374916.809  383075.9817  282417.27
 [4,]   33635.8807 6061169   22323.458   33578.9342  155071.53
 [5,]    4528.3030 6068762    9220.354    4853.0180  184233.88
 [6,]  334610.3284 5741660  333238.652  334210.4368  362697.07
 [7,]    2496.9431 6063405   10505.787    2213.9710  185941.91
 [8,] 7499090.6852 3154722 7506041.627 7498736.6408 7611323.90
 [9,]    3083.1890 6066155    8832.392    3226.2230  184401.52
[10,]  121113.1580 6036651  110133.983  121005.1084   86769.67
[11,]       0.0000 6065193   11834.521     410.6946  187429.35
[12,] 6065192.6107       0 6066609.694 6064782.0282 6091244.55
[13,]   11834.5210 6066610       0.000   11880.9217  175595.43
[14,]     410.6946 6064782   11880.922       0.0000  187473.14
[15,]  187429.3480 6091245  175595.426  187473.1380       0.00
[16,]    2054.2797 6065872    9838.574    2215.6090  185419.68
[17,]  338951.8755 5849919  331013.061  338679.2153  253887.03
            [,16]      [,17]
 [1,] 7499714.033 7436051.59
 [2,]    2987.428  339147.04
 [3,]  382202.119   48394.74
 [4,]   31806.814  311388.81
 [5,]    3012.334  339164.28
 [6,]  334788.237  203695.95
 [7,]    2474.609  336466.29
 [8,] 7500627.088 7437048.54
 [9,]    1030.070  337349.62
[10,]  119401.564  236505.79
[11,]    2054.280  338951.88
[12,] 6065871.627 5849918.57
[13,]    9838.574  331013.06
[14,]    2215.609  338679.22
[15,]  185419.682  253887.03
[16,]       0.000  337916.36
[17,]  337916.355       0.00

La cheatsheet del paquete sf describe operaciones geométricas disponibles, e incluye útiles esquemas.

👋 Resumen

• Hay muchas herramientas para trabajar con datos geoespaciales en R, en este taller introdujimos 2 paquetes considerando visualización: sf, leaflet.

• Con sf podemos cargar datos geoespaciales en R, vimos como leer datos en formato shapefile con la función st_read() y graficar con ggplot usando la geometría geom_sf().

• Para utilizar leaflet tenemos que reproyectar los datos al CRS 4326 usando st_transform() y con ello estamos listos para usar leaflet, que tiene una sintaxis similar a ggplot2.

leaflet(data = mi_capa_longlat) %>% 
  addTiles() %>% 
  addCircles() 

• Los detalles de que conforman un shapefile y las CRS dan contexto pero podemos empezar a hacer mapas mientras seguimos aprendiendo.

Licencia

Texto y gráficas con licencia Creative Commons Attribution CC BY 4.0, a menos que se indique lo contrario o tengan otra fuente (indicado en la nota). Código en https://github.com/tereom/taller-mapas-mv.