SpRaster

librerias

library(terra)

## terra version 1.0.10

library(raster)

## Loading required package: sp

library(sp)

entorno y procesamiento

input <- function(name,num,type,foldata,folder){
  ruta <- getwd()
  if (foldata == T){
    carpeta <- folder
    file <- paste0(ruta,'/',carpeta,'/',name, num, type)
  }else{
    file <- paste0(ruta,'/',name, num, type)
  }
}
bruto_landsat <- rast(input('LC08_L1TP_009057_20140820_20200911_02_T1_B',1:7,'.tif',T,'Datos'))
names(bruto_landsat) <- c('ultra-blue','blue','green','red','NIR','SWIR1','SWIR2')

Recorte

xmx <- 374838
xmn <- 358682
ymx <- 467578
ymn <- 454847
e <- ext(xmn, xmx, ymn, ymx)
# crop landsat by the extent
landsat <- crop(bruto_landsat, e) 
writeRaster(landsat, filename="mapa.tif", overwrite=TRUE)

Importar poligonos y el stack

LSstack <- brick(input('mapa','','.tif',F))
valle <- shapefile(input('valle','','.shp',T,'shapes'))

plotRGB(LSstack,r=5,g=5,b=5, axes = T)
plot(valle, add=T)

muestreo uniforme parecido a rasterisar

samples <- spsample(valle, 10000, type='regular')

## Warning in proj4string(obj): CRS object has comment, which is lost in output

numero <- length(samples)
df_samples <- as(samples,"SpatialPointsDataFrame")
df_samples@data = data.frame(ID=1:numero, size=1)
df_samples

## class       : SpatialPointsDataFrame 
## features    : 8821 
## extent      : 358747, 374782.4, 454862.2, 467527.1  (xmin, xmax, ymin, ymax)

## Warning in proj4string(x): CRS object has comment, which is lost in output

## crs         : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## variables   : 2
## names       :   ID, size 
## min values  :    1,    1 
## max values  : 8821,    1

plot

plot(LSstack$NIR)
plot(valle,  add=TRUE)
plot(df_samples, pch = 1, cex = (df_samples$size)/4, add=TRUE)

extraer valores del raster en los puntos de muestreo

df_samples$layer <- over(df_samples, valle)$layer
df <- raster::extract(LSstack, df_samples)
head(df)

##      ultra.blue  blue green   red   NIR SWIR1 SWIR2
## [1,]      10877 10271  9981 10478 16167 19614 14134
## [2,]      10951 10275 10005 10622 15323 19091 13957
## [3,]      10870 10164  9781 10119 16845 20308 14278
## [4,]      10348  9525  8902  8390 16399 14787 10704
## [5,]      10506  9809  9116  8883 15939 16301 11635
## [6,]      11228 10698 10413 11297 16321 21638 15886

perfiles espectrales :3

ms <- aggregate(df, list(df_samples$layer), mean)
rownames(ms) <- ms[,1]
ms <- ms[,-1]
ms

ms <- as.matrix(ms)

graficar

mycolor <- c('green', 'brown', 'blue', 'gray')

plot(0, ylim=c(0,18000), xlim = c(1,7), type='n', xlab="Bandas",
     ylab = "Reflectancia")

for (i in 1:nrow(ms)){
  lines(ms[i,], type = "l", lwd = 3, lty = 1, col = mycolor[i])
}
# Title
title(main="perfil espectral", font.main = 2)
# Legend
legend("topleft", rownames(ms),
       cex=0.8, col=mycolor, lty = 1, lwd =3, bty = "n")

reproducibilidad

sessionInfo()

## R version 4.0.4 (2021-02-15)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 17763)
## 
## Matrix products: default
## 
## locale:
## [1] LC_COLLATE=Spanish_Colombia.1252  LC_CTYPE=Spanish_Colombia.1252   
## [3] LC_MONETARY=Spanish_Colombia.1252 LC_NUMERIC=C                     
## [5] LC_TIME=Spanish_Colombia.1252    
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] raster_3.3-7 sp_1.4-2     terra_1.0-10
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_1.0.5        codetools_0.2-18  lattice_0.20-41   digest_0.6.25    
##  [5] grid_4.0.4        jsonlite_1.7.2    magrittr_1.5      evaluate_0.14    
##  [9] highr_0.8         rlang_0.4.10      stringi_1.4.6     vctrs_0.3.6      
## [13] rmarkdown_2.7     rgdal_1.5-23      tools_4.0.4       stringr_1.4.0    
## [17] xfun_0.21         yaml_2.2.1        compiler_4.0.4    htmltools_0.5.1.1
## [21] knitr_1.31