LABORATORIO 1: INTRODUCCIÓN AL MANEJO DE FICHEROS RÁSTER (VALORES CONTINUOS) CON R
G174: FOTOINTERPRETACIÓN Y TELEDETECCIÓN
Domingo F. Rasilla
2025-03-02
OBJETOS SpatRaster MONOBANDA
La actividad comenzará creando un objeto SpatRaster, que
denominaremos base. El procedimiento para crear este objeto
puede realizarse de varias maneras.
library(terra)## Warning: package 'terra' was built under R version 4.3.3## terra 1.8.21► Por ejemplo, este fichero puede crearse a partir de una matriz de datos, que luego se convierte en objeto ráster.
matriz <- matrix(1:36, 6, 6)
matriz## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 1 7 13 19 25 31
## [2,] 2 8 14 20 26 32
## [3,] 3 9 15 21 27 33
## [4,] 4 10 16 22 28 34
## [5,] 5 11 17 23 29 35
## [6,] 6 12 18 24 30 36base <- rast(matriz)Por si hemos olvidado a qué clase pertenece el objeto bas
class(base)## [1] "SpatRaster"
## attr(,"package")
## [1] "terra"La función crs() retorna el CRS de un objeto
SpatRaster.
crs(base)## [1] ""Este objeto ráster carece de proyección. Por ello, le añadiremos la proyección correspondiente al sistema WGS 84 / UTM zone 30N (“EPSG:32630”).
crs(base) <- "EPSG:4326"Cuando el objeto sí tiene proyección, la función
crs()[https://rspatial.github.io/terra/reference/crs.html]
devuelve una cadena de caracteres que contiene detalles sobre el sistema
de referencia de coordenadas (CRS) del ráster. Por defecto, la salida
está en formato de texto conocido (WKT).
crs(base)## [1] "GEOGCRS[\"WGS 84\",\n ENSEMBLE[\"World Geodetic System 1984 ensemble\",\n MEMBER[\"World Geodetic System 1984 (Transit)\"],\n MEMBER[\"World Geodetic System 1984 (G730)\"],\n MEMBER[\"World Geodetic System 1984 (G873)\"],\n MEMBER[\"World Geodetic System 1984 (G1150)\"],\n MEMBER[\"World Geodetic System 1984 (G1674)\"],\n MEMBER[\"World Geodetic System 1984 (G1762)\"],\n MEMBER[\"World Geodetic System 1984 (G2139)\"],\n ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n LENGTHUNIT[\"metre\",1]],\n ENSEMBLEACCURACY[2.0]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n CS[ellipsoidal,2],\n AXIS[\"geodetic latitude (Lat)\",north,\n ORDER[1],\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n AXIS[\"geodetic longitude (Lon)\",east,\n ORDER[2],\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n USAGE[\n SCOPE[\"Horizontal component of 3D system.\"],\n AREA[\"World.\"],\n BBOX[-90,-180,90,180]],\n ID[\"EPSG\",4326]]"El argumento parse = TRUE ofrece una salida más fácil
para su lectura en la pantalla.
crs(base, parse = TRUE)## [1] "GEOGCRS[\"WGS 84\","
## [2] " ENSEMBLE[\"World Geodetic System 1984 ensemble\","
## [3] " MEMBER[\"World Geodetic System 1984 (Transit)\"],"
## [4] " MEMBER[\"World Geodetic System 1984 (G730)\"],"
## [5] " MEMBER[\"World Geodetic System 1984 (G873)\"],"
## [6] " MEMBER[\"World Geodetic System 1984 (G1150)\"],"
## [7] " MEMBER[\"World Geodetic System 1984 (G1674)\"],"
## [8] " MEMBER[\"World Geodetic System 1984 (G1762)\"],"
## [9] " MEMBER[\"World Geodetic System 1984 (G2139)\"],"
## [10] " ELLIPSOID[\"WGS 84\",6378137,298.257223563,"
## [11] " LENGTHUNIT[\"metre\",1]],"
## [12] " ENSEMBLEACCURACY[2.0]],"
## [13] " PRIMEM[\"Greenwich\",0,"
## [14] " ANGLEUNIT[\"degree\",0.0174532925199433]],"
## [15] " CS[ellipsoidal,2],"
## [16] " AXIS[\"geodetic latitude (Lat)\",north,"
## [17] " ORDER[1],"
## [18] " ANGLEUNIT[\"degree\",0.0174532925199433]],"
## [19] " AXIS[\"geodetic longitude (Lon)\",east,"
## [20] " ORDER[2],"
## [21] " ANGLEUNIT[\"degree\",0.0174532925199433]],"
## [22] " USAGE["
## [23] " SCOPE[\"Horizontal component of 3D system.\"],"
## [24] " AREA[\"World.\"],"
## [25] " BBOX[-90,-180,90,180]],"
## [26] " ID[\"EPSG\",4326]]"El argumento `describe = TRUE devuelve un resumen abreviado que incluye el nombre de CRS y el código European Petroleum Specialty Group (EPSG)
crs(base, describe = TRUE)## name authority code area extent
## 1 WGS 84 EPSG 4326 World -180, 180, -90, 90► En otras ocasiones, el objeto SpatRáster puede ser creado diseñando primero un objeto vacío (sin datos) pero con dimensiones espaciales (número de filas y columnas), así como las coordenadas geográficas de los píxeles (con lo que ya se le adjunta un sistema de proyección).
base <- rast(ncols= 6, # Número de columnas del ráster
nrows= 6, # Número de filas
nlyrs=1, # Número de capas
xmax= -3, xmin= -4, # Coordenadas del eje X (considerémoslas la longitud)
ymax= 43, ymin= 42) # Coordenadas del eje y (idem latitud)
base## class : SpatRaster
## dimensions : 6, 6, 1 (nrow, ncol, nlyr)
## resolution : 0.1666667, 0.1666667 (x, y)
## extent : -4, -3, 42, 43 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (CRS84) (OGC:CRS84)A continuación, ese objeto ráster vacío se rellena con números
aleatorios generados por la función runiff dentro del rango
establecido por los valores min y max
values(base) <- as.integer(runif(ncell(base),
min = 1,
max=36))Para comprobar que la gestión ha sido correcta, se puede representar dicho objeto ráster.
plot(base, breaks=c(0, 8, 16, 24, 30, 36), col = terrain.colors(5))
text(base, digits=1)
Como hemos visto anteriormente, crs() también puede asignar un CRS a
un objeto SpatRaster, en este caso a una copia del objeto
base..
base_utm <- base
crs(base_utm) <- "EPSG:32630"A continuación debemos confirmar que el procedimiento ha tenido éxito
crs(base_utm)## [1] "PROJCRS[\"WGS 84 / UTM zone 30N\",\n BASEGEOGCRS[\"WGS 84\",\n ENSEMBLE[\"World Geodetic System 1984 ensemble\",\n MEMBER[\"World Geodetic System 1984 (Transit)\"],\n MEMBER[\"World Geodetic System 1984 (G730)\"],\n MEMBER[\"World Geodetic System 1984 (G873)\"],\n MEMBER[\"World Geodetic System 1984 (G1150)\"],\n MEMBER[\"World Geodetic System 1984 (G1674)\"],\n MEMBER[\"World Geodetic System 1984 (G1762)\"],\n MEMBER[\"World Geodetic System 1984 (G2139)\"],\n ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n LENGTHUNIT[\"metre\",1]],\n ENSEMBLEACCURACY[2.0]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n ID[\"EPSG\",4326]],\n CONVERSION[\"UTM zone 30N\",\n METHOD[\"Transverse Mercator\",\n ID[\"EPSG\",9807]],\n PARAMETER[\"Latitude of natural origin\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8801]],\n PARAMETER[\"Longitude of natural origin\",-3,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8802]],\n PARAMETER[\"Scale factor at natural origin\",0.9996,\n SCALEUNIT[\"unity\",1],\n ID[\"EPSG\",8805]],\n PARAMETER[\"False easting\",500000,\n LENGTHUNIT[\"metre\",1],\n ID[\"EPSG\",8806]],\n PARAMETER[\"False northing\",0,\n LENGTHUNIT[\"metre\",1],\n ID[\"EPSG\",8807]]],\n CS[Cartesian,2],\n AXIS[\"(E)\",east,\n ORDER[1],\n LENGTHUNIT[\"metre\",1]],\n AXIS[\"(N)\",north,\n ORDER[2],\n LENGTHUNIT[\"metre\",1]],\n USAGE[\n SCOPE[\"Navigation and medium accuracy spatial referencing.\"],\n AREA[\"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).\"],\n BBOX[0,-6,84,0]],\n ID[\"EPSG\",32630]]"La función project() reproyecta un
objeto SpatRaster a un nuevo CRS, en este caso ETRS89 / UTM zone
30N.
base_utm <- project(base, "EPSG:25830")
crs(base_utm)## [1] "PROJCRS[\"ETRS89 / UTM zone 30N\",\n BASEGEOGCRS[\"ETRS89\",\n ENSEMBLE[\"European Terrestrial Reference System 1989 ensemble\",\n MEMBER[\"European Terrestrial Reference Frame 1989\"],\n MEMBER[\"European Terrestrial Reference Frame 1990\"],\n MEMBER[\"European Terrestrial Reference Frame 1991\"],\n MEMBER[\"European Terrestrial Reference Frame 1992\"],\n MEMBER[\"European Terrestrial Reference Frame 1993\"],\n MEMBER[\"European Terrestrial Reference Frame 1994\"],\n MEMBER[\"European Terrestrial Reference Frame 1996\"],\n MEMBER[\"European Terrestrial Reference Frame 1997\"],\n MEMBER[\"European Terrestrial Reference Frame 2000\"],\n MEMBER[\"European Terrestrial Reference Frame 2005\"],\n MEMBER[\"European Terrestrial Reference Frame 2014\"],\n ELLIPSOID[\"GRS 1980\",6378137,298.257222101,\n LENGTHUNIT[\"metre\",1]],\n ENSEMBLEACCURACY[0.1]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n ID[\"EPSG\",4258]],\n CONVERSION[\"UTM zone 30N\",\n METHOD[\"Transverse Mercator\",\n ID[\"EPSG\",9807]],\n PARAMETER[\"Latitude of natural origin\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8801]],\n PARAMETER[\"Longitude of natural origin\",-3,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8802]],\n PARAMETER[\"Scale factor at natural origin\",0.9996,\n SCALEUNIT[\"unity\",1],\n ID[\"EPSG\",8805]],\n PARAMETER[\"False easting\",500000,\n LENGTHUNIT[\"metre\",1],\n ID[\"EPSG\",8806]],\n PARAMETER[\"False northing\",0,\n LENGTHUNIT[\"metre\",1],\n ID[\"EPSG\",8807]]],\n CS[Cartesian,2],\n AXIS[\"(E)\",east,\n ORDER[1],\n LENGTHUNIT[\"metre\",1]],\n AXIS[\"(N)\",north,\n ORDER[2],\n LENGTHUNIT[\"metre\",1]],\n USAGE[\n SCOPE[\"Engineering survey, topographic mapping.\"],\n AREA[\"Europe between 6°W and 0°W: Faroe Islands offshore; Ireland - offshore; Jan Mayen - offshore; Norway including Svalbard - offshore; Spain - onshore and offshore.\"],\n BBOX[35.26,-6,80.49,0.01]],\n ID[\"EPSG\",25830]]"Finalmente, puede representarse gráficamente.
plot(base_utm, legend = TRUE)
📝 Sistemas de Coordenadas de Referencia:
Los Sistemas de Coordenadas de Referencia serán tratados de manera más detallada más adelante.
`
► Otra de las posibilidades es importar un fichero en formato texto y
convertirlo en un fichero ráster. Este fichero debe contener al menos
dos columnas con las coordenadas x e y, y una
tercera con la variable. El procedimiento requiere importar el fichero
en formato *.csv y luego transformarlo. En este caso, se ha elegido un
fichero en el que la coma (,) es el separador entre columnas, y el punto
(.) separa los decimales (formato típico de los países
anglosajones).
lst_noche_csv <- read.csv("D:/G174_2025/LABORATORIO_1_Ficheros_raster_R/datos/lst_noche.csv",
header = TRUE,
sep = ",",
dec = ".")A continuación, el fichero se transforma en un
SpatRáster, que, como en el caso anterior, puede
representarse gráficamente.
lst_noche <- rast(lst_noche_csv, # Fichero de origen en formato *.csv
type = "xyz", # Estructura original del fichero
crs = "EPSG:32630") # Asignación del sistema proyección WGS84 30 UTM
plot(lst_noche)
► Sin embargo, la forma más habitual de trabajar los ficheros ráster en terra es importarlos en formato *.tif.
dem <- rast("D:/G174_2025/LABORATORIO_1_Ficheros_raster_R/datos/dem.tif")
plot(dem)
Información general sobre el objeto y metadatos
¿Cómo conocer la información que contiene un objeto SpatRaster? Fundamentalmente, llamando a dicho objeto ráster
base## class : SpatRaster
## dimensions : 6, 6, 1 (nrow, ncol, nlyr)
## resolution : 0.1666667, 0.1666667 (x, y)
## extent : -4, -3, 42, 43 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (CRS84) (OGC:CRS84)
## source(s) : memory
## name : lyr.1
## min value : 1
## max value : 35La información almacenada es la siguiente:
• class: tipo de objeto.
• dimensions: descripción geométrica del fichero (número de filas y columnas).
• resolution: tamaño (en este caso, en metros) de las celdas en eje x e y.
• extent: coordenadas extremas del área cubierta por el objeto ráster.
• coord.ref.: sistema de coordenadas de referencia (tipo de proyección).
• source(s): dónde reside la información (en archivo o memoria) del objeto ráster.
• name: nombre de la capa o variable.
• min value | max value: rango de los valores representados en el objeto ráster.
✅ EJERCICIO:
Llama a los objetos lst_noche y dem y
responde a las siguientes preguntas:
lst_noche► ¿Cuáles son las dimensiones del objeto lst_noche?
dim(lst_noche)► ¿Cuál es la resolución del objeto dem?
res(dem)► ¿Qué nombre recibe la única capa del objeto lst_noche?
names(lst_noche)► ¿En qué sistema de coordenadas está escrita la información del objeto dem?
crs(dem, describe = TRUE)Propiedades Geométricas
Las propiedades del objeto SpatRaster se pueden revisar mediante los mismos comandos con los que se trabajan las matrices:
dim(base); ncol(base); nrow(base)## [1] 6 6 1## [1] 6## [1] 6Otras características del objeto se pueden obtener con los siguientes
comandos. Por ejemplo, la resolución de los píxeles en
la dimensión X se consigue con terra::xres()
xres(base) # Tamaño eje x## [1] 0.1666667Y en la dimensión Y
yres(base) # Tamaño eje y## [1] 0.1666667En el caso de necesitar ambas dimensiones al mismo tiempo se utiliza
la función terra::res()
res(base) # Tamaño en ambos ejes## [1] 0.1666667 0.1666667El número total de píxeles
(terra::ncell())
ncell(base) ## [1] 36El número de capas terra::nlyr()
nlyr(base) # Número de capas## [1] 1La función terra::ext() devuelve la extensión (si no se
crea un objeto) o crea un objeto SpatExtent que contienen las
coordenadas geográficas del ráster, y podría servir más adelante para
seleccionar las dimensiones de un nuevo objeto ráster o incluso el marco
espacial para una selección geográfica.
ext(base) # Devuelve la extensión (objeto SpatExtent)## SpatExtent : -4, -3, 42, 43 (xmin, xmax, ymin, ymax)A continuación, se crea el objeto SpatExtent correspondiente al
objeto lst_noche y se confirma su clase.
base_ext <- ext(base)
class(base_ext)## [1] "SpatExtent"
## attr(,"package")
## [1] "terra"Las dimensiones concretas en las unidades del Sistema de Coordenadas de Referencia del objeto anterior se obtienen llamando a ese objeto SpatExtent
base_ext## SpatExtent : -4, -3, 42, 43 (xmin, xmax, ymin, ymax)# igual que
base_ext[1:4]## xmin xmax ymin ymax
## -4 -3 42 43✅ EJERCICIO:
► ¿Qué número de píxeles tienen los ejes X e Y del del objeto lst_noche?
dim(lst_noche)► ¿El número de columnas del objeto dem es superior a
1000?.
ncol(dem)► ¿Cuál de los dos objetos (lst_noche y
dem) tiene una mayor número de pixeles?
ncell(lst_noche)
ncell(dem)► ¿Cuál es la extensión de los objetos lst_noche y dem?
ext(lst_noche)
ext(dem)► ¿Cuál de los dos objetos tiene píxeles de menor resolución?
res(lst_noche)
res(dem)Resumen de la información numérica
La escritura del nombre de un objeto SpatRaster en la consola imprime información general sobre ese objeto. La función summary() proporciona algunas estadísticas descriptivas (mínimo, máximo, cuartiles, etc.). Otras estadísticas pueden ser calculadas con la función global().
La función terra::summary() proporciona una primera idea
de la distribución de los valores del objeto ráster. El argumento
size especifica el número de píxeles seleccionados
aleatoriamente para calcular los valores estadísticos.
summary(base, size = 1e6)## lyr.1
## Min. : 1.00
## 1st Qu.:11.25
## Median :21.00
## Mean :20.11
## 3rd Qu.:31.00
## Max. :35.00Esos mismos valores pueden ser calculados mediante estadísticas
globales, valores obtenidos a partir de todas las celdas de una banda,
mediante la función terra::global.
global(base, fun="min") # valor mínimo## min
## lyr.1 1global(base, fun= "max") # valor máximo## max
## lyr.1 35Otros estadísticos son la desviación típica (‘sd’), la suma de todos los valores (‘sum’), la diferencia entre los valores mínimo y máximo (“range”), raíz cuadrántica media (“rms”), el número de píxeles sin datos (“isNA”) o el de aquéllas que tienen datos (“notNA”).
Tampoco es necesario escribir el argumento fun.
global(base, 'mean') # Promedio de los valores## mean
## lyr.1 20.11111⚠️ ATENCIÓN:
En algunos casos si aplicamos directamente la función
terra::global a un objeto ráster, nos puede devolver
NaN (no hay resultado)
global(dem, 'sum') # Suma de todos los valores## sum
## dem NaNEsto puede deberse a la presencia de píxeles sin datos. Se puede verificar con la siguiente orden:
global(dem, "isNA")## isNA
## dem 3955Para evitarlos, simplemente se añade el argumento
na.rm=TRUE a la función
global(dem, 'sum', na.rm=TRUE) # Suma de todos los valores## sum
## dem 885547046La función acepta concatenar diferentes argumentos.
global(base, # objeto SpatRaster
c("sum", "mean", "sd"), # Estadísticos a calcular
na.rm=TRUE) # No tiene en cuenta los valores ausentes NA## sum mean sd
## lyr.1 724 20.11111 11.45869Dos funciones muy útiles para el manejo de objetos ráster son
terra::unique() y terra::freq().
unique
devuelven un vector con todos los elementos que no están duplicados.
unique(base)## lyr.1
## 1 1
## 2 2
## 3 3
## 4 6
## 5 8
## 6 9
## 7 12
## 8 14
## 9 16
## 10 17
## 11 18
## 12 20
## 13 22
## 14 24
## 15 26
## 16 28
## 17 31
## 18 32
## 19 34
## 20 35freq
construye una tabla de frecuencia para los valores contenidos (realiza
por defecto un redondeo a número entero y se puede controlar con el
parámetro digits).
freq(base, digits= 0)## layer value count
## 1 1 1 2
## 2 1 2 1
## 3 1 3 1
## 4 1 6 3
## 5 1 8 1
## 6 1 9 1
## 7 1 12 2
## 8 1 14 1
## 9 1 16 3
## 10 1 17 1
## 11 1 18 1
## 12 1 20 1
## 13 1 22 1
## 14 1 24 3
## 15 1 26 1
## 16 1 28 3
## 17 1 31 2
## 18 1 32 1
## 19 1 34 2
## 20 1 35 5Las estadísticas pueden ser visualizadas con funciones como
terra::hist() y terra::density().
hist(base) # Histograma
density(base) # Curva de densidad
✅ EJERCICIO:
El objeto lst_noche reproduce la temperatura de
superficie obtenida del sensor MODIS, correspondiente a un área del N de
España.
► ¿Cuál es la temperatura media de la imagen?
global(lst_noche, 'mean', na.rm=TRUE) ► ¿Cuáles son los valores máximo y mínimo?
global(lst_noche, c('min', 'max'), na.rm=TRUE) ► ¿Cuál es el valor (número entero, argumento
digits = 0) más frecuente?
freq(lst_noche,
digits= 0)Nombre de las capas
También podemos acceder a los nombres de las capas y modificarlas por claridad.
names(base)## [1] "lyr.1"A continuación se renombra la capa
names(base) <- 'capa_1'
base## class : SpatRaster
## dimensions : 6, 6, 1 (nrow, ncol, nlyr)
## resolution : 0.1666667, 0.1666667 (x, y)
## extent : -4, -3, 42, 43 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (CRS84) (OGC:CRS84)
## source(s) : memory
## name : capa_1
## min value : 1
## max value : 35✅ EJERCICIO:
► ¿Cómo se llama la única capa del objeto dem?.
names(dem)► Cambia su nombre por el de mdt y comprueba que se ha
realizado correctamente
names(dem) <- "mdt"
demIndexado y extracción de capas
Como se ha comentado anteriormente, los valores de los píxeles de un objeto ráster se almacenan internamente como un vector con sus elementos ordenados secuencialmente. El primer valor de la secuencia corresponde al valor situado en la esquina superior izquierda; a partir de él, la secuencia sigue hacia la derecha y hacia abajo. El último valor corresponde a la esquina inferior derecha, y es igual al número de celdas de todo el raster.
Atendiendo a este orden, se puede solicitar el valor de algunos píxeles concretos:
► La totalidad de los valores de un objeto SpatRaster puede consultarse de esta manera.
base[]## capa_1
## [1,] 31
## [2,] 16
## [3,] 8
## [4,] 2
## [5,] 18
## [6,] 1
## [7,] 6
## [8,] 9
## [9,] 35
## [10,] 24
## [11,] 3
## [12,] 24
## [13,] 28
## [14,] 16
## [15,] 35
## [16,] 14
## [17,] 1
## [18,] 34
## [19,] 24
## [20,] 35
## [21,] 22
## [22,] 35
## [23,] 34
## [24,] 12
## [25,] 6
## [26,] 20
## [27,] 35
## [28,] 28
## [29,] 16
## [30,] 31
## [31,] 6
## [32,] 32
## [33,] 28
## [34,] 26
## [35,] 12
## [36,] 17La función terra::values() produce el mismo
resultado
values(base)## capa_1
## [1,] 31
## [2,] 16
## [3,] 8
## [4,] 2
## [5,] 18
## [6,] 1
## [7,] 6
## [8,] 9
## [9,] 35
## [10,] 24
## [11,] 3
## [12,] 24
## [13,] 28
## [14,] 16
## [15,] 35
## [16,] 14
## [17,] 1
## [18,] 34
## [19,] 24
## [20,] 35
## [21,] 22
## [22,] 35
## [23,] 34
## [24,] 12
## [25,] 6
## [26,] 20
## [27,] 35
## [28,] 28
## [29,] 16
## [30,] 31
## [31,] 6
## [32,] 32
## [33,] 28
## [34,] 26
## [35,] 12
## [36,] 17► ¿Cúal es el valor del primer pixel (izquierda arriba; Celda con ID 1)
base[1]## capa_1
## 1 31► ¿Cuál es el valor del pixel 30 (5 fila, 6ª columna)
base[30]## capa_1
## 1 31► También podríamos solicitar varios píxeles, uno a continuación de otro
base[12:15]## capa_1
## 1 24
## 2 28
## 3 16
## 4 35Dado que el raster se organiza en filas y columnas, es mucho más
fácil extraer el valor de un pixel concreto señalando en qué fila y
columna se encuentra, mediante la sintaxis propia de R. En ella, se usan
los corchetes [ , ], teniendo en cuenta que el primer dígito antes de la
coma corresponde a las filas y el segundo, después de la coma, a las
columnas ([fila, columna]).
base[3, ] # Fila 3## capa_1
## 1 28
## 2 16
## 3 35
## 4 14
## 5 1
## 6 34base[ , 3] # Columna 3## capa_1
## 1 8
## 2 35
## 3 35
## 4 22
## 5 35
## 6 28base[ ,3:5] # Columnas 3 a 5## capa_1
## 1 8
## 2 2
## 3 18
## 4 35
## 5 24
## 6 3
## 7 35
## 8 14
## 9 1
## 10 22
## 11 35
## 12 34
## 13 35
## 14 28
## 15 16
## 16 28
## 17 26
## 18 12base[3, 5] # Píxel correspondiente a la 3ª fila y 5ª columna## capa_1
## 1 1base[3:5,3:5]## capa_1
## 1 35
## 2 14
## 3 1
## 4 22
## 5 35
## 6 34
## 7 35
## 8 28
## 9 16Los píxeles de cada uno de los vértices del ráster también se pueden consultar extrayendo primero el número de filas y columnas.
nr <- nrow(base)
nc <- ncol(base)
base[c(1, nr), c(1, nc)]## capa_1
## 1 31
## 2 1
## 3 6
## 4 17La sintaxis se expande para seleccionar todos los píxeles de la misma manera que podría abordar todos los valores de una matriz:
valores <- base[ , ]
valores## capa_1
## [1,] 31
## [2,] 16
## [3,] 8
## [4,] 2
## [5,] 18
## [6,] 1
## [7,] 6
## [8,] 9
## [9,] 35
## [10,] 24
## [11,] 3
## [12,] 24
## [13,] 28
## [14,] 16
## [15,] 35
## [16,] 14
## [17,] 1
## [18,] 34
## [19,] 24
## [20,] 35
## [21,] 22
## [22,] 35
## [23,] 34
## [24,] 12
## [25,] 6
## [26,] 20
## [27,] 35
## [28,] 28
## [29,] 16
## [30,] 31
## [31,] 6
## [32,] 32
## [33,] 28
## [34,] 26
## [35,] 12
## [36,] 17# equivalente a
valores <- base[]
valores## capa_1
## [1,] 31
## [2,] 16
## [3,] 8
## [4,] 2
## [5,] 18
## [6,] 1
## [7,] 6
## [8,] 9
## [9,] 35
## [10,] 24
## [11,] 3
## [12,] 24
## [13,] 28
## [14,] 16
## [15,] 35
## [16,] 14
## [17,] 1
## [18,] 34
## [19,] 24
## [20,] 35
## [21,] 22
## [22,] 35
## [23,] 34
## [24,] 12
## [25,] 6
## [26,] 20
## [27,] 35
## [28,] 28
## [29,] 16
## [30,] 31
## [31,] 6
## [32,] 32
## [33,] 28
## [34,] 26
## [35,] 12
## [36,] 17A continuación, se eliminan estos objetos
rm(valores)También se dispone de funciones para obtener las coordenadas a partir del número de celda
xFromCell(base,1) # Coordenada x de celda 1## x
## -3.916667yFromCell(base,1) # Coordenada y de celda 1## y
## 42.91667xyFromCell(base,1:3) # Vector de coordenadas de las celda 1 a 3## x y
## [1,] -3.916667 42.91667
## [2,] -3.750000 42.91667
## [3,] -3.583333 42.91667Y recíprocamente, a partir de coordenadas x e
yse obtiene el número de celda.
mat <- matrix(c(-3.75, 42.91667),
nrow = 1,
ncol = 2)
cellFromXY(base, mat) # Número de celda## [1] 2Con este tipo de numeración tenemos una serie de funciones que permiten obtener coordenadas reales a partir de la información de filas y columnas:
xFromCol(base,1) # Coordenada x de columna 1## [1] -3.916667yFromRow(base,1) # Coordenada y de fila 1## [1] 42.91667El inverso:
colFromX(base, -3.5) # Columna para coordenada x## [1] 4rowFromY(base, 42.5) # Fila para coordenada y## [1] 4Igualmente, se puede solicitar el pixel que corresponde a determinada fila y columna
colFromCell(base, 28) # A qué columna corresponde el pixel 28## [1] 4rowFromCell(base, c(1,5,15)) # A qué fila corresponen los píxeles 1, 5 y 15## [1] 1 1 3rowColFromCell(base, c(1,5,15)) # Cúales son las filas y columnas correspondientes a los píxeles## [,1] [,2]
## [1,] 1 1
## [2,] 1 5
## [3,] 3 3cellFromRowCol(base, 2, 2) # Cuál es el pixel correspondiente a la fila 2 y a la columna 2 ## [1] 8El intento de extraer un valor de una celda no existente retorna un NaN:
xFromCell(base,0) #id fuera del vector## x
## NaNxFromCell(base,82) #id fuera del vector## x
## NaN📝 ATENCIÓN:
Por ejemplo, utilizando el argumento drop=FALSE es
posible devolver un nuevo objeto raster con unas dimensiones
precisas
nuevo <- base[1:3,1:3,
drop=FALSE]
nuevo## class : SpatRaster
## dimensions : 3, 3, 1 (nrow, ncol, nlyr)
## resolution : 0.1666667, 0.1666667 (x, y)
## extent : -4, -3.5, 42.5, 43 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (CRS84) (OGC:CRS84)
## source(s) : memory
## name : capa_1
## min value : 6
## max value : 35✅ EJERCICIO:
► ¿Cuál es el valor del pixel 16 del objeto
lst_noche?
lst_noche[16]► ¿Y los valores de los píxeles 1142 a 1149 del objeto
dem?
dem[1142:1149]► ¿Y los valores correspondientes a los píxeles comprendidoes entre
las columnas 6-7 y las filas 2-3 del objeto lst_noche.
lst_noche[6:7,2:3]► ¿Cuáles son los píxeles de los vértices del objeto
lst_noche?
nr_lst_noche <- nrow(lst_noche)
nc_lst_noche <- ncol(lst_noche)
lst_noche[c(1, nr_lst_noche), c(1, nc_lst_noche)]► ¿Cuáles son las coordenadas correspondientes a las celdas 12 y 16
del objeto dem.
xyFromCell(dem, 12:16) Cálculos a partir de objetos SpatRaster.
Los valores de los píxeles de un objeto ráster se pueden modificar aplicando operaciones y funciones básicas.
q <- 100 * base^2
plot(q)
Por ejemplo, se puede crear una máscara binaria aplicando una expresión lógica
m <- q < 4000
plot(m)
También se pueden utilizar los corchetes para reemplazar los valores de un ráster. En el siguiente ejemplo todos los valores inferiores a 0.5 son reemplazados por el valor 0.5.
base[base < 20] <- 20
plot(base, range=c(1,36))
✅ EJERCICIO:
► Crea un nuevo objeto ráster a partir del objeto
lst_noche transformando los valores de temperatura de ºC a
ºKelvin (se debe sumar 273,15). Verifica que los resultados son
correctos
lst_noche_kelvin <- lst_noche +273.15
lst_noche_kelvinModificación del valor de uno o varios píxeles
El operador [] también puede utilizarse para modificar los valores de las celdas un objeto SpatRaster.
Por ejemplo, la modificación del valor de una celda
base_utm[1,1] <- -9999
base_utm[]## lyr.1
## [1,] -9999.000000
## [2,] 13.682685
## [3,] 5.416501
## [4,] 11.144663
## [5,] 5.887095
## [6,] 10.572392
## [7,] 16.121181
## [8,] 25.245052
## [9,] 11.780260
## [10,] 15.485473
## [11,] 20.038754
## [12,] 20.249407
## [13,] 27.151129
## [14,] 8.231881
## [15,] 22.441305
## [16,] 25.916086
## [17,] 27.027493
## [18,] 26.844051
## [19,] 21.972069
## [20,] 21.101221
## [21,] 18.127775
## [22,] 28.151371
## [23,] 29.611349
## [24,] 28.579096
## [25,] 21.913279
## [26,] 9.501218
## [27,] 26.771343
## [28,] 30.302029
## [29,] 19.953506
## [30,] 23.284414
## [31,] 11.290215
## [32,] 30.670319
## [33,] 27.077154
## [34,] 17.733528
## [35,] 15.597624Modificación de rangos de celdas
base_utm[3, c(1, 2)] <- -5555
base_utm[]## lyr.1
## [1,] -9999.000000
## [2,] 13.682685
## [3,] 5.416501
## [4,] 11.144663
## [5,] 5.887095
## [6,] 10.572392
## [7,] 16.121181
## [8,] 25.245052
## [9,] 11.780260
## [10,] 15.485473
## [11,] -5555.000000
## [12,] -5555.000000
## [13,] 27.151129
## [14,] 8.231881
## [15,] 22.441305
## [16,] 25.916086
## [17,] 27.027493
## [18,] 26.844051
## [19,] 21.972069
## [20,] 21.101221
## [21,] 18.127775
## [22,] 28.151371
## [23,] 29.611349
## [24,] 28.579096
## [25,] 21.913279
## [26,] 9.501218
## [27,] 26.771343
## [28,] 30.302029
## [29,] 19.953506
## [30,] 23.284414
## [31,] 11.290215
## [32,] 30.670319
## [33,] 27.077154
## [34,] 17.733528
## [35,] 15.597624plot(base_utm)
base_utm[1, 1:5] <- 0
base_utm[2, 1:5] <- 10
base_utm[3, 1:5] <- 15
base_utm[4, 1:5] <- 15
base_utm[5, 1:5] <- 20
base_utm[6, 1:5] <- 35Consulta de los valores
base[]## capa_1
## [1,] 31
## [2,] 20
## [3,] 20
## [4,] 20
## [5,] 20
## [6,] 20
## [7,] 20
## [8,] 20
## [9,] 35
## [10,] 24
## [11,] 20
## [12,] 24
## [13,] 28
## [14,] 20
## [15,] 35
## [16,] 20
## [17,] 20
## [18,] 34
## [19,] 24
## [20,] 35
## [21,] 22
## [22,] 35
## [23,] 34
## [24,] 20
## [25,] 20
## [26,] 20
## [27,] 35
## [28,] 28
## [29,] 20
## [30,] 31
## [31,] 20
## [32,] 32
## [33,] 28
## [34,] 26
## [35,] 20
## [36,] 20OBJETOS SpatRaster MULTIBANDA
A diferencia de lo que ocurría con el paquete precedente
(raster), terra() no requiere
diferenciar los archivos monobanda de los multibanda. Se puede crear un
objeto SpatRaster a partir de otros objetos o importando un
archivo (multibanda).
• Un ejemplo de la primera operación se realiza a continuación,
creando un objeto SpatRaster a partir de un
array tridimensional:
arr <- array(data = 1:6**3, # Datos a incluir en el array
dim = c(6, 6, 3)) # Dimensiones del arra, con 6 filas, 6 columnas y 3 "niveles".
arr## , , 1
##
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 1 7 13 19 25 31
## [2,] 2 8 14 20 26 32
## [3,] 3 9 15 21 27 33
## [4,] 4 10 16 22 28 34
## [5,] 5 11 17 23 29 35
## [6,] 6 12 18 24 30 36
##
## , , 2
##
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 37 43 49 55 61 67
## [2,] 38 44 50 56 62 68
## [3,] 39 45 51 57 63 69
## [4,] 40 46 52 58 64 70
## [5,] 41 47 53 59 65 71
## [6,] 42 48 54 60 66 72
##
## , , 3
##
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 73 79 85 91 97 103
## [2,] 74 80 86 92 98 104
## [3,] 75 81 87 93 99 105
## [4,] 76 82 88 94 100 106
## [5,] 77 83 89 95 101 107
## [6,] 78 84 90 96 102 108La función rast() convierte dicho objeto array en un
SpatRaster que contiene la misma información, pero formando
3 capas.
multi <- rast(arr)
multi## class : SpatRaster
## dimensions : 6, 6, 3 (nrow, ncol, nlyr)
## resolution : 1, 1 (x, y)
## extent : 0, 6, 0, 6 (xmin, xmax, ymin, ymax)
## coord. ref. :
## source(s) : memory
## names : lyr.1, lyr.2, lyr.3
## min values : 1, 37, 73
## max values : 36, 72, 108Para comprobar que se ha realizado correctamente la creación del objeto ráster multibanda se puede representar gráficamente.
plot(multi)
• No obstante, lo más habitual es añadir otro fichero en formato
ráster mediante la función c() (contatenar). En este caso,
se importará un nuevo fichero raster llamado lst_dia.tif,
transformándose en objeto SpatRáster, que se unirá al
objeto lst_noche, constituyendo un nuevo objeto que se
denominará modis_lst.
lst_dia <- rast("D:/G174_2025/LABORATORIO_1_Ficheros_raster_R/datos/lst_dia.tif")A continuación unimos ambos objetos en uno único usando la función
concatenar c().
lst <- c(lst_noche, lst_dia)
lst## class : SpatRaster
## dimensions : 27, 36, 2 (nrow, ncol, nlyr)
## resolution : 931.3632, 933.1531 (x, y)
## extent : 487488.3, 521017.4, 4779577, 4804772 (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 30N (EPSG:32630)
## sources : memory
## lst_dia.tif
## names : LST_noche, LST_dia
## min values : 6.956729, 19.14006
## max values : 12.535732, 29.10790Indexado de archivos multibanda.
La indexación de archivos multibanda es similar al de una lista en R.
Las capas se identifican con corchetes dobles [[]], que en
este caso devolverán una de las bandas del objeto ráster, no los valores
de los píxeles. En este caso, se va a seleccionar la segunda capa.
multi[[2]]## class : SpatRaster
## dimensions : 6, 6, 1 (nrow, ncol, nlyr)
## resolution : 1, 1 (x, y)
## extent : 0, 6, 0, 6 (xmin, xmax, ymin, ymax)
## coord. ref. :
## source(s) : memory
## name : lyr.2
## min value : 37
## max value : 72La consulta de más de una capa simultáneamente también se realiza con corchetes dobles:
► Con los números de las capas
multi[[1:2]]## class : SpatRaster
## dimensions : 6, 6, 2 (nrow, ncol, nlyr)
## resolution : 1, 1 (x, y)
## extent : 0, 6, 0, 6 (xmin, xmax, ymin, ymax)
## coord. ref. :
## source(s) : memory
## names : lyr.1, lyr.2
## min values : 1, 37
## max values : 36, 72► O con los nombres de cada capa.
multi[[c("lyr.2", "lyr.3")]]## class : SpatRaster
## dimensions : 6, 6, 2 (nrow, ncol, nlyr)
## resolution : 1, 1 (x, y)
## extent : 0, 6, 0, 6 (xmin, xmax, ymin, ymax)
## coord. ref. :
## source(s) : memory
## names : lyr.2, lyr.3
## min values : 37, 73
## max values : 72, 108Si las capas poseen un nombre, también se pueden seleccionar con nombre tras el signo $
multi$lyr.1## class : SpatRaster
## dimensions : 6, 6, 1 (nrow, ncol, nlyr)
## resolution : 1, 1 (x, y)
## extent : 0, 6, 0, 6 (xmin, xmax, ymin, ymax)
## coord. ref. :
## source(s) : memory
## name : lyr.1
## min value : 1
## max value : 36Para extraer capas del SpatRaster
original, puede utilizar el signo –c o bien la función
terra::subset(), como muestran los siguientes ejemplos.
multi.extraido <- multi[[-c(1,2)]]
multi.extraido## class : SpatRaster
## dimensions : 6, 6, 1 (nrow, ncol, nlyr)
## resolution : 1, 1 (x, y)
## extent : 0, 6, 0, 6 (xmin, xmax, ymin, ymax)
## coord. ref. :
## source(s) : memory
## name : lyr.3
## min value : 73
## max value : 108multi.extraido2 <- subset(multi, 3)Y también es equivalente a
multi.extraido3 <- subset(multi, "lyr.3")También se pueden añadir más capas a un objeto raster. En el
siguiente ejemplo, los valores de la capa layer.1 se
dividen entre dos, y esta nueva capa se agrega al objeto original:
nuevaCapa <- multi[[1]]/2
nuevaCapa## class : SpatRaster
## dimensions : 6, 6, 1 (nrow, ncol, nlyr)
## resolution : 1, 1 (x, y)
## extent : 0, 6, 0, 6 (xmin, xmax, ymin, ymax)
## coord. ref. :
## source(s) : memory
## name : lyr.1
## min value : 0.5
## max value : 18.0names(nuevaCapa) <- "nueva" # Cambio de nombre para que no coincida con la capa originalA continuación, se añadirá el objeto ráster multi la
capa nuevaCapa con la función terra::add()
add(multi) <- nuevaCapa
multi## class : SpatRaster
## dimensions : 6, 6, 4 (nrow, ncol, nlyr)
## resolution : 1, 1 (x, y)
## extent : 0, 6, 0, 6 (xmin, xmax, ymin, ymax)
## coord. ref. :
## source(s) : memory
## names : lyr.1, lyr.2, lyr.3, nueva
## min values : 1, 37, 73, 0.5
## max values : 36, 72, 108, 18.0La sintaxis anterior es equivalente a
multi[["miNuevaBanda"]] <- nuevaCapa
multi## class : SpatRaster
## dimensions : 6, 6, 5 (nrow, ncol, nlyr)
## resolution : 1, 1 (x, y)
## extent : 0, 6, 0, 6 (xmin, xmax, ymin, ymax)
## coord. ref. :
## source(s) : memory
## names : lyr.1, lyr.2, lyr.3, nueva, miNuevaBanda
## min values : 1, 37, 73, 0.5, 0.5
## max values : 36, 72, 108, 18.0, 18.0Al consultar los píxeles con corchetes simples, devolverá una matriz con una columna por capa y píxeles en filas:
multi[1:2, 3]## lyr.1 lyr.2 lyr.3 nueva miNuevaBanda
## 1 13 49 85 6.5 6.5
## 2 14 50 86 7.0 7.0Para consultar ciertos píxeles en una capa específica, se combinan los métodos descritos anteriormente, especificando primero la capa y luego la posición dentro de esa lista. Por ejemplo, para obtener la capa 2, las filas 2-4 y la columna 3:
multi[[2]][2:4, 3]## lyr.2
## 1 50
## 2 51
## 3 52Estos comandos no solo son útiles para consultar valores ráster; también se pueden usar al revés para establecer o cambiar valores.
borrame <- multi[[3]]
plot(borrame)
text(borrame, digits=1)
Podemos señalar en dicho objeto ráster los valores de un conjunto de filas y columnas determinadas para su posterior modificación.
borrame[2:3, 2:3] ## lyr.3
## 1 80
## 2 86
## 3 81
## 4 87Estos valores van a ser convertidos en NA, es decir, NoData.
borrame[2:3, 2:3] <- NA
plot(borrame)
text(borrame, digits=1)
Para cambiar a 0 todos los valores en borrame podríamos
hacer lo siguiente:
borrame[] <- 0
plot(borrame)
Para extraer una o varias capas del
SpatRaster original, se señala la(s) capa(s) entre
corchetes [[]], precedidas por el signo negativo
– (en el segundo caso, el signo negativo puede trasladarse
delante de -c), como muestran los siguientes ejemplos.
• Capa única
multi.extraido <- -multi[[2]]
multi.extraido## class : SpatRaster
## dimensions : 6, 6, 1 (nrow, ncol, nlyr)
## resolution : 1, 1 (x, y)
## extent : 0, 6, 0, 6 (xmin, xmax, ymin, ymax)
## coord. ref. :
## source(s) : memory
## name : lyr.2
## min value : -72
## max value : -37• Varias capas
multi.extraido <- -multi[[c(1,2)]]
multi.extraido## class : SpatRaster
## dimensions : 6, 6, 2 (nrow, ncol, nlyr)
## resolution : 1, 1 (x, y)
## extent : 0, 6, 0, 6 (xmin, xmax, ymin, ymax)
## coord. ref. :
## source(s) : memory
## names : lyr.1, lyr.2
## min values : -36, -72
## max values : -1, -37Si se desea eliminar capas de un objeto
multi.extraido <- multi.extraido[[2]]
multi.extraido## class : SpatRaster
## dimensions : 6, 6, 1 (nrow, ncol, nlyr)
## resolution : 1, 1 (x, y)
## extent : 0, 6, 0, 6 (xmin, xmax, ymin, ymax)
## coord. ref. :
## source(s) : memory
## name : lyr.2
## min value : -72
## max value : -37Uso de la función terra::subset
para extraer capas para un nuevo objeto.
multi.reducido <- subset(multi, c(1,2))
multi.reducido## class : SpatRaster
## dimensions : 6, 6, 2 (nrow, ncol, nlyr)
## resolution : 1, 1 (x, y)
## extent : 0, 6, 0, 6 (xmin, xmax, ymin, ymax)
## coord. ref. :
## source(s) : memory
## names : lyr.1, lyr.2
## min values : 1, 37
## max values : 36, 72Esto último es equivalente a
multi.reducido <- subset(multi, c("lyr.1","lyr.3"))
multi.reducido## class : SpatRaster
## dimensions : 6, 6, 2 (nrow, ncol, nlyr)
## resolution : 1, 1 (x, y)
## extent : 0, 6, 0, 6 (xmin, xmax, ymin, ymax)
## coord. ref. :
## source(s) : memory
## names : lyr.1, lyr.3
## min values : 1, 73
## max values : 36, 108Si se desea eliminar capas de ese objeto
multi.reducido <- subset(multi, -c(1,2))
multi.reducido## class : SpatRaster
## dimensions : 6, 6, 3 (nrow, ncol, nlyr)
## resolution : 1, 1 (x, y)
## extent : 0, 6, 0, 6 (xmin, xmax, ymin, ymax)
## coord. ref. :
## source(s) : memory
## names : lyr.3, nueva, miNuevaBanda
## min values : 73, 0.5, 0.5
## max values : 108, 18.0, 18.0También se pueden añadir más capas a un objeto
raster. En el siguiente ejemplo, los valores de la capa
layer.1 se dividen entre dos, y esta nueva capa se agrega
al objeto orignal con add(). Alternativamente, también se
puede especificar un nuevo nombre de capa entre corchetes.
nuevaCapa <- multi[[1]]/2
names(nuevaCapa) <- "nueva" # Cambio de nombre para que no coincida con la capa original
nuevaCapa## class : SpatRaster
## dimensions : 6, 6, 1 (nrow, ncol, nlyr)
## resolution : 1, 1 (x, y)
## extent : 0, 6, 0, 6 (xmin, xmax, ymin, ymax)
## coord. ref. :
## source(s) : memory
## name : nueva
## min value : 0.5
## max value : 18.0A continuación, se añadirá el objeto ráster multi la
capa nuevaCapa.
add(multi) <- nuevaCapa # Nueva banda que se añader
multi## class : SpatRaster
## dimensions : 6, 6, 6 (nrow, ncol, nlyr)
## resolution : 1, 1 (x, y)
## extent : 0, 6, 0, 6 (xmin, xmax, ymin, ymax)
## coord. ref. :
## source(s) : memory
## names : lyr.1, lyr.2, lyr.3, nueva, miNuevaBanda, nueva
## min values : 1, 37, 73, 0.5, 0.5, 0.5
## max values : 36, 72, 108, 18.0, 18.0, 18.0La sintaxis anterior es equivalent a
multi[["miNuevaBanda"]] <- nuevaCapa
multi## class : SpatRaster
## dimensions : 6, 6, 6 (nrow, ncol, nlyr)
## resolution : 1, 1 (x, y)
## extent : 0, 6, 0, 6 (xmin, xmax, ymin, ymax)
## coord. ref. :
## source(s) : memory
## names : lyr.1, lyr.2, lyr.3, nueva, miNuevaBanda, nueva
## min values : 1, 37, 73, 0.5, 0.5, 0.5
## max values : 36, 72, 108, 18.0, 18.0, 18.0Modificar valores de píxeles concretos
Al consultar los píxeles con corchetes simples, devolverá una matriz con una columna por capa y píxeles en filas:
multi[1:2, 3]## lyr.1 lyr.2 lyr.3 nueva miNuevaBanda nueva
## 1 13 49 85 6.5 6.5 6.5
## 2 14 50 86 7.0 7.0 7.0Para consultar ciertos píxeles en una capa específica, se combinan los métodos descritos anteriormente, especificando primero la capa y luego la posición dentro de esa lista. Por ejemplo, para obtener la capa 2, las filas 2-4 y la columna 3:
multi[[2]][2:4, 3]## lyr.2
## 1 50
## 2 51
## 3 52Estos comandos no solo son útiles para consultar valores ráster; también se pueden usar al revés para establecer o cambiar valores.
copia <- multi[[3]]
plot(copia)
text(copia, digits=1)
Podemos señalar en dicho objeto ráster los valores de un conjunto de filas y columnas determinadas para su posterior modificación.
copia[2:3, 2:3] ## lyr.3
## 1 80
## 2 86
## 3 81
## 4 87Estos valores van a ser convertidos en NA, es decir, NoData.
copia[2:3, 2:3] <- NA
plot(copia)
text(copia, digits=1)
Para cambiar a 0 todos los valores podríamos hacer lo siguiente:
copia[] <- 0
plot(copia)
Cálculos estadísticos entre varias bandas
De la misma manera que se realizó con los objetos monobanda, es posible realizar determinados cálculos estadísticos entre varias bandas, mediante la función terra::app.
A continuación se muestran algunos ejemplos. En el primer caso se
calcula la media de las dos bandas que componen el objeto
lst.
app(lst, fun=sum)## class : SpatRaster
## dimensions : 27, 36, 1 (nrow, ncol, nlyr)
## resolution : 931.3632, 933.1531 (x, y)
## extent : 487488.3, 521017.4, 4779577, 4804772 (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 30N (EPSG:32630)
## source(s) : memory
## name : sum
## min value : 29.05681
## max value : 39.62547sum(lst) # equivalente a## class : SpatRaster
## dimensions : 27, 36, 1 (nrow, ncol, nlyr)
## resolution : 931.3632, 933.1531 (x, y)
## extent : 487488.3, 521017.4, 4779577, 4804772 (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 30N (EPSG:32630)
## source(s) : memory
## name : sum
## min value : 29.05681
## max value : 39.62547Igualmente, se puede obtener el valor máximo de todos los píxeles de cada objeto ráster.
app(lst, max)## class : SpatRaster
## dimensions : 27, 36, 1 (nrow, ncol, nlyr)
## resolution : 931.3632, 933.1531 (x, y)
## extent : 487488.3, 521017.4, 4779577, 4804772 (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 30N (EPSG:32630)
## source(s) : memory
## name : max
## min value : 19.14006
## max value : 29.10790GRABACIÓN DE OBJETOs SpatRaster
Para grabar los objetos creatos o modificados por R existen dos posibilidades
Grabar en formatos geoespaciales habituales
Para guardar objetos raster en el disco duro es necesario utilizar la función writeRaster, siendo posible elegir diferentes opciones para almacenar los datos (más información en ?writeFormats).
► Si se intercambian frecuentemente los ficheros con otros software
SIG es recomendable el uso del formato GeoTIFF, pero leer y escribir
archivos GeoTiff en R es más lento que en el formato
*.grd.
► Este último formato crea un archivo de texto con un encabezado (.grd) que contiene información como nombres de proyección o de capa, y un archivo binario (.gri), que contiene los datos reales. Este formato es recomendable cuando la mayor parte del trabajo se realice sobre R, ya que no requiere ningún controlador externo y, lo que es más importante, mantiene el nombre de las capas originales. Su inconveniente es que no está comprimido, por lo tanto, los archivos “pesan” mucho, si bien, si no se van a utilizar con frecuencia, se pueden comprimir razonablemente con cualquier software de compresión de archivos. Por otro lado, no todos los programas son capaces de leer el formato .grd.
Otra precaución es la inclusión del argumento
overwrite=, que sobreescribe la última versión del objeto
ráster. Si deseamos sobreescribir, pero olvidamos incluir este argumento
en la sintaxis, proporcionará un mensaje de error. No obstante, hay que
tener cuidado, precisamente por su capacidad para reemplazar versiones
anteriores del mismo raster.
writeRaster(base, "D:/G174_2025/LABORATORIO_1_Ficheros_raster_R/resultados/base.tif", overwrite= TRUE)Guardar en otros formatos
writeRaster(base, "mi_raster.asc", overwrite=TRUE) # ASCII Grid
writeRaster(base, "mi_raster.img", overwrite=TRUE) # Erdas Imagine
writeRaster(base, "mi_raster.nc", overwrite=TRUE) # NetCDF
writeRaster(base, "mi_raster.gpkg", overwrite=TRUE) # GeoPackageIgualmente,
writeRaster(base,
"mi_raster.tif",
filetype="GTiff", # Tipo de formato
datatype="INT2S", # Define el tipo de dato (ej. FLT4S para flotantes).
gdal="COMPRESS=LZW", # Aplica compresión para reducir el tamaño del archivo.
overwrite=TRUE)Grabar objetos ráster con los formatos propios de R.
Las funciones habituales de R para guardar datos (como
save() o saveRDS ()) no son apropiadas para
trabajar con objetos en formato ráster. Esto se debe a que, por lo
general, los objetos ráster creados en R realmente no almacenan datos,
sólo algunos metadatos y la ruta a los archivos originales. Esto es un
“atajo” para el manejo de ficheros que “pesan” mucho y pueden causar
problemas si la memoria disponible en un ordenador personal no es
suficiente. Una vez que se realiza un cálculo sobre un objeto ráster, R
carga los trozos o mosaicos más pequeños en la memoria, realiza los
cálculos y vuelve a escribir los resultados en un archivo. Solo los
rásteres pequeños se leerán completamente en la memoria.
Este procedimiento tiene un inconveniente: si los archivos se cambian de directorio o se eliminan del directorio temporal, el enlace se romperá y el objeto raster quedará inutilizable.
Grabar un objeto en formato *.rds.
Para no destruir el objeto ráster original base lo
duplicamos.
base_bis <- baseGraba un único objeto.
saveRDS(base_bis, file = "D:/G174_2025/LABORATORIO_1_Ficheros_raster_R/resultados/datos_provisionales.rds")Carga en R ese objeto.
readRDS(file = "D:/G174_2025/LABORATORIO_1_Ficheros_raster_R/resultados/datos_provisionales.rds")## class : SpatRaster
## dimensions : 6, 6, 1 (nrow, ncol, nlyr)
## resolution : 0.1666667, 0.1666667 (x, y)
## extent : -4, -3, 42, 43 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (CRS84) (OGC:CRS84)
## source(s) : memory
## name : capa_1
## min value : 20
## max value : 35Grabar objetos en formato *.RData.
Graba un único objeto
save(base_bis, file = "D:/G174_2025/LABORATORIO_1_Ficheros_raster_R/resultados/datos_provisionales.RData")Graba múltiples objetos
save(base_bis, lst, file = "D:/G174_2025/LABORATORIO_1_Ficheros_raster_R/resultados/datos_provisionales_2.RData")Para cargar ese fichero
load("D:/G174_2025/LABORATORIO_1_Ficheros_raster_R/resultados/datos_provisionales_2.RData")Grabar todo el espacio de trabajo de manera simultánea
Es recomendable en sesiones de trabajo de gran longitud, para lo cual
se utiliza la función save.image().
save.image(file = "D:/G174_2025/LABORATORIO_1_Ficheros_raster_R/resultados/mi_espacio_trabajo.RData")Para restaurar todo este espacio de trabajo
load("D:/G174_2025/LABORATORIO_1_Ficheros_raster_R/resultados/mi_espacio_trabajo.RData")Si quisiéramos eliminar todo el Global Environment
rm(list = ls())