library(devtools)
## Loading required package: usethis
library(cwatmRutils)
library(raster)
## Loading required package: sp
## The legacy packages maptools, rgdal, and rgeos, underpinning this package
## will retire shortly. Please refer to R-spatial evolution reports on
## https://r-spatial.org/r/2023/05/15/evolution4.html for details.
## This package is now running under evolution status 0
library(ggplot2)
library(ncdf4)
library(lubridate)
## 
## Attaching package: 'lubridate'
## The following objects are masked from 'package:raster':
## 
##     intersect, union
## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union
library(mondate)
## 
## Attaching package: 'mondate'
## The following objects are masked from 'package:lubridate':
## 
##     as.difftime, day, month, quarter, year, ymd
## The following object is masked from 'package:base':
## 
##     as.difftime
point<- rbind(as.numeric(c('35.378', '-14.797')))
point<- as.data.frame(point)
colnames(point)<- c('x','y')
##############################################################

lake_outflow_lake_a_factor_3<- ncdf2raster(pth= "C:/Users/hinton/Documents/GitHub/MalawiLake/outputs/modflowrun6/lakeResOutflowM_daily.nc", transpose = TRUE, time = c(1, 12*3 ), spatial=point)
lake_outflow_lake_a_factor_3$date <- as.Date(lake_outflow_lake_a_factor_3$time, origin = "1901-01-01")

start<- mondate("1971-01-01")
end<- mondate("1973-12-31")
lake_outflow_lake_a_factor_3$date<- as.Date(mondate(start + seq(0, as.numeric(end - start), 1)))

ggplot(lake_outflow_lake_a_factor_3, aes(date, value, colour=var))+geom_line()

##############################################################

lake_outflow_lake_a_factor_1<- ncdf2raster(pth= "C:/Users/hinton/Documents/GitHub/MalawiLake/outputs/modflowrun7/lakeResOutflowM_daily.nc", transpose = TRUE, time = c(1, 12*3 ), spatial=point)
lake_outflow_lake_a_factor_1$date <- as.Date(lake_outflow_lake_a_factor_1$time, origin = "1901-01-01")

start<- mondate("1971-01-01")
end<- mondate("1973-12-31")
lake_outflow_lake_a_factor_1$date<- as.Date(mondate(start + seq(0, as.numeric(end - start), 1)))

ggplot(lake_outflow_lake_a_factor_1, aes(date, value, colour=var))+geom_line()

##############################################################
lake_outflow_lake_a_factor_1$Lake_A_factor<- "1"
lake_outflow_lake_a_factor_3$Lake_A_factor<- "3"

lake_outflow_lake_a_factor<- rbind(lake_outflow_lake_a_factor_1,lake_outflow_lake_a_factor_3)

ggplot(lake_outflow_lake_a_factor,aes(date, value, colour=Lake_A_factor))+geom_line()

Comparing the impact on lake discharge of “lake a factor” values. The higher lake a factor resulted in higher lake discharge

#Define path root for run looking to analyse
pth <- "C:/Users/hinton/Documents/GitHub/MalawiLake/outputs/modflowrun7"

fn <- "discharge_monthavg.nc"
# 
start<- mondate("1971-01-01")
end<- mondate("1973-12-31")

#Chiromo
pt<- rbind(as.numeric(c('35.125', '-16.375')))

pt<- as.data.frame(pt)
colnames(pt)<- c('x','y')

# extract data for the point(s)
point_data_chiromo <- ncdf2raster(pth = sprintf("%s/%s", pth, fn),
                 transpose = TRUE, time = c(1, 12*3 ), spatial = pt)

point_data_chiromo$date <- as.Date(point_data_chiromo$time, origin = "1901-01-01")

point_data_chiromo<- as.numeric(point_data_chiromo$value)
point_data_chiromo<- as.data.frame(point_data_chiromo)
point_data_chiromo$group<- "CWatM chiromo station"
point_data_chiromo$date<- as.Date(mondate(start + seq(0, as.numeric(end - start), 1)))
colnames(point_data_chiromo)<- c("Value", "group", "date")

######### 
#Discharge at Liwonde station

pt_liwonde<- rbind(as.numeric(c('35.29167', '-15.04167')))

pt_liwonde<- as.data.frame(pt_liwonde)
colnames(pt_liwonde)<- c('x','y')

# extract data for the point(s)
point_extract_liwonde <- ncdf2raster(pth = sprintf("%s/%s", pth, fn),
                 transpose = TRUE, time = c(1, 12*3 ), spatial = pt_liwonde)

point_extract_liwonde$date <- as.Date(point_extract_liwonde$time, origin = "1901-01-01")

point_extract_liwonde<- as.numeric(point_extract_liwonde$value)
point_extract_liwonde<- as.data.frame(point_extract_liwonde)
point_extract_liwonde$group<- "cwatm_liwonde"
point_extract_liwonde$date<- as.Date(mondate(start + seq(0, as.numeric(end - start), 1)))
colnames(point_extract_liwonde)<- c("Value", "group", "date")
Lake_flow_discharges<- as.data.frame(cbind(lake_outflow_lake_a_factor$value, lake_outflow_lake_a_factor$Lake_A_factor, lake_outflow_lake_a_factor$date))
Lake_flow_discharges$V3<- as.Date(mondate(start + seq(0, as.numeric(end - start), 1)))


colnames(Lake_flow_discharges)<- c("Value", "group", "date")
####
discharge_summary<- rbind(Lake_flow_discharges,point_extract_liwonde, point_data_chiromo)
discharge_summary$Value<- as.numeric(discharge_summary$Value)
####
ggplot(discharge_summary,aes(date, Value, colour=group))+geom_line()

####
discharge_summary_no_chiromo<- rbind(Lake_flow_discharges,point_extract_liwonde)
discharge_summary_no_chiromo$Value<- as.numeric(discharge_summary_no_chiromo$Value)
####
ggplot(discharge_summary_no_chiromo,aes(date, Value, colour=group))+geom_line()

comparing lake discharge with simulted downstream discharge. Lake discharge is lower than seen at the downstream simulated areas. Have checked the ncdf file and values for lake discharge are consistently low (not just missing the location of the main outflow)

Schwatke, C., Dettmering, D., Bosch, W., and Seitz, F.: DAHITI - an innovative approach for estimating water level time series over inland waters using multi-mission satellite altimetry: , Hydrol. Earth Syst. Sci., 19, 4345-4364, doi:10.5194/hess-19-4345-2015, 2015

#REFERENCE LAKE DATA

lake_malawi_water_levels<- nc_open("H:/MyDocuments/Validation data/DAHITI lake levels/lake_malawi_water_level_altimetry.nc")

lake_malawi_water_levels_data<- read.csv("H:/MyDocuments/Validation data/DAHITI lake levels/lake_malawi_water_level_data_2.csv")
colnames(lake_malawi_water_levels_data)<- c("date", "water_level", "error")
lake_malawi_water_levels_data$date<- as.Date(lake_malawi_water_levels_data$date)

ggplot(lake_malawi_water_levels_data, aes(date, water_level))+geom_point() +geom_line()