class: center, middle, inverse, title-slide # Datos espaciales con R ## Principios básicos ### Stephanie Orellana ### R Ladies Chile para RLadies La Paz ### 25-06-2020 --- # Hoy hablaremos de: - Introducción a R - Conceptos básicos de datos espaciales - Paquete sf - Paquete raster - Ejemplos prácticos --- # Sobre mí .right-column2[ <img src="img/ai.png" width="80%" style="display: block; margin: auto 0 auto auto;" /> <img src="img/ai2.png" width="60%" style="display: block; margin: auto;" /> ] -- - Stephanie Orellana (sporella@uc.cl) -- - Ingeniera Agrónoma y Magíster en Recursos Naturales -- - No aprendí a programar en ningún curso de la Universidad, sino que por necesidad. -- - Partí trabajando con modelos mixtos de variables climáticas espacializadas, usando: + Arcgis + Excel + R -- - Luego fui migrando completamente a R - Pero para "mirar" los rásters y vectores uso QGIS --- # Básicos de R -- - R es un lenguaje de **programación** y es un software libre -- - RStudio es la interfaz de usuario más conocida (hace que R sea amigable e incluye herramientas para facilitar el trabajo) -- - R es un lenguaje orientado a **objetos** en donde pueden aplicarse **funciones** a estos objetos. <img src="https://conceptosclaros.com/wp-content/uploads/2017/06/r-rstudio.jpg" width="80%" style="display: block; margin: auto 0 auto auto;" /> --- # Por qué usar R -- - Es un software libre y gratuito -- - Permite hacer investigación reproducible: -- - El flujo de trabajo queda escrito en código - Los análisis se hacen dentro del mismo programa -- - Es una herramienta potente y está en constante desarrollo de nuevas funcionalidades. -- - Existe una gran comunidad alrededor del mundo que siempre está dispuesta a enseñar, responder preguntas y desarrollar nuevos paquetes. -- - Aprender a **buscar en Google** es primordial - **Twitter** <div class="figure" style="text-align: right"> <img src="https://jules32.github.io/useR-2019-keynote/img/horst-welcome_to_rstats_twitter.png" alt="<a href='https://twitter.com/allison_horst'> @Allison Horst </a>" width="35%" /> <p class="caption"><a href='https://twitter.com/allison_horst'> @Allison Horst </a></p> </div> --- # R Base v/s tidyverse .pull-left[ <img src="https://i.ytimg.com/vi/MCB4gNVfysU/maxresdefault.jpg" width="90%" style="display: block; margin: auto;" /> - Programar con la sintaxis nativa de R - Indices de objetos se realizan con: ` $ , [ ], [,] o [[]]` - Funciones incluidas en el paquete base (base::) ] .pull-right[ <img src="https://miro.medium.com/max/4032/1*B-cwhqnFgGIbd9lWnzi_mQ.png" width="90%" style="display: block; margin: auto;" /> - Grupo de paquetes diseñados para hacer ciencia de datos - Busca una programación más intuitiva. Funciones como verbos. - Todos los paquetes comparten una filosofía de diseño, gramática y estructuras de datos subyacentes. - Operador "pipe" %>% ] --- # R Base v/s tidyverse ```r # Seleccionar columnas ## R base iris$Sepal.Length iris[, c("Petal.Length", "Sepal.Length", "Species")] # tidyverse library(dplyr) iris %>% select(Petal.Length, Sepal.Length, Species) ``` <img src="https://image.slidesharecdn.com/emilytidyversepresentation-180205195912/95/the-lesser-known-stars-of-the-tidyverse-5-638.jpg?cb=1517860826" width="60%" style="display: block; margin: auto;" /> -- .pull-center[ ¡Para procesar datos espaciales necesitaremos las dos formas! <img src="https://cdn160.picsart.com/upscale-235481347027212.png?type=webp&to=min&r=1024" width="25%" style="display: block; margin: auto;" /> ] --- class: middle # ¡Vamos a la parte espacial! <img src="https://encrypted-tbn0.gstatic.com/images?q=tbn%3AANd9GcRtbePUpnq_UKJKFI_yW88HotGcL3SHOdrXAZJCNtRyyNsUKss1&usqp=CAU" width="30%" style="display: block; margin: auto 0 auto auto;" /> --- # R vs GIS tradicionales |R | SIG Tradicionales --------|---------|--------- 1 | Código ilimitado, se pueden mezclar diferentes paquetes para diferentes operaciones y análisis (incluso puedes hacer tus reportes y presentaciones en R). | Algunos como ArcGis y QGis tienen una consola de código en Python. -- 2 | Visualización de datos espaciales no es tan detallada como en un SIG tradicional, aunque existen paquetes para explorar de forma dinámica. | Visualización rápida de grandes vectores y creación de estos con herramientas manuales especializadas. -- 3 | Se puede hacer todo tipo de análisis y gráficos de la información espacial (ej: `ggplot2`). | Sólo vistas de mapas e histogramas. -- 4 | Se pueden programar procesos iterativos para hacer la misma acción muchas veces (se puede programar todo). | Procesos en "batch" siguen siendo manuales. -- 5 | Fácil análisis de las tablas de atributos (se puede aplicar dplyr y contenidos del libro R for Data Science). | Análisis limitado de tablas de atributos. -- 6| Permite reproducibilidad, queda todo por escrito. | Poco reproducible, nada (o muy poco) queda por escrito. --- # R spatial [CRAN TASK VIEW](https://cran.r-project.org/web/views/Spatial.html) <iframe src="https://cran.r-project.org/web/views/Spatial.html" width="90%" height="400px" style="border: none;"></iframe> --- # Paquetes más descargados de #rspatial <img src="img/sppackages_stats.png" width="90%" style="display: block; margin: auto;" /> .footnote[ [*] 50 paquetes más descargados con la etiqueta RSpatial desde 2015 ] --- class: middle # Conceptos básicos a tomar en cuenta <img src="img/abolins180400032.jpg" width="30%" style="display: block; margin: auto 0 auto auto;" /> --- # 1) Tipo de dato/archivo .pull-left[ .h2[VECTORES] .p-caption[ <div class="figure" style="text-align: center"> <img src="img/vector.png" alt="Fuente: <a href='http://katiejolly.io/rnorth-19'> @Katie Jolly </a>" width="55%" /> <p class="caption">Fuente: <a href='http://katiejolly.io/rnorth-19'> @Katie Jolly </a></p> </div> ] <img src="img/vectorcat.png" width="40%" style="display: block; margin: auto;" /> .small[ - Mediciones de campo - Estaciones meteorológicas - Caminos - Rutas GPS - Comunas - Regiones - Áreas de estudio ] ] .pull-right[ .h2[RÁSTER] .p-caption[ <div class="figure" style="text-align: center"> <img src="img/raster.png" alt="Fuente: <a href='http://katiejolly.io/rnorth-19'> @Katie Jolly </a>" width="62%" /> <p class="caption">Fuente: <a href='http://katiejolly.io/rnorth-19'> @Katie Jolly </a></p> </div> ] <img src="img/pixelcat.png" width="40%" style="display: block; margin: auto;" /> .small[ - Imágenes satelitales - Imágenes aéreas - Modelos de elevación digital - Productos de tipos de usos de suelo - Datos (numéricos o categóricos) en grillas regulares ] ] --- # Vectores: paquete "sf" -- - Es el paquete "más moderno" para trabajar con vectores. -- - Actúa bajo la lógica de "simple features" en donde un objeto espacial se divide en dos componentes: - data.frame (tabla de atributos) - geometría **geom** (información espacial) -- - Es compatible con **tidyverse** <div class="figure" style="text-align: right"> <img src="https://user-images.githubusercontent.com/520851/50280460-e35c1880-044c-11e9-9ed7-cc46754e49db.jpg" alt="<a href='https://twitter.com/allison_horst'> @Allison Horst </a>" width="45%" /> <p class="caption"><a href='https://twitter.com/allison_horst'> @Allison Horst </a></p> </div> --- # Cargar datos con sf ```r library(sf) # Cargar datos ## Poligonos *provincias <- read_sf("data/provincias/provincias.shp") provincias ``` ``` ## Simple feature collection with 123 features and 5 fields ## geometry type: MULTIPOLYGON ## dimension: XY ## bbox: xmin: -69.64483 ymin: -22.90657 xmax: -57.45443 ymax: -9.669633 ## geographic CRS: WGS 84 ## # A tibble: 123 x 6 ## gml_id COD_DEP NOM_DEP COD_PROV NOM_PROV geometry ## <chr> <chr> <chr> <chr> <chr> <MULTIPOLYGON [°]> ## 1 provinc~ <NA> Potosi 14 Daniel C~ (((-67.77305 -19.95029, -67.7731~ ## 2 provinc~ 0 LAGO <NA> Lago Poo~ (((-67.06019 -18.44845, -67.0602~ ## 3 provinc~ 0 LAGO <NA> Lago Tit~ (((-69.32165 -15.54775, -69.3212~ ## 4 provinc~ 0 LAGO <NA> Lago Uru~ (((-67.10087 -18.02564, -67.1007~ ## 5 provinc~ 0 SALAR <NA> Salar de~ (((-67.98348 -19.03682, -67.9834~ ## 6 provinc~ 0 SALAR 0 Salar de~ (((-67.48285 -19.53102, -67.4830~ ## 7 provinc~ 01 Chuqui~ 07 Nor Cinti (((-64.71348 -19.7808, -64.71298~ ## 8 provinc~ 01 Chuqui~ 09 Sud Cinti (((-64.26963 -20.52255, -64.2672~ ## 9 provinc~ 02 La Paz 02 Omasuyos (((-68.53956 -16.10705, -68.5297~ ## 10 provinc~ 02 La Paz 12 Los Andes (((-68.39181 -16.01124, -68.3901~ ## # ... with 113 more rows ``` --- # Cargar datos con sf ```r ## Puntos *aeropuetos <- read_sf("data/aero1/aero1.shp") aeropuetos ``` ``` ## Simple feature collection with 33 features and 4 fields ## geometry type: POINT ## dimension: XY ## bbox: xmin: -69.40331 ymin: -22.6977 xmax: -57.84501 ymax: -10.82314 ## geographic CRS: WGS 84 ## # A tibble: 33 x 5 ## gml_id ID AERO_INTER AEROPUERTO geometry ## <chr> <dbl> <chr> <chr> <POINT [°]> ## 1 aero1.1 0 INTER VIRU VIRU (-63.10903 -17.80115) ## 2 aero1.2 0 NAL ROBORE (-59.69179 -18.31267) ## 3 aero1.3 0 NAL SALVADOR OGAYA (Pto. Suarez) (-57.84501 -18.93599) ## 4 aero1.4 0 NAL YACUIBA (-63.63405 -21.97789) ## 5 aero1.5 0 NAL CAMIRI (-63.46358 -20.04498) ## 6 aero1.6 0 INTER JORGE WILSTERMAN (-66.12517 -17.40998) ## 7 aero1.7 0 INTER J. .F. KENNEDY (-68.22983 -16.51672) ## 8 aero1.8 0 NAL <NA> (-68.09211 -13.74058) ## 9 aero1.9 0 NAL RURRENABAQUE (-67.44159 -14.43207) ## 10 aero1.10 0 NAL APOLO (-68.38563 -14.712) ## # ... with 23 more rows ``` --- # Rásters: paquete "raster" -- - Es el paquete más utilizado para el procesamiento de rásters, aunque se está desarrollando el paquete `stars` (del mismo creador de `sf`) el cual permite el procesamiento de ráster en tidyverse -- - El paquete ráster actúa bajo la lógica de una matriz de datos en donde el inicio se encuentra en el extremo superior izquierdo. -- - Si sabes cómo procesar matrices con R Base, es fácil acceder a los valores de un ráster. -- - El paquete tiene funciones para leer, transformar, analizar y guardar rásters. <img src="https://i.stack.imgur.com/xeqec.png" width="30%" style="display: block; margin: auto;" /> --- # Cargar datos con raster ```r library(raster) # Cargar datos ## Ráster individual *lst_01 <- raster("data/raster/lst_01.tif") lst_01 ``` ``` ## class : RasterLayer ## dimensions : 1546, 1420, 2195320 (nrow, ncol, ncell) ## resolution : 0.008983153, 0.008983153 (x, y) ## extent : -69.87096, -57.11489, -23.28433, -9.396378 (xmin, xmax, ymin, ymax) ## crs : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 ## source : C:/Users/PEZ/Desktop/git/datos_espaciales_lapaz/data/raster/lst_01.tif ## names : lst_01 ``` --- # Cargar datos con raster ```r library(raster) # Cargar datos ## Grupos de rásters l <- list.files("data/raster/", full.names = T) *lst <- stack(l) lst ``` ``` ## class : RasterStack ## dimensions : 1546, 1420, 2195320, 12 (nrow, ncol, ncell, nlayers) ## resolution : 0.008983153, 0.008983153 (x, y) ## extent : -69.87096, -57.11489, -23.28433, -9.396378 (xmin, xmax, ymin, ymax) ## crs : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 ## names : lst_01, lst_02, lst_03, lst_04, lst_05, lst_06, lst_07, lst_08, lst_09, lst_10, lst_11, lst_12 ``` --- # 2) Proyección Espacial (CRS) - Es fundamental conocer en qué proyección espacial se encuentran nuestros datos. Las proyecciones están determinadas por un CRS (coordinate reference system). -- - Cada CRS puede estar definido por un EPSG (sigla de European Petroleum Survey Group) el cual es un código que se encuentra registrado en un listado de parámetros geodésicos (https://spatialreference.org/ref/) -- - Si quiero hacer operaciones espaciales, todas las capas deben tener el mismo sistema de proyección. -- - En Chile, los sistemas de proyección más usados son: - EPSG:4326: WGS 84 - EPSG:32719: WGS 84 / UTM zone 19S - EPSG:3857: WGS84 Web (Pseudo)Mercator (Auxiliary Sphere) - EPSG:24879: PSAD56 / UTM zone 19S --- # 2) Proyección espacial <iframe src="img/proj_face1.mp4?showcase=0" width="100%" height="400px"></iframe> [Andy Woodruff: Projection Face!](http://bl.ocks.org/awoodruff/9216081?utm_campaign) --- # Cambiar proyección de un vector ```r # Cambiar proyección ## Proyección actual st_crs(provincias)$epsg st_crs(provincias)$proj4string ``` ``` ## [1] 4326 ## [1] "+proj=longlat +datum=WGS84 +no_defs" ``` -- ```r ## Cambiar a Pseudo Mercator *provincias_psm <- st_transform(provincias, crs = 3857) st_crs(provincias_psm)$epsg st_crs(provincias_psm)$proj4string ``` ``` ## [1] 3857 ## [1] "+proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs" ``` --- # Cambiar proyección de un ráster ```r # Conocer CRS actual crs(lst_01) ``` ``` ## CRS arguments: ## +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 ``` ```r # Proyectar a UTM *lst_01_psm <- projectRaster(lst_01, crs = crs("+init=epsg:3857")) crs(lst_01_psm) ``` ``` ## CRS arguments: ## +init=epsg:3857 +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 ## +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs ``` --- # 3) Extensión y resolución espacial -- - Muy importante cuando trabajo con rásters. -- - Si mis rásters no tienen la misma extensión (bbox) espacial, no pueden tratarse como stacks. - Tampoco si es que tienen diferente resolución -- .pull-left[ ```r extent(lst_01) ``` ``` ## class : Extent ## xmin : -69.87096 ## xmax : -57.11489 ## ymin : -23.28433 ## ymax : -9.396378 ``` ```r extent(lst_01_psm) ``` ``` ## class : Extent ## xmin : -7783000 ## xmax : -6353000 ## ymin : -2672081 ## ymax : -1045521 ``` ] -- .pull-right[ ```r res(lst_01) ``` ``` ## [1] 0.008983153 0.008983153 ``` ```r res(lst_01_psm) ``` ``` ## [1] 1000 1040 ``` ] --- class:: middle # Visualización --- # Cortar Vectores por atributos y visualizar ```r library(sf) library(tidyverse) # Seleccionar sólo comunas en el departamento de Cochabamba provincias_coch <- provincias %>% * filter(NOM_DEP == "Cochabamba") ``` -- ```r ## Gráfico rápido plot(provincias_coch) ``` -- <img src="presentacion_files/figure-html/unnamed-chunk-32-1.png" style="display: block; margin: auto;" /> --- # Vectores y ggplot2 ```r ## Gráfico con ggplot library(ggplot2) ggplot(provincias_coch) + * geom_sf(aes(fill = NOM_PROV)) ``` -- <img src="presentacion_files/figure-html/unnamed-chunk-34-1.png" style="display: block; margin: auto;" /> --- # Cortar rásters por máscaras y visualizar ```r # Creo máscara provincias_lapaz_psm <- provincias %>% filter(NOM_DEP == "La Paz") %>% st_transform(crs = 3857) *m_lst_01 <- mask(lst_01_psm, provincias_lapaz_psm) ``` ```r plot(m_lst_01) ``` <img src="presentacion_files/figure-html/unnamed-chunk-36-1.png" width="40%" style="display: block; margin: auto;" /> --- # Rásters y ggplot2 ```r # Convertir a tabla rt <- data.frame(rasterToPoints(m_lst_01)) ggplot() + geom_raster(data = rt , aes(x = x, y = y, fill=lst_01))+ scale_fill_gradientn(colours=rev(terrain.colors(15)))+ coord_equal() ``` <img src="presentacion_files/figure-html/unnamed-chunk-37-1.png" width="40%" style="display: block; margin: auto;" /> --- # leaflet
--- class: middle, invert # Recomendaciones --- # Libros .pull-left[ <div class="figure" style="text-align: center"> <img src="https://geocompr.robinlovelace.net/images/cover.png" alt=" <a href='https://geocompr.robinlovelace.net/'>Geocomputation with R</a>" /> <p class="caption"> <a href='https://geocompr.robinlovelace.net/'>Geocomputation with R</a></p> </div> ] .pull-right[ <div class="figure" style="text-align: center"> <img src="https://d33wubrfki0l68.cloudfront.net/b88ef926a004b0fce72b2526b0b5c4413666a4cb/24a30/cover.png" alt=" <a href='https://r4ds.had.co.nz/'>R for data science</a>" /> <p class="caption"> <a href='https://r4ds.had.co.nz/'>R for data science</a></p> </div> ] --- # A quién seguir <img src="https://ohmygeek.net/wp-content/uploads/2012/06/twitter-logo-jun-2012.jpg" width="20%" style="display: block; margin: auto 0 auto auto;" /> - https://twitter.com/CivicAngela `#rspatialchat` - https://twitter.com/Paula_Moraga_ - https://twitter.com/hadleywickham - https://twitter.com/edzerpebesma - https://twitter.com/RogerBivand - https://twitter.com/jakub_nowosad - `#rspatial` - `#rstats` --- # Recomendaciones generales -- - Si se demoró mucho en cargar, mejor no visualizar. -- - Cuando se trabaja con datos muy grandes, mejor hacer todas las pruebas con una selección pequeña. -- - R tampoco es mágico, si le damos las intrucciones incorrectas, no obtendremos el resultado esperado. -- - R es rápido, pero hay procesos que se demoran. Esperar. <img src="https://image.freepik.com/vector-gratis/lego-isometrico_71884-8.jpg" width="40%" style="display: block; margin: auto;" />