#ALUMNOS: #jhon jairo cuenca campos #angie dannielle cediel zabala

#install.packages("raster")
#install.packages("rgdal")
#install.packages("sp")

library(raster)
## Loading required package: sp
library(rgdal)
## rgdal: version: 1.5-23, (SVN revision 1121)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.2.1, released 2020/12/29
## Path to GDAL shared files: C:/Users/Jairo Cuenca/Documents/R/win-library/4.0/rgdal/gdal
## GDAL binary built with GEOS: TRUE 
## Loaded PROJ runtime: Rel. 7.2.1, January 1st, 2021, [PJ_VERSION: 721]
## Path to PROJ shared files: C:/Users/Jairo Cuenca/Documents/R/win-library/4.0/rgdal/proj
## PROJ CDN enabled: FALSE
## Linking to sp version:1.4-5
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading rgdal.
## Overwritten PROJ_LIB was C:/Users/Jairo Cuenca/Documents/R/win-library/4.0/rgdal/proj
library(sp)
#blue
b2 <- raster("C:/Users/Jairo Cuenca/Desktop/10 - DECIMO SEMESTRE/GEOMATICA/DATOS 1/rs/LC08_044034_20170614_B2.tif")
#green
b3 <- raster("C:/Users/Jairo Cuenca/Desktop/10 - DECIMO SEMESTRE/GEOMATICA/DATOS 1/rs/LC08_044034_20170614_B3.tif")
#red
b4 <- raster("C:/Users/Jairo Cuenca/Desktop/10 - DECIMO SEMESTRE/GEOMATICA/DATOS 1/rs/LC08_044034_20170614_B4.tif")
# near infrared (NIR)
b5 <- raster("C:/Users/Jairo Cuenca/Desktop/10 - DECIMO SEMESTRE/GEOMATICA/DATOS 1/rs/LC08_044034_20170614_B5.tif")
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 
## source     : C:/Users/Jairo Cuenca/Desktop/10 - DECIMO SEMESTRE/GEOMATICA/DATOS 1/rs/LC08_044034_20170614_B2.tif 
## names      : LC08_044034_20170614_B2 
## values     : 0.0748399, 0.7177562  (min, max)
b3
## 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 
## source     : C:/Users/Jairo Cuenca/Desktop/10 - DECIMO SEMESTRE/GEOMATICA/DATOS 1/rs/LC08_044034_20170614_B3.tif 
## names      : LC08_044034_20170614_B3 
## values     : 0.04259216, 0.6924697  (min, max)
b4
## 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 
## source     : C:/Users/Jairo Cuenca/Desktop/10 - DECIMO SEMESTRE/GEOMATICA/DATOS 1/rs/LC08_044034_20170614_B4.tif 
## names      : LC08_044034_20170614_B4 
## values     : 0.02084067, 0.7861769  (min, max)
b5
## 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 
## source     : C:/Users/Jairo Cuenca/Desktop/10 - DECIMO SEMESTRE/GEOMATICA/DATOS 1/rs/LC08_044034_20170614_B5.tif 
## names      : LC08_044034_20170614_B5 
## values     : 0.0008457669, 1.012432  (min, max)
#coordinate reference system (CRS)
crs(b2)
## CRS arguments:
##  +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs
#number of cells, rows, colums
ncell(b2)
## [1] 1863765
dim(b2)
## [1] 1245 1497    1
#spatial resolucion
res(b2)
## [1] 30 30
#number of nads
nlayers(b2)
## [1] 1
compareRaster(b2,b3)
## [1] TRUE
s <- stack(b5, b4,b3)
#check the properties of the resterstack
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 
## 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
par(mfrow = c(2,2))
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))

landsatRGB <- stack(b4, b3, b2)
plotRGB(landsatRGB, axes=TRUE, stretch="lin", main="landsat true color composite")

par(mfrow=c(1,2))
plotRGB(landsatRGB, axes=TRUE, stretch="lin", main="landsat true color composite")
landsatFCC <- stack(b5, b4, b3)
plotRGB(landsatFCC, axes=TRUE, stretch="lin", main="landsat false color composite")

filenames <- paste0('C:/Users/Jairo Cuenca/Desktop/10 - DECIMO SEMESTRE/GEOMATICA/DATOS 1/rs/LC08_044034_20170614_B', 1:11, ".tif") 
filenames
##  [1] "C:/Users/Jairo Cuenca/Desktop/10 - DECIMO SEMESTRE/GEOMATICA/DATOS 1/rs/LC08_044034_20170614_B1.tif" 
##  [2] "C:/Users/Jairo Cuenca/Desktop/10 - DECIMO SEMESTRE/GEOMATICA/DATOS 1/rs/LC08_044034_20170614_B2.tif" 
##  [3] "C:/Users/Jairo Cuenca/Desktop/10 - DECIMO SEMESTRE/GEOMATICA/DATOS 1/rs/LC08_044034_20170614_B3.tif" 
##  [4] "C:/Users/Jairo Cuenca/Desktop/10 - DECIMO SEMESTRE/GEOMATICA/DATOS 1/rs/LC08_044034_20170614_B4.tif" 
##  [5] "C:/Users/Jairo Cuenca/Desktop/10 - DECIMO SEMESTRE/GEOMATICA/DATOS 1/rs/LC08_044034_20170614_B5.tif" 
##  [6] "C:/Users/Jairo Cuenca/Desktop/10 - DECIMO SEMESTRE/GEOMATICA/DATOS 1/rs/LC08_044034_20170614_B6.tif" 
##  [7] "C:/Users/Jairo Cuenca/Desktop/10 - DECIMO SEMESTRE/GEOMATICA/DATOS 1/rs/LC08_044034_20170614_B7.tif" 
##  [8] "C:/Users/Jairo Cuenca/Desktop/10 - DECIMO SEMESTRE/GEOMATICA/DATOS 1/rs/LC08_044034_20170614_B8.tif" 
##  [9] "C:/Users/Jairo Cuenca/Desktop/10 - DECIMO SEMESTRE/GEOMATICA/DATOS 1/rs/LC08_044034_20170614_B9.tif" 
## [10] "C:/Users/Jairo Cuenca/Desktop/10 - DECIMO SEMESTRE/GEOMATICA/DATOS 1/rs/LC08_044034_20170614_B10.tif"
## [11] "C:/Users/Jairo Cuenca/Desktop/10 - DECIMO SEMESTRE/GEOMATICA/DATOS 1/rs/LC08_044034_20170614_B11.tif"
landsat <- stack(filenames)
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 
## 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
landsat <- subset(landsat, 1:7)
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"
#using extent
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)
x <- writeRaster(landsatcrop, filename="cropped-landsat.tif", overwrite=TRUE)
pairs(landsatcrop[[1:2]], main="ultra-blue versus blue")

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

pairs(landsatcrop[[5:6]], main="NIR versus SWIR1")

pairs(landsatcrop[[2:5]], main="blue versus NIR")