R/sen_clip_stack.R

In envima/GEDItools:

#' Clip and stack bands from Sentinel 2a
#'
#' This function allows you to clip and stack several bands from sentinel 2a data
#' @param inpath characterstring to unzipped folder from sentinel download
#' @param xmin numeric, minimal longitude of study area
#' @param xmax numeric, maximal longitude of study area
#' @param ymin numeric, minimal latitude of study area
#' @param ymax numeric, maximal latitude of study area
#' @param write_out boolean write out sentinel stack as tif and rds
#'
#' @keywords clip stack sentinel
#' @export
#' @examples
#' sen_clip_stack
#'
#' @import raster
#' @import sf
#'

##university forest
# xmin = 8.661830
# xmax = 8.694962
# ymin = 50.831518
# ymax = 50.848372

# ##larger MR area
# ymax <- 50.848933
# ymin <- 50.740749
# xmax <- 8.820696
# xmin <- 8.516494

##Madagaskar
# xmin = 43.645483310
# xmax = 43.778625076
# ymin = -24.229014539
# ymax = -23.726632957
#
#
#
#
#

# inpath <- "C:/Users/Alice/Uni/Projekte/GEDI/data/Sentinel/S2A_MSIL2A_20190420T103031_N0211_R108_T32UMB_20190420T111631"
# inpath <- "C:/Users/Alice/Uni/Diverses/Madagaskar/S2A_MSIL2A_20190613T071211_N0212_R020_T38KLU_20190613T103107"

sen_clip_stack <- function(inpath, xmin = 0, xmax = 0, ymin = 0, ymax = 0, crop = T, write_out = T){

########################################################################################
### settings
########################################################################################
  #####
  ###folder and name extraction
  #####
  print(paste0("Processing: ", inpath))
  basenm <- basename(inpath)
  granule_dir <- list.dirs(paste0(inpath, "/", basenm, ".SAFE/GRANULE/"), recursive = F, full.names = F)
  inpath_sen <- paste0(inpath, "/", basenm, ".SAFE/GRANULE/", granule_dir, "/IMG_DATA/")
  setwd(inpath_sen)

  granule_splt <- strsplit(granule_dir, "_")
  basenm_splt <- strsplit(basenm, "_")
  granule_nm <- paste0(granule_splt[[1]][2], "_", basenm_splt[[1]][3])


  #####
  ###study area as sf utm
  #####
  sen1 <- raster(paste0("R60m/", granule_nm, "_B01_60m.jp2"))
  area_utm <- extent(sen1)
  crs(sen1)
  if(crop ==T){
  study_area <- data.frame(xmin = xmin, xmax = xmax, ymin = ymin, ymax = ymax) # #xmin, xmax, ymin, ymax
  mat <- matrix(c(study_area$xmax, study_area$ymax,
                  study_area$xmin, study_area$ymax,
                  study_area$xmin, study_area$ymin,
                  study_area$xmax, study_area$ymin,
                  study_area$xmax, study_area$ymax),  ## need to close the polygon
                ncol =2, byrow = T)
  ## create polygon objects
  pol <- st_polygon(list(mat))
  ## create sf object
  area_wgs <- st_sf(study_area, st_sfc(pol), crs = "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")
  ## transofrm sf object
  area_utm <- st_transform(area_wgs, crs = crs(sen1))
  }
  #####
  ###read files
  #####

  #load and crop sentinel bands
  bd1 <- crop(raster(paste0("R60m/", granule_nm, "_B01_60m.jp2")), area_utm)
  bd2 <- crop(raster(paste0("R10m/", granule_nm, "_B02_10m.jp2")), area_utm)
  bd3 <- crop(raster(paste0("R10m/", granule_nm, "_B03_10m.jp2")), area_utm)
  bd4 <- crop(raster(paste0("R10m/", granule_nm, "_B04_10m.jp2")), area_utm)
  bd5 <- crop(raster(paste0("R20m/", granule_nm, "_B05_20m.jp2")), area_utm)
  bd6 <- crop(raster(paste0("R20m/", granule_nm, "_B06_20m.jp2")), area_utm)
  bd7 <- crop(raster(paste0("R20m/", granule_nm, "_B07_20m.jp2")), area_utm)
  bd8 <- crop(raster(paste0("R10m/", granule_nm, "_B08_10m.jp2")), area_utm)
  bd8a<- crop(raster(paste0("R20m/", granule_nm, "_B8A_20m.jp2")), area_utm)
  bd9 <- crop(raster(paste0("R60m/", granule_nm, "_B09_60m.jp2")), area_utm)
  bd11<- crop(raster(paste0("R20m/", granule_nm, "_B11_20m.jp2")), area_utm)
  bd12<- crop(raster(paste0("R20m/", granule_nm, "_B12_20m.jp2")), area_utm)


  ########################################################################################
  ########################################################################################
  ########################################################################################
  ###Do it (Don't change anything past this point except you know what you are doing!) ###
  ########################################################################################
  ########################################################################################
  ########################################################################################

  #calc NDVI
  NDVI <- (bd8 - bd4)/(bd8 + bd4)
  names(NDVI) <- "NDVI"

  # resample
  bd1 <- resample(bd1, bd2)
  bd5 <- resample(bd5, bd2)
  bd6 <- resample(bd6, bd2)
  bd7 <- resample(bd7, bd2)
  bd8a <- resample(bd8a, bd2)
  bd9 <- resample(bd9, bd2)
  bd11 <- resample(bd11, bd2)
  bd12 <- resample(bd12, bd2)

  #stack
  sen_stack <- stack(bd1, bd2, bd3, bd4, bd5, bd6, bd7, bd8, bd8a, bd9, bd11, bd12, NDVI)
  rm(list = c("bd1", "bd2", "bd3", "bd4", "bd5", "bd6", "bd7", "bd8", "bd8a", "bd9", "bd11", "bd12"))

  if(write_out ==T){
    if (file.exists(paste0(inpath, "/clipped"))==F){
      dir.create(paste0(inpath, "/clipped"), recursive = T)
    }
    #write RDS
    saveRDS(sen_stack, file = paste0(inpath, "/clipped/sen_stack_", basenm_splt[[1]][3], "_123456788a91112NDVI.rds"))
    writeRaster(sen_stack, file = paste0(inpath, "/clipped/sen_stack_", basenm_splt[[1]][3], "_123456788a91112NDVI.tif"), overwrite = T)
  }


  return(sen_stack)
}