A continuación se presentan los comandos básicos en R para explorar la imágenes satelitales:

-Propiedades de imagen -Información de imagen y estadísticas -Banda única y mapas compuestos -Subconjunto y renombrar bandas -Subconjunto espacial o recorte -Guardando resultados en el disco -Relación entre bandas -Extraer valores de píxeles

-Propiedades de imagen Inicialmente se debe crear un objeto por cada una de las bandas de la imagen.

library(raster)
## Loading required package: sp
##Un RasterLayerobjeto representa datos ráster de una sola capa (variable)
# Banda Azul
b2 <- raster('data/rs/LC08_044034_20170614_B2.tif')
# Banda Verde
b3 <- raster('data/rs/LC08_044034_20170614_B3.tif')
# Banda roja
b4 <- raster('data/rs/LC08_044034_20170614_B4.tif')
# Infrarojo cercano (NIR)
b5 <- raster('data/rs/LC08_044034_20170614_B5.tif')

Propiedades de la banda Azul

b2
## class      : RasterLayer 
## dimensions : 1245, 1497, 1863765  (nrow, ncol, ncell)
## resolution : 30, 30  (x, y)
## extent     : 594090, 639000, 4190190, 4227540  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
## source     : C:/UNAL_MAESTRIA_UNAL/PERCEPCION_REMOTA/PRACTICAS/data/rs/LC08_044034_20170614_B2.tif 
## names      : LC08_044034_20170614_B2 
## values     : 0.0748399, 0.7177562  (min, max)

Información de imagen y estadísticas A continuación se muestra cómo puede acceder a varias propiedades desde un objeto Ráster * (esto es lo mismo para cualquier conjunto de datos ráster).

# coordinate reference system (CRS) - Referencia del Sistema de Coordenadas
crs(b2)
## CRS arguments:
##  +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs +ellps=WGS84
## +towgs84=0,0,0
# Número de celdas (columnas * filas)
ncell(b2)
## [1] 1863765
## Número de columnas y filas
dim(b2)
## [1] 1245 1497    1
# Resolución Espacial
res(b2)
## [1] 30 30
# Número de Bandas
nlayers(b2)
## [1] 1
# Compara con otra banda para saber si tienen las mismas características
compareRaster(b2,b3)
## [1] TRUE

Puede crear un RasterStack (un objeto con varias capas) a partir de los objetos RasterLayer (banda única) existentes.

# Crear un RasterStackes
# RasterStackes una colección de RasterLayerobjetos con la misma extensión y resolución espacial.
# se puede formar a partir de archivos separados y / o de unas pocas capas ('bandas') de un solo archivo.
s <- stack(b5, b4, b3)
# Comprueba las propiedades del RasterStack
s
## class      : RasterStack 
## dimensions : 1245, 1497, 1863765, 3  (nrow, ncol, ncell, nlayers)
## resolution : 30, 30  (x, y)
## extent     : 594090, 639000, 4190190, 4227540  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
## names      : LC08_044034_20170614_B5, LC08_044034_20170614_B4, LC08_044034_20170614_B3 
## min values :            0.0008457669,            0.0208406653,            0.0425921641 
## max values :               1.0124315,               0.7861769,               0.6924697

También puede crear el RasterStack usando los nombres de archivo.

# crear una lista de capas rasterizadas para usar
filenames <- paste0('data/rs/LC08_044034_20170614_B', 1:11, ".tif")
# Muestra los elemnetos de la lista
filenames
##  [1] "data/rs/LC08_044034_20170614_B1.tif" 
##  [2] "data/rs/LC08_044034_20170614_B2.tif" 
##  [3] "data/rs/LC08_044034_20170614_B3.tif" 
##  [4] "data/rs/LC08_044034_20170614_B4.tif" 
##  [5] "data/rs/LC08_044034_20170614_B5.tif" 
##  [6] "data/rs/LC08_044034_20170614_B6.tif" 
##  [7] "data/rs/LC08_044034_20170614_B7.tif" 
##  [8] "data/rs/LC08_044034_20170614_B8.tif" 
##  [9] "data/rs/LC08_044034_20170614_B9.tif" 
## [10] "data/rs/LC08_044034_20170614_B10.tif"
## [11] "data/rs/LC08_044034_20170614_B11.tif"
# Crear un RasterStackes con los elemento de la lista.
landsat <- stack(filenames)
# Comprueba las propiedades del RasterStack
landsat
## class      : RasterStack 
## dimensions : 1245, 1497, 1863765, 11  (nrow, ncol, ncell, nlayers)
## resolution : 30, 30  (x, y)
## extent     : 594090, 639000, 4190190, 4227540  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
## names      : LC08_044034_20170614_B1, LC08_044034_20170614_B2, LC08_044034_20170614_B3, LC08_044034_20170614_B4, LC08_044034_20170614_B5, LC08_044034_20170614_B6, LC08_044034_20170614_B7, LC08_044034_20170614_B8, LC08_044034_20170614_B9, LC08_044034_20170614_B10, LC08_044034_20170614_B11 
## min values :            9.641791e-02,            7.483990e-02,            4.259216e-02,            2.084067e-02,            8.457669e-04,           -7.872183e-03,           -5.052945e-03,            3.931751e-02,           -4.337332e-04,             2.897978e+02,             2.885000e+02 
## max values :              0.73462820,              0.71775615,              0.69246972,              0.78617686,              1.01243150,              1.04320455,              1.11793602,              0.82673049,              0.03547901,             322.43139648,             317.99530029

Banda única y mapas compuesto

Puede mostrar las capas individuales de un RasterStack de una imagen multiespectral.

# Matriz de 2*2
par(mfrow = c(2,2))
plot(b2, main = "Azul", col = gray(0:100 / 100))
plot(b3, main = "Verde", col = gray(0:100 / 100))
plot(b4, main = "Rojo", col = gray(0:100 / 100))
plot(b5, main = "NIR", col = gray(0:100 / 100))

Echa un vistazo a las leyendas de los mapas creados anteriormente. Pueden variar entre 0 y 1. Observe la diferencia en el sombreado y el rango de leyendas entre las diferentes bandas. Esto se debe a que las diferentes características de la superficie reflejan la radiación solar incidente de manera diferente. Cada capa representa la cantidad de radiación solar incidente que se refleja para un rango de longitud de onda particular. Por ejemplo, la vegetación refleja más energía en NIR que otras longitudes de onda y, por lo tanto, parece más brillante. Por el contrario, el agua absorbe la mayor parte de la energía en la longitud de onda NIR y parece oscura.

No obtenemos tanta información de estas parcelas en escala de grises; a menudo se combinan en un “compuesto” para crear tramas más interesantes. Puede obtener más información sobre los compuestos de color en la teledetección aquí y también en la sección a continuación.

Para hacer una imagen de “color verdadero (o natural)”, es decir, algo que parece una fotografía normal (vegetación en verde, azul agua, etc.), necesitamos bandas en las regiones roja, verde y azul. Para esta imagen Landsat, se pueden usar las bandas 4 (rojo), 3 (verde) y 2 (azul). El plotRGBmétodo se puede utilizar para combinarlos en un solo compuesto. También puede proporcionar argumentos adicionales plotRGBpara mejorar la visualización (por ejemplo, un estiramiento lineal de los valores, utilizando ).strecth = “lin”

# Crear un RasterStackes con los layers de las bandas Rojo, Verde y azul
landsatRGB <- stack(b4, b3, b2)
#Muestra RasterStackes -composición a color
plotRGB(landsatRGB, axes=TRUE, stretch="lin", main="Landsat True Color Composite")

La composición de color verdadero revela mucho más sobre el paisaje que las anteriores imágenes grises. Otro método popular de visualización de imágenes en la teledetección es la conocida imagen de “color falso” en la que se combinan bandas NIR, rojas y verdes. Esta representación es popular ya que facilita la visualización de la vegetación (en rojo).

# Matriz de 1*2
par(mfrow = c(1,2))
#Muestra RasterStackes -composición a color
plotRGB(landsatRGB, axes=TRUE, stretch="lin", main="Landsat True Color Composite")
# Crear un RasterStackes con los layers de las bandas NIR. Rojo y Verde
landsatFCC <- stack(b5, b4, b3)
#Muestra RasterStackes -composición en falso color
plotRGB(landsatFCC, axes=TRUE, stretch="lin", main="Landsat False Color Composite")

Nota: Revise siempre la documentación del paquete (help(plotRGB)) para otros argumentos que se pueden añadir (como la escala) para mejorar o modificar la imagen.

Pregunta 1: Utiliza la función plotRGB con RasterStack landsat para crear una composición de color verdadera y falsa (recuerda la posición de las bandas en la pila).

Subconjunto y renombrar las bandas

Puedes seleccionar capas específicas (bandas) usando la función de subconjunto, o a través de la indexación.

# Selecciona sólo las 3 primeras bandas
landsatsub1 <- subset(landsat, 1:3)
# igual 
landsatsub2 <- landsat[[1:3]]
# Número de bandas en los datos originales y nuevos
nlayers(landsat)
## [1] 11
nlayers(landsatsub1)
## [1] 3
nlayers(landsatsub2)
## [1] 3

No usaremos las últimas cuatro bandas en landsat. Puedes quitarlas usando

# Selecciona sólo las 7 primeras bandas
landsat <- subset(landsat, 1:7)

Para mayor claridad, es útil establecer los nombres de las bandas.

names(landsat)
## [1] "LC08_044034_20170614_B1" "LC08_044034_20170614_B2"
## [3] "LC08_044034_20170614_B3" "LC08_044034_20170614_B4"
## [5] "LC08_044034_20170614_B5" "LC08_044034_20170614_B6"
## [7] "LC08_044034_20170614_B7"
names(landsat) <- c('ultra-blue', 'blue', 'green', 'red', 'NIR', 'SWIR1', 'SWIR2')
names(landsat)
## [1] "ultra.blue" "blue"       "green"      "red"        "NIR"       
## [6] "SWIR1"      "SWIR2"

Subconjunto espacial o cultivo

El subconjunto espacial puede utilizarse para limitar el análisis a un subconjunto geográfico de la imagen. Los subconjuntos espaciales pueden crearse con la función de cultivo, utilizando un objeto de extensión, u otro objeto espacial del que se pueda extraer una extensión.

# Usando la extensión
extent(landsat)
## class      : Extent 
## xmin       : 594090 
## xmax       : 639000 
## ymin       : 4190190 
## ymax       : 4227540
e <- extent(624387, 635752, 4200047, 4210939)
# crop landsat by the extent
landsatcrop <- crop(landsat, e)

Pregunta 2: La selección interactiva de la imagen también es posible. Use “dibujar el texto” y “dibujar el polietileno” para seleccionar un área de interés.

Pregunta 3: Usar el “Landsatcrop” de RasterStack para crear un compuesto de color verdadero y falso.

Guardar los resultados en el disco

En esta etapa podemos querer guardar el raster en el disco usando la función writeRaster. Se admiten varios tipos de archivos. Usaremos el formato GeoTiff comúnmente usado. Mientras que el orden de las capas se conserva, los nombres de las capas se pierden lamentablemente en el formato GeoTiff.

x <- writeRaster(landsatcrop, filename="cropped-landsat.tif", overwrite=TRUE)

Alternativamente puedes usar el formato ‘raster-grd’.

writeRaster(landsatcrop, filename="cropped-landsat.grd", overwrite=TRUE)

Una ventaja de este formato es que guarda los nombres de las capas. La desventaja de este formato es que no muchos otros programas pueden leer los datos, en contraste con el formato GeoTiff.

Nota: Revise la documentación del paquete (help(writeRaster)) para ver si hay argumentos útiles adicionales que se puedan añadir.

Relación entre las bandas

Una matriz de dispersión puede ser útil para explorar las relaciones entre las capas rasterizadas. Esto puede hacerse con la función pairs() del paquete raster.

Gráfica de reflexión en la longitud de onda ultra azul contra la reflexión en la longitud de onda azul.

pairs(landsatcrop[[1:2]], main = "Ultra-blue versus Blue")

Trama de la reflexión en la longitud de onda roja contra la reflexión en la longitud de onda NIR.

pairs(landsatcrop[[4:5]], main = "Red versus NIR")

La primera trama revela altas correlaciones entre las regiones de longitud de onda del azul. Debido a la alta correlación, podemos usar una de las bandas azules sin perder mucha información.

Esta distribución de puntos en el segundo gráfico (entre el NIR y el rojo) es única debido a su forma triangular. La vegetación se refleja muy bien en el rango NIR que el rojo y crea la esquina superior cerca del eje NIR (y). El agua absorbe la energía de todas las bandas y ocupa el lugar cercano al origen. El ángulo más lejano se crea debido a las características de la superficie altamente reflectante como el suelo brillante o el hormigón.

Extraer los valores de los píxeles

A menudo queremos obtener los valores de las células raster para lugares o áreas geográficas específicas. La función de extracción se utiliza para obtener valores raster en las ubicaciones de otros datos espaciales. Se pueden utilizar puntos, líneas, polígonos o un objeto Extent (rectángulo). También puede utilizar números de celda para extraer valores. Cuando se utilizan puntos, la extracción devuelve los valores de un objeto Raster* para las celdas en las que caen un conjunto de puntos.

# Cargar los polígonos con información sobre la cobertura de la tierra
samp <- readRDS('data/rs/samples.rds')
# Generar muestras de 300 puntos de los polígonos
ptsamp <- spsample(samp, 300, type='regular')
# Agregar la clase de cobertura de tierra a los puntos
ptsamp$class <- over(ptsamp, samp)$class
# Extraer valores con puntos
df <- extract(landsat, ptsamp)
# Para ver algunos de los valores de reflectancia
head(df)
##      ultra.blue      blue     green       red       NIR     SWIR1     SWIR2
## [1,]  0.1334367 0.1163694 0.1022516 0.1009504 0.1749445 0.2161054 0.1812119
## [2,]  0.1361692 0.1180610 0.1025335 0.1017962 0.1769830 0.2105320 0.1744024
## [3,]  0.1290343 0.1210754 0.1267139 0.1703470 0.3089017 0.3323447 0.1927708
## [4,]  0.1326776 0.1288608 0.1400077 0.1872624 0.3399783 0.3316074 0.1949394
## [5,]  0.1470774 0.1529545 0.1812336 0.2523217 0.4037582 0.4015462 0.2377918
## [6,]  0.1382294 0.1213139 0.1076298 0.1093213 0.1753782 0.2318280 0.2001008

Perfiles espectrales

Un trazado del espectro (todas las bandas) para los píxeles que representan ciertas características de la superficie terrestre (por ejemplo, el agua) se conoce como perfil espectral. Esos perfiles demuestran las diferencias en las propiedades espectrales de diversos rasgos de la superficie terrestre y constituyen la base del análisis de las imágenes. Los valores espectrales pueden extraerse de cualquier conjunto de datos multiespectrales utilizando la función de extracción. En el ejemplo anterior, se extrajeron los valores de los datos del Landsat para las muestras. Estas muestras incluyen: tierra de cultivo, agua, barbecho, construida y abierta. Primero calculamos los valores medios de reflectancia para cada clase y cada banda.

ms <- aggregate(df, list(ptsamp$class), mean)
# en lugar de la primera columna, usamos nombres de fila
rownames(ms) <- ms[,1]
ms <- ms[,-1]
ms

Ahora trazamos los espectros medios de estas características.

# Crear un vector de color para las clases de cobertura terrestre para su uso en el trazado
mycolor <- c('darkred', 'yellow', 'burlywood', 'cyan', 'blue')
#Transformar ms de un data.frame a una matriz
ms <- as.matrix(ms)
# Primero crea una trama vacía
plot(0, ylim=c(0,0.6), xlim = c(1,7), type='n', xlab="Bands", ylab = "Reflectance")
# Agregar las diferentes clases
for (i in 1:nrow(ms)){
  lines(ms[i,], type = "l", lwd = 3, lty = 1, col = mycolor[i])
}
# Título
title(main="Spectral Profile from Landsat", font.main = 2)
# Leyenda
legend("topleft", rownames(ms),
       cex=0.8, col=mycolor, lty = 1, lwd =3, bty = "n")