library(unikn)         #For fancy colormaps
library(raster)        # For raster operations
library(RColorBrewer)  #Another library for fancy color palettes 

# Import data from the downloaded folder
sm=raster("./SampleData-master/raster_files/SMAP_SM.tif") # SMAP soil moisture data

# Make custom color palette
mypal=unikn::seecol(pal = usecol(c("brown","orange","yellow","darkgreen","cyan","darkblue","black"), n = 20)) #User defined color palette using unikn library

#Alternatively, use "brewer.pal" for color palettes 
mypal=brewer.pal(11,"Spectral") #We can use color palettes from RColorBrewer

# Arrange the plots in user defined grid layout
layout(matrix(c(1,2,3,4,5,6), 2,3, byrow = TRUE)) #Use layout to arrange the plots

#~~ Projection 1: NAD83 (EPSG: 4269)
sm_proj1 = projectRaster(sm, crs = CRS("+init=epsg:4267"));
plot(sm_proj1, 
     main="NAD83",          # Title of the plot
     col=mypal,             # Colormap for the plot
     axes=FALSE,            # Disable axes
     box=FALSE,             # Disable box around the plots
     asp = NA,              # No fixed aspect ratio
     legend=FALSE,          # Disable legend
     interpolate=TRUE)      # Interpolate output raster   

#~~ Projection 2: Mercator (EPSG: 3857): Used in Google Maps, Open Street Maps, etc.
sm_proj2 = projectRaster(sm, crs = CRS("+init=epsg:3857"));
plot(sm_proj2, main="Mercator", col=mypal, axes=FALSE, 
     box=FALSE,asp = NA,legend=FALSE,interpolate=TRUE)

#~~ Projection 3: World Robinson Projection system
test_proj="+proj=robin +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"
sm_proj3 = projectRaster(sm, crs = test_proj);
plot(sm_proj3, main="Robinson", col=mypal, axes=FALSE, 
     box=FALSE,asp = NA,legend=FALSE,interpolate=TRUE)

#~~ Projection 4: Azimuthal equidistant projection centered on Lat=0 and Lon=0
test_proj="+proj=aeqd +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"
sm_proj4 = projectRaster(sm, crs = test_proj);
plot(sm_proj4, main="Azimuthal", col=mypal, axes=FALSE, 
     box=FALSE,asp = NA,legend=FALSE,interpolate=TRUE)

#~~ Projection 5:  World Mercator
test_proj="+proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"
sm_proj5 = projectRaster(sm, crs = test_proj);
plot(sm_proj5, main="World Mercator", col=mypal, axes=FALSE, 
     box=FALSE,asp = NA,legend=FALSE,interpolate=TRUE)

library(ggplot2)

#~~~ 1.2.1. Transformation of raster to x-y-z format. 
sm_df_smap = as.data.frame(as(sm, "SpatialPixelsDataFrame")) # Spatial dataframe of SM data
colnames(sm_df_smap) = c("SM", "Long", "Lat")        # Rename the columns of the dataframe
head(sm_df_smap)                                     # Checkout the spatial dataframe

sm_world=ggplot() +                             # Initialize the plot
  geom_tile(data=sm_df_smap,                    # Add soil moisture data
            aes(x=Long, y=Lat, fill=SM))+       # Provide X, Y and Z data
   scale_fill_gradientn(colors=mypal,           # Use user-defined colormap
   name = "SM",                                 # Name of the colorbar
   na.value = "transparent",                    # transparent NA cells
   labels=(c("0", "0.2", "0.4", "0.6", "0.8")), # Labels of colorbar
   breaks=seq(0,0.8,by=0.2),                    # Set breaks of colorbar
   limits=c(0,0.8))+                            # Set limits of colorbar
  coord_equal()+                                # Aspect ratio of plot
  xlab("Longitude")+                            # X-axis label
  ylab("Latitude")+                             # Y-axis label
  borders("world",                              # Add global landmass boundaries
          colour="gray43",                      # Fill light-gray color to the landmass
          fill="transparent")+                  # Transparent background
  borders("state",                              # Add US state borders
          colour="gray43",                      # Use light-gray color
          fill="transparent")                   # Use transparent background      
print(sm_world)

#~~~ 1.2.2. Plot a fancier re-projected maps in Robinson projection
source("./SampleData-master/functions/robin_ras.R") # User-defined function for plotting
robin_plot=robin_ras(sm,                            # Raster to plot
                     map_col=rev(mypal),            # Colormap
                     lim=c(0,0.8),                  # Set min and max z-limits
                     brk=seq(0,0.8,0.2),            # Set breaks in colormap
                     name="Global Surface Soil Moisture from SMAP")   # Name of the plot
robin_plot                          # print the reprojected plot

sm_aus=ggplot() +                               # Initialize the plot
  geom_tile(data=sm_df_smap,                    # Add soil moisture data
            aes(x=Long, y=Lat, fill=SM))+       # Provide X, Y and Z data
   scale_fill_gradientn(colors=mypal,           # Use user-defined colormap
   name = "SM",                                 # Name of the colorbar
   na.value = "transparent",                    # transparent NA cells
   labels=(c("0", "0.2", "0.4", "0.6", "0.8")), # Labels of colorbar
   breaks=seq(0,0.8,by=0.2),                    # Set breaks of colorbar
   limits=c(0,0.8))+                            # Set limits of colorbar
  coord_equal()+                                # Aspect ratio of plot
  xlab("Longitude")+                            # X-axis label
  ylab("Latitude")+                             # Y-axis label
  borders("world",                              # Add global landmass boundaries
          colour="gray43",                      # Fill light-gray color to the landmass
          fill="transparent")+                  # Transparent background
  borders("state",                              # Add US state borders
          colour="gray43",                      # Use light-gray color
          fill="transparent")+                  # Use transparent background
  coord_fixed(xlim = c(114,153),                # Add extent for Australia
              ylim = c(-43,-11))+               
  theme_bw()                                    # Try black-and-white theme. 

print(sm_aus)

###############################################################
#~~~ PART 2.1.1: Importing and visualizing shapefiles
library(sf)
library(raster)
library(rgdal)
library(fields)

# Import the shapefile of global IPCC climate reference regions (only for land) 
IPCC_shp = readOGR("./SampleData-master/CMIP_land/CMIP_land.shp",verbose = FALSE)
# View attribute table of the shapefile
head(IPCC_shp@data)

# Download global coastlines
download.file("http://www.naturalearthdata.com/http//www.naturalearthdata.com/download/10m/physical/ne_10m_coastline.zip",destfile = 'coastlines.zip')
# Unzip the downloaded file
unzip(zipfile = "coastlines.zip",exdir = 'ne-coastlines-10m')
# Load global coastline shapefile 
coastlines= readOGR("ne-coastlines-10m/ne_10m_coastline.shp",verbose = FALSE)

# Plot a raster with overlaying shapefiles
plot(sm, col=brewer.pal(10,"Spectral"))  # Global SMAP soil moisture
plot(IPCC_shp,add=TRUE)                  # Add IPCC land regions in blue color
plot(coastlines, add=TRUE)               # Add global coastline

#Subset shapefile for Eastern North-America
plot(IPCC_shp[4,], main="Polygon for Eastern North-America")

###############################################################
#~~~ PART 2.2.1: Using shapefile to summarize a raster
sm_IPCC_df=raster::extract(sm,       # Raster to be summarized
                          IPCC_shp,  # Shapefile/ polygon to summarize the raster
                          df=TRUE,   # Gives the summary statistics as a dataframe
                          fun=mean,  # Desired statistic
                          na.rm=TRUE)# Ignore NA values? yes! 

head(sm_IPCC_df)

###############################################################
#~~~ PART 2.2.2: Convert a raster to a shapefile
aridity=raster("./SampleData-master/raster_files/aridity_36km.tif") #Global aridity
plot(aridity, 
     col=mypal, 
     interpolate=TRUE, 
     lab.breaks=c(0:5), 
     breaks=c(0:5),
     main="Global aridity map [1=driest, 5= wettest]") 

# Convert raster to shapefile
arid_poly=raster::rasterToPolygons(aridity, # Input raster
                           na.rm=TRUE,      # Ignore NA values
                           digits=12,       # No. of digits to round the coordinates
                           dissolve=TRUE)   # FALSE will give boundaries for each cell

plot(arid_poly, col=arid_poly$aridity_36km) # Plot aridity polygon

sm_arid_df=raster::extract(sm,       # Raster to be summarized
                          arid_poly, # Shapefile/ polygon to summarize the raster
                          df=TRUE,   # Gives the summary statistics as a dataframe
                          fun=mean,  # Desired statistic
                          na.rm=TRUE)# Ignore NA values? yes! 

# Lets plot the climate-wise mean of surface soil moisture
plot(sm_arid_df,xaxt = "n",
     xlab="Aridity", ylab="Soil moisture",
     type="b", 
     col="blue", 
     main="Climate-wise mean of surface soil moisture")
axis(1, at=1:5, labels=c("Hyper-arid", "Arid", "Semi-Arid","Sub-humid","Humid"))

###############################################################
#~~~ PART 2.3.1: Converting shapefile to raster
samp_ras = raster(ncol=1000,nrow=1000)     # Create an empty sample raster
extent(samp_ras) = extent(IPCC_shp)        # Set raster extent to the shapefile 

# Convert polygon to raster
IPCC_ras=raster::rasterize(IPCC_shp,          # Input shapefile
                    samp_ras,                 # Sample raster
                    as.numeric(IPCC_shp$V4))  # Rasterize using field ("V4")

###############################################################
#~~~ PART 2.3.2: Resampling for common extent and resolution
IPCC_ras_resamp=resample(IPCC_ras,sm,method='bilinear') # Resample using bilinear or ngb (nearest neighbor) interpolation
print(IPCC_ras_resamp) #Notice resolution, extent and coordiate. ref.

plot(IPCC_ras_resamp, 
     asp=NA,      # Fits plot to window
     main="Raster of IPCC6 land regions")

# Convert resampled raster to matrix
sp_mat=as.matrix(IPCC_ras_resamp)

# Notice the transpose and inversion used in sp_mat
fields::image.plot(t(sp_mat[nrow(sp_mat):1,]),main="IPCC land regions") 

###############################################################
#~~~ PART 3.1.1: Create and plot NDVI RasterBrick/ Stack 
library(raster)
library(fields)
library(unikn)

# Location of the NDVI raster files
ndvi_path="./SampleData-master/raster_files/NDVI/" #Specify location of NDVI rasters

# List of all NDVI rasters
ras_list=list.files(ndvi_path,pattern='*.tif',full.names=TRUE)
head(ras_list)

# Use lapply to create raster layers from the raster location
ras_brick=brick(lapply(paste(ras_list, sep=""), raster)) 
dim(ras_brick) #23 raster layers with 456x964 cells

# Subsetting raster stack/brick (notice the double [[]] bracket and similarity to lists)
sub_ras_brick=ras_brick[[c(1,3,5,10,12)]] #Select 1st, 3rd, 5th, 10th and 12th layers

#Plot first 4 elements of NDVI raster brick
plot(sub_ras_brick[[1:4]],         # Select raster layers to plot
     col=brewer.pal(10,"Spectral"),# Specify colormap
     interpolate=TRUE,             # Interpolate the output map
     asp=NA)                       # Aspect ratio= NA

mapTheme = rasterTheme(region = brewer.pal(10,"Spectral")) # Set raster theme (color etc.)
my.at=seq(-1,1,0.1)                               # Set colorkey/ colorbar sequence
myColorkey = list(at=my.at, space="right")        # Set colorkey/colorbar to the right

plt= rasterVis::levelplot(sub_ras_brick[[1:4]],   # Select first 4 raster layers to plot
               margin = F,          
               pretty=TRUE,             
               layout=c(2,2),           # Layout of the plots
               row.values=2,            # Number of rows
               par.settings = mapTheme, # Add user define theme
               colorkey=myColorkey,     # Add user defined 
               main="NDVI")
plt=plt + layer(sp.lines(coastlines, col="black", lwd=0.1)) # Add coastline to the plots
print(plt)

###############################################################
#~~~ PART 3.2.1: Layer-wise cell-statistics 

# Layer-wise Mean
cellStats(sub_ras_brick, mean) # Mean of each raster layer. Try modal, median, min etc. 

# Layer-wise quantiles
quantile(sub_ras_brick,probs=c(.25,.75))

# User-defined statistics by defining own function
quant_fun = function(x, na.rm=TRUE) # Remember to add "na.rm" option
                    {quantile(x, probs = c(0.25, 0.75), na.rm=TRUE)} 
cellStats(sub_ras_brick, quant_fun) # 25th, and 75th percentile of each layer

###############################################################
#~~~ PART 3.2.2: Global cell-statistics

min_val=min(raster::minValue(sub_ras_brick))
max_val=max(raster::maxValue(sub_ras_brick))
print(c(min_val,max_val))

# Or by vectorizing the RasterBrick and finding statistics
min_val=min(as.vector(sub_ras_brick), na.rm=T) 
max_val=max(as.vector(sub_ras_brick), na.rm=T)
print(c(min_val,max_val))

# Another example of statistics of vectorized RasterBrick
quant=quantile(as.vector(sub_ras_brick), c(0.01,0.99), na.rm=TRUE) #1st and 99th percentiles
print(quant)

###############################################################
#~~~ PART 3.2.3: Layer-wise operations on RasterStack/RasterBrick

# Arithmetic operations on RasterBrick/ RatserStacks are same as lists
add=sub_ras_brick+10                  # Add a number to raster layers
mult=sub_ras_brick*5                  # Multiply a number to raster layers
subset_mult=sub_ras_brick[[1:3]]*10   # Multiply a number to a subset of raster layers

# Data filtering based on cell-value
filter_brick=sub_ras_brick            # Create a RasterBrick to operate on
filter_brick[filter_brick<0.5]=NA     # Assign NA to any value less than 0.5

# Let's plot the filtered rasters
plot(filter_brick, 
     asp=NA,                         # Aspect ratio. NA= No white space
     colNA="gray95",                 # Background color
     nr=2)                           # Number of rows to arrange the plots

# Layer-wise Log-transformation
log_ras_brick=log(sub_ras_brick) 

# Normalize raster layers to [0,1] based on min and max of raster brick/stack
norm_brk=(sub_ras_brick-min_val)/(max_val-min_val) # Notice that the values are between [0,1]

# Plot in Robinson projection
source("./SampleData-master/functions/robin_ras.R")  # User defined function
ndvi_robin_plot=robin_ras(norm_brk[[1]],             # Select first raster
                     map_col=rev(mypal),             # Colormap
                     lim=c(0,1),                     # Set min and max z-limits
                     brk=c(0, 0.2, 0.4, 0.6, 0.8, 1),# Set breaks in colormap
                     name="Normalized NDVI")         # Name of the plot
ndvi_robin_plot # print the reprojected plot

###############################################################
#~~~ PART 3.2.4: Cell-wise operations on raster stack/brick

# Mean of all cells accross layers
mean_ras=mean(sub_ras_brick, na.rm=TRUE)  #Try modal, median, min etc.

# Approximate NA values by interpolating between Layers
approx_brick=approxNA(sub_ras_brick,method="linear") #Try "constant"

# Using calc function for user defined operations 
quant_fun = function(x, na.rm=TRUE) # Remember to add "na.rm" option
                    {quantile(x, probs = c(0.25), na.rm=TRUE)} 

# Apply the function "quant_fun" on rasterBrick using "calc" and "stackApply" 
quant_ras=calc(sub_ras_brick, 
               fun = quant_fun)
quant_ras=raster::stackApply(sub_ras_brick, 
                             fun=quant_fun,
                             indices=c(1,1,1,1,1)) 

# Plot output raster
plot(quant_ras, main="25th percentile of time series for each cell")

# "Indices" specifies the group/subset of layers on which to apply the "fun"
quant_ras = raster::stackApply(sub_ras_brick, 
                             fun=quant_fun,
                             indices=c(1,1,2,2,2)) 
names(quant_ras) = c("Output of rasters 1 to 2", "Output of rasters 3 to 5")
# Two layers are formed, one for each group of indices

# Lets plot the two output rasters
plot(quant_ras[[1:2]], # Layers to plot
     asp=NA,           # Aspect ratio. NA= No white space
     colNA="gray95",   # Background color
     nr=2)             # Number of rows to arrange the plots

###############################################################
#~~~~ "NetCDF" file
library(ncdf4) 
file_name="./SampleData-master/SMOS_nc/SM_RE06_MIR_CDF3SA_20181225.nc"
nc= nc_open(file_name)
print(nc) # Notice the structure and attributes of the file

# Extract variables from NC file
sm=ncvar_get(nc,"Soil_Moisture")  # Soil moisture variable
Latitude=ncvar_get(nc, 'lat')     # Latitude
Longitude=ncvar_get(nc, 'lon')    # Longitude
flag=ncvar_get(nc,"Quality_Flag") # 0: data OK, 1: data not recommended
nc_close(nc)  #Close .nc file. Imp if multiple files are opened recursively

# Filter dataset
sm[!(flag==0)]=NA # Filter-out flagged data. 
sm[sm<0]=NA       # Negative values are floored at 0

# Convert the matrix to raster using lat-long
source("./SampleData-master/functions/mat2ras.R")
library(reshape2)
sm_ras=mat2ras(sm, Longitude,Latitude, FALSE)

plot(sm_ras, main="Swath of global SM retrieval (2018-12-25)")
plot(coastlines,add=TRUE)

###############################################################
#~~~~ ".H5" files 
#list all SMAP-L4 .H5 files
smapl4_files=list.files("./SampleData-master/SMAPL4_H5", full.names = TRUE) 
# Open the first file as an example
nc= nc_open(smapl4_files[1]) 
print(nc) # Notice the variables and the structure of the file

# Extract fields/ variables
sm=ncvar_get(nc,"Geophysical_Data/sm_surface")  # Surface soil moisture
Latitude=ncvar_get(nc, 'cell_lat')     # Latitude
Longitude=ncvar_get(nc, 'cell_lon')    # Longitude
nc_close(nc)                          # Close connection to .nc file

print(Longitude[1:5,1:5]) #Notice how the values are repeated in each row

print(Latitude[1:5,1:5]) #Notice how the values are repeated in each column

# Subset the Lat-long fields to retrieve unique values
Latitude=Latitude[1,]
Longitude=Longitude[,1]

# Conver data to raster using mat2ras.R function
source("./SampleData-master/functions/mat2ras.R")
library(raster)
library(reshape2)
library(RColorBrewer)

sm_ras=mat2ras(sm, Longitude,Latitude, FALSE)
plot(sm_ras, main="SMAP-L4 Surface SM retrieval", col=brewer.pal(11, "Spectral"))
plot(coastlines,add=TRUE)

source("./SampleData-master/functions/mat2ras.R")
library(reshape2) # Required for mat2ras.R

#list all SMAP-L4 .H5 files in the folder 
smapl4_files=list.files("./SampleData-master/SMAPL4_H5", full.names = TRUE) 
# List of ".H5" files 
print(smapl4_files)  

# Create and empty brick to save rasters
collate_ras=brick() 
# Create an empty list to store SM matrix
collate_mat=list() 

for (file_num in 1:length(smapl4_files)){
#print(file_num) #Print loop count
  
# Extract data fields
nc= nc_open(smapl4_files[file_num]) 
# Extract fields/ variables
sm=ncvar_get(nc,"Geophysical_Data/sm_surface")  # Soil moisture
Latitude=ncvar_get(nc, 'cell_lat')              # Latitude
Longitude=ncvar_get(nc, 'cell_lon')             # Longitude
nc_close(nc)                        # Close connection to .nc file

# Subset the Lat-long fields to retrieve unique values
Latitude=Latitude[1,]
Longitude=Longitude[,1]

# Convert soil moisture data matrix to raster
sm_ras=mat2ras(sm, Longitude,Latitude, # Inputs
               plot_ind=FALSE)         # Disable output plot

# Store layer data to RaterBrick or list
collate_mat[[file_num]]=sm
collate_ras=addLayer(collate_ras,sm_ras)
}

# Plot using RasterBrick
library(RColorBrewer)
collate_ras[collate_ras<0]=NA   #Flooring negative values to zero
plot(collate_ras,
     col=brewer.pal(11, "Spectral"))

# Plot first element of list (2D matrix)
library(fields)
collate_mat[[1]][collate_mat[[1]]<0]=NA  #Flooring negative values to zero
image.plot(Longitude,
           rev(Latitude),
           collate_mat[[1]][,ncol(collate_mat[[1]]):1],
           col=brewer.pal(11, "Spectral"), xlab="Longitude",ylab="Latitude")

# ~~~ Un-comment and run the following 5 lines if the package is not available on CRAN

# if (!requireNamespace("BiocManager", quietly = TRUE))
# install.packages("BiocManager")
# library(BiocManager)
# BiocManager::install(version = "3.11")
# BiocManager::install(c("rhdf5"))

#~~~~ Using rhdf5 for ".h5" files in R
# Import necessary functions and load libraries
source("./SampleData-master/functions/mat2ras.R")
library(reshape2) # Required for mat2ras.R
library(rhdf5)
library(raster)

smapl4_files=list.files("./SampleData-master/SMAPL4_H5", full.names = TRUE)
#View attributes of ".H5" files using h5ls function
h5ls(smapl4_files[1])

# Extract data fields
sm=h5read(smapl4_files[1],"Geophysical_Data/sm_surface")
Longitude = h5read(smapl4_files[1],"/cell_lon")
Latitude = h5read(smapl4_files[1],"/cell_lat")
# Subset the Lat-long fields to retrieve unique values
Latitude=Latitude[1,]
Longitude=Longitude[,1]
# Conver data to raster
sm_ras=mat2ras(sm, Longitude,Latitude, FALSE)
plot(sm_ras, main="SMAP-L4 Surface SM retrieval")
plot(coastlines,add=TRUE)