##In this files, some examples on how you can process and compare the data of
##the output files when you have the netCDF output from TARGET and the GeoTIFF files ready
##In this scrip, I use "sq" for status quo simulation, and "s1" for greening scenarrio 1
##These are ideas on what you could do, but you don't need to follow this scrip

##set the directory of your files
setwd("E:/GBI/Scenarios/output") 

##these two packages a useful to explore the netCDF (.nc) output files, you only need to install them once
#install.packages("raster")
#install.packages("ncdf4")
library(raster)
library(ncdf4)

##this provides you with a list of all the files within your working directory
list.files()

##here, you can read in and explore a bit your netCDF output file
nc = nc_open("epfl_.nc") #put the directory and name of your output file
print(nc)
names(nc)
names(nc$var) #here you can see the names of the variables of your output file

## you can for example get the time series of your simulation
#install.packages("ncdf4.helpers")
library(ncdf4.helpers)
time <- nc.get.time.series(nc,"TEMP_AIR")


##now, you can read in your output data in the GeoTIFF format (this one is used for the plotting)
sq <- stack("sq.tiff") #or .tif
s1 <- stack("s1.tiff")



## If you don't have them already, install and library the following packages:
library(reshape2)
library(tidyverse)
library(scales)

## convert your data stacks into data frames
##If you open he data frames, you can see that the first two columns are
##the x and y coordinates of each grid cell (rows), and the other columns are the time steps
stack_sq <- as.data.frame(sq, xy = TRUE)
stack_s1 <- as.data.frame(s1, xy = TRUE)
stack_t <- as.data.frame(test, xy = TRUE)

##here, you can calculate the overall mean of your scenarios
q<-stack_sq[,3:241]
mean(as.matrix(q), na.rm=T)

a<-stack_s1[,3:241]
mean(as.matrix(a), na.rm=T)


##here, you can visualize the daily temperature progression over your grid for the first day of your simulation
#the numbers do correspond to the columns, thus the numbers depend on the temporal resolution of your simulation (I have 30min)
day1 <- stack_sq[,c(1:3,9,15,31,37,43,49,55,61)]
#I chose to show the development by representing the temperature every 3 hours
colnames(day1)<-c("Lat","Lon","0h00","3h00","6h00","9h00","12h00","15h00","18h00","21h00","24h00")
#this rearranges your data to allow easier plotting
day1_re <- as.data.frame(day1, xy = TRUE) %>%
  melt(id.vars = c('Lat','Lon'))
#the variable corresponds now to the 3h steps (the bands/columns you chose above)
colnames(day1_re) <- c("Lat","Lon","variable","Temp.")

#create the plot
ggplot() +
  geom_raster(data = day1_re , aes(x = Lat, y = Lon, fill = Temp.)) +
  theme_classic()+
  facet_wrap(~ variable)


##you can also plot the time series of specific cells to see differences between different
#landcover classes or to compare the same cell but in different scenarios (as done below)
cell_sq <- stack_sq[c(4127),] #chose the cell you would like to plot
cell_sq$id <- row_number(cell_sq)
cell_sq_re <- as.data.frame(cell_sq, xy = TRUE) %>%
  melt(id.vars = c("id","x","y"))
cell_sq_re$datetime <- as.POSIXct(time)# here you can use the time series from the beginning
str(cell_sq_re) #this gives you a small summary of your data structure

cell_s1 <- stack_s1[c(4127),]
cell_s1$id <- row_number(cell_s1)
cell_s1_re <- as.data.frame(cell_s1, xy = TRUE) %>%
  melt(id.vars = c("id","x","y"))
cell_s1_re$datetime <- as.POSIXct(time)

#creates the plot
ggplot()+
  geom_line(data=cell_sq_re, aes(x=datetime, y=value, color="SQ"))+
  geom_line(data=cell_s1_re, aes(x=datetime, y=value, color="S1"))+
  scale_color_manual(name = "Scenarios", values = c("SQ" = "black", "S1" = "green3"))+ #chose your colors
  theme_classic()+
  labs(y="Temperature [°C]", x = "datetime") #chose the labels of you x- and y-axis

#if the lines of your plot are overlapping, this is because the landcover class wasn't changed in your
#scenario and thus TARGET modelled the same temperature as in the status quo


#here, you can calculate the mean temperature over the whole grid for each timestep
mean_sq <- as.data.frame(apply(stack_sq, 2, mean, na.rm=T))
mean_sq <- as.data.frame(mean_sq[-c(1,2),])
mean_sq$datetime <- as.POSIXct(time)
colnames(mean_sq) <- c("Ta","datetime")

mean_s1 <- as.data.frame(apply(stack_s1, 2, mean, na.rm=T))
mean_s1 <- as.data.frame(mean_s1[-c(1,2),])
mean_s1$datetime <- as.POSIXct(time)
colnames(mean_s1) <- c("Ta","datetime")

##you can also include the measured temperature to compare your modeled data with it
# Ta<- read.csv("temp_laus.txt", header = T, skip=1, sep = ";")
# Ta$datetime <- ymd_hm(Ta$time)
# Ta$Ta <- as.numeric(Ta$tre200s0)
# Ta_30 = Ta[seq(1, nrow(Ta), 3), ]
# Ta = Ta[6913:7633,4:5]
# Ta_30 = Ta[seq(1, nrow(Ta), 3), ]


ggplot()+
  geom_line(data=mean_sq, aes(x=datetime, y=Ta, color="SQ"))+
  geom_line(data=mean_s1, aes(x=datetime, y=Ta, color="S1"))+
  # geom_line(data=Ta, aes(x=datetime, y=Ta,color="Measured"))+ #measured temperature
  scale_color_manual(name = "Scenarios", values = c("Measured" = "black","SQ"="red3", "S1" = "green3"))+
  theme_classic()+
  scale_x_datetime(breaks = "1 day", date_labels = ("%d %b"))+ #adjust dateformat on x-axis
  labs(y=expression(Delta*"T [°C]"), x = "Date")


##Or you can plot the difference of the mean air temperature (delta T) between your status quo and scenario
mean_s1$dif <- mean_sq$Ta - mean_s1$Ta

ggplot()+
  geom_line(data=mean_s1, aes(x=datetime, y=dif, color="S1"))+
  scale_color_manual(name = "Scenarios", values = c("S1" = "green3"))+
  theme_classic()+
  geom_abline(slope=0,intercept = 0, color="black",linetype="dashed")+ #draw a line
  scale_x_datetime(breaks = "1 day", date_labels = ("%d %b"))+
  labs(y=expression(Delta*"T [°C]"), x = "Date")


##Another analysis you could do is to calculate the mean delta t
#of each grid cell over the whole simulation period
mw_sq <- as.data.frame(rowMeans(stack_sq[,3:241]))
mw_sq$x <- stack_sq$x
mw_sq$y <- stack_sq$y
colnames(mw_sq) <- c("Temp.","x","y")

mw_s1 <- as.data.frame(rowMeans(stack_s1[,3:241]))
mw_s1$x <- stack_s1$x
mw_s1$y <- stack_s1$y
colnames(mw_s1) <- c("Temp.","x","y")
mw_s1$dif <- (mw_sq$Temp.- mw_s1$Temp.) #delta T


##This packages helps to create nice color gradients for your plots
#install.packages("RColorBrewer")
library(RColorBrewer)
display.brewer.all() #check these out!

ggplot() +
  geom_raster(data = mw_sq , aes(x = x, y = y, fill = Temp.)) +
  labs(y="Northing [m]", x = "Easting [m]")+
  scale_fill_gradientn(name = "T [°C]",colors = rev(brewer.pal(10,"RdYlBu")), na.value = "white")+
  theme_classic()

ggplot() +
  geom_raster(data = mw_s1 , aes(x = x, y = y, fill = dif)) +
  labs(y="Northing [m]", x = "Easting [m]")+
  scale_fill_gradientn(name = expression(Delta*"T [°C]"), colors = brewer.pal(10,"RdBu"), na.value = "white",limits=c(-2.1,2.1))+
  theme_classic()


##there are a lot of other ways on how to analyse the data, for example only during the day
# or only at night. Or you only want to look at the grid cells that were affected by your scenarios