LABORATORIO 5. ÁLGEBRA ESPACIAL: OPERACIONES ZONALES
G174: FOTOINTERPRETACIÓN Y TELEDETECCIÓN
Domingo F. Rasilla
2025-03-31
💡 MATERIALES PARA LA ACTIVIDAD:
Los materiales para el desarrollo de esta actividad son los siguientes.
INTRODUCCIÓN
Las operaciones zonales aplican una función matemática a un grupo de celdas dentro de un área o zona específica. Esta(s) zona(s) puede tener un formato vectorial o ráster. Por ejemplo, si deseamos estudiar un determinado fenómeno dentro de los espacio protegidos de España, dispondremos de una primera capa en formato ráster en la que aparece ese fenómeno (por ejemplo, píxeles con el número de cuevas dentro de ese píxel), sobre la que superpondremos una segunda capa con los límites de los espacios protegidos (LIC, ZEPAS etc…). A partir de ahí es posible contar cuántas veces aparece ese fenómeno en cada espacio protegido, etc. Constituyen una herramienta muy adecuada para obtener información, tanto en bruto como a través del cálculo de algún estadístico, como el valor de la pendiente media de una parcela, la lluvia recogida en un valle etc…
ACTIVIDADES INICIALES
Activación de paquetes de R
Para esta actividad es necesario activar los siguientes paquetes:
library(terra)## Warning: package 'terra' was built under R version 4.3.3## terra 1.8.21library(tidyverse)## Warning: package 'tidyverse' was built under R version 4.3.2## Warning: package 'ggplot2' was built under R version 4.3.3## Warning: package 'tidyr' was built under R version 4.3.2## Warning: package 'readr' was built under R version 4.3.2## Warning: package 'purrr' was built under R version 4.3.2## Warning: package 'dplyr' was built under R version 4.3.2## Warning: package 'stringr' was built under R version 4.3.2## Warning: package 'lubridate' was built under R version 4.3.3## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.1 ✔ tibble 3.2.1
## ✔ lubridate 1.9.4 ✔ tidyr 1.3.1
## ✔ purrr 1.0.2## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ tidyr::extract() masks terra::extract()
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorsPreparación de los datos
Deberá establecerse el directorio de trabajo.
setwd("D:/G174_2025/LABORATORIO_5_Algebra_espacial/")Desde ese directorio, se importará un fichero multibanda en formato GeoTIFF, solicitándo a R que muestre sus atributos.
imagen <- rast("D:/G174_2025/LABORATORIO_5_Algebra_espacial/datos/operaciones_zonales/L8_estadisticos_zonales.tif")
imagen## class : SpatRaster
## dimensions : 1395, 1848, 4 (nrow, ncol, nlyr)
## resolution : 30, 30 (x, y)
## extent : 423405, 478845, 4452735, 4494585 (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 30N (EPSG:32630)
## source : L8_estadisticos_zonales.tif
## names : L8_esta~nales_1, L8_esta~nales_2, L8_esta~nales_3, L8_esta~nales_4
## min values : 9287, 8398, 7472, 7127
## max values : 54511, 43380, 54694, 65535Corresponde a un fichero multibanda pero sólo con 4 bandas. Para facilitar el trabajo con este fichero, cada banda se identificará con su nombre.
names(imagen) <- c("blue", "green", "red", "NIR")
names(imagen)## [1] "blue" "green" "red" "NIR"En las escenas Landsat originales, alrededor de la imagen aparece unas regiones con valores ausentes. Estos valores ausentes se codifican como 0. Si no se avisa al programa de esta circunstancia, considerará ese 0 como un valor numérico real. Por ello es necesario especificar esta circunstancia. Además, dado que se realizarán diferentes cálculos estadísticos, conviene evitar la conversión de los valores originales a notación científica.
NAflag(imagen) <- 0 # Convertir los píxeles con 0 en NA.
options(scipen = 999) # Para que no aparezca la notación científicaA continuación se representará gráficamente el objeto importado mediante una imagen en falso color:
plotRGB(imagen,
r=4, g=3, b=2,
stretch = "lin",
axes = TRUE)
✅ DIGITALIZACIÓN DE OBJETOS VECTORIALES CON TERRA:
Después de llamar a la función, comience a hacer clic en el mapa. Cuando haya terminado, pulse ESC. También puede preestablecer el número máximo de clics. En ocasiones, esta función no funciona bien en el dispositivo de trazado por defecto de RStudio. Para evitarlo, primero puede ejecutar dev.new(noRStudioGD = TRUE) que creará una ventana separada para el trazado, a continuación, utilice plot() seguido de click() y haga clic en el mapa.
draw(x=“extent”, # Tipo de geometría a dibujar: “extent”, “polygon”, “line”, or “points” col=“red”, # Color usado lwd=2, # Anchura de la líneas id TRUE, # Si TRUE, un ID numérico se dibuja sobre el mapa n=1000) # Número máximo de clicks (no se aplica cuando x==“extent” -en este caso n es 2)
Respecto al tipo de geometría
Puntos: “points” Captura coordenadas individuales
Líneas: “lines” Crea líneas conectadas
Polígonos: “polygons” Genera áreas cerradas
¿Cómo funciona? Es conveniente iniciar la actividad representando gráficamente la imagen sobre la que se van a digitalizar las geometrías:
plotRGB(imagen[[4:2]],
stretch = "lin",
axes = TRUE)A continuación se ejecutará el código, apareciendo una ventana
gráfica. Sobre esa ventana gráfica se digitalizarán las geometrías
haciendo clic en el mapa para añadir puntos, líneas, polígonos etc. Para
finalizar la digitalización deberá presionarse en “Esc”. El objeto
resultante aparecerá en el Global Envornment. 🔹
Digitalización de puntos:
puntos <- draw("points") # Iniciar la digitalización de puntos
print(puntos) #Mostrar el resultado en pantalla
plot(puntos, col = "red", pch = 16, cex = 2) # Mostrar el resultado de forma gráfica🔹 Digitalización de líneas (se dibuja una línea conectando los puntos en el orden en que se seleccionan)
lineas <- draw("lines")
print(lineas)
plot(lineas, col = "blue", lwd = 2)🔹 Digitalización de polígonos. Una vez se llama a la función, comienza la digitalización de los vértices de los polígonos. Éstos se digitalizan de manera consecutiva, repitiendo la operación tantas veces como polígonos se desee crear. El último punto se conecta automáticamente con el primero para cerrar el polígono). Para finalizar, se hace click con el botón de la derecha del ratón y se selecciona stop
poligono1 <- draw("polygon")
print(poligono1)
plot(poligono1, col = "green", border = "black")En caso de líneas o polígonos, es conveniente agruparlos en un único fichero
zonas <- rbind(poligono1,poligono2, poligono3, poligono4, poligono5, poligono6)
zonas <- project(poligono, imagen)
writeVector(zonas, ".shp")CONSULTAS ESPACIALES: EXTRACCIÓN DE VALORES DE LAS IMÁGENES
En el caso de las imágenes de satélite, se pueden obtener los valores
correspondientes a una determinada función, por ejemplo, la media,
superponiendo un fichero vectorial con polígonos a la imagen en formato
raster, también mediante la función terra::extract
Las funciones disponibles deben incluirse con el argumento
fun. Como en casos anteriores, tanto la imagen raster como
el fichero vectorial deberían tener el mismo CRS. A continuación se van
a plantear varios ejemplos.
Extracción de los valores correspondientes a puntos.
En caso de no disponer de esos puntos, se pueden crear esos puntos utilizando diferentes procedimientos.
● Un único punto
En primer lugar, se crear un dataframe
df_punto <- data.frame(x = 453405, # Dimensión horizontal
y = 4472735) # Dimension verticalDado que df_punto es un dataframe, y no un objeto
espacial, es necesaria su transformación, adquiriendo además la misma
proyección del objeto ráster con el que se está trabajando.
crs(imagen, describe = TRUE)## name authority code
## 1 WGS 84 / UTM zone 30N EPSG 32630
## area
## 1 Between 6°W and 0°W, northern hemisphere between equator and 84°N, onshore and offshore. Algeria. Burkina Faso. Côte' Ivoire (Ivory Coast). Faroe Islands - offshore. France. Ghana. Gibraltar. Ireland - offshore Irish Sea. Mali. Mauritania. Morocco. Spain. United Kingdom (UK)
## extent
## 1 -6, 0, 0, 84Conversión en un SpatVector de puntos
punto <- vect(df_punto,
geom = c("x", "y"),
crs = imagen)Tras ello, se deben superponer los datos vectoriales a la imagen para visualizar la ubicación del punto.
plotRGB(imagen, r=4, g=3, b=2, stretch = "lin", axes = TRUE)
plot(punto,
col = "yellow", # Color de los puntos
cex = 3, # Tamaño del símbolo
pch = 19, # Tipo de puntos
add = TRUE)
A continuación, con la función terra::extract se extraen los valores correspondientes a cada banda bajo ese punto
terra::extract(imagen, # Objeto ráster
punto, # Puntos de los que se va a extraer los valores de un pixel
xy = TRUE) # Incluir coordenadas del punto ## ID blue green red NIR x y
## 1 1 12923 13080 13872 16965 453420 4472730● Múltiples puntos, en este caso, se representarán 10 puntos
distribuidos aleatoriamente por la zona de análisis. Estos puntos se
situarán dentro de las coordenadas del objeto ráster (argumentos
xmin, xmax, ymin, ymin`). El
procedimiento es muy similar al caso anterior. Primero se creará un
dataframe con las coordenadas de los puntos
df_puntos <- data.frame(ID = 1:10, # 10 puntos, cada uno con su Id
x = runif(10, xmin(imagen)+5000, xmax(imagen)-5000), # Dimensiones de los puntos en el eje X
y = runif(10, ymin(imagen)+5000, ymax(imagen)-5000), # Dimensiones de los puntos en el eje Y
tipo = c("A", "B", "A", "C", "B")) # Atributo categóricoEn segundo lugar, ese dataframe se convertirá en un objeto espacial
(SpatVector) de puntos, con el mismo sistema de proyección
de coordenadas el objeto ráster imagen.
puntos <- vect(df_puntos,
geom = c("x", "y"),
crs = imagen)
puntos## class : SpatVector
## geometry : points
## dimensions : 10, 2 (geometries, attributes)
## extent : 429449.7, 472466.2, 4460370, 4478902 (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 30N (EPSG:32630)
## names : ID tipo
## type : <int> <chr>
## values : 1 A
## 2 B
## 3 AIgualmente, se representa gráficamente ese conjunto:
plotRGB(imagen, r=4, g=3, b=2, stretch = "lin", axes = TRUE)
plot(puntos,
col = "yellow", # Color de los puntos
cex = 3, # Tamaño del símbolo
pch = 19, # Tipo de puntos
add = TRUE)
A partir de este momento ya podemos extraer los valores de todas las
bandas de la imagen en las localizaciones especificadas mediante la
función terra::extract, almacenando además estos datos como
un dataframe.
puntos_df_extraidos <- terra::extract(imagen, # Objeto ráster
puntos, # Objeto vectorial
xy = TRUE) # Inclusión en el dataframe las coordenadas de los puntos.
puntos_df_extraidos## ID blue green red NIR x y
## 1 1 11060 11267 12292 14684 469200 4474530
## 2 2 11498 11529 11846 16680 443130 4469220
## 3 3 11871 12216 13344 20520 472470 4468530
## 4 4 10685 10584 11048 16234 438570 4478910
## 5 5 11315 10089 9885 11779 454230 4460370
## 6 6 11499 11480 12049 17437 456180 4473150
## 7 7 10627 10308 10411 16855 458970 4474710
## 8 8 14785 14759 15501 17815 451770 4477020
## 9 9 12774 12417 12781 15321 440790 4461300
## 10 10 11292 10710 11001 13086 429450 4474920Estos valores extraidos de la imagen raster pueden ser mantenidos
como un dataframe puntos_df y exportados (por ejemplo, como
.csv). En caso de querer exportarlo como fichero *csv, el
código sería el siguiente:
write.csv2(puntos_df_extraidos,
file = "D:/G174_2025/LABORATORIO_5_Algebra_espacial/datos/operaciones_zonales/fichero_puntos_extraido.csv",
sep = ";",
row.names = FALSE,
col.names = TRUE)Como poseen las coordenadas proyectadas de los puntos, también podría ser convertido en un fichero espacial (Spatvector) de la siguiente manera:
valores_espectrales <- vect(puntos_df_extraidos,
geom = c("x", "y"),
crs = imagen)
valores_espectrales## class : SpatVector
## geometry : points
## dimensions : 10, 5 (geometries, attributes)
## extent : 429450, 472470, 4460370, 4478910 (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 30N (EPSG:32630)
## names : ID blue green red NIR
## type : <num> <num> <num> <num> <num>
## values : 1 1.106e+04 1.127e+04 1.229e+04 1.468e+04
## 2 1.15e+04 1.153e+04 1.185e+04 1.668e+04
## 3 1.187e+04 1.222e+04 1.334e+04 2.052e+04plot(valores_espectrales, "NIR")
Los datos extraídos pueden ser representados gráficamente para visualizar la ditribución de valores de una banda (en este caso la banda azul).
plot(puntos_df_extraidos$ID,
puntos_df_extraidos$blue, # Indicación para la primera columa
type = "o",
axes = TRUE,
ylim = c(10000,15000), # Valores límite del eje x
main = "Representacion espectral de los puntos extraidos")
Se puede mejorar esta representacion gráfica incluyendo más detalles.
plot(puntos_df_extraidos$ID,
puntos_df_extraidos$blue,
type = "o", # Tipo de símbolo
axes = TRUE,
ylim = c(10000,22000), # Valores límite del eje Y
ylab = "Valores DN",
xlab = "ID de los puntos",
col = "blue",
main = "Representacion espectral de los puntos extraidos")
lines(puntos_df_extraidos$green, type = "o", col = "green")
lines(puntos_df_extraidos$red, type = "o", col = "red")
lines(puntos_df_extraidos$NIR, type = "o", col = "purple")
# Leyenda
legend("topleft",
legend = c("blue", "green", "red", "NIR"),
lty = 1,
bty = "n",
col = c("blue", "green", "red", "purple"),
y.intersp = 0.9,
cex = 1,
title = "PUNTOS")
PERFIL ESPECTRAL
Un gráfico representando los valores de los píxeles de una
determinada superfie en todas las bandas de una imagen multiespectral se
conoce como perfil espectral. Estos perfiles muestran
las diferencias en las propiedades espectrales de los distintos objetos
sobre la superficie terrestre y constituyen la base para el análisis
posterior de las imágenes de satélite. Los valores espectrales pueden
extraerse de cualquier conjunto de datos multiespectrales utilizando esa
misma función extract().
En el ejemplo anterior, hemos extraído valores de datos Landsat para las muestras. Estas muestras incluyen: tierras de cultivo, agua, barbecho, construido y abierto. Primero calculamos los valores medios de reflectancia para cada clase y cada banda.
imagen2 <- rast('D:/G174_2025/LABORATORIO_5_Algebra_espacial/perfil_espectral/imagen2.tif')
imagen2## class : SpatRaster
## dimensions : 1245, 1497, 7 (nrow, ncol, nlyr)
## resolution : 30, 30 (x, y)
## extent : 594090, 639000, 4190190, 4227540 (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 10N (EPSG:32610)
## source : imagen2.tif
## names : ultra-blue, blue, green, red, NIR, SWIR1, ...
## min values : 0.09641791, 0.0748399, 0.04259216, 0.02084067, 0.0008457669, -0.007872183, ...
## max values : 0.73462820, 0.7177562, 0.69246972, 0.78617686, 1.0124315023, 1.043204546, ...Por su parte, el fichero usos es un fichero vectorial en formato ESRI shp.
usos <- vect('D:/G174_2025/LABORATORIO_5_Algebra_espacial/perfil_espectral/usos.shp')
usos## class : SpatVector
## geometry : points
## dimensions : 48, 2 (geometries, attributes)
## extent : 594910.1, 638328.8, 4190718, 4221800 (xmin, xmax, ymin, ymax)
## source : usos.shp
## coord. ref. : +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs
## names : id clase_uso
## type : <chr> <chr>
## values : 1 cultivado
## 1 cultivado
## 1 cultivadoEste fichero vectorial contiene un total de 48 puntos que se han
clasificado en varios grupos de usos de suelo. Mediante la función
extract() se extran los valores espectrales de los puntos
anteriores y, al mismo tiempo, se comprueba los valores de reflectancia
en las diferentes bandas correspondientes a una serie de puntos.
df <- terra::extract(imagen2, usos)
head(df)## ID ultra-blue blue green red NIR SWIR1
## 1 1 0.1154152 0.09286134 0.08362291 0.05464983 0.43908536 0.16941448
## 2 2 0.1109912 0.08893609 0.08464217 0.05553897 0.38280904 0.09902029
## 3 3 0.1083455 0.08470724 0.07037250 0.04621380 0.33631334 0.16342902
## 4 4 0.1132900 0.09097461 0.08505422 0.05432453 0.45357192 0.20305014
## 5 5 0.1414823 0.12404644 0.10318408 0.08483735 0.06360633 0.04842582
## 6 6 0.1240898 0.10487562 0.07991453 0.05924735 0.03840668 0.01986478
## SWIR2
## 1 0.08999872
## 2 0.04517285
## 3 0.07896033
## 4 0.08722286
## 5 0.03894885
## 6 0.01570098A continuación se van a calcular los valores medios de cada uno de
los usos de suelo en todos los puntos correspondientes a su categoría,
mediante la función aggregate(). Esta función tiene la
siguiente sintaxis:
aggregate(formula, data, FUN)
formula: Forma variable ~ grupo, indicando qué variable se debe agrupar.
data: El conjunto de datos donde están las variables que servirán para establecer las categorías resultantes
FUN: La función que se aplicará a cada grupo.
En este caso, la sintaxis es la siguiente
ms <- aggregate(df[ ,-1], # El citado dataframe, pero sin el encabezado (primera fila).
list(usos$clase_uso), # Reemplazo es primera columna con las superficies según usos
mean) # Calculo el valor medio de cada uso
ms## Group.1 ultra-blue blue green red NIR SWIR1
## 1 agua 0.1327735 0.11559101 0.09775448 0.07715008 0.04756178 0.0319213
## 2 artificial 0.1786493 0.16950844 0.17206383 0.18493112 0.22927266 0.2335666
## 3 barbecho 0.1304385 0.11181529 0.09206436 0.09047041 0.12907769 0.2044923
## 4 cultivado 0.1120105 0.08936982 0.08092295 0.05268178 0.40294491 0.1587285
## 5 sin cultivar 0.1390433 0.13788239 0.15237471 0.20631468 0.34316624 0.3646445
## SWIR2
## 1 0.02590160
## 2 0.19642493
## 3 0.19113344
## 4 0.07533869
## 5 0.21842583En vez de la primera columna, se usan los nombres de las filas.
rownames(ms) <- ms[ ,1] # Creo una nueva columna sin encabezado
ms <- ms[ ,-1] # Elimino la 1ª columna antigua
ms## ultra-blue blue green red NIR SWIR1
## agua 0.1327735 0.11559101 0.09775448 0.07715008 0.04756178 0.0319213
## artificial 0.1786493 0.16950844 0.17206383 0.18493112 0.22927266 0.2335666
## barbecho 0.1304385 0.11181529 0.09206436 0.09047041 0.12907769 0.2044923
## cultivado 0.1120105 0.08936982 0.08092295 0.05268178 0.40294491 0.1587285
## sin cultivar 0.1390433 0.13788239 0.15237471 0.20631468 0.34316624 0.3646445
## SWIR2
## agua 0.02590160
## artificial 0.19642493
## barbecho 0.19113344
## cultivado 0.07533869
## sin cultivar 0.21842583A continuación se representará gráficamente el perfil espectral de las siguientes superficies en la imagen landsat anterior.
# Vector de colores que representará cada una de las clases
mycolor <- c('darkred', 'yellow', 'burlywood', 'cyan', 'blue')
# Transformación del dataframe ms en una matriz
ms <- as.matrix(ms)
# Se crea un dibujo vacío que se rellena a continuación
plot(0, ylim=c(0,0.6), xlim = c(1,7), type='n', xlab="Bandas", ylab = "Reflectancia")
# Añadimos al dibujo vacío las diferentes categorías clases_usos
for (i in 1:nrow(ms)){
lines(ms[i,], type = "l", lwd = 3, lty = 1, col = mycolor[i])
}
# Se añade el título
title(main="Firma espectral", font.main = 2)
# Se añade la leyenda
legend("topleft", rownames(ms),
cex=0.8, col=mycolor, lty = 1, lwd =3, bty = "n")
Extracción de valores correspondientes a una línea (creación de un transecto)
También es posible extraer los valores a lo largo de una línea,
creando un transecto. Para ello, se recurrirá a la función
terra::draw() antes señalada. Es conveniente que la
digitalización se realice sobre una imagen lo más real posible, para lo
cual se utilizará la función terra::plotRGB(). El resultado
es un objeto en formato vectorial.
plotRGB(imagen, r=4, g=3, b=2, stretch = "lin", axes = TRUE)
linea <- draw('line')A continuación se extraen los valores con la función
terra::extract().
transecto <- extract(imagen, line)Finalmente, se puede dibujar el transecto, en este caso correspondiente a la banda NIR.
plot(transecto$NIR, type = 'l', ylab = "Reflectancia superficial")Extracción de valores correspondientes a polígonos.
La opción más habitual es que estemos interesados en extraer los valores correspondientes a una o varias zonas de interés (ROI), que habitualmente conforman objetos espaciales en formato vectorial (polígonos). Estos polígonos pueden ser digitalizados u obtenidos de una fuente de información externa.
En las siguientes lineas se importará un fichero vectorial que incluye varios poligonos para calcular posteriormente la moda de todos los píxeles situados dentro de cada polígono, en cada una de las bandas.
zonas <- vect("D:/G174_2025/LABORATORIO_5_Algebra_espacial/datos/operaciones_zonales/zonas.shp")
zonas## class : SpatVector
## geometry : polygons
## dimensions : 4, 1 (geometries, attributes)
## extent : 427981.7, 476274.3, 4456098, 4492225 (xmin, xmax, ymin, ymax)
## source : zonas.shp
## coord. ref. : Undefined Cartesian SRS with unknown unit
## names : FID
## type : <num>
## values : 0
## 1
## 2Dado que el nuevo fichero vectorial no comparte el mismo CRS que la imagen, es requisito ineludible asignarle el mismo CRS.
crs(zonas) <- imagen
zonas## class : SpatVector
## geometry : polygons
## dimensions : 4, 1 (geometries, attributes)
## extent : 427981.7, 476274.3, 4456098, 4492225 (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 30N (EPSG:32630)
## names : FID
## type : <num>
## values : 0
## 1
## 2A continuación, se superpone la imagen raster y el fichero vectorial con los polígonos. Obsérvese que se representa una imagen en falso color (4,3,2 o imagen [[4:2]])
plotRGB(imagen[[4:2]], stretch = "lin", axes = TRUE)
plot(zonas,
col = "yellow",
add =TRUE)
Para obtener todos los valores de los píxeles ubicados dentro de los
polígonos, se usa la función terra::extract(). Esta función
requiere:
El ráster del que deseamos extraer valores
La capa vectorial que contiene los polígonos que deseamos utilizar como límite o límites
Por ejemplo, para la extracción de los valores promedios de los píxeles localizados dentro de los polígonos
promedios <- terra::extract(imagen, # Objeto raster
zonas, # Objeto vectorial
fun = "mean", # Función a calcular: el promedio.
na.rm = TRUE) # Omitir NA en caso de existir.Para ver los primeros datos del dataframe
head(promedios)## ID blue green red NIR
## 1 1 10646.36 10350.68 10758.91 15819.12
## 2 2 11858.66 11697.60 12214.97 16219.51
## 3 3 11857.08 12422.25 14107.29 18270.96
## 4 4 11925.88 12245.39 13494.40 18048.09También podemos usar la función summary() para ver
estadísticas descriptivas que incluyen valores mínimo, 1er cuartil,
mediana, promedio, 3er cuartil y máximo de la banda analizada.
summary(promedios$blue)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 10646 11554 11858 11572 11875 11926También pueden calcularse los 5 números de Tukey (valor minimo, percentil 25, mediana (percentil 50), percentil 75 y valor máximo)
fivenum(promedios$red)## [1] 10758.91 11486.94 12854.68 13800.84 14107.29Si se qusiera calcular los 5 números de Tukey para todas las bandas habría que recurrir a la siguiente sintaxis.
sapply(promedios[c('blue', 'green', 'red', 'NIR')], fivenum)## blue green red NIR
## [1,] 10646.36 10350.68 10758.91 15819.12
## [2,] 11251.72 11024.14 11486.94 16019.32
## [3,] 11857.87 11971.49 12854.68 17133.80
## [4,] 11892.27 12333.82 13800.84 18159.52
## [5,] 11925.88 12422.25 14107.29 18270.96Operaciones zonales a partir de un ráster
Aunque no tan frecuentemente, las mismas operaciones zonales se pueden aplicar a las zonas delimitadas por los valores de un segundo ráster, estando este último conformado por categorías, como por ejemplo los usos de suelo. Debe tenerse además la precaución de que ambos raster tengan la misma extensión y resolución espacial.
ndvi <- rast("D:/G174_2025/LABORATORIO_5_Algebra_espacial/datos/operaciones_zonales/NDVI_zonal.tif")
ndvi## class : SpatRaster
## dimensions : 80, 84, 1 (nrow, ncol, nlyr)
## resolution : 100, 100 (x, y)
## extent : 3635900, 3644300, 2407500, 2415500 (xmin, xmax, ymin, ymax)
## coord. ref. : ETRS89-extended / LAEA Europe (EPSG:3035)
## source : NDVI_zonal.tif
## name : T31TCK_20211012T105011_B08_10m
## min value : 0.1080967
## max value : 0.8858864plot(ndvi)
A continuación se carga el objeto ráster con los usos del suelo
usos <- rast("D:/G174_2025/LABORATORIO_5_Algebra_espacial/datos/operaciones_zonales/usos_zonal.tif")
plot(usos)
¿Tiene etiquetas? No, sólo categorías numéricas.
levels(usos)## [[1]]
## [1] ""A continuación se asignarán etiquetas a las categorías numéricas
levels(usos) <- data.frame(ID = 1:4, etiqueta = c("Matorral", "Cultivos", "Bosque", "Río"))
levels(usos) # Verificar las categorías asignadas## [[1]]
## ID etiqueta
## 1 1 Matorral
## 2 2 Cultivos
## 3 3 Bosque
## 4 4 RíoEl resultado sería el siguiente:
promedio_ndvi <- zonal(ndvi,
usos,
"mean",
na.rm = TRUE)
names(promedio_ndvi) <- c("Usos", "Promedio NDVI")📝 ACTIVIDAD DE EVALUACIÓN: OPERACIONES ZONALES:
Los materiales para esta actividad están incluidos en el mismo fichero que incluye las actividades locales y las globales.
ACTIVIDAD 1. Esta actividad se llevará a cabo a partir de los siguientes materiales:
Imagen_Reflectancia_superficial.tif: imágen multiespectral procedente de la plataforma Landsat 8. Para reducir su tamaño y volumen, sólo incluye las bandas blue, green, red y NIR.puntos.shp: puntos que representan localizaciones.poligonos.shp: límites una serie de polígonos.
A continuación se enumeran los pasos que conforman las actividades:
Convierte el fichero
Imagen_Reflectancia_superficial.tifen un objeto SpatRaster y cambia el nombre original de las bandas por su nombre real.Convierte el fichero
puntos.shpen un objeto SpatVector, comprobando que su crs es el mismo que el del objeto SpatRaster importado anteriormente, y representándolo gráficamente superpuesto a la banda NIR del SpatRaster.Haz lo mismo con el fichero
poligonos.shp.Extrae los valores medios de cada una de las bandas sobre los puntos, pero sin incluir las coordenadas de éstos últimos.
Representa gráficamente los datos extraidos para visualizar la ditribución de valores de las bandas.
ACTIVIDAD 2. Esta actividad se llevará a cabo a partir de los siguientes materiales:
lst.tif: imagen con la temperatura de la superficie terrestre correspondiente al día 18-07-2022 centrado en la ciudad de Madrid.spots_lst: imagen conteniendo las áreas definidas como “puntos” fríos o calientes en el ejercicio sobre operaciones locales.ndvi.tif: imagen conteniendo el índice NDVI correspondiente al día 18-07-2022 centrado en la ciudad de Madrid.evi.tif: imagen conteniendo el índice EVI correspondiente al día 18-07-2022 centrado en la ciudad de Madrid.savi.tif: imagen conteniendo el índice SAVI correspondiente al día 18-07-2022 centrado en la ciudad de Madrid.distritos.geojson: distritos en los que se divide la ciudad de Madrid (se ha omitido el distrito de El Pardo-Fuencarral por sus características singulares.
A continuación, cada alumno deberá elaborar un script respondiendo a las siguientes cuestiones:
Crea un fichero ráster multibanda que incluya todos los ficheros ráster señalados anteriormente:
Carga el fichero vectorial
distritos.geojsonPosteriormente, calcula el valor medio y la desviación típica de cada uno de los objetos ráster en cada uno de los distritos de la ciudad de Madrid
Importa el ráster
spots_lst.Extraer los valores del raster dentro de cada distritos
:::