Practica 2 percepción remota

getwd()

list.files(path = ".")

library(raster)

library(sp)

# Blue
b2 <- raster('C:\\Users\\USS\\Desktop\\Pedro\\2019 - 2020\\Percepcion remota\\P2\\Practica\\rsdata\\rs\\LC08_044034_20170614_B2.TIF') #con \\ en vez de /, importante poner ' en la ruta
b2

# Green
b3 <- raster("C:/Users/USS/Desktop/Pedro/2019 - 2020/Percepcion remota/P2/Practica/rsdata/rs/LC08_044034_20170614_B3.TIF") #con / en vez de \\
b3

# Red
b4 <- raster("C:/Users/USS/Desktop/Pedro/2019 - 2020/Percepcion remota/P2/Practica/rsdata/rs/LC08_044034_20170614_B4.TIF")
b4

# Near Infrared (NIR)
b5 <- raster("C:/Users/USS/Desktop/Pedro/2019 - 2020/Percepcion remota/P2/Practica/rsdata/rs/LC08_044034_20170614_B5.TIF")
b5

crs(b2) #Para ver sistema de coordenadas de referencia

ncell(b2) #Para ver el número de celdas

nrow(b2) #Para ver el número de filas

ncol(b2) #Para ver el número de columnas

dim(b2) #Para ver la dimension

res(b2) #Para ver la resolución espacial

nlayers(b2) #Número de bandas

compareRaster(b2,b3) #Compara si las bandas tienen la misma extensión, número de filas y columnas, proyección, resolución y origen

s <- stack(b5, b4, b3) #Crear RasterStack con múltiples bandas a través de RasterLayers (bandas simples) existentes 
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 

#Otra forma de crear RasterStack

filenames <- paste0("C:/Users/USS/Desktop/Pedro/2019 - 2020/Percepcion remota/P2/Practica/rsdata/rs/LC08_044034_20170614_B", 1:11, ".TIF") #crear una lista de RastersLayers para usar 
filenames

 [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(filenames) #Crear RasterStack a través de los nombres de los archivos
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 

compareRaster(s,landsat)

[1] TRUE

par(mfrow = c(2,2)) #Para dividir la pantalla plot en 2 filas y 2 columnas
plot(b2, main = "Blue", col = gray(0:100 / 100))
plot(b3, main = "Green", col = gray(0:100 / 100))

plot(b4, main = "Red", col = gray(0:100 / 100))
plot(b5, main = "NIR", col = gray(0:100 / 100))

plot(b2, main = "Blue", col = gray(0:100 / 100)) #Se ve mejor si se grafica solo una banda

landsatRGB <- stack(b4, b3, b2)
plotRGB(landsatRGB, axes = TRUE, stretch = "lin", main = "Landsat True Color Composite")

landsatFCC <- stack(b5, b4, b3) #Imagen de "color falso" en la que se combinan las bandas NIR, rojo y verde, se aprecia bien la vegetación en rojo
plotRGB(landsatFCC, axes=TRUE, stretch="lin", main="Landsat False Color Composite")

landsatsub1 <- subset(landsat, 1:3) # Para crear un objeto con 3 bandas de otro objeto (landsat)
landsatsub1

landsatsub2 <- landsat[[1:3]] # Otro metodo
landsatsub2

# Comprobacion de numero de bandas de los objetos
nlayers(landsat) 
nlayers(landsatsub1)
nlayers(landsatsub2)

landsat <- subset(landsat, 1:7) # Se crea un objeto con las primeras siete bandas de landsat

nlayers(landsat)

names(landsat) # Para ver los nombres de las bandas de landsat

[1] "LC08_044034_20170614_B1" "LC08_044034_20170614_B2" "LC08_044034_20170614_B3"
[4] "LC08_044034_20170614_B4" "LC08_044034_20170614_B5" "LC08_044034_20170614_B6"
[7] "LC08_044034_20170614_B7"

names(landsat) <- c('ultra-blue', 'blue', 'green', 'red', 'NIR', 'SWIR1', 'SWIR2') #Para cambiar los nombres de las bandas
names(landsat)

[1] "ultra.blue" "blue"       "green"      "red"        "NIR"        "SWIR1"      "SWIR2"     

extent(landsat) # Funcion extent

class      : Extent 
xmin       : 594090 
xmax       : 639000 
ymin       : 4190190 
ymax       : 4227540 

e <- extent(624387, 635752, 4200047, 4210939) #Se crea un objeto extent para cortar
e

class      : Extent 
xmin       : 624387 
xmax       : 635752 
ymin       : 4200047 
ymax       : 4210939 

landsatcrop <- crop(landsat, e) # Se corta landsat con la zona delimitada por e
landsatcrop

class      : RasterBrick 
dimensions : 363, 379, 137577, 7  (nrow, ncol, ncell, nlayers)
resolution : 30, 30  (x, y)
extent     : 624390, 635760, 4200060, 4210950  (xmin, xmax, ymin, ymax)
crs        : +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
source     : memory
names      :  ultra.blue,        blue,       green,         red,         NIR,       SWIR1,       SWIR2 
min values : 0.101796150, 0.075989284, 0.044045158, 0.030556191, 0.007243267, 0.001865030, 0.003144530 
max values :   0.4548080,   0.4746728,   0.5034941,   0.5592065,   0.7013395,   0.9038041,   0.9750224 

x <- writeRaster(landsatcrop, filename="cropped-landsat.tif", overwrite=TRUE) # Para guardar el raster recortado en formato .tif

y <-writeRaster(landsatcrop, filename="cropped-landsat.grd", overwrite=TRUE) # Para guardar el raster recortado en formato .grd (formato raster-grd)

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

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

# Se cargan los poligonos con información sobre uso del suelo y cobertura del suelo
samp <- readRDS("C:/Users/USS/Desktop/Pedro/2019 - 2020/Percepcion remota/P2/Practica/rsdata/rs/samples.rds")
samp

class       : SpatialPolygonsDataFrame 
features    : 49 
extent      : 594416.7, 638822.2, 4190232, 4224228  (xmin, xmax, ymin, ymax)
crs         : +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
variables   : 2
names       : id, class 
min values  :  1, built 
max values  :  5, water 

plot(samp)

# Se generan muestras de 300 puntos de los polígonos
ptsamp <- spsample(samp, 300, type='regular')
ptsamp

class       : SpatialPoints 
features    : 286 
extent      : 594996.8, 638799.8, 4190289, 4223952  (xmin, xmax, ymin, ymax)
crs         : +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 

# Para agregar la clase de samp a ptsamp
ptsamp$class <- over(ptsamp, samp)$class
list(ptsamp$class)

[[1]]
  [1] "fallow"   "fallow"   "fallow"   "open"     "open"     "open"     "fallow"   "fallow"  
  [9] "fallow"   "open"     "open"     "open"     "open"     "cropland" "fallow"   "fallow"  
 [17] "open"     "open"     "open"     "open"     "open"     "cropland" "cropland" "fallow"  
 [25] "fallow"   "open"     "open"     "open"     "open"     "open"     "open"     "open"    
 [33] "open"     "open"     "open"     "open"     "open"     "open"     "open"     "open"    
 [41] "open"     "open"     "open"     "open"     "open"     "open"     "open"     "open"    
 [49] "cropland" "cropland" "open"     "open"     "open"     "open"     "open"     "open"    
 [57] "open"     "open"     "open"     "cropland" "cropland" "open"     "open"     "open"    
 [65] "open"     "open"     "open"     "open"     "open"     "open"     "open"     "cropland"
 [73] "cropland" "cropland" "fallow"   "open"     "cropland" "cropland" "cropland" "fallow"  
 [81] "fallow"   "open"     "open"     "open"     "open"     "cropland" "fallow"   "open"    
 [89] "open"     "open"     "open"     "open"     "fallow"   "fallow"   "open"     "open"    
 [97] "fallow"   "fallow"   "fallow"   "fallow"   "fallow"   "fallow"   "fallow"   "fallow"  
[105] "built"    "built"    "built"    "built"    "open"     "built"    "fallow"   "open"    
[113] "fallow"   "fallow"   "fallow"   "fallow"   "open"     "fallow"   "fallow"   "open"    
[121] "open"     "fallow"   "open"     "open"     "open"     "built"    "built"    "built"   
[129] "built"    "built"    "built"    "built"    "built"    "cropland" "cropland" "cropland"
[137] "cropland" "built"    "built"    "built"    "built"    "built"    "built"    "built"   
[145] "cropland" "water"    "cropland" "cropland" "cropland" "water"    "water"    "water"   
[153] "water"    "built"    "built"    "water"    "water"    "water"    "water"    "water"   
[161] "water"    "water"    "water"    "water"    "water"    "water"    "water"    "cropland"
[169] "water"    "water"    "water"    "water"    "water"    "water"    "water"    "water"   
[177] "water"    "cropland" "cropland" "water"    "water"    "water"    "water"    "water"   
[185] "water"    "water"    "water"    "water"    "water"    "water"    "water"    "water"   
[193] "water"    "cropland" "cropland" "water"    "water"    "water"    "water"    "water"   
[201] "water"    "water"    "water"    "water"    "water"    "water"    "water"    "water"   
[209] "cropland" "water"    "water"    "water"    "water"    "water"    "water"    "water"   
[217] "water"    "water"    "water"    "water"    "water"    "water"    "water"    "water"   
[225] "water"    "water"    "water"    "water"    "water"    "water"    "water"    "water"   
[233] "water"    "water"    "water"    "water"    "water"    "water"    "water"    "water"   
[241] "water"    "water"    "water"    "water"    "water"    "water"    "cropland" "cropland"
[249] "cropland" "cropland" "water"    "water"    "cropland" "water"    "water"    "water"   
[257] "water"    "water"    "water"    "water"    "water"    "cropland" "cropland" "water"   
[265] "water"    "water"    "water"    "water"    "water"    "water"    "water"    "water"   
[273] "water"    "water"    "water"    "water"    "water"    "water"    "water"    "water"   
[281] "water"    "water"    "water"    "water"    "water"    "water"   

# Para extraer los valores donde hay puntos
df <- extract(landsat, ptsamp)
head(df) #Para visualizar los primeros valores

     ultra.blue      blue     green       red       NIR     SWIR1     SWIR2
[1,]  0.1410052 0.1248705 0.1113816 0.1148080 0.1793469 0.2322618 0.1928575
[2,]  0.1346511 0.1186031 0.1097117 0.1124225 0.2006646 0.2405243 0.2032236
[3,]  0.1372535 0.1197959 0.1036178 0.1024034 0.1818842 0.2088187 0.1686338
[4,]  0.1298801 0.1218344 0.1257163 0.1667470 0.3061041 0.3189208 0.1924455
[5,]  0.1289693 0.1201212 0.1211838 0.1598074 0.2862176 0.3166437 0.1924021
[6,]  0.1480100 0.1449739 0.1573568 0.2008598 0.3047379 0.3821368 0.2652685

ms <- aggregate(df, list(ptsamp$class), mean) # Se calculan los valores medios de reflectancia para cada clase y cada banda.
ms

# Ponemos los nombres de fila de ms con los valores de la primera columna de ms
rownames(ms) <- ms[,1]
rownames(ms)

[1] "built"    "cropland" "fallow"   "open"     "water"   

ms <- ms[,-1] # Quitamos la primera columna de ms
ms

# Se crea un vector de color para las clases de cubierta de tierra para usar en el plot
mycolor <- c('darkred', 'yellow', 'burlywood', 'cyan', 'blue')

# Trasformamos ms de data.frame a matriz
ms <- as.matrix(ms)

# Primero se crea un plot vacío
plot(0, ylim=c(0,0.6), xlim = c(1,7), type='n', xlab="Bands", ylab = "Reflectance")

# Se añaden 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")

#Si no se ejecuta todo el código junto da error