Operaciones con geometrias

Preparativos

Carga de paquetes

library(sf)
library(raster)
library(dplyr)
library(rmapshaper)
library(spData)

Conjuntos de datos utilizados

# Carga de la capa de provincias
provincias <-
  st_read(
    "https://raw.githubusercontent.com/gf0604-procesamientodatosgeograficos/2021i-datos/main/ign/delimitacion-territorial-administrativa/cr_provincias_simp_wgs84.geojson",
    quiet = TRUE
  )

# Carga de la capa de cantones
cantones <-
  st_read(
    "https://raw.githubusercontent.com/gf0604-procesamientodatosgeograficos/2021i-datos/main/ign/delimitacion-territorial-administrativa/cr_cantones_simp_wgs84.geojson",
    quiet = TRUE
  )

# Carga de la capa de áreas silvestres protegidas (ASP)
asp <-
  st_read(
    "https://raw.githubusercontent.com/gf0604-procesamientodatosgeograficos/2021i-datos/main/sinac/asp/asp-wgs84.geojson",
    quiet = TRUE
  )

# Carga de la capa de red vial
red_vial <-
  st_read(
    "https://raw.githubusercontent.com/gf0604-procesamientodatosgeograficos/2021i-datos/main/ign/red-vial/cr-redvial-simp-wgs84.geojson",
    quiet = TRUE
  )

# Carga de la capa de mamíferos
mamiferos <-
  st_read(
    "https://raw.githubusercontent.com/gf0604-procesamientodatosgeograficos/2021i-datos/main/gbif/mammalia-cr-registros.csv",
    options = c(
      "X_POSSIBLE_NAMES=decimalLongitude",
      "Y_POSSIBLE_NAMES=decimalLatitude"
    ),
    quiet = TRUE
  )

# Asignación del sistema de coordenadas
st_crs(mamiferos) <- 4326

Introducción

Esta lección brinda una visión general de las operaciones con geometrías en datos vectoriales implementadas en el paquete sf y en datos raster implementadas en el paquete raster. Estas operaciones trabajan con la columna de geometrías (ej. geometry, geom) del paquete sf, para el caso de los datos vectoriales, y con la localización geográfica de los pixeles para el caso de los datos raster. En la sección final, se muestran varias operaciones de interacción entre los modelos raster y vectorial.

Esta lección está basada en el capítulo 5 del libro Lovelace, R., Nowosad, J., & Muenchow, J. (2019). Geocomputation with R..

Datos vectoriales

Las operaciones con geometrías en datos vectoriales incluyen:

Operaciones con geometrías con el paquete sf

Estas operaciones modifican las geometrías de objetos vectoriales (sf).

Simplificación

La simplificación puede realizarse en geometrías de líneas y polígonos. Reduce la cantidad de memoria, disco y ancho de banda que utilizan las geometrías. Para simplificar geometrías, sf incluye el método st_simplify, basado en el algoritmo de Douglas-Peucker, el cual recibe el argumento dTolerance para controlar el nivel de generalización de las unidades del mapa. Este argumento se expresa en las unidades de medida del CRS de la capa, por lo que es conveniente utilizar un CRS con unidades de medida de distancias (ej. metros).

# Conversión de la capa de provincias a CRTM05 (unidad de medida = metros)
provincias <-
  provincias %>%
  st_transform(crs = 5367)

# Mapa de la capa de provincias sin simplificación
plot(
  provincias$geometry,
  extent = extent(280000, 660000, 880000, 1250000),  
  main = "Provincias con geometrías no simplificadas",
  axes = TRUE,
  graticule = TRUE)


# Simplificación sin preservación de topología
provincias_simp <-
  provincias %>%
  st_simplify(dTolerance = 5000, preserveTopology = FALSE)

# Mapa de la capa de provincias con simplificación y sin preservación de topología
plot(
  provincias_simp$geometry,
  extent = extent(280000, 660000, 880000, 1250000),  
  main = "Provincias simplificadas sin preservación de topología",
  axes = TRUE,
  graticule = TRUE)


# Simplificación con preservación de topología
provincias_simp_topo <-
  provincias %>%
  st_simplify(dTolerance = 5000, preserveTopology = TRUE)

# Mapa de la capa de provincias con simplificación y con preservación de topología
plot(
  provincias_simp_topo$geometry,
  extent = extent(280000, 660000, 880000, 1250000),  
  main = "Provincias simplificadas con preservación de topología",
  axes = TRUE,
  graticule = TRUE)


# Tamaño de la capa original
object.size(provincias)

840296 bytes

# Tamaño de la capa simplificada sin preservación de topología
object.size(provincias_simp)

18720 bytes

# Tamaño de la capa simplificada con preservación de topología
object.size(provincias_simp_topo)

70400 bytes

La función rmapshaper::ms_simplify() proporciona un método alternativo para la simplificación de geometrías, el cual preserva la topología.

# Simplificación con rmapshaper::ms_simplify()
provincias_simp <-
  provincias %>%
  rmapshaper::ms_simplify(keep = 0.1, keep_shapes = TRUE)

# Mapa de la capa de provincias con simplificación mediante rmapshaper::ms_simplify()
plot(
  provincias_simp$geometry,
  extent = extent(280000, 660000, 880000, 1250000),  
  main = "Provincias simplificadas con rmapshaper::ms_simplify()",
  axes = TRUE,
  graticule = TRUE)


# Tamaño de la capa simplificada con rmapshaper::ms_simplify()
object.size(provincias_simp)

115016 bytes

Centroides

Un centroide es un punto que identifica el centro de un objeto geográfico. Puede calcularse para geometrías de líneas y de polígonos y se utilizan para brindar una representación simplificada de geometrías más complejas. Existen varios métodos para calcularlos.

El paquete sf incluye la función st_centroid() la cual calcula el centroide geográfico (comúnmente llamado “el centroide”). Es posible que el centroide geográfico se ubique fuera de la geometría “padre” (ej. en una con forma de rosca). Para evitar este resultado, la función st_point_on_surface() se asegura de que el centroide esté siempre dentro de la geometría “padre”.

# Costa Rica y sus centroides calculados con st_centroid() y st_point_on_surface()
plot(
  st_union(provincias),
  main = "Costa Rica centroides st_centroid (rojo) st_point_on_surface (verde)",
  axes = TRUE,
  graticule = TRUE)

plot(st_centroid(st_union(provincias)),
     add = TRUE,
     pch = 16,
     col = "red")

plot(
  st_point_on_surface(st_union(provincias)),
  add = TRUE,
  pch = 16,
  col = "green")


# Coordenadas del centroide calculado con st_centroid()
# CRTM05
st_coordinates(st_centroid(st_union(provincias)))

       X       Y
1 478674 1102734
# WGS84
st_coordinates(st_transform(st_centroid(st_union(provincias)), crs = 4326))

          X        Y
1 -84.19452 9.972735

# Coordenadas del centroide calculado con st_point_on_surface()
# CRTM05
st_coordinates(st_point_on_surface(st_union(provincias)))

         X       Y
1 539376.2 1065149
# WGS84
st_coordinates(st_transform(st_point_on_surface(st_union(provincias)), crs = 4326))

          X        Y
1 -83.64122 9.632758


# Provincias de Costa Rica y sus centroides calculados con st_centroid() y st_point_on_surface()
plot(
  provincias$geometry,
  extent = extent(280000, 660000, 880000, 1250000),    
  main = "Provincias centroides st_centroid (rojo) st_point_on_surface (verde)",
  axes = TRUE,
  graticule = TRUE)

plot(st_centroid(provincias),
     add = TRUE,
     pch = 16,
     col = "red")

plot(
  st_point_on_surface(provincias),
  add = TRUE,
  pch = 16,
  col = "green")



# Ruta 32 y sus centroides calculados con st_centroid() y st_point_on_surface()
plot(
  provincias$geometry,
  extent = extent(280000, 660000, 880000, 1250000),    
  main = "Ruta 32 y centroides st_centroid (rojo) y st_point_on_surface (verde)",
  axes = TRUE,
  graticule = TRUE)

ruta_32 <-
  red_vial %>%
  filter(num_ruta == "32") %>%
  st_transform(crs = 5367)

plot(
  ruta_32$geometry,
  add = TRUE,
  col = "blue")

plot(
  st_centroid(st_union(ruta_32)),
  add = TRUE,
  pch = 16,
  col = "red")

plot(
  st_point_on_surface(st_union(ruta_32)),
  add = TRUE,
  pch = 16,
  col = "green")

Áreas de amortiguamiento (buffers)

Los buffers son polígonos creados alrededor de otra geometría, ya sea otro polígono, una línea o un punto. El paquete sf incluye la función st_buffer() para la generación de buffers.

# Buffer alrededor de la ruta 32
plot(
  st_buffer(st_union(ruta_32), 5000),
  main = "Buffer alrededor de la ruta 32",
  axes = TRUE,
  graticule = TRUE)

plot(
  ruta_32$geometry,
  col = "blue",
  add = TRUE
)

Es común el uso de buffers en análisis de datos, para responder preguntas como, por ejemplo, “¿cuántos puntos hay alrededor de una línea?” o “¿cuáles especies pueden encontrarse en las márgenes de un río?”.

# Registros de presencia de mamíferos terrestres ubicados alrededor de la ruta 32
mamiferos <-
  mamiferos %>%
  filter(order != "Chiroptera") %>% # se excluye el orden de los murciélagos
  st_transform(crs = 5367)

buffer_ruta_32 <-
  ruta_32 %>%
  st_buffer(dist = 5000)

mamiferos_buffer_ruta_32 <-
  st_join(mamiferos, buffer_ruta_32) %>%
  filter(!is.na(codigo))

plot(
  st_union(buffer_ruta_32),
  main = "Mamíferos terrestres alrededor de la ruta 32",
  axes = TRUE,
  graticule = TRUE
)

plot(ruta_32$geometry,
     col = "blue",
     add = TRUE)

plot(
  mamiferos_buffer_ruta_32,
  pch = 16,
  col = "orange",
  add = TRUE
)


# 10 especies con mayor cantidad de registros
mamiferos_buffer_ruta_32 %>% 
  st_drop_geometry() %>%
  filter(!is.na(species) & species != "") %>%
  group_by(species) %>% 
  summarise(registros = n()) %>%
  arrange(desc(registros)) %>%
  slice(1:10)

# A tibble: 10 x 2
   species                  registros
   <chr>                        <int>
 1 Sciurus variegatoides           84
 2 Bradypus variegatus             81
 3 Heteromys desmarestianus        48
 4 Choloepus hoffmanni             47
 5 Peromyscus mexicanus            38
 6 Melanomys caliginosus           30
 7 Alouatta palliata               20
 8 Orthogeomys cherriei            20
 9 Didelphis marsupialis           16
10 Cebus capucinus                 15

Recortes (clipping)

El recorte de una geometría con base en la forma de otra puede realizarse con el método st_intersection(), el cual retorna la intersección entre dos geometrías.

# Recorte de la sección del Parque Internacional La Amistad (PILA) ubicada en Puntarenas
pila <-
  asp %>%
  filter(nombre_asp == "Internacional La Amistad") %>%
  st_transform(crs = 5367)

puntarenas_limon <-
  provincias %>%
  filter(provincia == "Puntarenas" | provincia == "Limón")

# Mapa de Puntarenas, Limón y el PILA
plot(
  puntarenas_limon$geometry,
  main = "Puntarenas, Limón y el PILA",
  extent = extent(350000, 660000, 880000, 1200000),
  axes = TRUE,
  graticule = TRUE)

plot(
  pila$geometry,
  border = "green",
  add = TRUE)


puntarenas <-
  provincias %>%
  filter(provincia == "Puntarenas")

# Recorte de la sección del PILA ubicada en Puntarenas
pila_puntarenas <- st_intersection(pila, puntarenas)

# Mapa de la sección recortada
plot(
  pila_puntarenas$geometry,
  main = "Sección del PILA ubicada en Puntarenas",
  col = "red",
  axes = TRUE,
  graticule = TRUE)

plot(
  puntarenas_limon$geometry,
  add = TRUE
)

Uniones de geometrías

En lecciones anteriores, se ha mostrado como agregar geometrías mediante los métodos agregate() y summarize(). Internamente, ambos métodos utilizan la función st_union() para combinar las geometrías y disolver sus límites.

# Cantones de la provincia de San José
cantones_sanjose <-
  cantones %>%
  filter(provincia == "San José")

plot(
  cantones_sanjose$geometry,
  main = "Cantones de San José",
  axes = TRUE,
  graticule = TRUE)  


# Cantones de la provincia de San José unificados
cantones_sanjose_unificados <- 
  st_union(cantones_sanjose)

plot(
  cantones_sanjose_unificados,
  main = "Cantones de San José unificados",
  axes = TRUE,
  graticule = TRUE)

Datos raster

Las operaciones con geometrías en datos raster incluyen:

Operaciones con geometrías con el paquete raster

Intersecciones geométricas

En lecciones anteriores, se mostró como extraer valores de un raster al que se sobreponen otros objetos espaciales. Para recuperar otro objeto espacial, solamente debe especificarse el argumento drop = FALSE.

El siguiente bloque de código retornará un objeto raster en el cual los puntos centrales de la capa de altitud se sobreponen con el de otro objetos raster que lo interseca (o lo recorta).

# Especificación del directorio de trabajo (debe ser una ruta existente)
setwd("C:/Users/mfvargas/Downloads")

# Obtención de la capa de altitud
alt <-
  getData(
    "worldclim",
    var = "alt",
    res = .5,
    lon = -84,
    lat = 10
  )

# Reproyección de la capa de altitud a CRTM05
alt <-
  alt %>%
  projectRaster(crs = 5367)

# Recorte de la capa de altitud con base en la capa vectorial de provincias
altitud <-
  alt %>%
  crop(provincias) %>%
  mask(provincias)

plot(
  altitud,
  ext = extent(280000, 660000, 880000, 1250000),
  main = "Altitud de Costa Rica",
  axes = TRUE
)

# Objeto raster para intersecar con el de altitud
recorte <-
  raster(
    xmn = -84.5,
    xmx = -83.5,
    ymn = 9.5,
    ymx = 10.5,
    res = 0.5
  )

recorte <-
  recorte %>%
  projectRaster(crs = 5367)

# Intersección de altitud y recorte
altitud_recorte <- altitud[recorte, drop = FALSE]

plot(
  provincias$geometry,
  extent = extent(280000, 660000, 880000, 1250000),
  main = "Recorte de altitud de Costa Rica",
  axes = TRUE,
  graticule = TRUE
)

plot(altitud_recorte,
     add = TRUE)

El mismo resultado puede obtenerse mediante el uso de las funciones intersect() y crop().

Cambio de extensión y de origen

Muchas operaciones que se realizan con múltiples objetos raster requieren que estos tengan la misma resolución, proyección, origen y extensión. Como ejemplos de estas operaciones pueden mencionarse el álgebra raster y la combinación de imágenes satelitales con diferentes proyecciones y resoluciones. Si hay diferencias en las características de los objetos raster que impiden realizar estas operaciones, deben “alinearse” primero.

Cambio de extensión

A continuación, se ilustra el caso de objetos raster que difieren solamente en su extensión (i.e. cantidad de filas y columnas). El siguiente bloque de código agrega una fila y dos columnas a los bordes del objeto raster elev del paquete spData, por medio del método extend(). A las nuevas celdas, se les asigna un valor de 1000.

# Objeto raster de elevación
plot(elev)

# Modificación de la extensión del objeto raster
elev_2 <- extend(elev, c(1, 2), value = 1000)
plot(elev_2)

Al intentar realizar una operación algebraica que requiere que los raster tengan las mismas dimensiones, la operación se ejecuta en la intersección y, además, se muestra una advertencia.

# Suma de dos objetos raster no alineados
elev_3 <- elev + elev_2
plot(elev_3)

Ambos objetos raster pueden alinearse con extent(). En este caso, no se especifica la cantidad de filas y columnas que se desan añadir a un raster, sino otro raster cuyas dimensiones se utilizan como referencia. Con el argumento value, se especificar el valor de las nuevas celdas (su valor por defecto es NA).

# Alineación de la extensión de dos objetos raster
elev_4 <- extend(elev, elev_2)
plot(elev_4)

# Suma de los dos objetos raster alineados
plot(elev + elev_4)

Cambio de origen

El origen de un objeto raster es la celda esquinera más cercana a las coordenas (0, 0). La función origin() retorna las coordenadas del origen.

En el ejemplo del siguiente bloque de código, existe una celda con las coordenadas (0, 0), pero ese no siempre es el caso.

# Impresión del origen de un objeto raster
origin(elev)

[1] 0 0

Seguidamente, se muestra el efecto de cambiar el origen de un objeto raster.

# Cambio del origen de un objeto raster
origin(elev_4) <- c(0.25, 0.25)

plot(elev_4)
plot(elev, add = TRUE)

Para realizar operaciones aritméticas con álgebra raster, los orígenes de ambos objetos deben coincidir. Algunas operaciones (ej. cambio de resolución) pueden provocar un cambio en el origen.

Agregación y desagregación

La resolución de un raster puede disminuirse con la función aggregate() o aumentarse con la función disaggregate().

En el siguiente bloque de código, se utiliza la función aggregate() para disminuir la resolución del raster de altitud de Costa Rica por un factor de 4, como se especifica con el argumentofactor = 4. Así, se generará una celda en el raster agrupado por cada ocho celdas del raster original. El argumento fun = mean indica que en el raster agrupado, cada celda será el promedio de las cuatro celdas correspondientes en el raster original.

# Capa de altitud de Costa Rica
plot(
  altitud,
  ext = extent(280000, 660000, 880000, 1250000),
  main = "Capa de altitud de Costa Rica",
  axes = TRUE)

# Agrupación de la capa de altitud de Costa Rica
altitud_agregada <- 
  altitud %>%
  aggregate(fact = 8, fun = mean)

# Capa de altitud agrupada de Costa Rica
plot(
  altitud_agregada,
  ext = extent(280000, 660000, 880000, 1250000),
  main = "Capa agrupada de altitud de Costa Rica",
  axes = TRUE)

# Metadatos de la capa original de altitud
altitud

class      : RasterLayer 
dimensions : 679, 555, 376845  (nrow, ncol, ncell)
resolution : 906, 932  (x, y)
extent     : 156171.6, 659001.6, 608735.2, 1241563  (xmin, xmax, ymin, ymax)
crs        : +proj=tmerc +lat_0=0 +lon_0=-84 +k=0.9999 +x_0=500000 +y_0=0 +ellps=WGS84 +units=m +no_defs 
source     : memory
names      : alt_23 
values     : -6.604836, 3717.273  (min, max)
# Metadatos de la capa agrupada de altitud
altitud_agregada

class      : RasterLayer 
dimensions : 85, 70, 5950  (nrow, ncol, ncell)
resolution : 7248, 7456  (x, y)
extent     : 156171.6, 663531.6, 607803.2, 1241563  (xmin, xmax, ymin, ymax)
crs        : +proj=tmerc +lat_0=0 +lon_0=-84 +k=0.9999 +x_0=500000 +y_0=0 +ellps=WGS84 +units=m +no_defs 
source     : memory
names      : alt_23 
values     : 3.816522, 3174.02  (min, max)
# Tamaño de la capa original de altitud
object.size(altitud)

3028992 bytes
# Tamaño de la capa agrupada de altitud
object.size(altitud_agregada)

61832 bytes

La función disaggregate() genera varias celdas por cada celda del raster original. Debe especificarse un método con el argumento method. Su valor puede ser "", para copiar los valores de la celda de entrada o bilinear, que es un método de interpolación.

En el siguiente bloque de código, se desagrega el raster de elevación por un factor de dos.

# Desagregación de la capa de elevación
elev_desagregada <- 
  elev %>%
  disaggregate(fact = 2, method = "bilinear")

plot(elev_desagregada)

# Metadatos de la capa original de elevación
elev

class      : RasterLayer 
dimensions : 6, 6, 36  (nrow, ncol, ncell)
resolution : 0.5, 0.5  (x, y)
extent     : -1.5, 1.5, -1.5, 1.5  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
source     : memory
names      : layer 
values     : 1, 36  (min, max)
# Metadatos de la capa desagrupada de elevación
elev_desagregada

class      : RasterLayer 
dimensions : 12, 12, 144  (nrow, ncol, ncell)
resolution : 0.25, 0.25  (x, y)
extent     : -1.5, 1.5, -1.5, 1.5  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
source     : memory
names      : layer 
values     : -0.75, 37.75  (min, max)
# Tamaño de la capa original de elevación
object.size(elev)

12776 bytes
# Metadatos de la capa desagrupada de elevación
object.size(elev_desagregada)

13784 bytes

El proceso de calcular valores para nuevos pixeles se denomina “remuestreo” (resampling). El paquete raster incluye una función llamada resample(), el cual alinea simultáneamente varias propiedades de objetos raster: extensión, origen y resolución.

Interacciones raster-vector

Recorte (cropping) de datos raster

Los métodos crop() y mask() pueden utilizarse para recortar (crop) un objeto raster con base en el contorno de un objeto vectorial.

# crop() a capa de altitud
altitud <-
  alt %>%
  crop(provincias)

plot(
  altitud,
  ext = extent(280000, 660000, 880000, 1250000),
  main = "crop() a capa de altitud de Costa Rica",
  axes = TRUE)

# crop() + mask() a capa de altitud
altitud <-
  alt %>%
  crop(provincias) %>%
  mask(provincias)

plot(
  altitud,
  ext = extent(280000, 660000, 880000, 1250000),
  main = "crop() + mask() a capa de altitud de Costa Rica",
  axes = TRUE)

Extracción de valores raster

La extracción de valores raster es el proceso de identificar y retornar los valores asociados con una capa raster en localizaciones específicas, generalmente determinadas por objetos vectoriales, ya sean puntos, líneas o polígonos. En el paquete raster, la extracción se implementa a través de la función extract().

El siguiente bloque de código extrae el valor del objeto raster de altitud para cada punto de la capa vectorial de mamíferos de Costa Rica.

# Columna con la altitud de cada registro de mamíferos con base en la capa raster de altitud
mamiferos$altitud <- raster::extract(altitud, mamiferos)

mamiferos %>%
  select(order, family, genus, species, altitud) %>%
  arrange(desc(altitud)) %>%
  slice(1:10)

Simple feature collection with 10 features and 5 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 554896.9 ymin: 1050485 xmax: 554896.9 ymax: 1050485
Projected CRS: CR05 / CRTM05
            order    family      genus           species  altitud
1  Perissodactyla Tapiridae  Tapirella Tapirella bairdii 3573.324
2  Perissodactyla Tapiridae  Tapirella Tapirella bairdii 3573.324
3  Perissodactyla Tapiridae  Tapirella Tapirella bairdii 3573.324
4  Perissodactyla Tapiridae  Tapirella Tapirella bairdii 3573.324
5      Lagomorpha Leporidae Sylvilagus  Sylvilagus dicei 3573.324
6      Lagomorpha Leporidae Sylvilagus  Sylvilagus dicei 3573.324
7      Lagomorpha Leporidae Sylvilagus  Sylvilagus dicei 3573.324
8      Lagomorpha Leporidae Sylvilagus  Sylvilagus dicei 3573.324
9      Lagomorpha Leporidae Sylvilagus  Sylvilagus dicei 3573.324
10     Lagomorpha Leporidae Sylvilagus  Sylvilagus dicei 3573.324
                   geometry
1  POINT (554896.9 1050485)
2  POINT (554896.9 1050485)
3  POINT (554896.9 1050485)
4  POINT (554896.9 1050485)
5  POINT (554896.9 1050485)
6  POINT (554896.9 1050485)
7  POINT (554896.9 1050485)
8  POINT (554896.9 1050485)
9  POINT (554896.9 1050485)
10 POINT (554896.9 1050485)

Ejercicio: genere un data frame con la altitud promedio para cada especie de felinos (familia Felidae).

Rasterización

La rasterización es la conversión de objetos vectoriales en su representación raster. En el paquete raster, la función rasterize() realiza esta labor.

El siguiente bloque de código asigna a cada celda de un objeto raster un valor de 1 si contiene puntos de una capa vectorial y 0 en caso contrario.

# Generación de un raster con celda = 1 si hay presencia, celda = 0 si no hay presencia

# Plantilla de raster
raster_plantilla <-
  altitud %>%
  aggregate(fact = 32)

# Rasterización
mamiferos_raster_presencia <-
  rasterize(filter(mamiferos, species == "Puma concolor"),
            raster_plantilla,
            field = 1)

# Mapeo
plot(
  mamiferos_raster_presencia,
  ext = extent(280000, 660000, 880000, 1250000),
  main = "Presencia de Puma concolor",
  axes = TRUE
)

plot(provincias$geometry,
     add = TRUE)

El siguiente bloque de código asigna a cada celda de un objeto raster un valor correspondiente a la cantidad de puntos ubicados en la celda.

# Generación de un raster con celda = cantidad de registros

# Plantilla de raster
raster_plantilla <-
  altitud %>%
  aggregate(fact = 32)

# Rasterización
mamiferos_raster_registros <-
  rasterize(mamiferos,
            raster_plantilla,
            field = 1,
            fun = "count")

# Mapeo
plot(
  mamiferos_raster_registros,
  ext = extent(280000, 660000, 880000, 1250000),
  main = "Cantidad de registros de mamíferos",
  axes = TRUE
)

plot(provincias$geometry,
     add = TRUE)

El siguiente bloque de código asigna a cada celda de un objeto raster un valor correspondiente a la cantidad de valores diferentes de un atributo de los puntos ubicados en la celda.

# Generación de un raster con celda = cantidad de especies

# Plantilla de raster
raster_plantilla <-
  altitud %>%
  aggregate(fact = 32)

# Rasterización
mamiferos_raster_especies <-
  rasterize(filter(mamiferos, family == "Felidae"),
            raster_plantilla,
            field = "species",
            fun = function(x, ...) {length(unique(na.omit(x)))})

# Mapeo
plot(
  mamiferos_raster_especies,
  ext = extent(280000, 660000, 880000, 1250000),
  main = "Cantidad de especies de felinos",
  axes = TRUE
)

plot(provincias$geometry,
     add = TRUE)

Vectorización espacial

La vectorización es la operación opuesta a la rasterización. Convierte datos raster en vectoriales. En el paquete raster, se implementa a través de operaciones como rasterToPoints(), rasterToContour() y rasterToPolygons().

En el siguiente bloque de código se realiza una rasterización a puntos.

# Vectorización a puntos de la capa de altitud de Costa Rica
altitud_puntos = rasterToPoints(aggregate(altitud, fact = 8), spatial = TRUE) %>%
  st_as_sf()

plot(
  altitud_puntos,
  ext = extent(280000, 660000, 880000, 1250000),
  main = "Centroides de celdas raster",
  axes = TRUE)

En el siguiente bloque de código se realiza una rasterización a líneas de contorno.

# Vectorización a líneas de contorno de la capa de altitud de Costa Rica
altitud_contorno = rasterToContour(altitud)

plot(
  altitud,
  ext = extent(280000, 660000, 880000, 1250000),
  main = "Líneas de contorno de altitud",
  axes = TRUE
)
plot(altitud_contorno, add = TRUE)

Seguidamente, se utiliza hillShade() para generar un efecto de “sombras de montañas” (hill shading). Esta función recibe como argumentos slope y aspect, los cuales pueden calcularse con la función terrain().

# Creación de una sombra de montañas
sombra = hillShade(
  slope = terrain(altitud_recorte, "slope"),
  aspect = terrain(altitud_recorte, "aspect")
)

plot(
  sombra, 
  main = "Sombra de montañas y líneas de contorno",
  col = gray(0:100 / 100), 
  legend = FALSE)
plot(
  altitud_recorte,
  col = terrain.colors(25),
  alpha = 0.5,
  legend = FALSE,
  add = TRUE
)
contour(altitud_recorte, col = "white", add = TRUE)

Finalmente, se ejemplifica el uso de la función rasterToPolygons(), la cual convierte cada celda raster a un polígono.

# Capa raster de tipos de granos en el suelo
plot(
  grain,
  main="Raster")


# Celdas raster convertidas a polígonos individuales
grain_poligonos <-
  rasterToPolygons(grain) %>% 
  st_as_sf()
plot(
  grain_poligonos,
  main="Polígonos")


# Celdas raster convertidas a polígonos agrupados
grain_poligonos_agrupados <-
  grain_poligonos %>% 
  group_by(layer) %>%
  summarize()
plot(
  grain_poligonos_agrupados,
  main="Polígonos agrupados")


# Tamaño del objeto raster
object.size(grain)

14248 bytes

# Tamaño del objeto vectorial de polígonos individuales
object.size(grain_poligonos)

34224 bytes

# Tamaño del objeto vectorial de polígonos agrupados
object.size(grain_poligonos_agrupados)

16480 bytes

Corrections

If you see mistakes or want to suggest changes, please create an issue on the source repository.

Reuse

Text and figures are licensed under Creative Commons Attribution CC BY-SA 4.0. Source code is available at https://github.com/gf0604-procesamientodatosgeograficos/2021i-leccion-15/, unless otherwise noted. The figures that have been reused from other sources don't fall under this license and can be recognized by a note in their caption: "Figure from ...".