OBJETIVOS

MATERIALES

PAQUETES DE R NECESARIOS PARA LA ACTIVIDAD

1 INTRODUCCION

El formato habitual de las imágenes de satélite es el formato ráster, uno de los más comunes empleados en los Sistemas de Información Geográfica.

Este formato divide la imagen en celdas regulares (cuadrículas o píxeles), que constituyen la unidad de análisis espacial básica.

Imagen raster
Imagen raster

Aunque hay un gran número de formatos que soportan la información ráster, los más habituales son los JPG, TIFF, ECW, MrSID, JP2, IMG, BIL, GRID o ASC. Esta diversidad tiene su explicación por la diferente capacidad de almacenar datos en los píxeles, también conocida como profundidad de bits por píxel, que influye a la hora de representar la información. Cuanto mayor sea, los archivos incluirán una mayor gama de valores. Cada uno de esos formatos puede almacenar un número determinado de valores en cada píxel, de acuerdo con la siguiente fórmula: \(2^n\).

Para ilustrar este concepto, debe tenerse en cuenta que 1 bit almacena dos colores y, por tanto, dos valores; 8 bits almacenan 256 valores de colores. En niveles superiores, 16 bits albergan más de 65000 valores, y 24 bits permiten trabajar con 16 millones de colores. Un archivo JPG alberga píxeles de 8 bits mientras que archivos TIFF permiten trabajar hasta 24 bits.

Además, algunos formatos tampoco incluyen valores NoData, lo cual puede complicar el análisis en los márgenes de las imágenes, particularmente cuando están formadas por mosaicos superpuestos.

Por último, algunos formatos no incorporan información espacial en las cabeceras de los archivos, mientras que otros sí lo hacen, obligando a que la imagen esté acompañada de ficheros adicionales, lo que complica la estructura de archivos. Los formatos que no contienen en su cabecera información espacial acompañan el archivo de imagen por uno o varios archivos complementarios encargados de localizar espacialmente la imagen. Este es el caso de los archivos auxiliares como TFW (para algunos archivos TIF) o JGW (para los archivos JPG).

Otros formatos, como los archivos GRID, aparecen con una estructura de archivos basados en carpetas. De esta forma, un archivo GRID está formado por dos carpetas, una de ellas denominada siempre como INFO y otra denominada con el nombre de referencia del archivo ráster. Ambas carpetas han de trabajarse en conjunto y son dependientes entre sí. Separarlas y compartirlas excluyendo la carpeta INFO es uno de los fallos más comunes cuando no tenemos conocimiento del manejo de este tipo de archivos, especialmente cuando la carpeta INFO es compartida por múltiples archivos.

Dependiendo del número de bandas incluidas en un archivo ráster, se pueden distinguir

2 FICHEROS RASTER CON R

A continuación se debe activar el paquete raster. Puede ocurrir que no se recuerde si se ha instalado previamente dicho paquete. En ese caso, se puede preguntar a R que verifique si está instalado con la función if (!require()); en caso afirmativo, no lo vuelve a instalar.

if (!require("raster")) install.packages("raster")
## Loading required package: raster
## Warning: package 'raster' was built under R version 4.3.2
## Loading required package: sp
## Warning: package 'sp' was built under R version 4.3.2
library(raster)

Aunque existen paquetes más actuales y potentes para el tratamiento de objetos en formato ráster (pe. “Terra” y “Stars”), su grado de compatibilidad con otros paquetes dedicados al análisis geoespacial no es tan grande como el veterano paquete raster. Además, la mayoría de las páginas web sobre esta temática, que pueden servir de complemento a las actividades desarrolladas aquí aún usan dicho paquete.

Precaución

En R pueden utilizarse dos notaciones para determinar la dirección de un directorio:

  • Mediante dos barras inclinadas a la derecha (\) .
  • Mediante una única barra, inclinada a la izquierda (/).

La barra única inclinada a la derecha () se usa en R como tecla de escape.

2.1 Objetos ráster monobada (RasterLayer)

El paquete raster proporciona varias funciones para leer archivos raster.

  • Para cargar una única capa (archivo monobanda) se emplea la función raster::raster()`.

  • Para cargar simultáneamente varias capas (archivo multibanda se usan las funciones raster::stack() y raster::brick(). Volveremos a ellas más adelantes.

A continuación se van a mostrar algunas operaciones básicas para el manejo de ficheros ráster en R, que luego serán utilizadas para el análisis de imágenes de satélite.

2.1.1 Operaciones básicas

El primer paso es la creación de un objeto ráster, denominado base. El procedimiento para crear este objeto puede realizarse de varias maneras.

Una de ellas, por ejemplo implica crear inicialmente una matriz, que luego se convierte en objeto ráster.

mat <- matrix(1:36, 6, 6)
mat
##      [,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
base <- raster(mat)
base
## class      : RasterLayer 
## dimensions : 6, 6, 36  (nrow, ncol, ncell)
## resolution : 0.1666667, 0.1666667  (x, y)
## extent     : 0, 1, 0, 1  (xmin, xmax, ymin, ymax)
## crs        : NA 
## source     : memory
## names      : layer 
## values     : 1, 36  (min, max)

En otras ocasiones, el ráster puede ser creado estableciendo primero el marco espacial: sus dimensiones (número de filas y columnas), así como las coordenadas geográficas de los píxeles.

base <- raster(ncol= 6,                   # Número de columnas del ráster
               nrow= 6,                   # Número de filas
               xmx= -3,xmn= -4,           # Coordenadas del eje X (considerémoslas la longitud)
               ymx= 43,ymn= 42)           # Coordenadas del eje y (idem latitud)
base
## class      : RasterLayer 
## dimensions : 6, 6, 36  (nrow, ncol, ncell)
## resolution : 0.1666667, 0.1666667  (x, y)
## extent     : -4, -3, 42, 43  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs

En segundo lugar, dado que el ráster no contiene ningún dato, ese marco espacial 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)

¿Qué información contiene dicho objeto?

base
## class      : RasterLayer 
## dimensions : 6, 6, 36  (nrow, ncol, ncell)
## resolution : 0.1666667, 0.1666667  (x, y)
## extent     : -4, -3, 42, 43  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : memory
## names      : layer 
## values     : 1, 35  (min, max)

La información almacenada que describe a continuación:

• class: tipo de objeto.

• dimensions: descripción geométrica del modelo con su 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.

• crs: sistema de coordenadas de referencia (tipo de proyección).

• source: dónde reside la información (en archivo o memoria) del ráster.

• names: nombre de la capa o variable.

• values: rango de valores del objeto.

2.1.2 Estructura Interna slots

El objeto RasterLayer creado en el punto anterior, internamente asocia una serie de atributos y propiedades que en R se denominan slots (ranuras). Cada slot tiene un nombre y una clase asociada. La extracción de un slot devuelve un objeto de esa clase.

Debe usarse el operador @ (arroba) para extraer o modificar el contenido de un slot (función que reside en el paquete base). Por ejemplo, para listar los slots que componen un objeto debemos usar:

slotNames(base)
##  [1] "file"     "data"     "legend"   "title"    "extent"   "rotated" 
##  [7] "rotation" "ncols"    "nrows"    "crs"      "srs"      "history" 
## [13] "z"

Y para revisar el contenido de uno de esos slots basta usar su nombre:

slot(base, 
     name = 'extent')
## class      : Extent 
## xmin       : -4 
## xmax       : -3 
## ymin       : 42 
## ymax       : 43

Alternativamente, se usa la sintaxis objeto@nombre

base@file
## An object of class ".RasterFile"
## Slot "name":
## [1] ""
## 
## Slot "datanotation":
## [1] "FLT4S"
## 
## Slot "byteorder":
## [1] "little"
## 
## Slot "nodatavalue":
## [1] -Inf
## 
## Slot "NAchanged":
## [1] FALSE
## 
## Slot "nbands":
## [1] 1
## 
## Slot "bandorder":
## [1] "BIL"
## 
## Slot "offset":
## [1] 0
## 
## Slot "toptobottom":
## [1] TRUE
## 
## Slot "blockrows":
## [1] 0
## 
## Slot "blockcols":
## [1] 0
## 
## Slot "driver":
## [1] ""
## 
## Slot "open":
## [1] FALSE

2.1.3 Propiedades Geométricas

El objeto rasterlayer posee una serie de propiedades que se pueden revisar mediante los mismos comandos con los que se trabaja con matrices:

dim(base); ncol(base); nrow(base)
## [1] 6 6 1
## [1] 6
## [1] 6

Y algunos especialmente diseñado para él:

xres(base)                       # Tamaño eje x
## [1] 0.1666667
yres(base) #Tamaño eje y
## [1] 0.1666667
res(base) #Tamaño en ambos ejes
## [1] 0.1666667 0.1666667
ncell(base) #Número total de elementos
## [1] 36
extent(base) #retorna objeto extensión
## class      : Extent 
## xmin       : -4 
## xmax       : -3 
## ymin       : 42 
## ymax       : 43

2.1.4 Estadísticas Generales

Los valores contenidos en el objeto ráster pueden ser explorados utilizando estadísticos globales mediante la función raster::cellStats:

cellStats(base,stat=min) #valor mínimo
## [1] 1
cellStats(base,stat=max) #valor máximo
## [1] 35

La función acepta diferentes argumentos que funcionan más rápido en ráster grandes.

cellStats(base, 'mean') #Promedio de los valores
## [1] 20.30556
cellStats(base,'sd') #Desviación estándar
## [1] 10.44665
cellStats(base, 'sum') #Sumatoria
## [1] 731

Para extraer el contenido y su frecuencia existen dos funciones, unique y freq. La primera devuelven un vector con todos los elementos que no están duplicados.

unique(base)
##  [1]  1  2  3  5  8  9 11 12 13 14 15 16 17 20 21 24 26 27 28 29 30 31 35

La segunda 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)
##       value count
##  [1,]     1     1
##  [2,]     2     1
##  [3,]     3     1
##  [4,]     5     2
##  [5,]     8     1
##  [6,]     9     1
##  [7,]    11     1
##  [8,]    12     2
##  [9,]    13     1
## [10,]    14     1
## [11,]    15     1
## [12,]    16     1
## [13,]    17     2
## [14,]    20     1
## [15,]    21     1
## [16,]    24     1
## [17,]    26     3
## [18,]    27     1
## [19,]    28     1
## [20,]    29     5
## [21,]    30     1
## [22,]    31     3
## [23,]    35     3

2.1.5 Nombre de las capas

También podemos acceder a los nombres de las capas y modificarlas por claridad.

base
## class      : RasterLayer 
## dimensions : 6, 6, 36  (nrow, ncol, ncell)
## resolution : 0.1666667, 0.1666667  (x, y)
## extent     : -4, -3, 42, 43  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : memory
## names      : layer 
## values     : 1, 35  (min, max)
names(base) <- 'mi_variable'
base
## class      : RasterLayer 
## dimensions : 6, 6, 36  (nrow, ncol, ncell)
## resolution : 0.1666667, 0.1666667  (x, y)
## extent     : -4, -3, 42, 43  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : memory
## names      : mi_variable 
## values     : 1, 35  (min, max)

2.1.6 Indexado de un objeto RasterLayer (monobanda)

Un objeto RasterLayer es un objeto ráster con una única banda. 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 número es igual al número de celdas de todo el raster.

¿Cúal es el valor del primer pixel (izquierda arriba)

base[1]
## [1] 13

¿Cuál es el valor del pixel 30 (5 fila, sexta columna)

base[30]
## [1] 28

También podríamos solicitar varios píxeles, uno a continuación de otro

base[12:15]
## [1] 20 31 35 31

Dado 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.

base[3,  ]                                                    # Fila 3
## [1] 31 35 31 31 30 27
base[  , 3]                                                   # Columna 3
## [1] 16 11 31  3  5  8
base[ ,3:5]                                                   # Columnas 3 a 5
##  [1] 16 11 31  3  5  8 15 21 31 17 26 26 17 14 30 35  1  5
base[3, 5]                                                    # Píxel correspondiente a la 3ª fila y 5ª columna
## [1] 30
base[3:5,3:5]
## [1] 31 31 30  3 17 35  5 26  1

Los 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. Los resultados, sin embargo, ofrecen valores NA, como se ha señalado anteriormente.

nr <- nrow(base)
nc <- ncol(base) 

base[c(1, nr), c(1, nc)]
## [1] 13 12 35 12

La 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
##  [1] 13 24 16 15 17 12 29  2 11 21 14 20 31 35 31 31 30 27 29 26  3 17 35 29  9
## [26] 29  5 26  1 28 35 29  8 26  5 12
# equivalente a
valores <- base[]
valores
##  [1] 13 24 16 15 17 12 29  2 11 21 14 20 31 35 31 31 30 27 29 26  3 17 35 29  9
## [26] 29  5 26  1 28 35 29  8 26  5 12

A continuación, se eliminan estos objetos

rm(nc, nr, valores, valores)
## Warning in rm(nc, nr, valores, valores): objeto 'valores' no encontrado

También se dispone de funciones para obtener las coordenadas coordenadas reales a partir del número de celda

xFromCell(base,1)                   # Coordenada x de celda 1
## [1] -3.916667
yFromCell(base,1)                   # Coordenada y de celda 1
## [1] 42.91667
xyFromCell(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.91667

Y recíprocamente, a partir de coordenadas x e yse obtiene el número de celda.

cellFromXY(base, c(-3.75, 42.91667))       # Número de celda
## [1] 2

Igualmente, se puede solicitar el pixel que corresponde a determinada fila y columna

colFromCell(base, 28)                          # A qué columna corresponde el pixel 28
## [1] 4
rowFromCell(base, c(1,5,15))                   # A qué fila corresponen los píxeles 1, 5 y 15
## [1] 1 1 3
rowColFromCell(base, c(1,5,15))                # Cúales son las filas y columnas correspondientes a los píxeles
##      row col
## [1,]   1   1
## [2,]   1   5
## [3,]   3   3
cellFromRowCol(base, 2, 2)                     # Cuál es el pixel correspondiente a la fila 2 y a la columna 2 
## [1] 8

2.2 Objetos ráster multibanda

Los tipos de objetos R ráster que maneja el paquete raster son un RasterLayer, un RasterStack y un RasterBrick. Como ha sido señalado anteriormente,

  • Para crear un objeto ráster monobanda (RasterLayer) se emplea la función raster().

  • Para crear un objeto ráster multibanda podemos usar dos opciones, la función raster::stack() y la función raster::brick().

Los objetos RasterStack y RasterBrick sólo pueden crearse si todas las bandas tienen la misma resolución espacial y extensión. Tanto un RasterStackcomo un RasterBrick son archivos multibanda, pero almacenan cada banda de manera diferente.

  • Las bandas en un RasterStackson almacenadas como vínculos (links) a los ficheros ráster individuales almacenados en el ordenador, es decir, en realidad no conforman físicamente un único archivo.

  • Un RasterBrick almacena realmente todos los ficheros en un único objeto. Dado que el RasterBrick se aloja en la memoria del ordenador, conlleva más tiempo crearlo. Sin embargo, a la hora de realizar operaciones, los cálculos son más rápidos que con un RasterStack.

Raster stack vs raster brick
Raster stack vs raster brick

Se puede crear un RasterBrick a partir de objetos tipo RasterLayer, RasterStack o desde un archivo (multicapa). A continuación se creará uno 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  108

La función brick convierte el objeto array en un RasterBrick que contendrá la misma información, pero formando 3 capas.

multi <- brick(arr)
multi
## class      : RasterBrick 
## dimensions : 6, 6, 36, 3  (nrow, ncol, ncell, nlayers)
## resolution : 0.1666667, 0.1666667  (x, y)
## extent     : 0, 1, 0, 1  (xmin, xmax, ymin, ymax)
## crs        : NA 
## source     : memory
## names      : layer.1, layer.2, layer.3 
## min values :       1,      37,      73 
## max values :      36,      72,     108
plot(multi)

2.2.1 Indexado de archivos multibanda (Raster stack y Rasterbricks).

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      : RasterLayer 
## dimensions : 6, 6, 36  (nrow, ncol, ncell)
## resolution : 0.1666667, 0.1666667  (x, y)
## extent     : 0, 1, 0, 1  (xmin, xmax, ymin, ymax)
## crs        : NA 
## source     : memory
## names      : layer.2 
## values     : 37, 72  (min, max)

La 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      : RasterBrick 
## dimensions : 6, 6, 36, 2  (nrow, ncol, ncell, nlayers)
## resolution : 0.1666667, 0.1666667  (x, y)
## extent     : 0, 1, 0, 1  (xmin, xmax, ymin, ymax)
## crs        : NA 
## source     : memory
## names      : layer.1, layer.2 
## min values :       1,      37 
## max values :      36,      72
  • O con los nombres de cada capa.
multi[[c("layer.2", "layer.3")]]
## class      : RasterBrick 
## dimensions : 6, 6, 36, 2  (nrow, ncol, ncell, nlayers)
## resolution : 0.1666667, 0.1666667  (x, y)
## extent     : 0, 1, 0, 1  (xmin, xmax, ymin, ymax)
## crs        : NA 
## source     : memory
## names      : layer.2, layer.3 
## min values :      37,      73 
## max values :      72,     108

Si las capas poseen un nombre, también se pueden seleccionar con nombre tras el signo $

multi$layer.1
## class      : RasterLayer 
## dimensions : 6, 6, 36  (nrow, ncol, ncell)
## resolution : 0.1666667, 0.1666667  (x, y)
## extent     : 0, 1, 0, 1  (xmin, xmax, ymin, ymax)
## crs        : NA 
## source     : memory
## names      : layer.1 
## values     : 1, 36  (min, max)

Para extraer capas del RasterBrick original, puede utilizar el signo –c o bien la función raster::dropLayer(), como muestran los siguientes ejemplos.

multi.extraido <- multi[[-c(1,2)]]
multi.extraido
## class      : RasterLayer 
## dimensions : 6, 6, 36  (nrow, ncol, ncell)
## resolution : 0.1666667, 0.1666667  (x, y)
## extent     : 0, 1, 0, 1  (xmin, xmax, ymin, ymax)
## crs        : NA 
## source     : memory
## names      : layer.3 
## values     : 73, 108  (min, max)

Si se desea eliminar capas de un objeto, se puede realizar también con la función raster::dropLayer.

multi.reducido <- dropLayer(multi, c(1,3))
multi.reducido
## class      : RasterLayer 
## dimensions : 6, 6, 36  (nrow, ncol, ncell)
## resolution : 0.1666667, 0.1666667  (x, y)
## extent     : 0, 1, 0, 1  (xmin, xmax, ymin, ymax)
## crs        : NA 
## source     : memory
## names      : layer.2 
## values     : 37, 72  (min, max)

Y también es equivalente a

multi.reducido <- dropLayer(multi, c("layer.1", "layer.3"))
multi.reducido
## class      : RasterLayer 
## dimensions : 6, 6, 36  (nrow, ncol, ncell)
## resolution : 0.1666667, 0.1666667  (x, y)
## extent     : 0, 1, 0, 1  (xmin, xmax, ymin, ymax)
## crs        : NA 
## source     : memory
## names      : layer.2 
## values     : 37, 72  (min, max)

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 orignal con raster::addLayer(). Alternativamente, también se puede especificar un nuevo nombre de capa entre corchetes:

nuevaCapa <- multi[[1]]/2
nuevaCapa
## class      : RasterLayer 
## dimensions : 6, 6, 36  (nrow, ncol, ncell)
## resolution : 0.1666667, 0.1666667  (x, y)
## extent     : 0, 1, 0, 1  (xmin, xmax, ymin, ymax)
## crs        : NA 
## source     : memory
## names      : layer.1 
## values     : 0.5, 18  (min, max)
names(nuevaCapa) <- "nueva"                           # Cambio de nombre para que no coincida con la capa original

A continuación, se añadirá el objeto ráster multi la capa nuevaCapa.

multi.add <- addLayer(multi,          # Objeto raster al que se añadirá una nueva banda
                      nuevaCapa)             # Nueva banda que se añader
multi.add
## class      : RasterStack 
## dimensions : 6, 6, 36, 4  (nrow, ncol, ncell, nlayers)
## resolution : 0.1666667, 0.1666667  (x, y)
## extent     : 0, 1, 0, 1  (xmin, xmax, ymin, ymax)
## crs        : NA 
## names      : layer.1, layer.2, layer.3, nueva 
## min values :     1.0,    37.0,    73.0,   0.5 
## max values :      36,      72,     108,    18

La sintaxis anterior es equivalent a

multi.add[["miNuevaBanda"]] <- nuevaCapa
multi.add
## class      : RasterStack 
## dimensions : 6, 6, 36, 5  (nrow, ncol, ncell, nlayers)
## resolution : 0.1666667, 0.1666667  (x, y)
## extent     : 0, 1, 0, 1  (xmin, xmax, ymin, ymax)
## crs        : NA 
## names      : layer.1, layer.2, layer.3, nueva, miNuevaBanda 
## min values :     1.0,    37.0,    73.0,   0.5,          0.5 
## max values :      36,      72,     108,    18,           18

Al consultar los píxeles de un RasterStacko de un RasterBrick con corchetes simples, devolverá una matriz con una columna por capa y píxeles en filas:

multi[1:2, 3]
##      layer.1 layer.2 layer.3
## [1,]      13      49      85
## [2,]      14      50      86

Para 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]
## [1] 50 51 52

Estos comandos no solo son útiles para consultar valores ráster; también se pueden usar al revés para establecer o cambiar valores.

layerCopy <- multi[[3]]
plot(layerCopy)
text(layerCopy, 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.

layerCopy[2:3, 2:3]                  
## [1] 80 86 81 87

Estos valores van a ser convertidos en NA, es decir, NoData.

layerCopy[2:3, 2:3] <- NA
plot(layerCopy)
text(layerCopy, digits=1)

Para cambiar a 0 todos los valores en layerCopy podríamos hacer lo siguiente:

layerCopy[] <- 0
plot(layerCopy)

rm(list = ls())

ACTIVIDAD DE EVALUACIÓN CONTINUA

Con el fin de que el alumno ponga en práctica los conocimientos adquiridos, cada alumno deberá elaborar un script que contenga los procedimientos necesarios para responder a las cuestiones planteadas en Introduccion al manejo de ficheros raster con R. Una vez concluido, el script será enviado posteriormente al profesor a través de correo electrónico.

3 DESCARGA DE LA IMAGEN RÁSTER

La imagen que será analizada en este capítulo se encuentra disponible en:

Para comenzar a trabajar con ellos, es necesario crear una carpeta donde se alojarán los ficheros descargados. Esto puede hacerse directamente desde R. Como sugerencia, en primer lugar, se creará una carpeta en el disco D:/ denominada teledeteccion, y se establecerá esta carpeta como carpeta de trabajo por defecto.

#if(!dir.exists("D:/G174_OCW/")) dir.create("D:/G174_OCW/")
#setwd("D:/G174_OCW")

A continuación, el siguiente fragmento de código buscará si existe una carpeta denominada /Lab_2_manejo_ficheros_raster_R/datos en el directorio /G174_OCW/y si no la encuentra, la creará.

#if(!dir.exists("./Lab_2_manejo_ficheros_raster_R/datos/")) dir.create("./Lab_2_manejo_ficheros_raster_R/datos/")

Para cargarlos en R es necesario establecer dicha carpeta como carpeta de trabajo.

## [1] "D:/G174_OCW/Lab_2_manejo_ficheros_raster_R/datos"

Una vez creado esta carpeta, el siguiente fragmento descomprimirá el fichero *.tar con las imágenes de landsat en dicha carpeta

#untar("./datos/LC08_L1TP_201032_20210723_20210729_02_T1.tar", 
#      exdir = "./datos/")

Se puede comprobar esta relación de listando todos los ficheros del directorio de trabajo.

list.files("./datos/")
##  [1] "LC08_L1TP_201032_20210723_20210729_01_T1.rar"             
##  [2] "LC08_L1TP_201032_20210723_20210729_02_T1.tar"             
##  [3] "LC08_L1TP_201032_20210723_20210729_02_T1_ANG.txt"         
##  [4] "LC08_L1TP_201032_20210723_20210729_02_T1_B1.TIF"          
##  [5] "LC08_L1TP_201032_20210723_20210729_02_T1_B10.TIF"         
##  [6] "LC08_L1TP_201032_20210723_20210729_02_T1_B11.TIF"         
##  [7] "LC08_L1TP_201032_20210723_20210729_02_T1_B2.TIF"          
##  [8] "LC08_L1TP_201032_20210723_20210729_02_T1_B3.TIF"          
##  [9] "LC08_L1TP_201032_20210723_20210729_02_T1_B4.TIF"          
## [10] "LC08_L1TP_201032_20210723_20210729_02_T1_B5.TIF"          
## [11] "LC08_L1TP_201032_20210723_20210729_02_T1_B6.TIF"          
## [12] "LC08_L1TP_201032_20210723_20210729_02_T1_B7.TIF"          
## [13] "LC08_L1TP_201032_20210723_20210729_02_T1_B8.TIF"          
## [14] "LC08_L1TP_201032_20210723_20210729_02_T1_B9.TIF"          
## [15] "LC08_L1TP_201032_20210723_20210729_02_T1_MTL.json"        
## [16] "LC08_L1TP_201032_20210723_20210729_02_T1_MTL.txt"         
## [17] "LC08_L1TP_201032_20210723_20210729_02_T1_MTL.xml"         
## [18] "LC08_L1TP_201032_20210723_20210729_02_T1_QA_PIXEL.TIF"    
## [19] "LC08_L1TP_201032_20210723_20210729_02_T1_QA_RADSAT.TIF"   
## [20] "LC08_L1TP_201032_20210723_20210729_02_T1_SAA.TIF"         
## [21] "LC08_L1TP_201032_20210723_20210729_02_T1_stac.json"       
## [22] "LC08_L1TP_201032_20210723_20210729_02_T1_SZA.TIF"         
## [23] "LC08_L1TP_201032_20210723_20210729_02_T1_thumb_large.jpeg"
## [24] "LC08_L1TP_201032_20210723_20210729_02_T1_thumb_small.jpeg"
## [25] "LC08_L1TP_201032_20210723_20210729_02_T1_VAA.TIF"         
## [26] "LC08_L1TP_201032_20210723_20210729_02_T1_VZA.TIF"

El fichero descomprimido contiene tres tipos de archivos:

• txt: Metadatos

• xml: Metadatos en formato xml

• TIF: Archivos de imágenes

Podemos revisar el contenido total:

dir(path = "D:/G174_OCW/Lab_2_manejo_ficheros_raster_R/datos")
##  [1] "LC08_L1TP_201032_20210723_20210729_01_T1.rar"             
##  [2] "LC08_L1TP_201032_20210723_20210729_02_T1.tar"             
##  [3] "LC08_L1TP_201032_20210723_20210729_02_T1_ANG.txt"         
##  [4] "LC08_L1TP_201032_20210723_20210729_02_T1_B1.TIF"          
##  [5] "LC08_L1TP_201032_20210723_20210729_02_T1_B10.TIF"         
##  [6] "LC08_L1TP_201032_20210723_20210729_02_T1_B11.TIF"         
##  [7] "LC08_L1TP_201032_20210723_20210729_02_T1_B2.TIF"          
##  [8] "LC08_L1TP_201032_20210723_20210729_02_T1_B3.TIF"          
##  [9] "LC08_L1TP_201032_20210723_20210729_02_T1_B4.TIF"          
## [10] "LC08_L1TP_201032_20210723_20210729_02_T1_B5.TIF"          
## [11] "LC08_L1TP_201032_20210723_20210729_02_T1_B6.TIF"          
## [12] "LC08_L1TP_201032_20210723_20210729_02_T1_B7.TIF"          
## [13] "LC08_L1TP_201032_20210723_20210729_02_T1_B8.TIF"          
## [14] "LC08_L1TP_201032_20210723_20210729_02_T1_B9.TIF"          
## [15] "LC08_L1TP_201032_20210723_20210729_02_T1_MTL.json"        
## [16] "LC08_L1TP_201032_20210723_20210729_02_T1_MTL.txt"         
## [17] "LC08_L1TP_201032_20210723_20210729_02_T1_MTL.xml"         
## [18] "LC08_L1TP_201032_20210723_20210729_02_T1_QA_PIXEL.TIF"    
## [19] "LC08_L1TP_201032_20210723_20210729_02_T1_QA_RADSAT.TIF"   
## [20] "LC08_L1TP_201032_20210723_20210729_02_T1_SAA.TIF"         
## [21] "LC08_L1TP_201032_20210723_20210729_02_T1_stac.json"       
## [22] "LC08_L1TP_201032_20210723_20210729_02_T1_SZA.TIF"         
## [23] "LC08_L1TP_201032_20210723_20210729_02_T1_thumb_large.jpeg"
## [24] "LC08_L1TP_201032_20210723_20210729_02_T1_thumb_small.jpeg"
## [25] "LC08_L1TP_201032_20210723_20210729_02_T1_VAA.TIF"         
## [26] "LC08_L1TP_201032_20210723_20210729_02_T1_VZA.TIF"

También podemos crear listas por tipo de archivo:

tifs <- dir(path = "./datos/", pattern ="TIF")
tifs
##  [1] "LC08_L1TP_201032_20210723_20210729_02_T1_B1.TIF"       
##  [2] "LC08_L1TP_201032_20210723_20210729_02_T1_B10.TIF"      
##  [3] "LC08_L1TP_201032_20210723_20210729_02_T1_B11.TIF"      
##  [4] "LC08_L1TP_201032_20210723_20210729_02_T1_B2.TIF"       
##  [5] "LC08_L1TP_201032_20210723_20210729_02_T1_B3.TIF"       
##  [6] "LC08_L1TP_201032_20210723_20210729_02_T1_B4.TIF"       
##  [7] "LC08_L1TP_201032_20210723_20210729_02_T1_B5.TIF"       
##  [8] "LC08_L1TP_201032_20210723_20210729_02_T1_B6.TIF"       
##  [9] "LC08_L1TP_201032_20210723_20210729_02_T1_B7.TIF"       
## [10] "LC08_L1TP_201032_20210723_20210729_02_T1_B8.TIF"       
## [11] "LC08_L1TP_201032_20210723_20210729_02_T1_B9.TIF"       
## [12] "LC08_L1TP_201032_20210723_20210729_02_T1_QA_PIXEL.TIF" 
## [13] "LC08_L1TP_201032_20210723_20210729_02_T1_QA_RADSAT.TIF"
## [14] "LC08_L1TP_201032_20210723_20210729_02_T1_SAA.TIF"      
## [15] "LC08_L1TP_201032_20210723_20210729_02_T1_SZA.TIF"      
## [16] "LC08_L1TP_201032_20210723_20210729_02_T1_VAA.TIF"      
## [17] "LC08_L1TP_201032_20210723_20210729_02_T1_VZA.TIF"
txts <- dir(path = "./datos/", pattern ="txt")
txts
## [1] "LC08_L1TP_201032_20210723_20210729_02_T1_ANG.txt"
## [2] "LC08_L1TP_201032_20210723_20210729_02_T1_MTL.txt"
xmls <- dir(path = "./datos/", pattern ="xml")
xmls
## [1] "LC08_L1TP_201032_20210723_20210729_02_T1_MTL.xml"

3.1 Los metadatos

El volumen de información que contiene es elevado. Por un lado contiene:

  • Nombre del satélite y del sensor.

  • Identificador (ID) de la escena.

  • Fecha de adquisición y procesado de la escena.

  • Número de pasada y de fila.

  • Proyección.

Además, también incluye la información necesaria para llevar a cabo el proceso de corrección, como es:

  • Ángulo del azimuth del sol.

  • Elevación del sol.

  • Distancia de la Tierra al Sol.

  • Bandas disponibles.

  • Resolución espacial y radiométrica de cada banda.

Es importante tener acceso a los metadatos de la imagen para extraer información de la escena. R cuenta con una familia de funciones para la carga de archivos ASCII. La siguiente configuración nos permite cargar los datos existente en el fichero *_MTL.txt`

parametros <- read.csv(file = "./datos/LC08_L1TP_201032_20210723_20210729_02_T1_MTL.txt",
                       header = FALSE,
                       sep = "=", 
                       stringsAsFactors = FALSE, 
                       strip.white = TRUE)

names(parametros) <- c('id','valor')

Para visualizar el contenido podemos utilizar el mismo entorno abriendo un nuevo panel y desplegar los datos en forma de planilla de datos

View(parametros)

Del conjunto de datos disponibles, podemos revisar algunos interesantes para el análisis posterior:

• Porcentaje de cobertura Nubosa de la imagen:

cloud <- parametros[which(parametros$id == "CLOUD_COVER"), 2]
cloud
## [1] "2.97"

• Porcentaje de cobertura nubosa sobre superficie terrestre:

cloud_land <- parametros[which(parametros$id == "CLOUD_COVER_LAND"), 2]
cloud_land
## [1] "2.97"

• Fecha de toma de la imagen:

fecha <- parametros[which(parametros$id == "DATE_ACQUIRED"), 2]
fecha
## [1] "2021-07-23"

• Nombre oficial de la imagen; se usa el índice [1] ya que existe varias veces en el archivo:

id <- parametros[which(parametros$id == "LANDSAT_PRODUCT_ID"), 2][1]
id
## [1] "LC08_L1TP_201032_20210723_20210729_02_T1"

3.2 Procesando las Bandas en R

3.2.1 QA_pixel

Uno de los ficheros más importantes del conjunto es la imagen de control de calidad de los píxeles. Esta banda, en formato *.tif, presenta la etiqueta QA_PIXEL.

qa_px <- raster("./datos/LC08_L1TP_201032_20210723_20210729_02_T1_QA_PIXEL.tif")

En relación al contenido de la banda podemos revisar sus valores y revisar su significado

unique(qa_px)
##  [1]     1 21762 21824 21890 21952 22018 22080 22146 22280 23826 23888 24082
## [13] 24144 54534 54596 54790 54852 55052 56598 56660 56854 56916
Raster stack vs raster brick
Raster stack vs raster brick

Esos datos los podríamos utilizar como máscara para futuros procesos (Figura 8.31):

mascara <- qa_px
mascara[mascara != 21824] <- NA               # Si el valor no es igual a  21824 se transforma en valor NA
plot(mascara)

rm(imagen, mascara, qa_px)
## Warning in rm(imagen, mascara, qa_px): objeto 'imagen' no encontrado

3.3 Manejo de imágenes de satélite en formato ráster.

A partir de este momento, nos centraremos en el manejo de las imágenes de satélite en formato raster. Se va a trabajar fundamentalmente con los ficheros que corresponden a las bandas de Landsat 8.

Bandas de Landsat 8
Bandas de Landsat 8

Para importar cada banda de manera individual (en este caso, la banda NIR) la sintaxis es la siguiente:

imagen <- raster("./datos/LC08_L1TP_201032_20210723_20210729_02_T1_B5.TIF")

Al citar el nombre del objeto, R devuelve un resumen con información sobre dicho objeto.

imagen
## class      : RasterLayer 
## dimensions : 7911, 7791, 61634601  (nrow, ncol, ncell)
## resolution : 30, 30  (x, y)
## extent     : 316485, 550215, 4345785, 4583115  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=30 +datum=WGS84 +units=m +no_defs 
## source     : LC08_L1TP_201032_20210723_20210729_02_T1_B5.TIF 
## names      : LC08_L1TP_201032_20210723_20210729_02_T1_B5

La información proporcionada (visto en un apartado anterior) es la siguiente:

  • class: el tipo de datos (en lenguaje R).

  • dimensions: las dimensiones, según el número filas (nrow), columnas (ncol) y de celdas (ncell).

  • spatial resolution: la resolución espacial, en este caso en metros (30).

  • extent: los límites (según el sistema de proyección; dado que es UTM, serían metros) de la imagen.

  • source: ubicación del fichero.

  • names: nombre(s) de la banda importada. Este nombre puede ser cambiado por el investigador.

  • values: valores mínimo y máximo de los píxeles (en este caso, números digitales).

La importanción individual de diferentes bandas de la escena se puede realizar varias maneras. Una de ellas es utilizar la función raster() de manera consecutiva.

temp1 <- raster("./datos/LC08_L1TP_201032_20210723_20210729_02_T1_B1.TIF") 
temp2 <- raster("./datos/LC08_L1TP_201032_20210723_20210729_02_T1_B2.TIF") 
temp3 <- raster("./datos/LC08_L1TP_201032_20210723_20210729_02_T1_B3.TIF")

No obstante, lo habitual es convertir esos ficheros en un RasterStack, mediante la función raster::stack.

imagen.stack <- stack(temp1, temp2, temp3)
imagen.stack
## class      : RasterStack 
## dimensions : 7911, 7791, 61634601, 3  (nrow, ncol, ncell, nlayers)
## resolution : 30, 30  (x, y)
## extent     : 316485, 550215, 4345785, 4583115  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=30 +datum=WGS84 +units=m +no_defs 
## names      : LC08_L1TP_201032_20210723_20210729_02_T1_B1, LC08_L1TP_201032_20210723_20210729_02_T1_B2, LC08_L1TP_201032_20210723_20210729_02_T1_B3

Si se desea convertir todas las bandas en una única imagen multibanda, debe utilizarse la función raster::brick, que puede ser aplicada tanto a imágenes ya cargadas como a las almacenadas en el disco duro.

imagen.brick <- brick(imagen.stack)
imagen.brick
## class      : RasterBrick 
## dimensions : 7911, 7791, 61634601, 3  (nrow, ncol, ncell, nlayers)
## resolution : 30, 30  (x, y)
## extent     : 316485, 550215, 4345785, 4583115  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=30 +datum=WGS84 +units=m +no_defs 
## source     : r_tmp_2024-03-20_170249.58233_23188_81499.grd 
## names      : LC08_L1TP_201032_20210723_20210729_02_T1_B1, LC08_L1TP_201032_20210723_20210729_02_T1_B2, LC08_L1TP_201032_20210723_20210729_02_T1_B3 
## min values :                                        8486,                                        7594,                                        7003 
## max values :                                       65535,                                       57121,                                       65535

Una vez ilustrado este ejemplo, y dado que no necesitaremos más los objetos anteriores, se eliminan del Global Environment.

rm(temp1, temp2, temp3, imagen.stack, imagen.brick)

El procedimiento de importanción anterior es farrogoso y propenso a errores, dado que es fácil equivocarse por la complejidad de los nombres de los ficheros originales. Una opción alternativa es crear una lista con las bandas, aprovechando que comparten la misma denominación (salvo el número de banda) y que estan alojadas en la misma carpeta. Una vez creada la lista, será cargados con la función stack().

escena.completa <- list.files("./datos/",                                 # Directorio de los datos
                              pattern = ".*B[1234567]\\.tif$",            # El signo $ lista todos los ficheros que acaban igual
                              ignore.case=TRUE,                           # Sensible a mayúsculas y minúsculas
                              full.names = TRUE)
stack.escena.completa <- stack(escena.completa)
stack.escena.completa
## class      : RasterStack 
## dimensions : 7911, 7791, 61634601, 7  (nrow, ncol, ncell, nlayers)
## resolution : 30, 30  (x, y)
## extent     : 316485, 550215, 4345785, 4583115  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=30 +datum=WGS84 +units=m +no_defs 
## names      : LC08_L1TP//9_02_T1_B1, LC08_L1TP//9_02_T1_B2, LC08_L1TP//9_02_T1_B3, LC08_L1TP//9_02_T1_B4, LC08_L1TP//9_02_T1_B5, LC08_L1TP//9_02_T1_B6, LC08_L1TP//9_02_T1_B7

Otro procedimiento alternativo para crear esa lista es el siguiente

escena.completa2 <- paste0("./datos/LC08_L1TP_201032_20210723_20210729_02_T1_B", 
                  1:7,
                  pattern = ".TIF")
stack.escena.completa2 <- stack(escena.completa2)
stack.escena.completa2
## class      : RasterStack 
## dimensions : 7911, 7791, 61634601, 7  (nrow, ncol, ncell, nlayers)
## resolution : 30, 30  (x, y)
## extent     : 316485, 550215, 4345785, 4583115  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=30 +datum=WGS84 +units=m +no_defs 
## names      : LC08_L1TP//9_02_T1_B1, LC08_L1TP//9_02_T1_B2, LC08_L1TP//9_02_T1_B3, LC08_L1TP//9_02_T1_B4, LC08_L1TP//9_02_T1_B5, LC08_L1TP//9_02_T1_B6, LC08_L1TP//9_02_T1_B7

En ambos casos, sólo se han cargado las primeras 7 bandas, dado que las bandas 8, 9, 10 y 11 tienen una resolución espacial diferente al resto.

Dado que los nombres de banda son bastante largos y redundantes, una alternativa es acortarlos, pero sin que pierdan significado. El uso de nombres simples facilita el acceso a las capas individuales. Existen varias posibilidades para realizar el cambio. Por un lado, se pueden cambiar utilizando la función names() de acuerdo con la siguiente sintaxis:

names(stack.escena.completa) <- gsub(pattern = "LC08_L1TP_201032_20210723_20210729_02_T1_",   # conjunto de caracteres que se desea sustituir
                                     replace ="",                                             # caracteres que sustituirán a los anteriores
                                     x <- names(stack.escena.completa))                       # vector o dataframe que sustituirán a los anteriores
stack.escena.completa
## class      : RasterStack 
## dimensions : 7911, 7791, 61634601, 7  (nrow, ncol, ncell, nlayers)
## resolution : 30, 30  (x, y)
## extent     : 316485, 550215, 4345785, 4583115  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=30 +datum=WGS84 +units=m +no_defs 
## names      : B1, B2, B3, B4, B5, B6, B7

Con una sintaxis más fácil, es posible también reducir el nombre de las bandas así.

names(stack.escena.completa) <- paste("b",
                    1:7,
                    sep="")
stack.escena.completa
## class      : RasterStack 
## dimensions : 7911, 7791, 61634601, 7  (nrow, ncol, ncell, nlayers)
## resolution : 30, 30  (x, y)
## extent     : 316485, 550215, 4345785, 4583115  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=30 +datum=WGS84 +units=m +no_defs 
## names      : b1, b2, b3, b4, b5, b6, b7

3.4 Exportación e importación de imágenes en formato ráster.

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. Pero este diseño 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.

Para guardar objetos raster en el disco duro es necesario utilizar la función raster::writeRaster(), con la que es posible elegir diferentes opciones para almacenar los datos (más información en ?writeFormats).

  • El paquete raster tiene su propio formato de archivo: raster::writeRaster() salva el raster como 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.

  • Si se intercambian frecuentemente los ficheros con otros software SIG es recomendable el uso del formato GeoTIFF. Sin embargo, téngase en cuenta que los nombres de las bandas no se guardan bajo este formato, por lo que, en caso de archivos multicapa tipo RasterBrick se deberá recordar el orden de las bandas. Por otro lado, leer y escribir archivos GeoTiff en R es más lento que en el formato ráster .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(stack.escena.completa, 
            filename = "./resultados/landsat_multibanda.tif",
            overwrite = TRUE)

Podemos volver a importar este fichero para comprobar que todo se ha realizado correctamente. Si múltiples bandas de la misma escena ya están integradas en una única imagen, se puede cargar con la función brick.

landsat_brick <- brick("./resultados/landsat_multibanda.tif")
landsat_brick
## class      : RasterBrick 
## dimensions : 7911, 7791, 61634601, 7  (nrow, ncol, ncell, nlayers)
## resolution : 30, 30  (x, y)
## extent     : 316485, 550215, 4345785, 4583115  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=30 +datum=WGS84 +units=m +no_defs 
## source     : landsat_multibanda.tif 
## names      : landsat_multibanda_1, landsat_multibanda_2, landsat_multibanda_3, landsat_multibanda_4, landsat_multibanda_5, landsat_multibanda_6, landsat_multibanda_7 
## min values :                 8486,                 7594,                 7003,                 5998,                 5207,                 4755,                 4954 
## max values :                65535,                57121,                65535,                65535,                65535,                65535,                65535

Obsérvese que en este caso los nombres de las capas ya no son los originales, ya que como hemos señalado, al salvar un fichero en formato geotiff se pierde el nombre original de la capa.

En el caso de preferir el formato .grd

writeRaster(stack.escena.completa, 
            filename = "./resultados/landsat_multibanda.grd",
            overwrite = TRUE)

En el caso de importar el fichero *.grd se mantiene el nombre original de las capas, lo cual es recomendable.

landsat_brick <- brick("./resultados/landsat_multibanda.grd")
landsat_brick
## class      : RasterBrick 
## dimensions : 7911, 7791, 61634601, 7  (nrow, ncol, ncell, nlayers)
## resolution : 30, 30  (x, y)
## extent     : 316485, 550215, 4345785, 4583115  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=30 +datum=WGS84 +units=m +no_defs 
## source     : landsat_multibanda.grd 
## names      :    b1,    b2,    b3,    b4,    b5,    b6,    b7 
## min values :  8486,  7594,  7003,  5998,  5207,  4755,  4954 
## max values : 65535, 57121, 65535, 65535, 65535, 65535, 65535

Atención Si se utiliza la función raster para importar un fichero multibanda, sólo se carga la primera banda, e ignora las demás.

Ejemplo del comentario anterior

borrame <- raster("./resultados/landsat_multibanda.grd")
borrame
## class      : RasterLayer 
## band       : 1  (of  7  bands)
## dimensions : 7911, 7791, 61634601  (nrow, ncol, ncell)
## resolution : 30, 30  (x, y)
## extent     : 316485, 550215, 4345785, 4583115  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=30 +datum=WGS84 +units=m +no_defs 
## source     : landsat_multibanda.grd 
## names      : b1 
## values     : 8486, 65535  (min, max)

Además de los objetos raster, RasterStack y RasterBrick, es probable que haya que trabajar en algún momento con objetos de las clases SpatialPixelsDataFrame y SpatialGridDataFrames, generadas por el paquete sp. Ambas clases están estrechamente relacionadas con los objetos ráster anteriores, pero requieren que todos los datos se carguen en la memoria. Esta circunstancia es su gran inconveniente: el mensaje “Error: cannot allocate vector of size …” es una señal de que está alcanzando los límites de memoria de la máquina. Estos objetos se pueden importar con la función rgdal::readGDAL() del paquete rgdal, pero también desde el paquete raster con las funciones de conversión raster::as(), raster::raster(), raster::stack() y raster::brick().

Otros formatos son similares al formato (.grd); de hecho, la mayoría sólo difiere en el archivo de encabezado. Se puede utilizar esta característica para generar, por ejemplo, un archivo de encabezado ENVI adicional (.hdr), mediante la función raster::hdr().

hdr(landsat_brick, 
    format = "ENVI")

Si se establece como opción ráster global addheader = "ENVI" la función writeRaster() produce un archivo .hdr adicional cada vez que se guarda un archivo .grd, lo que facilita el intercambio de ficheros entre R y otros SIG clásico como QGIS:

rasterOptions(addheader = "ENVI")

A continuación, se eliminan todos los objetos ráster innecesarios

rm(borrame, brick.escena.completa, brick.escena.completa2, stack.escena.completa, stack.escena.completa2, escena.completa, escena.completa2, x)
## Warning in rm(borrame, brick.escena.completa, brick.escena.completa2,
## stack.escena.completa, : objeto 'brick.escena.completa' no encontrado
## Warning in rm(borrame, brick.escena.completa, brick.escena.completa2,
## stack.escena.completa, : objeto 'brick.escena.completa2' no encontrado

3.5 Atributos e indexado de imágenes en formato ráster.

El paquete ráster incluye también un conjunto de herramientas para abordar el análisis, bien de forma individual, bien conjunta, de los píxeles y de las capas que conforman los objetos raster. A este respecto, un RasterLayer es equivalente a una matriz de datos (con filas y columnas), mientras que un RasterStack o RasterBrick equivalen a una sucesión de matrices unidas entre sí.

Capa única vs raster multibanda
Capa única vs raster multibanda
https://www.neonscience.org/resources/learning-hub/tutorials/dc-multiband-rasters-r

3.5.1 Atributos

En ocasiones, puede ser interesante conocer los atributos específicos de una imagen. Si no interesan todos en bloque, es posible acceder por separado.

ncell(landsat_brick)                   # Número de celdas.
## [1] 61634601
nrow(landsat_brick)                    # Número de filas.
## [1] 7911
ncol(landsat_brick)                    # Número de columnas.
## [1] 7791
nlayers(landsat_brick)                 # Número de bandas.
## [1] 7
dim(landsat_brick)                     # Dimensiones.
## [1] 7911 7791    7
xres(landsat_brick)                    # Resolución espacial (dimensión x).
## [1] 30
yres(landsat_brick)                    # Resolución espacial (dimensión y).
## [1] 30
res(landsat_brick)                     # Resolución espacial.
## [1] 30 30
crs(landsat_brick)                     # Sistema de Coordenadas de Referencia (CRS).
## Coordinate Reference System:
## Deprecated Proj.4 representation:
##  +proj=utm +zone=30 +datum=WGS84 +units=m +no_defs 
## WKT2 2019 representation:
## PROJCRS["WGS 84 / UTM zone 30N",
##     BASEGEOGCRS["WGS 84",
##         DATUM["World Geodetic System 1984",
##             ELLIPSOID["WGS 84",6378137,298.257223563,
##                 LENGTHUNIT["metre",1]]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433]],
##         ID["EPSG",4326]],
##     CONVERSION["UTM zone 30N",
##         METHOD["Transverse Mercator",
##             ID["EPSG",9807]],
##         PARAMETER["Latitude of natural origin",0,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8801]],
##         PARAMETER["Longitude of natural origin",-3,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8802]],
##         PARAMETER["Scale factor at natural origin",0.9996,
##             SCALEUNIT["unity",1],
##             ID["EPSG",8805]],
##         PARAMETER["False easting",500000,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8806]],
##         PARAMETER["False northing",0,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8807]]],
##     CS[Cartesian,2],
##         AXIS["(E)",east,
##             ORDER[1],
##             LENGTHUNIT["metre",1]],
##         AXIS["(N)",north,
##             ORDER[2],
##             LENGTHUNIT["metre",1]],
##     USAGE[
##         SCOPE["Engineering survey, topographic mapping."],
##         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)."],
##         BBOX[0,-6,84,0]],
##     ID["EPSG",32630]]
landsat_brick@crs                      # Similar al anterior
## Coordinate Reference System:
## Deprecated Proj.4 representation: NA
extent(landsat_brick)                  # Límites (coordenadas).
## class      : Extent 
## xmin       : 316485 
## xmax       : 550215 
## ymin       : 4345785 
## ymax       : 4583115
landsat_brick@extent                   # Similar al anterior
## class      : Extent 
## xmin       : 316485 
## xmax       : 550215 
## ymin       : 4345785 
## ymax       : 4583115
cellStats(landsat_brick, 
          stat = "min")                # Valores mínimos de cada banda.
##   b1   b2   b3   b4   b5   b6   b7 
## 8486 7594 7003 5998 5207 4755 4954
object.size(landsat_brick)             # Tamaño de la imagen en el directorio
## 15944 bytes

La denominación de cada una de las bandas también puede modificarse para que sea más fácil recordarlas. En las líneas de código a continuación, se sustituye el nombre de las bandas por una etiqueta con la descripción real de cada banda. Primero, se comprobarán los nombres actuales:

names(landsat_brick)
## [1] "b1" "b2" "b3" "b4" "b5" "b6" "b7"

Posteriormente, se crea un vector con los nuevos nombres de cada una de las bandas, que serán asignados a cada banda que compone el RasterBrick.

names(landsat_brick) <- c("coastal_aerosol", "blue", "green", "red", "NIR", "SWIR_1", "SWIR_2")

Conviene comprobar si el proceso ha tenido éxito.

names(landsat_brick)
## [1] "coastal_aerosol" "blue"            "green"           "red"            
## [5] "NIR"             "SWIR_1"          "SWIR_2"

3.6 Indexado.

La indexación de archivos multibanda es similar al de una lista. 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 banda 2, que corresponde a la longitud de onda blue.

landsat_brick[[2]]
## class      : RasterLayer 
## band       : 2  (of  7  bands)
## dimensions : 7911, 7791, 61634601  (nrow, ncol, ncell)
## resolution : 30, 30  (x, y)
## extent     : 316485, 550215, 4345785, 4583115  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=30 +datum=WGS84 +units=m +no_defs 
## source     : landsat_multibanda.grd 
## names      : blue 
## values     : 7594, 57121  (min, max)

La consulta de más de una capa también se realiza con corchetes dobles:

  • Con los números de la banda
landsat_brick[[1:5]]
## class      : RasterStack 
## dimensions : 7911, 7791, 61634601, 5  (nrow, ncol, ncell, nlayers)
## resolution : 30, 30  (x, y)
## extent     : 316485, 550215, 4345785, 4583115  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=30 +datum=WGS84 +units=m +no_defs 
## names      : coastal_aerosol,  blue, green,   red,   NIR 
## min values :            8486,  7594,  7003,  5998,  5207 
## max values :           65535, 57121, 65535, 65535, 65535
  • O con los nombres de cada banda.
landsat_brick[[c("blue", "NIR")]]
## class      : RasterStack 
## dimensions : 7911, 7791, 61634601, 2  (nrow, ncol, ncell, nlayers)
## resolution : 30, 30  (x, y)
## extent     : 316485, 550215, 4345785, 4583115  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=30 +datum=WGS84 +units=m +no_defs 
## names      :  blue,   NIR 
## min values :  7594,  5207 
## max values : 57121, 65535

Si las capas poseen un nombre, también se pueden seleccionar con nombre tras el signo $

landsat_brick$blue
## class      : RasterLayer 
## band       : 2  (of  7  bands)
## dimensions : 7911, 7791, 61634601  (nrow, ncol, ncell)
## resolution : 30, 30  (x, y)
## extent     : 316485, 550215, 4345785, 4583115  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=30 +datum=WGS84 +units=m +no_defs 
## source     : landsat_multibanda.grd 
## names      : blue 
## values     : 7594, 57121  (min, max)

Para extraer capas del RasterBrick original, puede utilizar el signo –c o bien la función raster::dropLayer(), como muestran los siguientes ejemplos.

landsat.reducido <- landsat_brick[[-c(1,5)]]
landsat.reducido
## class      : RasterStack 
## dimensions : 7911, 7791, 61634601, 5  (nrow, ncol, ncell, nlayers)
## resolution : 30, 30  (x, y)
## extent     : 316485, 550215, 4345785, 4583115  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=30 +datum=WGS84 +units=m +no_defs 
## names      :  blue, green,   red, SWIR_1, SWIR_2 
## min values :  7594,  7003,  5998,   4755,   4954 
## max values : 57121, 65535, 65535,  65535,  65535

Si se desea eliminar capas de un objeto, se puede realizar con la función raster::dropLayer.

landsat.reducido <- dropLayer(landsat.reducido, c(1,3))
landsat.reducido
## class      : RasterStack 
## dimensions : 7911, 7791, 61634601, 3  (nrow, ncol, ncell, nlayers)
## resolution : 30, 30  (x, y)
## extent     : 316485, 550215, 4345785, 4583115  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=30 +datum=WGS84 +units=m +no_defs 
## names      : green, SWIR_1, SWIR_2 
## min values :  7003,   4755,   4954 
## max values : 65535,  65535,  65535

Y también es equivalente a

landsat.reducido <- dropLayer(landsat_brick, c("green", "SWIR_1"))
landsat.reducido
## class      : RasterStack 
## dimensions : 7911, 7791, 61634601, 5  (nrow, ncol, ncell, nlayers)
## resolution : 30, 30  (x, y)
## extent     : 316485, 550215, 4345785, 4583115  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=30 +datum=WGS84 +units=m +no_defs 
## names      : coastal_aerosol,  blue,   red,   NIR, SWIR_2 
## min values :            8486,  7594,  5998,  5207,   4954 
## max values :           65535, 57121, 65535, 65535,  65535

También se pueden añadir más capas a un objeto raster. En el siguiente ejemplo, los valores de la banda 1 (aerosol costero) se dividen entre dos, y esta nueva capa se agrega al objeto oriignal con raster::addLayer(), o especificar un nuevo nombre de capa entre corchetes:

nuevaBanda <- landsat_brick[[1]]/2
nuevaBanda
## class      : RasterLayer 
## dimensions : 7911, 7791, 61634601  (nrow, ncol, ncell)
## resolution : 30, 30  (x, y)
## extent     : 316485, 550215, 4345785, 4583115  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=30 +datum=WGS84 +units=m +no_defs 
## source     : memory
## names      : coastal_aerosol 
## values     : 4243, 32767.5  (min, max)
names(nuevaBanda) <- "nueva"                                             # Cambio de nombre para que no coincida con la capa original
landsat_brick.add <- addLayer(landsat_brick,                         # Objeto raster al que se añadirá una nueva banda
                              nuevaBanda)                              # Nueva banda que se añader

La sintaxis anterior es equivalent a

landsat_brick.add[["miNuevaBanda"]] <- nuevaBanda
landsat_brick.add
## class      : RasterStack 
## dimensions : 7911, 7791, 61634601, 9  (nrow, ncol, ncell, nlayers)
## resolution : 30, 30  (x, y)
## extent     : 316485, 550215, 4345785, 4583115  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=30 +datum=WGS84 +units=m +no_defs 
## names      : coastal_aerosol,    blue,   green,     red,     NIR,  SWIR_1,  SWIR_2,   nueva, miNuevaBanda 
## min values :            8486,    7594,    7003,    5998,    5207,    4755,    4954,    4243,         4243 
## max values :         65535.0, 57121.0, 65535.0, 65535.0, 65535.0, 65535.0, 65535.0, 32767.5,      32767.5
rm(list = ls())

ACTIVIDADES DE EVALUACIÓN CONTINUA

Las actividades de este bloque son dos:

Cada alumno deberá elaborar un script que contenga los procedimientos necesarios para responder a las cuestiones planteadas en dichas actividades. Una vez concluidas, el script será enviado posteriormente al profesor a través de correo electrónico.