practica 4

Including Plots

You can also embed plots, for example:

Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot.

library(raster)

## Loading required package: sp

raslist <- paste0('D:/documentos/universidad/perc_remota/practicas en R/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 razón (vegetación). En la función a continuación, imges un objeto Raster * de múltiples capas iy kson los índices de las capas (números de capa) utilizados para calcular el índice de vegetación.

FUNCION Y CODIGO PARA NDVI

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

# 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 en el verdor desde la trama.

Una forma alternativa de lograr esto es así

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")

Adapte el código que se muestra arriba para calcular los índices para identificar i) agua y ii) acumulados. Sugerencia: Use el diagrama del perfil espectral para encontrar las bandas que tienen reflectancia máxima y mínima para estas dos clases

aGUA

ndvi3 <- overlay(landsat[[1]], landsat[[7]], fun=vi2)
plot(ndvi3, col=rev(terrain.colors(10)), main="Landsat-NDVI AGUA")

ndvi4 <- overlay(landsat[[4]], landsat[[6]], fun=vi2)
plot(ndvi4, col=rev(terrain.colors(10)), main="Landsat-NDVI BARBECHO")

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

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))

Pregunta 2 : Haga histogramas de los valores que los índices de vegetación desarrollaron en la pregunta 1.

hist(ndvi3,
     main = "Distribution of NDVI values -AGUA",
     xlab = "NDVI AGUA",
     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))

BARBECHO

hist(ndvi4,
     main = "Distribution of NDVI values -BARBECHO",
     xlab = "NDVI BARBECHO",
     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))

Umbralización Podemos aplicar reglas básicas para obtener una estimación de la extensión espacial de diferentes características de la superficie terrestre. Tenga en cuenta que los valores NDVI están estandarizados y varían entre -1 y +1. Los valores más altos indican más cobertura verde.

Las celdas con valores de NDVI superiores a 0.4 son definitivamente 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')

Vamos a mapear el á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 = 'LAND?')

Estas pueden ser las áreas abiertas. Puede trazar landsobre el landsatFCCráster original para obtener más información.

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 puede 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')

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

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 otras PC explican la varianza adicional en orden decreciente.

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

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

La función prcomp() es una de las múltiples funciones en R que realizan PCA. Por defecto, prcomp() centra las variables para que tengan media cero, pero si se quiere además que su desviación estándar sea de uno, hay que indicar scale = TRUE.

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 destaca el segundo componente principal. Vamos a intentar umbral de nuevo:

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)

practica 4

Julieth Perez

24/3/2020

R Markdown

Including Plots