class: center, middle, inverse, title-slide # Datos espaciales con R ## Principios básicos ### Stephanie Orellana ### R Ladies Chile ### 09-05-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** <img src="https://jules32.github.io/useR-2019-keynote/img/horst-welcome_to_rstats_twitter.png" width="40%" style="display: block; margin: auto 0 auto auto;" /> --- # 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 Existen algunas diferencias entre R (izquierda) y los GIS tradicionales (derecha): .pull-left[ 1. Código ilimitado, se pueden mezclar diferentes paquetes para diferentes operaciones y análisis (incluso puedes hacer tus reportes y presentaciones en R). 2. Visualización de datos espaciales no es tan detallada como en un SIG tradicional, aunque existen paquetes para explorar de forma dinámica. 3. Se puede hacer todo tipo de análisis y gráficos de la información espacial (ej: `ggplot2`). ] .pull-right[ 1. Algunos como ArcGis y QGis tienen una consola de código en Python. 2. Visualización rápida de grandes vectores y creación de estos con herramientas manuales especializadas. 3. Sólo vistas de mapas e histogramas. ] --- # R vs GIS tradicionales Existen algunas diferencias entre R (izquierda) y los GIS tradicionales (derecha): .pull-left[ <ol start=4> <li> Se pueden programar procesos iterativos para hacer la misma acción muchas veces (se puede programar todo). <br></br> <li> Fácil análisis de las tablas de atributos (se puede aplicar dplyr y R4DS). <br></br> <li> Permite reproducibilidad, queda todo por escrito. </ol> ] .pull-right[ <ol start=4> <li> Procesos en "batch" siguen siendo manuales. <br></br> <li> Análisis limitado de tablas de atributos. <br></br> <li> Poco reproducible, nada (o muy poco) queda por escrito. </ol> ] --- # 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;" /> - 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;" /> - 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 ] --- # Diferentes tipos, diferentes paquetes. ### 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** <img src="https://user-images.githubusercontent.com/520851/50280460-e35c1880-044c-11e9-9ed7-cc46754e49db.jpg" width="55%" style="display: block; margin: auto;" /> --- # Cargar datos con sf ```r library(sf) # Cargar datos ## Poligonos comunas <- read_sf("data/comunas.shp") comunas ``` ``` ## Simple feature collection with 346 features and 11 fields ## geometry type: MULTIPOLYGON ## dimension: XY ## bbox: xmin: -12184470 ymin: -7554436 xmax: -7393642 ymax: -1978920 ## CRS: 3857 ## # A tibble: 346 x 12 ## objectid shape_leng dis_elec cir_sena cod_comuna codregion st_area_sh ## <dbl> <dbl> <int> <int> <int> <int> <dbl> ## 1 48 170039. 16 8 6204 6 9.69e8 ## 2 29 125730. 15 8 6102 6 4.16e8 ## 3 30 63026. 15 8 6103 6 1.45e8 ## 4 31 89841. 15 8 6104 6 3.26e8 ## 5 78 122626. 23 11 9121 9 6.99e8 ## 6 79 279936. 23 11 9103 9 3.13e9 ## 7 81 213152. 23 11 9105 9 1.45e9 ## 8 82 142386. 22 11 9106 9 9.26e8 ## 9 34 56157. 15 8 6106 6 1.65e8 ## 10 37 109701. 15 8 6109 6 3.24e8 ## # ... with 336 more rows, and 5 more variables: st_length_ <dbl>, Region <chr>, ## # Comuna <chr>, Provincia <chr>, geometry <MULTIPOLYGON [m]> ``` --- # Cargar datos con sf ```r ## Puntos observatorios <- read_sf("data/observatorios.shp") observatorios ``` ``` ## Simple feature collection with 11 features and 10 fields ## geometry type: POINT ## dimension: XY ## bbox: xmin: 301486.6 ymin: 6546377 xmax: 364430 ymax: 6785395 ## projected CRS: SIRGAS 2000 / UTM zone 19S ## # A tibble: 11 x 11 ## OBJECTID_1 OBJECTID_2 OBJECTID NOMBRE DESCRIP Descp POINT_X POINT_Y COMUNA ## <dbl> <dbl> <dbl> <chr> <chr> <chr> <dbl> <dbl> <chr> ## 1 1 1 1 Cerro~ Observ~ Cien~ 331767. 6.76e6 La Hi~ ## 2 2 2 2 Cerro~ Observ~ Cien~ 326033. 6.66e6 Vicuña ## 3 3 3 3 Cerro~ Observ~ Cien~ 335575. 6.79e6 La Hi~ ## 4 4 4 4 Collo~ Observ~ Turi~ 301487. 6.65e6 Andac~ ## 5 5 5 5 Cerro~ Observ~ Cien~ 332896 6.65e6 Vicuña ## 6 6 6 6 Cerro~ Observ~ Turi~ 337459. 6.68e6 Vicuña ## 7 7 7 7 Mayu Observ~ Turi~ 304778 6.68e6 La Se~ ## 8 8 8 8 Cruz ~ Observ~ Turi~ 309392 6.55e6 Comba~ ## 9 9 9 9 del P~ Observ~ Turi~ 336805 6.66e6 Vicuña ## 10 10 10 10 Canca~ Observ~ Turi~ 364430 6.66e6 Paihu~ ## 11 11 11 11 Hacie~ Observ~ Turi~ 327489. 6.62e6 Monte~ ## # ... with 2 more variables: PROVINCIA <chr>, geometry <POINT [m]> ``` --- # Diferentes tipos, diferentes paquetes. ### 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 matrix 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/lst_01.tif") lst_01 ``` ``` ## class : RasterLayer ## dimensions : 266, 295, 78470 (nrow, ncol, ncell) ## resolution : 0.008983153, 0.008983153 (x, y) ## extent : -72.12573, -69.4757, -34.33361, -31.94409 (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_basic/data/raster/LST/lst_01.tif ## names : lst_01 ``` --- # Cargar datos con raster ```r library(raster) # Cargar datos ## Grupos de rásters l <- list.files("data/raster/LST/", full.names = T) lst <- stack(l) lst ``` ``` ## class : RasterStack ## dimensions : 266, 295, 78470, 12 (nrow, ncol, ncell, nlayers) ## resolution : 0.008983153, 0.008983153 (x, y) ## extent : -72.12573, -69.4757, -34.33361, -31.94409 (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(comunas)$epsg st_crs(comunas)$proj4string ``` ``` ## [1] 3857 ## [1] "+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 +wktext +no_defs" ``` -- ```r ## Cambiar a UTM comunas_utm <- st_transform(comunas, crs = 32719) st_crs(comunas_utm)$epsg st_crs(comunas_utm)$proj4string ``` ``` ## [1] 32719 ## [1] "+proj=utm +zone=19 +south +datum=WGS84 +units=m +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_utm <- projectRaster(lst_01, crs = crs("+init=epsg:32719")) crs(lst_01_utm) ``` ``` ## CRS arguments: ## +init=epsg:32719 +proj=utm +zone=19 +south +datum=WGS84 +units=m ## +no_defs +ellps=WGS84 +towgs84=0,0,0 ``` --- # 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 : -72.12573 ## xmax : -69.4757 ## ymin : -34.33361 ## ymax : -31.94409 ``` ```r extent(lst_01_utm) ``` ``` ## class : Extent ## xmin : 199485.7 ## xmax : 460941.7 ## ymin : 6191763 ## ymax : 6470643 ``` ] -- ```r res(lst_01) ``` ``` ## [1] 0.008983153 0.008983153 ``` ```r res(lst_01_utm) ``` ``` ## [1] 838 996 ``` --- class:: middle # Visualización --- # Cortar Vectores por atributos y visualizar ```r # Seleccionar sólo comunas en la región de Coquimbo comunas_coquimbo <- comunas %>% filter(Region == "Región de Coquimbo") ``` -- ```r ## Gráfico rápido plot(comunas_coquimbo) ``` -- <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(comunas_coquimbo)+ geom_sf(aes(fill=Comuna)) ``` -- <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 comunas_valparaiso_utm <- comunas %>% filter(Region == "Región de Valparaíso") %>% st_transform( crs = 32719) m_lst_01 <- mask(lst_01_utm, comunas_valparaiso_utm) ``` ```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/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;" />