GEDItools:

#' Clip and stack bands from Sentinel 1 grd
#'
#' This function allows you to clip and stack several bands from sentinel 1 data
#' @param inpath characterstring to unzipped folder from sentinel download .SAFE
#' @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 sentinel1
#' @export
#' @examples
#' sen1_clip_stack
#'
#' @import raster
#' @import sf
#' @import gdalUtils
#'
#'
# 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


###Needs to be implemented and run with snapgraphsvia termial

after_snap <- raster("C:/Users/Alice/Uni/Projekte/GEDI/data/Sentinel1/snap/S1_GRDH_20190628.tif")
# stack_after_snap <- stack("C:/Users/Alice/Uni/Projekte/GEDI/data/Sentinel1/snap/S1_GRDH_20190628.tif")
after_snap <- raster("C:/Users/Alice/Uni/Projekte/GEDI/data/Sentinel1/snap/Subset_S1A_IW_GRDH_1SDV_20190625T170813_20190625T170838_027839_032493_8A5E_Orb_NR_Cal_TF_TC.tif")
mapview(after_snap)
# -> snap in terminal works ==> needs to be controlled via R and

# #
# #
# #
# # a <- raster(inpath)
# sen1_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)
#   inpath_sen <- paste0(list.dirs(inpath, recursive = F), "/measurement/")
#   # granule_dir <- list.dirs(paste0(inpath, "/", basenm, ".SAFE/GRANULE/"), recursive = F, full.names = F)
#   # inpath_sen <- paste0(inpath, "/", basenm, ".SAFE/measurement/")
#   setwd(inpath_sen)
#
#   # granule_splt <- strsplit(granule_dir, "_")
#   # basenm_splt <- strsplit(basenm, "_")
#   # granule_nm <- paste0(granule_splt[[1]][2], "_", basenm_splt[[1]][3])
#
#   sen1_files <- list.files(inpath_sen)
#   #get names for vv and vh
#   nm_vv <- sen1_files[grepl("vv", sen1_files)]
#   nm_vh <- sen1_files[grepl("vh", sen1_files)]
# # setwd("C:/Users/Alice/Desktop/temp/")
#   # raster("s1a-iw-grd-vh-20190522t054201-20190522t05.tiff")
#   #####
#   ###study area as sf utm
#   #####
#   sen1 <- raster(paste0(inpath_sen, nm_vv), crs = "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0" )
#   sen1 <- raster(paste0(inpath_sen, nm_vv))
#
#   # 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))
#
#     wkt <- st_as_text(pol)
#     ## 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
#   bd_vv <- crop(raster(paste0(inpath_sen, nm_vv)), area_wgs)
#   bd_vh <- crop(raster(paste0(inpath_sen, nm_vh)), 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)
# }
#