Exploracion de imagenes satelitales

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

## Loading required package: sp

library(sp)
library(rgdal)

## Please note that rgdal will be retired by the end of 2023,
## plan transition to sf/stars/terra functions using GDAL and PROJ
## at your earliest convenience.
## 
## rgdal: version: 1.5-25, (SVN revision 1143)
## 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/NICOLAS/Documents/R/win-library/4.1/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/NICOLAS/Documents/R/win-library/4.1/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 sp or rgdal.
## Overwritten PROJ_LIB was C:/Users/NICOLAS/Documents/R/win-library/4.1/rgdal/proj

setwd("D:/mis documentos/ingenieria civil/semestre 10/Electiva 2/Taller 2")
# Llamar imagines satelitales
# Azul banda 2
b2 <- raster("D:/mis documentos/ingenieria civil/semestre 10/Electiva 2/Taller 2/LC08_L1TP_008058_20210731_20210810_02_T1_B2.tif")
# Verde banda 3
b3 <- raster("D:/mis documentos/ingenieria civil/semestre 10/Electiva 2/Taller 2/LC08_L1TP_008058_20210731_20210810_02_T1_B3.tif")
# Rojo banda 4
b4 <- raster("D:/mis documentos/ingenieria civil/semestre 10/Electiva 2/Taller 2/LC08_L1TP_008058_20210731_20210810_02_T1_B4.tif")
b5 <- raster("D:/mis documentos/ingenieria civil/semestre 10/Electiva 2/Taller 2/LC08_L1TP_008058_20210731_20210810_02_T1_B5.tif")

# Informacion de la imagen
b2

## class      : RasterLayer 
## dimensions : 7741, 7581, 58684521  (nrow, ncol, ncell)
## resolution : 30, 30  (x, y)
## extent     : 411585, 639015, 203385, 435615  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## source     : LC08_L1TP_008058_20210731_20210810_02_T1_B2.TIF 
## names      : LC08_L1TP_008058_20210731_20210810_02_T1_B2 
## values     : 0, 65535  (min, max)

crs(b2) # consultar el sistema de referencia

## CRS arguments:
##  +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs

ncell(b2) # Consultar numero de celdas

## [1] 58684521

dim(b2) # Dimensiones, es decir alto x ancho

## [1] 7741 7581    1

res(b2) # Resolucion

## [1] 30 30

nlayers(b2) # Numero de capas

## [1] 1

compareRaster(b2,b3) # Comparar 2 imagenes

## [1] TRUE

s <- stack(b5,b4,b3)
s

## class      : RasterStack 
## dimensions : 7741, 7581, 58684521, 3  (nrow, ncol, ncell, nlayers)
## resolution : 30, 30  (x, y)
## extent     : 411585, 639015, 203385, 435615  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## names      : LC08_L1TP_008058_20210731_20210810_02_T1_B5, LC08_L1TP_008058_20210731_20210810_02_T1_B4, LC08_L1TP_008058_20210731_20210810_02_T1_B3 
## min values :                                           0,                                           0,                                           0 
## max values :                                       65535,                                       65535,                                       65535

par(mfrow = c(2,2))
plot(b2, main = "azul", col = gray(0:100 / 100))
plot(b3, main = "verde", col = gray(0:100 / 100))
plot(b4, main = "rojo", col = gray(0:100 / 100))
plot(b5, main = "infra rojo", col = gray(0:100 / 100))

# Combinacion para lograr el color verdadero
landsatRGB <- stack(b4,b3,b2)
plotRGB(landsatRGB, axes = TRUE, stretch = "lin", main = "Landsat color verdadero")