library(terra)
library(magrittr)
###Define Inputs----------------------------------------------------------------
#Image Directory Folder
Image_Directory = "C:/temp/S2A_MSIL1C_20180609T161901_N0206_R040_T16SGJ_20180609T194342.SAFE/GRANULE/L1C_T16SGJ_A015481_20180609T162519/IMG_DATA"
### Extacts all raster files from Image Directory with extention .jp2
Rasters = dir(Image_Directory, pattern = "*.jp2$", full.names = TRUE)
#Import Shapefile of Area of Interest
AOI = terra::vect("C:/temp/S2/eastfork_lake.shp")
#Reproject AOI if needed (Example - UTM Zone 16)
Projection = "+proj=utm +zone=16 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"
AOI = terra:project(AOI, Projection)
### Create Template for Final Raster Stack--------------------------------------
#For Sentinel-2 we want our final output to be in 20m so use Band 5
raster_B5 = raster(Rasters[[5]])
raster_template = terra::rast(raster_B5)
### Resample, Crop, & Stack All Images -----------------------------------------
raster_stack = Rasters %>%
lapply(terra::rast) %>%
lapply(terra::resample, raster_template) %>%
lapply(terra::crop, AOI) %>%
stack()
### Mask cropped image----------------------------------------------------------
#Further reducing size of stacked image
raster_mask <- terra::mask(raster_stack,AOI)
### Band Subset-----------------------------------------------------------------
S2_Harsha <- raster_mask[1:9] #Sentinel-2 Algorithms Only Use Bands 1-8A
### Save Final Stacked Image as Tiff -------------------------------------------
writeRaster(x = S2_Harsha,
filename= "C:/temp/S2_Harsha.tif", # save as a tif
datatype ="FLT4S", # save as a float 4 significant digits
overwrite = FALSE) #Overwrites same named file
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.