Sistemas de Información Geográfica con R

Introducción

knitr::opts_chunk$set(fig.align = "center",
                      warning = FALSE,
                      message = FALSE,
                      fig.width = 9)

• Notas tomadas del curso “Sistemas de Información Geográfica en R (SIG en R)” con Derek Corcoran.

Generalidades

• Datos
  • Datos de captura remota (satelites lidar)
  • Datos modelados (SDM)
  • Datos de expertos
  • Planes futuros
• Formatos
  • Datos vectoriales (.shp –> Shapefiles)
    • Puntos –> Ciudades, eventos, sitio de muestreo
    • Líneas
    • Polígonos
  • Datos en grilla (raster)

Shapefiles

Data countriesHigh

library(sf)
library(tidyverse)
library(rworldxtra)
data("countriesHigh")
class(countriesHigh)

[1] "SpatialPolygonsDataFrame"
attr(,"package")
[1] "sp"

Objeto sf

• Transformación a objeto sf (simple feature): este tipo de objeto es más moderno que el tipo sp para manejar datos espaciales. También facilita la conversión a data.frame, de tal manera que se facilita el manejo a través del tidyverse.

class(mundo)

[1] "sf"         "data.frame"

Gráficos

• Gráfico con gggplot2:

mundo %>% 
  ggplot(aes(fill = POP_EST)) +
  geom_sf() +
  labs(title = "Población mundial")

• Filtrando sólo África:

library(viridis)
mundo %>% 
  filter(continent == "Africa") %>% 
  ggplot(aes(fill = POP_EST)) +
  geom_sf() +
  scale_fill_viridis_c()

Exportando shape

• Exportando subconjunto de datos como .shp:

mundo %>% 
  filter(continent == "Africa") %>% 
  dplyr::select(NAME, POP_EST, GDP_MD_EST, GLOCAF) ->
  africa
write_sf(africa, "ejemplo_africa.shp")

• Lectura de archivo shape:

# Función biblioteca raster
africa2 <- shapefile("ejemplo_africa.shp")
class(africa2)

[1] "SpatialPolygonsDataFrame"
attr(,"package")
[1] "sp"

#Función biblioteca sf
africa3 <- read_sf("ejemplo_africa.shp") 
class(africa3)

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

• Gráficos por defecto con la función plot():

plot(africa2)

plot(africa3)

Obtención de shapes con R

• Los nombres ISO-3 de cada país se pueden obtener de la siguiente manera:

getData(name = "ISO3")

• Obtención de datos con función getData() para Colombia:

  • level = 0: división a nivel de país.
  • level = 1: división regional (departamentos)
  • level = 2: división municipal (municipios)

colombia <- getData(name = "GADM", country = "COL", level = 0)
colombia_sf <- st_as_sf(colombia)
colombia_sf %>% 
  ggplot() + 
  geom_sf()

• Proyección del sistema de coordenadas de objetos “colombia” y “colombia_sf”:

crs(colombia)

CRS arguments:
 +proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs 

crs(colombia_sf)

[1] "+proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs"

• Añadiendo punto al mapa. Con el formato sf es posible hacer esto facilmente con ggplot2. Es necesario ingresar a la función st_as_sf() la posición de las columnas longitud y latitud, respectivamente. Además es necesario incoporar la proyección del sistema de coordenadas.

prueba <- data.frame(lon = -70, lat = 5, punto = "ejemplo")

prueba %>% 
  st_as_sf(coords = c(1, 2),
           crs = "+proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs") ->
  prueba_sf
prueba_sf

Simple feature collection with 1 feature and 1 field
geometry type:  POINT
dimension:      XY
bbox:           xmin: -70 ymin: 5 xmax: -70 ymax: 5
CRS:            +proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs
    punto      geometry
1 ejemplo POINT (-70 5)

• Gráfico de shape + punto creado previamete:

ggplot() + 
  geom_sf(data = colombia_sf)  +
  geom_sf(data = prueba_sf, color = "red", size = 4)

Raster

• Datos en grilla, similares a una imagen (pixeles). Poseen características de resolución, proyección, unidades (m, km, …)

Datos CHELSA

• Individuales:

raster1 <- raster("images-tif/CHELSA_Bio1_4km2.tif")
raster2 <- raster("images-tif/CHELSA_Bio6_4km2.tif")
raster1

class      : RasterLayer 
dimensions : 3900, 4980, 19422000  (nrow, ncol, ncell)
resolution : 0.016666, 0.016666  (x, y)
extent     : -112.0029, -29.00618, -39.99518, 25.00222  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs 
source     : D:/DocumentosEdimer/Github/Spatial-Data-Science/SIG_R/images-tif/CHELSA_Bio1_4km2.tif 
names      : CHELSA_Bio1_4km2 
values     : -169.4425, 293.2272  (min, max)

raster2

class      : RasterLayer 
dimensions : 3900, 4980, 19422000  (nrow, ncol, ncell)
resolution : 0.016666, 0.016666  (x, y)
extent     : -112.0029, -29.00618, -39.99518, 25.00222  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs 
source     : D:/DocumentosEdimer/Github/Spatial-Data-Science/SIG_R/images-tif/CHELSA_Bio6_4km2.tif 
names      : CHELSA_Bio6_4km2 
values     : -286.5303, 264  (min, max)

• Stack de raster: es posible almacenar varios raster en un sólo objeto.

mi_stack <- stack(raster1, raster2)
mi_stack

class      : RasterStack 
dimensions : 3900, 4980, 19422000, 2  (nrow, ncol, ncell, nlayers)
resolution : 0.016666, 0.016666  (x, y)
extent     : -112.0029, -29.00618, -39.99518, 25.00222  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs 
names      : CHELSA_Bio1_4km2, CHELSA_Bio6_4km2 
min values :        -169.4425,        -286.5303 
max values :         293.2272,         264.0000 

Gráficos

plot(mi_stack, colNA = "black")

• Escalas equivalentes: a través de la biblioteca rasterVis es posible generar mapas que compartan la misma leyenda.

library(rasterVis)
levelplot(mi_stack)

Operaciones simples

• Suma de raster:

total <- mi_stack[[1]] + mi_stack[[2]]
levelplot(total)

• Promedios:

promedios <- mean(mi_stack)
levelplot(promedios)

Corte de raster

• Es posible cortar el raster para una región específica, en este caso sólo para Colombia. La función crop() permitirá extraer el cuadrante donde se encuenta la región de interés.

raster_colombia <- promedios %>% crop(colombia_sf) # también funciona: crop(colombia)
levelplot(raster_colombia)

• Gráfico con la biblioteca base:

plot(raster_colombia, colNA = "black")

• Corte de sólo Colombia: además de la función crop() es necesario utilizar la función mask(), ambas de la biblioteca raster.

solo_colombia <- promedios %>% 
  crop(colombia_sf) %>% 
  mask(colombia_sf)

Discarded datum Unknown based on WGS84 ellipsoid in CRS definition,
 but +towgs84= values preserved

levelplot(solo_colombia)

• Gráfico con la biblioteca base:

plot(solo_colombia, colNA = "black")

Extracción de valores

• Con los datos creados previamente para representar un punto dentro del mapa, se ejemplifica la extracción de valores del raster (temperatura en este caso) correspondiente a la respectiva longitud y latitud. En este caso equivale a 236. Nota: recordar que los valores de las variables bioclimáticas están multiplicados por 10, de tal manera que la temperatura para el punto ficticio es 236/10 = 23.6°C.

extract(solo_colombia, prueba_sf)

[1] 236

• Se podría añadir este valor a la base de datos:

prueba_sf

Simple feature collection with 1 feature and 2 fields
geometry type:  POINT
dimension:      XY
bbox:           xmin: -70 ymin: 5 xmax: -70 ymax: 5
CRS:            +proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs
    punto      geometry temperatura
1 ejemplo POINT (-70 5)         236

• También es posible aplicar el proceso anterior sobre un stack de raster:

extract(mi_stack, prueba_sf)

     CHELSA_Bio1_4km2 CHELSA_Bio6_4km2
[1,]              251              221

cbind(prueba_sf, extract(mi_stack, prueba_sf)) %>% 
  as.data.frame()

Exportar raster

writeRaster(solo_colombia, "colombia.grd", overwrite = TRUE)

class      : RasterLayer 
dimensions : 1209, 898, 1085682  (nrow, ncol, ncell)
resolution : 0.016666, 0.016666  (x, y)
extent     : -81.8374, -66.87133, -4.229944, 15.91925  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs 
source     : D:/DocumentosEdimer/Github/Spatial-Data-Science/SIG_R/colombia.grd 
names      : layer 
values     : -47.5102, 269.0904  (min, max)

• Lectura del raster:

nuevo_col <- raster("colombia.grd")
plot(nuevo_col, colNA = "black")

Proyecciones

• Importante cuando el área del pixel ingresa en los calculos del análisis, por ejemplo, cuando se está obteniendo abundancia (relativa) de especies. Lugares cercanos al ecuador podrían tener tamaños más grandes de pixel y ocasionar sesgo si no se controla la proyección. En la página web de Projection Wizard es posible obtener proyecciones para lugares específicos, garantizando así que tenemos un mapa de igual área. En este caso para colombia la proyección se muestra en la siguiente imagen:

• Para más información acerca de proyecciones geográficas, recomiendo ver el siguiente video.

• Proyección inicial:

proj4string(solo_colombia)

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

• Cambiando proyección:

colombia_igual <- projectRaster(
  solo_colombia,
  crs = "+proj=laea +lon_0=-74.0917969 +lat_0=0 +datum=WGS84 +units=m +no_defs"
)
plot(colombia_igual, colNA = "black")

Shape + Raster

• Conversión de raster a base de datos (data.frame).

colombia_data <- solo_colombia %>% 
  as("SpatialPixelsDataFrame") %>% 
  as.data.frame()
head(colombia_data)

• Mapa con gggplot2:

# Opcional: cambiando nombre de variable layer
colombia_data <- colombia_data %>% dplyr::rename(Temp = layer)

# Gráfico
ggplot() +
  geom_tile(data = colombia_data, aes(x = x, y = y, fill = Temp)) +
  geom_sf(data = colombia_sf, alpha = 0) +
  scale_fill_viridis_c() +
  theme_bw()

Recursos recomendados

• Diva-GIS
• Spatial Data Science with R
• Worldclim
• CHELSA
• Projection Wizard