RASTER

LIBRERIA

library(raster)

## Loading required package: sp

library(terra)

## terra version 1.0.10

DESCARGA

dir.create('data', showWarnings = FALSE)
if (!file.exists('data/rs/samples.rds')) {
    download.file('https://biogeo.ucdavis.edu/data/rspatial/rsdata.zip', dest = 'data/rsdata.zip')
    unzip('data/rsdata.zip', exdir='data')
} 

ENTORNO DE TRABAJO E IMPORTAR DATOS

ruta <- getwd()
carpeta <- 'Datos'
banda1 <- 'LC08_044034_20170614_B1.tif'
banda2 <- 'LC08_044034_20170614_B2.tif'
banda3 <- 'LC08_044034_20170614_B3.tif'
banda4 <- 'LC08_044034_20170614_B4.tif'
banda5 <- 'LC08_044034_20170614_B5.tif'
banda6 <- 'LC08_044034_20170614_B6.tif'
banda7 <- 'LC08_044034_20170614_B7.tif'
banda8 <- 'LC08_044034_20170614_B8.tif'
banda9 <- 'LC08_044034_20170614_B9.tif'
banda10 <- 'LC08_044034_20170614_B10.tif'
banda11 <- 'LC08_044034_20170614_B11.tif'
 
B1 <- raster(paste0(ruta,'/',carpeta,'/',banda1))
B2 <- raster(paste0(ruta,'/',carpeta,'/',banda2))
B3 <- raster(paste0(ruta,'/',carpeta,'/',banda3))
B4 <- raster(paste0(ruta,'/',carpeta,'/',banda4))
B5 <- raster(paste0(ruta,'/',carpeta,'/',banda5))
B6 <- raster(paste0(ruta,'/',carpeta,'/',banda6))
B7 <- raster(paste0(ruta,'/',carpeta,'/',banda7))
B8 <- raster(paste0(ruta,'/',carpeta,'/',banda8))
B9 <- raster(paste0(ruta,'/',carpeta,'/',banda9))
B10 <- raster(paste0(ruta,'/',carpeta,'/',banda10))
B11 <- raster(paste0(ruta,'/',carpeta,'/',banda11))

ACCEDER A PROPIEDADES

crs(B2)

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

ncell(B2)

## [1] 1863765

dim(B2)

## [1] 1245 1497    1

res(B2)

## [1] 30 30

nlayers(B2)

## [1] 1

compareRaster(B2, B3)

## [1] TRUE

RASTERSTACK (VARIAS CAPAS)

s <- stack(B5,B4,B3)
nlayers(s)

## [1] 3

#############UN SOLO RASTER DE 11 BANDAS##########
file <- paste0(ruta,'/',carpeta,'/','LC08_044034_20170614_B', 1:11, ".tif")
landsat <- stack(file)
nlayers(landsat)

## [1] 11

landsat <- subset(landsat,1:7)#eliminar banda del 8 al 11
names(landsat) <- c('ultra-blue','blue','green','red','NIR','SWIR1','SWIR2')

MAPAS DE UNA CAPA

par(mfrow = c(2,2))
plot(B2, main = "Blue", col = rgb((0:100)/100, red = 0, green = 0))
plot(B3, main = "Green", col = rgb((0:100 / 100), red = 0, blue = 0))
plot(B4, main = "Red", col = rgb((0:100 / 100),green = 0, blue = 0))
plot(B5, main = "NIR", col = rgb((0:100 / 100), red = 1, green=0,blue=0))

MAPA COLOR NATURAL

par(mfrow = c(1,2))
bandasRGB <- stack(B4, B3, B2)
plotRGB(bandasRGB, axes = TRUE, stretch = "lin", main = "LandSat natural Bandas")
########################
landsatRGB <- stack(landsat$red,landsat$green,landsat$blue)
plotRGB(landsatRGB, axes = TRUE, stretch = "lin", main = "Landsat natural Capas")

FALSO COLOR

par(mfrow = c(1,2))
bandasFCC <- stack(B5, B4, B3)
plotRGB(bandasFCC, axes=TRUE, stretch="lin", main="Landsat Falso Color Bandas")
##################################
landsatFCC <- stack(landsat$NIR,landsat$red,landsat$green)
plotRGB(landsatFCC, axes=TRUE, stretch="lin", main="Landsat Falso Color Capas")

FUNCION EXTENT

e <- extent(624387, 635752, 4200047, 4210939)
landsatcrop <- crop(landsat, e)
landsatcropRGB <- stack(landsatcrop$red,landsatcrop$green,landsatcrop$blue)
####################
longitude <- c(624387,624387,635752,635752)
latitude <- c(4200047,4210939,4210939,4200047)
lonlat <- cbind(longitude, latitude)
pols <- vect(lonlat, type='polygons')
####################
par(mfrow = c(1,2))
plotRGB(landsatRGB, axes = TRUE, stretch = "lin", main = "Landsat natural")
plot(pols,las=1,borde = 'red' , add=T)
plotRGB(landsatcropRGB, axes = TRUE, stretch = "lin", main = "Landsat recorte")