Operaciones espaciales con datos rasters

Dr. Francisco Zambrano

Contenidos

Subconjunto espaciales con datos raster
Algebra raster
Operaciones locales
Operaciones globales
Union de rasters (merge)

1. Subconjunto espacial raster

Subconjunto espacial

library(spData)
library(terra)
elev <- rast(system.file("raster/elev.tif", package = "spData"))
grain <- rast(system.file("raster/grain.tif", package = "spData"))
par(mfrow=c(1,2))
plot(elev,main = 'elev');text(elev);
plot(grain,main = 'grain');text(grain)

Subconjunto espacial

Extrae la posición de la celda con cellFromXY y luego con indexación extrae el valor de esa celda con elev[id]

id = cellFromXY(elev, xy = matrix(c(0.1, 0.1), ncol = 2))
elev[id]

  elev
1   16

Extrae el valor que se encuentra en las coordenadas indicadas con la función extract

terra::extract(elev, matrix(c(0.1, 0.1), ncol = 2))

  elev
1   16

Subconjunto espacial

También se puede extraer utilizando otro raster, como rater de indexación.

clip <- rast(xmin = 0.9, xmax = 1.8, ymin = -0.45, ymax = 0.45,
            resolution = 0.3, vals = rep(1, 9))
par(mfrow = c(1,2))
plot(elev);text(elev)

Subconjunto espacial

También se puede extraer utilizando otro raster, como rater de indexación.

clip <- rast(xmin = 0.9, xmax = 1.8, ymin = -0.45, ymax = 0.45,
            resolution = 0.3, vals = rep(1, 9),crs="EPSG:4326")
par(mfrow = c(1,2))
plot(elev);text(elev)
plot(clip,add = TRUE)

Subconjunto espacial

Extraer desde el objeto elev los valores que se encuentran en la extensión del objeto clip

elev[clip] #como valores númericos

  elev
1   18
2   24

elev_clip <- elev[clip,drop = FALSE]
elev_clip

class       : SpatRaster 
size        : 2, 1, 1  (nrow, ncol, nlyr)
resolution  : 0.5, 0.5  (x, y)
extent      : 1, 1.5, -0.5, 0.5  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326) 
source(s)   : memory
varname     : elev 
name        : elev 
min value   :   18 
max value   :   24

Subconjunto espacial

La función crop sirve para realizar elev[clip,drop = FALSE]

elev_crop <- crop(elev,clip)
par(mfrow = c(1,2))
plot(elev_clip)
plot(elev_crop)

Subconjunto espacial

Con indexación

elev[1:2,]

elev[1:2,drop = FALSE]

class       : SpatRaster 
size        : 1, 2, 1  (nrow, ncol, nlyr)
resolution  : 0.5, 0.5  (x, y)
extent      : -1.5, -0.5, 1, 1.5  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326) 
source(s)   : memory
varname     : elev 
name        : elev 
min value   :    1 
max value   :    2

Subconjunto espacial

Utilizando un raster con valores lógicos para ser utilizado como mascara

rmask <- elev
values(rmask) <- sample(c(FALSE, TRUE), 36, replace = TRUE)
par(mfrow = c(1,2))
plot(elev);text(elev)
plot(rmask)

Subconjunto espacial

em <- elev[rmask,drop = FALSE]
par(mfrow = c(1,2))
plot(elev);text(elev)
plot(em);text(em)

Utilizando un raster con valores NA para ser utilizado como mascara con la función mask

par(mfrow = c(1,2))
plot(elev);text(elev)
plot(rmask,add = TRUE)

Subconjunto espacial

Utilizando un raster con valores NA para ser utilizado como mascara.

Se puede utilizar la función mask

rmask[isFALSE(rmask)] <- NA #transforma FALSE to NA para aplicar mask
em <- mask(elev,rmask) 
par(mfrow = c(1,2))
plot(elev);text(elev)
plot(em);text(em)

Subconjunto espacial

Se pueden crear mascaras con operaciones lógicas.

m <- elev
m[m >30] <- NA
em2 <- mask(elev,m)
par(mfrow = c(1,2))
plot(elev);text(elev)
plot(em2);text(em2)

2. Algebra raster

Algunas consideraciones

Si usas multiples raster, estos deben tener el mismo punto origen¹ y la misma resolución
Generalmenet tienen la misma extensión, pero si no, devuelve la operación dentro de la intersección de los dos
El reciclado de raster en operaciones algebraicas funciona de las misma forma que en el caso de los atomic vectors

2. Algebra raster

Muchas de las operaciones genéricas han sido implementadas

Operaciones: +, -, *, /,
Operadores lógicos: >, >=, <, ==, !
Funciones como: abs, round, ceiling, floor, trunc, sqrt, log, log10, exp, cos, sin, atan, tan, max, min, range, prod, sum, any, all

3. Operaciones locales

e1 <- elev + elev
e2<- elev^2
e3 <- log(elev)
e4 <- elev > 5
par(mfrow=c(1,4))
plot(e1,main = 'elev+elev')
plot(e2,main = 'elev^2')
plot(e3,main = 'log(elev)')
plot(e4,main = 'elev > 5')

3. Operaciones locales

Reclasificación

Agrupamos los valores por un rango que define una clase y se le asigna un valor

rcl <- matrix(c(0, 12, 1, 12, 24, 2, 24, 36, 3), ncol = 3, byrow = TRUE)
rcl

     [,1] [,2] [,3]
[1,]    0   12    1
[2,]   12   24    2
[3,]   24   36    3

3. Operaciones locales

Reclasificación

recl <- classify(elev, rcl = rcl)
par(mfrow = c(1,2))
plot(elev);text(elev)
plot(recl);text(recl)

3. Operaciones locales

Funciones de resumen: min, max, mean, prod, sum, Median, cv, range, any, all, siempre devuelven un SpatRaster

mr <- c(e1,e2,e3,e4) #raster multicapa
names(mr) <- c('elev+elev','elev^2','log(elev)','elev > 5')
plot(mr)

3. Operaciones locales

Funciones de resumen: min, max, mean, prod, sum, Median, cv, range, any, all, siempre devuelven un SpatRaster

mr_mean <- mean(mr)
mr_max <- max(mr)
mr_min <- min(mr)
mr_sum <- sum(mr)
par(mfrow= c(1,4))
plot(mr_mean);plot(mr_max);plot(mr_min);plot(mr_sum)

4. Operaciones globales

Usar la función global para obtener un valor resumen de todo el raster. En general operaciones de estadística descriptiva.

global(e1, 'mean')

     mean
elev   37

global(e1, 'max')

     max
elev  72

global(e1, 'min')

     min
elev   2

global(e1, 'sum')

      sum
elev 1332

4. Operaciones globales

Cuando se aplica sobre un raster multicapa, devuleve el resultado de la función sobre cada capa del raster.

global(mr, 'mean')

                 mean
elev+elev  37.0000000
elev^2    450.1666667
log(elev)   2.6588804
elev > 5    0.8611111

global(mr, 'max')

                  max
elev+elev   72.000000
elev^2    1296.000000
log(elev)    3.583519
elev > 5     1.000000

global(mr, 'min')

          min
elev+elev   2
elev^2      1
log(elev)   0
elev > 5    0

global(mr, 'sum')

                  sum
elev+elev  1332.00000
elev^2    16206.00000
log(elev)    95.71969
elev > 5     31.00000

5. Unión de rasters (merge)

Cuando tenemos dos o más raster de una misma área geográfica y queremos unirlos.

r1 <- rast(xmin = 0, xmax=5,ymin=0,ymax=5,vals = 1:25,ncol=5,nrow = 5)
r2 <- rast(xmin = 3, xmax=8,ymin=0,ymax=5,vals = 25:1,ncol=5, nrow = 5)
par(mfrow = c(1,2))
plot(r1);text(r1)
plot(r2);text(r2)

5. Unión de rasters (merge)

Cuando tenemos dos o más raster de una misma área geográfica y queremos unirlos.

plot(r1,xlim = c(0,8));text(r1)
plot(r2,add = TRUE);text(r2)

5. Unión de rasters (merge)

Para unirlos utilizamos la función merge.

Los raster deben tener el mismo origen y resolución.
Los valores del primer raster se mantienen por defecto.

m <- merge(r1,r2)
plot(m);text(m)

5. Unión de rasters (merge)

También podemos utilizar la función mosaic.
mosaic permite definir una función para resolver en la zona de intersección, por defecto utiliza el promedio, pero se puede cambiar la función.
mosaic es más lento que merge,

mos <- mosaic(r1,r2)
plot(mos);text(mos)

5. Unión de rasters (merge)

mos2 <- mosaic(r1,r2,fun='max')
plot(mos2);text(mos2)