################################################################################################################################################
################################################ Introduction to R raster
################################################################################################################################################
# load required library
library(raster)
library(rgdal)

### Handling a raster

# Here below we will see how to open, visualize and modify a raster layer in R.

ras = raster('MAR_cov/MAR_cov.grd') # open Landcover layer

# Raster layers are imported as the object type RasterLayer. 
ras # Shows the different properties of a raster layer object

res(ras) # shows the raster resolution 
extent(ras) # Shows the raster extent
projection(ras) #Shows the raster layer projection system

plot(ras) # visualize the raster file as a plot
ras[] # visualize the raster file as a vector (each number is the value of a pixel)

plot(ras==22) # visualize only pixel with value equal to 22 (i.e. artifical areas)

# Note that raster with same extent and resolution can be added, subtracted, multiplied, divided
plot((ras==22)+(ras==16)) # Visualize the areas with landcover classes 22 or 16. 

### Stacks

## Now we will use again the Cloud Fraction February means. We will calculate the standard deviation from 2001 to 2016.

filenames <- list.files('Cloud Fraction/', pattern="*.tif$", full.names=TRUE) # list all the .tif in the Cloud Fraction folder

CLD = stack(filenames) # Load all the raster in R as a stack of raster
CLD # shows the properties of stack object. 

# A stack is an object that contains several raster with same extent, resolution and projection.

CLD$MODAL2_M_CLD_OT_2000.08 # You can access any raster included in the stack by searching it by name after the $ separator...
CLD[[1]] # ... or by indexing it as a list.

CLD[1,1] # indexing as a matrix shows the value for a given pixel (in this case, the top-left pixel in the map) for all the raster in the stack (i.e. the 17 years)

# The goal of a stack is to perform fast calculations on different snapshots of the same area. 

plot(mean(CLD)) # Applying the mean function to the stack calculates the mean for the 17 years contained in it, as we did in QGIS.

CLDsd <- calc(CLD,sd) # The calc function allows to apply any function to a stack. Here we compute the standard deviation. 

plot(CLDsd)

# Finally we can write a raster output so that we can visualize the new layer on QGIS.

writeRaster(CLDsd, 'CLDsd.tif', format='GTiff')




