Operaciones matemáticas básicas

Operaciones matemáticas básicas

El paquete rasterizado soporta muchas operaciones matemáticas. Las operaciones matemáticas se realizan generalmente por píxel (célula de la cuadrícula). Primero haremos algunas operaciones aritméticas básicas para combinar las bandas. En el primer ejemplo escribimos una función matemática personalizada para calcular el Índice de Diferencia Normalizada de Vegetación (NDVI). Aprende más sobre los índices de vegetación aquí y el NDVI.

Utilizamos los mismos datos de Landsat que en el capítulo 2.

library(raster)

## Loading required package: sp

library(rgdal)

## rgdal: version: 1.4-8, (SVN revision 845)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 2.2.3, released 2017/11/20
##  Path to GDAL shared files: C:/Users/YEIMY MONTILLA/Documents/R/win-library/3.6/rgdal/gdal
##  GDAL binary built with GEOS: TRUE 
##  Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
##  Path to PROJ.4 shared files: C:/Users/YEIMY MONTILLA/Documents/R/win-library/3.6/rgdal/proj
##  Linking to sp version: 1.4-1

raslist <- paste0('data/rs/LC08_044034_20170614_B', 1:11, ".tif")
landsat <- stack(raslist)
landsatRGB <- landsat[[c(4,3,2)]]
landsatFCC <- landsat[[c(5,4,3)]]

Índices de vegetación

Definamos una función general para un índice basado en la proporción (vegetación). En la siguiente función, img es un objeto Raster* de varias capas y i y k son los índices de las capas (números de capa) utilizados para calcular el índice de vegetación.

vi <- function(img, k, i) {
  bk <- img[[k]]
  bi <- img[[i]]
  vi <- (bk - bi) / (bk + bi)
  return(vi)
}

Aplicacion la funcion NDVI y ploteo de resultados

# For Landsat NIR = 5, red = 4.
ndvi <- vi(landsat, 5, 4)
plot(ndvi, col = rev(terrain.colors(10)), main = "Landsat-NDVI")

Puedes ver la variación del verde en el gráfico.

Otra forma de operar las capas y generar el NDVI

vi2 <- function(x, y) {
    (x - y) / (x + y)
}
ndvi2 <- overlay(landsat[[5]], landsat[[4]], fun=vi2)
plot(ndvi2, col=rev(terrain.colors(10)), main="Landsat-NDVI")

Pregunta 1: Adapte el código que se muestra arriba para calcular los índices para identificar i) el agua y ii) las acumulaciones. Sugerencia: Utilice la gráfica del perfil espectral para encontrar las bandas que tienen una reflectancia máxima y mínima para estas dos clases.

Histograma Podemos explorar la distribución de los valores contenidos en nuestro raster usando la función hist() que produce un histograma. Los histogramas son a menudo útiles para identificar valores atípicos y malos datos en nuestros datos de raster.

# ver el histograma de datos
hist(ndvi,
     main = "Distribution of NDVI values",
     xlab = "NDVI",
     ylab= "Frequency",
     col = "wheat",
     xlim = c(-0.5, 1),
     breaks = 30,
     xaxt = 'n')
axis(side=1, at = seq(-0.5,1, 0.05), labels = seq(-0.5,1, 0.05))

Umbrales

Podemos aplicar reglas básicas para obtener una estimación de la extensión espacial de las diferentes características de la superficie de la Tierra. Nótese que los valores del NDVI están estandarizados y varían entre -1 y +1. Valores más altos indican más cobertura verde.

Las células con valores de NDVI superiores a 0,4 son definitivamente vegetación. La siguiente operación enmascara todas las células que tal vez no son vegetación.

La siguiente operación enmascara todas las celdas que tal vez no sean vegetación

veg <- reclassify(ndvi, cbind(-Inf, 0.4, NA))
plot(veg, main='Vegetation')

Cartografiemos el área que corresponde al pico entre 0,25 y 0,3 en el histograma del NDVI.

área que corresponde al pico entre 0.25 y 0.3 en el histograma NDVI.

land <- reclassify(ndvi, c(-Inf, 0.25, NA,  0.25, 0.3, 1,  0.3, Inf, NA))
plot(land, main = 'What is it?')

Estas pueden ser las áreas abiertas. Puedes trazar la tierra sobre la trama original de la FCC para saber más.

Composicion en felso Color

plotRGB(landsatRGB, r=1, g=2, b=3, axes=TRUE, stretch="lin", main="Landsat False Color Composite")
plot(land, add=TRUE, legend=FALSE)

También puedes crear clases para diferentes cantidades de presencia de vegetación.

vegc <- reclassify(ndvi, c(-Inf,0.25,1, 0.25,0.3,2, 0.3,0.4,3, 0.4,0.5,4, 0.5,Inf, 5))
plot(vegc,col = rev(terrain.colors(4)), main = 'NDVI based thresholding')

Pregunta 3: ¿Es posible encontrar agua utilizando el umbral del NDVI o cualquier otro índice?

Análisis de componentes principales

Análisis de los componentes principales Los datos multiespectrales a veces se transforman para ayudar a reducir la dimensionalidad y el ruido de los datos. La transformación de los componentes principales es un método genérico de reducción de datos que puede utilizarse para crear unas pocas bandas no correlacionadas a partir de un conjunto mayor de bandas correlacionadas.

Se puede calcular el mismo número de componentes principales que el número de bandas de entrada. El primer componente principal (PC) explica el mayor porcentaje de varianza y otros PC explican adicionalmente la varianza en orden decreciente.

set.seed(1)
sr <- sampleRandom(landsat, 10000)
plot(sr[,c(4,5)], main = "NIR-Red plot")

Esto se conoce como gráfico de vegetación y línea de suelo (igual que el gráfico de dispersión de la sección anterior).

pca <- prcomp(sr, scale = TRUE)
pca

## Standard deviations (1, .., p=11):
##  [1] 2.52687022 1.40533428 1.09901821 0.92507423 0.53731672 0.42657919
##  [7] 0.28002273 0.12466139 0.09031384 0.04761419 0.03609857
## 
## Rotation (n x k) = (11 x 11):
##                                PC1        PC2         PC3         PC4
## LC08_044034_20170614_B1  0.2973469 -0.3516438 -0.29454767 -0.06983456
## LC08_044034_20170614_B2  0.3387629 -0.3301194 -0.16407835 -0.03803678
## LC08_044034_20170614_B3  0.3621173 -0.2641152  0.07373240  0.04090884
## LC08_044034_20170614_B4  0.3684797 -0.1659582  0.10260552  0.03680852
## LC08_044034_20170614_B5  0.1546322  0.1816252  0.71354112  0.32017620
## LC08_044034_20170614_B6  0.3480230  0.2108982  0.23064060  0.16598851
## LC08_044034_20170614_B7  0.3496281  0.2384417 -0.11662258  0.07600209
## LC08_044034_20170614_B8  0.3490464 -0.2007305  0.08765521  0.02303421
## LC08_044034_20170614_B9  0.1314827  0.1047365  0.33741447 -0.92325315
## LC08_044034_20170614_B10 0.2497611  0.4912132 -0.29286315 -0.02950655
## LC08_044034_20170614_B11 0.2472765  0.4931489 -0.29515754 -0.03176549
##                                  PC5          PC6         PC7          PC8
## LC08_044034_20170614_B1   0.49474685  0.175510232 -0.23948553  0.215745065
## LC08_044034_20170614_B2   0.22121122  0.094184121  0.06447037  0.216537517
## LC08_044034_20170614_B3   0.08482031  0.009040232  0.30511210 -0.518233675
## LC08_044034_20170614_B4  -0.33835490 -0.066844529  0.60174786  0.012437959
## LC08_044034_20170614_B5   0.51960822 -0.059561476 -0.07348455 -0.083217504
## LC08_044034_20170614_B6  -0.29437062  0.317984598 -0.02106132  0.632178645
## LC08_044034_20170614_B7  -0.25404931  0.525411720 -0.40543545 -0.478543437
## LC08_044034_20170614_B8  -0.31407992 -0.673584139 -0.52642131  0.003527306
## LC08_044034_20170614_B9   0.03040161  0.059642905 -0.03152221 -0.002775800
## LC08_044034_20170614_B10  0.16317572 -0.243735973  0.14341520  0.041736319
## LC08_044034_20170614_B11  0.19294569 -0.241611777  0.11997475 -0.022446494
##                                   PC9         PC10          PC11
## LC08_044034_20170614_B1   0.122812108  0.535959306  0.1203473847
## LC08_044034_20170614_B2   0.091063964 -0.773112627 -0.1817872036
## LC08_044034_20170614_B3  -0.644305383  0.070860458  0.0540175730
## LC08_044034_20170614_B4   0.543822097  0.218324141  0.0135097158
## LC08_044034_20170614_B5   0.209682702 -0.040186292 -0.0004965182
## LC08_044034_20170614_B6  -0.397210135  0.089423617 -0.0045305608
## LC08_044034_20170614_B7   0.248871690 -0.074907393  0.0003958921
## LC08_044034_20170614_B8  -0.012642132  0.002411975  0.0007749022
## LC08_044034_20170614_B9   0.003176077 -0.000832654  0.0019986985
## LC08_044034_20170614_B10 -0.003574004 -0.158727672  0.6900281980
## LC08_044034_20170614_B11 -0.043475408  0.148102091 -0.6878990264

screeplot(pca)

pci <- predict(landsat, pca, index = 1:2)
plot(pci[[1]])

El primer componente principal resalta los límites entre las clases de uso de la tierra o los detalles espaciales, que es la información más común entre todas las longitudes de onda. Es difícil entender lo que el segundo componente principal está resaltando. Intentemos de nuevo con la composición :

pc2 <- reclassify(pci[[2]], c(-Inf,0,1,0,Inf,NA))
par(mfrow = c(1,2))
plotRGB(landsatFCC, r = 1, g = 2, b = 3, axes = TRUE, stretch = "lin", main = "Landsat False Color Composite")
plotRGB(landsatFCC, r = 1, g = 2, b = 3, axes = TRUE, stretch = "lin", main = "Landsat False Color Composite")
plot(pc2, legend = FALSE, add = TRUE)