######################### HOMEWORK 1 CI71F: USEFUL COMMANDS #############################
#
# Created by Pablo Mendoza on July 8 2024
#
# This script contains commands to read and analyze data from the catchment dataset
# used by Munoz-Castro et al. (2023)
#
#  References:
#
#
#
#
############################ PART 1: Define 'hard-coded' parameters #####################
if(!require(airGR))  {install.packages("airGR"); require(airGR)}
if(!require(sfsmisc)){install.packages("sfsmisc"); require(sfsmisc)}

##### Define directories
UsernamePC    <- "pablo"
ForcingsDir   <- paste0('/Users/',UsernamePC, # Directory with meteorological forcings
                       '/Dropbox/Docencia/CI71F/GR_models/01_Forcings')
HypsoDir      <- paste0('/Users/',UsernamePC, # Directory with hypsometric curves
                       '/Dropbox/Docencia/CI71F/GR_models/02_Hypsometry')
MetadataDir   <- paste0('/Users/',UsernamePC, # Directory with metadata
                       '/Dropbox/Docencia/CI71F/GR_models/03_Catchment_Attributes')

#####  Define BNA basin code
BasinID     <- 6027001  # Note that you can change this!

#####  Load data
BasinObsData <- read.csv(paste0(ForcingsDir,'/',BasinID,'.csv')) # Hydromet time series
BasinHypso  <- read.csv(paste0(HypsoDir,'/',BasinID,'.csv')) # Hypsometric curve
BasinObsData$DatesR <- as.POSIXlt(BasinObsData$DatesR,tz="UTC")
BasinMetadata   <- read.csv(paste0(MetadataDir,'/Catchment_attributes.csv'),
                            fileEncoding = "Latin1", check.names = F) # Metadata
BasinIndex  <- which(BasinMetadata[,1]==BasinID)
BasinName   <- as.vector(BasinMetadata[BasinIndex,2])

#####  Define months
Months      <- c('Jan','Feb','Mar','Apr','May','Jun',
                'Jul','Aug','Sep','Oct','Nov','Dec')
Months2     <- c('J','F','M','A','M','J','J','A','S','O','N','D')
Nmonths     <- length(Months)

######## Declare useful functions
# Funcion para estimar Curva de duracion (FDC)
FDC <- function(data){
  
  Qaux  <- data[!is.na(data)]
  Pexc  <- (1:length(Qaux))/(1+length(Qaux))
  Qfdc  <- sort(Qaux, decreasing = T)
  
  Result <- list(Pexc = Pexc, Qfdc = Qfdc)
  
  return(Result)
  
}




############################# PART 2: Explore observed data #############################

#### Time series with hydrometeorological variables
par(mfrow=c(2,1)) # plot with 2 rows and 1 column (i.e., 2 panels)
### Daily precipitation
plot(BasinObsData$DatesR, BasinObsData$P, xlab = "Time", ylab = "Precipitation (mm/d)",
     type = "l", col = "blue", las = 1)
mtext(BasinName,      side = 3, line = 2, cex = 1.1)
### Daily runoff
plot(BasinObsData$DatesR, BasinObsData$Qmm, xlab = "Time", ylab = "Runoff (mm/d)",
     type = "l", col = "red", las = 1)
### Daily air temperature
plot(BasinObsData$DatesR, BasinObsData$T, xlab = "Time", ylab = "Temperature (C)",
     type = "l", col = "purple", las = 1)
mtext(BasinName,      side = 3, line = 2, cex = 1.1)
### Daily potential evapotranspiration
plot(BasinObsData$DatesR, BasinObsData$E, xlab = "Time", ylab = "PET (mm/d)",
     type = "l", col = "darkgreen", las = 1)

### Monthly cycles
# First, compute monthly amounts (P,PET,R) or averages (T)
Tyears    <- as.numeric(format(BasinObsData$DatesR,"%Y"))   # Time series with years only 
Tmonths   <- as.numeric(format(BasinObsData$DatesR,"%m"))   # Time series with months only
Tdays     <- as.numeric(format(BasinObsData$DatesR,"%d"))   # Time series with days only

# Define variables for number of years and number of months
Nyears    <- length(unique(Tyears)) 

# Initialize matrices with monthly values
MTemp         <- array(NA, dim=c(Nyears,Nmonths))  # Temperature
MPrecip       <- array(NA, dim=c(Nyears,Nmonths))  # Precipitation
MPET          <- array(NA, dim=c(Nyears,Nmonths))  # Potential ET
MRobs         <- array(NA, dim=c(Nyears,Nmonths))  # Monthly runoff

# Before computing averages, examine the number of days with missing flow data
# If the data is complete, outputs should be zero!
print(paste0('The number of days with NAs is ', length(which(is.na(BasinObsData$Qmm)))))
print(paste0('The number of days with NaNs is ', length(which(is.nan(BasinObsData$Qmm))))) 

# Compute monthly values per year
for (iyear in 1:Nyears){ # Start loop over years
  for (imonth in 1:Nmonths) { # Start loop over months
    mindex = intersect(which(Tyears == unique(Tyears)[iyear]),
                       which(Tmonths == imonth))
    MTemp[iyear,imonth]     <- mean(BasinObsData$T[mindex])  # Mean monthly temperature
    MPrecip[iyear,imonth]   <- sum(BasinObsData$P[mindex])   # Montly precipitation
    MPET[iyear,imonth]      <- sum(BasinObsData$E[mindex])   # Monthly PET
    MRobs[iyear,imonth]     <- sum(BasinObsData$Qmm[mindex]) # Monthly Runoff
  
  } # End loop over months
} # End loop over years

# Compute mean monthly values for the entire period
Monthly_T     <- apply(MTemp,    2, mean)
Monthly_P     <- apply(MPrecip,  2, mean)
Monthly_PET   <- apply(MPET,     2, mean)
Monthly_Robs  <- apply(MRobs,    2, mean, na.rm=T)

# Now plot the results!
# Analyze the annual cycles of the variables
par(mfrow=c(1,1), mar=c(5,4,4,5)) # plot with 2 rows and 1 column (i.e., 2 panels)
plot(1:Nmonths,Monthly_P,type="l",xlab="",ylab="Variable (mm)",
     ylim=range(Monthly_P,Monthly_PET,Monthly_Robs),las=1,col="blue",xaxt="n",lwd=2)
lines(1:Nmonths,Monthly_PET,xlab="",ylab="Variable (mm)",
     las=1,col="red",xaxt="n",lwd=2)
lines(1:Nmonths,Monthly_Robs,xlab="",ylab="Variable (mm)",
     las=1,col="green",xaxt="n",lwd=2)
axis(1,at=1:Nmonths,Months2)
legend("topright",c("P","PET","Q","T"),horiz="TRUE",
       lty=c(rep(1,3),2),lwd=rep(2,4),col=c("blue","red","green","purple"))

par(new=TRUE)
plot(1:Nmonths,Monthly_T,type="l",xlab="",ylab="",axes=FALSE,
     ylim=range(Monthly_T),col="purple",xaxt="n",lwd=2,lty=2)
axis(4,las=1)#, ylim=range(Monthly_T),col="purple",las=1)  ## las=1 makes horizontal labels
text(par("usr")[2]*1.15,mean(par("usr")[3:4])+2,"Temperature (C)",srt = -90, xpd = TRUE, pos = 4)
box()


######## Compute and plot a flow duration curve
par(mfrow=c(1,1))
QFDC     <- FDC(BasinObsData$Qmm)
plot(QFDC$Pexc, QFDC$Qfdc, type="l",log="y",yaxt = "n", xlab="Exceedance probability",
     ylab="Runoff (mm/d)",las = 1)
eaxis(2, n.axp =1)




############################# PART 3: Run model simulations #############################

## Create input forcing files
InputsModelBasGR4J <- CreateInputsModel(FUN_MOD = RunModel_GR4J, 
                                         DatesR = BasinObsData$DatesR,
                                         Precip = BasinObsData$P,
                                         PotEvap = BasinObsData$E)

## Define simulation period
Ind_Run <- seq(which(format(BasinObsData$DatesR, format = "%d/%m/%Y %H:%M")=="01/01/1990 00:00"), 
               which(format(BasinObsData$DatesR, format = "%d/%m/%Y %H:%M")=="31/12/1999 00:00"))

## preparation of the RunOptions object for both basins, GR4J model
RunOptionsGR4JBas <- CreateRunOptions(FUN_MOD       = RunModel_GR4J, # Specify model
                                      InputsModel   = InputsModelBasGR4J, # Forcing dataset
                                      IndPeriod_Run = Ind_Run)   # Simulation period

## Run the model GR4J with a pre-defined set of parameter values
Param <- c(257.238, 1.012, 88.235, 2.208)
OutputsGR4JBas  <- RunModel_GR4J(InputsModel = InputsModelBasGR4J,
                                  RunOptions = RunOptionsGR4JBas,
                                  Param = Param) # Run for basin 1

str(OutputsGR4JBas)

# Scatter plot of simulated versus observed flows
maxQ = max(c(BasinObsData$Qmm[Ind_Run],
             OutputsGR4JBas$Qsim),na.rm=TRUE)
plot(BasinObsData$Qmm[Ind_Run],OutputsGR4JBas$Qsim,las = 1,
     xlab="Qobs (mm/d)", ylab="Qsim (mm/d)",xlim=c(0,maxQ),ylim=c(0,maxQ))
lines(c(0,maxQ),c(0,maxQ))


######## Include CemaNeige in simulations!
ParamSnow <- c(257.238, 1.012, 88.235, 2.208, 0.962, 2.249) # Vector with parameters
InputsModelBasGR4JCemma <- CreateInputsModel(FUN_MOD = RunModel_CemaNeigeGR4J,
                                 DatesR = BasinObsData$DatesR,
                                 Precip = BasinObsData$P,
                                 PotEvap = BasinObsData$E,
                                 TempMean = BasinObsData$T,
                                 ZInputs = median(BasinHypso$HypsoData),
                                 HypsoData = BasinHypso$HypsoData, NLayers = 5)
RunOptionsGR4JBasCemma <- CreateRunOptions(FUN_MOD = RunModel_CemaNeigeGR4J,
                                     InputsModel = InputsModelBasGR4JCemma,
                                     IndPeriod_Run = Ind_Run)
OutputsGR4JBasCemma  <- RunModel_CemaNeigeGR4J(InputsModel = InputsModelBasGR4JCemma,
                                                  RunOptions = RunOptionsGR4JBasCemma,
                                                  Param = ParamSnow)
str(OutputsGR4JBasCemma)

# Scatter plot of simulated versus observed flows
maxQ = max(c(BasinObsData$Qmm[Ind_Run],
             OutputsGR4JBasCemma$Qsim),na.rm=TRUE)
plot(BasinObsData$Qmm[Ind_Run],OutputsGR4JBasCemma$Qsim,las = 1,
     xlab="Qobs (mm/d)", ylab="Qsim (mm/d)",xlim=c(0,maxQ),ylim=c(0,maxQ))
lines(c(0,maxQ),c(0,maxQ))




##################################### END OF SCRIPT #####################################

