Práctica 3

Basic mathematical operations

library(sp)

library(raster)

raslist <- paste0("C:/Users/USS/Desktop/Pedro/2019 - 2020/Percepcion remota/P2/Practica/rsdata/rs/LC08_044034_20170614_B", 1:11, ".TIF")
raslist

 [1] "C:/Users/USS/Desktop/Pedro/2019 - 2020/Percepcion remota/P2/Practica/rsdata/rs/LC08_044034_20170614_B1.TIF" 
 [2] "C:/Users/USS/Desktop/Pedro/2019 - 2020/Percepcion remota/P2/Practica/rsdata/rs/LC08_044034_20170614_B2.TIF" 
 [3] "C:/Users/USS/Desktop/Pedro/2019 - 2020/Percepcion remota/P2/Practica/rsdata/rs/LC08_044034_20170614_B3.TIF" 
 [4] "C:/Users/USS/Desktop/Pedro/2019 - 2020/Percepcion remota/P2/Practica/rsdata/rs/LC08_044034_20170614_B4.TIF" 
 [5] "C:/Users/USS/Desktop/Pedro/2019 - 2020/Percepcion remota/P2/Practica/rsdata/rs/LC08_044034_20170614_B5.TIF" 
 [6] "C:/Users/USS/Desktop/Pedro/2019 - 2020/Percepcion remota/P2/Practica/rsdata/rs/LC08_044034_20170614_B6.TIF" 
 [7] "C:/Users/USS/Desktop/Pedro/2019 - 2020/Percepcion remota/P2/Practica/rsdata/rs/LC08_044034_20170614_B7.TIF" 
 [8] "C:/Users/USS/Desktop/Pedro/2019 - 2020/Percepcion remota/P2/Practica/rsdata/rs/LC08_044034_20170614_B8.TIF" 
 [9] "C:/Users/USS/Desktop/Pedro/2019 - 2020/Percepcion remota/P2/Practica/rsdata/rs/LC08_044034_20170614_B9.TIF" 
[10] "C:/Users/USS/Desktop/Pedro/2019 - 2020/Percepcion remota/P2/Practica/rsdata/rs/LC08_044034_20170614_B10.TIF"
[11] "C:/Users/USS/Desktop/Pedro/2019 - 2020/Percepcion remota/P2/Practica/rsdata/rs/LC08_044034_20170614_B11.TIF"

landsat <- stack(raslist) # Creamos un objeto con las imagenes raster
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 

landsatRGB <- landsat[[c(4,3,2)]] # Creamos un objeto con las capas RGB de landsat
landsatRGB

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_B4, LC08_044034_20170614_B3, LC08_044034_20170614_B2 
min values :              0.02084067,              0.04259216,              0.07483990 
max values :               0.7861769,               0.6924697,               0.7177562 

landsatFCC <- landsat[[c(5,4,3)]] # Creamos un objeto con las capas FCC (Falso color) de landsat
landsatFCC

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 

Vegetation indices

# Se crea la funcion "vi" (indice de vegetacion) img es un  objeto mutilayer Raster y k e i son los indices de las capas

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

# Para Landsat NIR = 5, red = 4.

ndvi <- vi(landsat, 5, 4)

plot(ndvi, col = rev(terrain.colors(10)), main = "Landsat-NDVI")

# Otra forma de obtener 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")

Question 1

# Para insertar una imagen

knitr::include_graphics("C:/Users/USS/Desktop/Pedro/2019 - 2020/Percepcion remota/P3/perfil.png")  

# Segun el perfil espectral, utilizando la banda 1 y 7, el agua es la única superficie con el indice entre 0.5 y 1 (la reflectancia es mayor en la banda 1 que en la 7), ya que las tierras de cultivo tienen un indice cercano a 0 (la reflectancia es prácticamente igual en la banda 1 que en la 7) y el resto tienen un indice negativo (la reflectancia es mayor en la banda 7 que en la 1)

ndwi <- vi(landsat, 1, 7)

plot(ndwi, col = rev(topo.colors(255)), main = "NDWI")

# Otra forma de obtener el NDWI es NDWI = (green-NIR)/ (green+NIR). En landsat 8 estas bandas son la 3 (green) y la 5 (NIR)

ndwi2 <- vi(landsat, 3, 5)
plot(ndwi2, col = bpy.colors(255), main = "NDWI")

# Para el terreno construido se utiliza el NDBI, NDBI= (SWIR-NIR)/ (SWIR+NIR). Para landsat 8 estas bandas son la 6 SWIR y 5 la NIR

ndbi <- vi(landsat, 6, 5)
plot(ndbi, col = bpy.colors(255), main = "NDBI")

# Debe haber algún error en algún dato de esa capa, ya que da un rango de valores por encima del esperado (-1,1), impidiendo su identificación visual mediante NDBI

Histogram

# Para Visualizar el histograma de los 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))

Question 2

# Visualización del histograma de el índice NDWI para identificar agua

hist(ndwi,
     main = "Distribution of NDWI (water, layers 1 and 7) values",
     xlab = "NDWI",
     ylab= "Frequency",
     col = "blue",
     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))

# Visualización del histograma de el índice NDBI para identificar terreno construido

hist(ndbi,
     main = "Distribution of NDBI (Built-up, layers 6 and 5) values",
     xlab = "NDBI",
     ylab= "Frequency",
     col = "red",
     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))

Thresholding

# Las celdas con valores NDVI superiores a 0.4 son definitivamente vegetación, por lo que eliminamos los valores por debajo de 0.4 NDVI, para visualizar mejor la vegetación

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

# Representamos el pico entre 0.25 y 0.3 en el histograma

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

# Puede tratarse de tierras en barbecho. Podemos plotear el raster land sobre landsatFCC para visualizarlo

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

Question 3

water <- reclassify(ndwi, cbind(-Inf, 0.5, NA))
colors <- colorRampPalette(c("grey","deepskyblue1","blue"))
plot(water,col = colors(50), main = 'water')

ndbi2 <- reclassify(ndbi, c(-Inf, -1, NA, 1, Inf, NA))
plot(ndbi2, col = rev(bpy.colors(255)), main='ndbi2')

# Los terrenos construidos pueden tratarse de las zonas rojas

colors <- colorRampPalette(c("white","grey","wheat","red"))
plot(ndbi2,col = colors(255), main = 'water')

Principal component analysis

set.seed(1) # set.seed() garantiza que obtenga el mismo resultado si comienza con la misma semilla (seed) cada vez que ejecuta el mismo proceso 
sr <- sampleRandom(landsat, 10000) # Se crean puntos aleatorios
plot(sr[,c(4,5)], main = "NIR-Red plot")

# prcomp realiza un análisis de componentes principales en la matriz de datos dada y devuelve los resultados como un objeto de clase prcomp

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

Standard deviations (1, .., p=11):
 [1] 2.52687022 1.40533428 1.09901821 0.92507423 0.53731672 0.42657919 0.28002273 0.12466139
 [9] 0.09031384 0.04761419 0.03609857

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

plot(pca)

# Se puede apreciar que la mayor varianza está en los dos primeros componentes principales, por lo que son los que contienen mayor información

# Predecimos un nuevo raster a traves del análisis de componentes principales, basandose en landsat 

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

pci

class      : RasterBrick 
dimensions : 1245, 1497, 1863765, 2  (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 
source     : memory
names      :        PC1,        PC2 
min values :  -6.334633, -24.223508 
max values :   41.63416,   10.23197 

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 

# Para entender lo que destaca el segundo componente principal se utiliza nuevamente thresholding

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)