###############################################################################################################
##############################  Environmental Sampling
###############################################################################################################

## Load required packages

library(raster)

## load data

filelist <- c(list.files('.', '*tif$'), list.files('.', '*grd$')) # list of raster files path, detected by extension (.tif or .grd).

# Create an object list containing all raster layers
raslist <- list()
for (i in filelist) {raslist <- c(raslist, raster(i)) }
names(raslist) <- filelist
# All the raster layers aree contained in the raslist list. Stack can not be used since rasters carry different extent and resolution. 

samples <- shapefile('./MOOA_ENV.shp')

# The shapefile function allows to import vector data. 
samples # shows the different properties of the object of the class 'spatialpointsdataframe'


plot(raslist$MAR_alt.grd) # plot altitude
points(samples) # add samples coordinates to the plot. 

### Sample Environmental Data

## The raster package has a function called 'extract' that can sample values from a raster based on a georefrenced object. 
## The extraction can be point (at a determined location) or area based (around a determined location).
## We will use both methods, the first one will perform the exact same operation that we performed in QGIS with the Point Sampling Tool.
## The second one performs what we did using the python scripts calling the Zonal Statistics function. 

# Point Sampling
env <- c()
for (i in 1:length(raslist)) { env <- cbind(env, extract(raslist[[i]], samples))} # fills matrix env with point sampling values for each variable
colnames(env) <- names(raslist)
point_samples <- cbind(samples, env) # add environmental values to sample coordinates

# Area Sampling, 10 km : this might take some time to run (5-10 minutes)
env <- c()
for (i in 1:length(raslist)) { env <- cbind(env, extract(raslist[[i]], samples,buffer=10000,fun=mean,na.rm=T))} # fills matrix env with the mean value of 10 km of buffer around each point. 
colnames(env) <- names(raslist)
buff10k_samples <- cbind(samples, env) # add environmental values to sample coordinates

# Save the Result of the Area Sampling tool. 

shapefile(buff10k_samples, './MOOA_ENV_buf10k_R.shp')

# A warning will prompt: ESRI shapefiles are strict in concern of the length of table headers. 
# If you load the MOOA_ENV_buf10k_R.shp in QGIS you will notice that headers have been abbreviated.
# It is usually a good practice to foresee this issue and use abbreviations that are set by the user. 

