---
title: "ESSD2021_Klingleretal_reviewcheck"
output:
  pdf_document: default
  html_document:
    df_print: paged
  word_document: default
---
## Purpose
This document shows basic quality checks applied to LamaH | Large-Sample Data for Hydrology and Environmental
Sciences for Central Europe, submitted by Christoph Klingler, Karsten Schulz, and Mathew Herrnegger to
ESSD in March 2021.

## Set-up
```{r setup, include=T, echo=F, include=F}
knitr::opts_chunk$set(echo = TRUE)

library(broom)
library(data.table)
library(plotly)
library(ggplot2)
library(magrittr)
library(rgdal)
library(rprojroot)
library(sp)

rootdir <- find_root(has_dir("src"))
resdir <- file.path(rootdir, "results")
datadir <- file.path(rootdir, "data")
```

## Basic checks

# Check download link
Note: directory is called "CAMELS_AT". Please correct it to "LamaH".
```{r download}
path_lamah <- file.path(datadir, 'CAMELS_AT')
path_gz <- file.path(datadir, '2_LamaH_daily.tar.gz')

if (!dir.exists(path_lamah)) {
  dlv <- download.file(url="https://zenodo.org/record/4609826/files/2_LamaH_daily.tar.gz?download=1",
                       destfile=path_gz)
  print(dlv)
  
  untar(path_gz, exdir = datadir)
}

```

# Check folder structure
```{r structure, echo=T}
filelist <- list.files(path = path_lamah, pattern = NULL, all.files = T,
                       full.names = T, recursive = T, include.dirs = T)
mapply(function(u, v) print(filelist[u:v]), 
       u=seq(1, length(filelist), 1000),
       v=c(seq(1000, length(filelist), 1000), length(filelist))
)
```

# Check number of daily record files
```{r dailycheck, echo=T}
maindirs <- c('A_basins_total_upstrm',
              'B_basins_diff_upstrm_all',
              'C_basins_diff_upstrm_lowimp',
              'D_gauges')

lapply(maindirs, function(dir) {
  print(dir)
  length(list.files(path=file.path(path_lamah, dir, "2_timeseries/daily/"),
                    pattern = 'ID.*',
                    all.files = T)
  )
})
```

# Check A_basins_total_upstrm
## Check catchment attributes
```{r Astructure, echo=T}
path_A <- file.path(path_lamah, 'A_basins_total_upstrm')


cattri <- fread(file.path(path_A, '1_Attributes', 'Catchment_attributes.csv'))
str(cattri)
summary(cattri)
nrow(cattri)
which(duplicated(cattri$ID))
```
###Check that max. sum of LCs doesn't exceed 1 beyond rounding.
```{r Alc, echo=T}
cattri[, sum(c(agr_fra, bare_fra, forest_fra, glac_fra, lake_fra, urban_fra)), 
       by=ID][, range(V1)]
```

###Check that max. sum of geology doesn't exceed 1 beyond rounding.
```{r Ageo, echo=T}
range(
  cattri[,rowSums(.SD, na.rm=T), 
         .SDcols= grep('gc.*', names(cattri), value=T)[-1]]
)
```

###Check that max. sum of texture doesn't exceed 1 -> issue here
```{r Asoil, echo=T}
cattri[, sum(sand_fra, silt_fra, clay_fra, grav_fra, oc_fra), 
       by=ID][, range(V1)]

cattri[, sum(sand_fra, silt_fra, clay_fra, grav_fra), 
       by=ID][, range(V1)]
```
In this case, there seems to be an issue. Please adjust.

## Check stream dist
```{r Adist, echo=T}
sdist <- fread(file.path(path_A, '1_Attributes','Stream_dist.csv'))
nrow(sdist)
summary(sdist)
which(duplicated(sdist$ID))
which(!sdist$ID %in% cattri$ID)
```

## Check water balance
```{r Awb, echo=T}
wb <- fread(file.path(path_A, '1_Attributes','Water_balance.csv'))
nrow(wb)
summary(wb)
which(duplicated(wb$ID))
which(!wb$ID %in% cattri$ID)

ggplot(wb, aes(x=P_1989_2009 - Q_1989_2009,
               y=ETA_1989_2009)) + 
  geom_point()
```
There are 165 NAs in Q_1989_2009. I may be missing something, but this appears to be a glitch. Check with gauges whether these are the ones which do not have data 1989-2009


## Check daily time series
Inspect 3 time series at random
```{r Adaily}
Adpath <- file.path(path_A, '2_timeseries', 'daily')

Adlist <- list.files(Adpath) %>%
  sample(3) %>%
  file.path(Adpath, .)


Adtab1 <- fread(Adlist[1])
str(Adtab1)
summary(Adtab1)

Adtab2 <- fread(Adlist[2])
str(Adtab2)
summary(Adtab2)

Adtab3 <- fread(Adlist[3])
str(Adtab3)
summary(Adtab3)

```

## Inspect shapefiles
```{r A shapefile}
UAT <- readOGR(dsn = file.path(path_A, '3_shapefiles'), 
               layer = 'Upstrm_area_total')
names(cattri) %in% names(UAT) #Check that all attributes are there

#Check CRS
UAT@proj4string

# Plot it
UATformat <- fortify(UAT)

ggplotly(
  ggplot(UATformat, aes(x = long, y = lat, group=group, fill=group)) + 
    geom_polygon(alpha=0.2, color='black') + 
    theme_bw() +
    theme(legend.position='none')
)
  

```

# Check B_basins_diff_upstrm_all
## Check catchment attributes
```{r Bstructure, echo=T}
path_B <- file.path(path_lamah, 'B_basins_diff_upstrm_all')


cattriB <- fread(file.path(path_B, '1_Attributes', 'Catchment_attributes.csv'))
nrow(cattriB)
str(cattriB)
names(cattriB) %in% names(cattri)
summary(cattriB)
```
###Check that max. sum of LCs doesn't exceed 1 beyond rounding.
```{r Blc, echo=T}
cattriB[, sum(c(agr_fra, bare_fra, forest_fra, glac_fra, lake_fra, urban_fra)), 
       by=ID][, range(V1)]
```

###Check that max. sum of geology doesn't exceed 1 beyond rounding.
```{r Bgeo, echo=T}
range(
  cattriB[,rowSums(.SD, na.rm=T), 
         .SDcols= grep('gc.*', names(cattri), value=T)[-1]]
)
```

###Check that max. sum of texture doesn't exceed 1 -> issue here
```{r Bsoil, echo=T}
cattriB[, sum(sand_fra, silt_fra, clay_fra, grav_fra, oc_fra), 
       by=ID][, range(V1)]

cattriB[, sum(sand_fra, silt_fra, clay_fra, grav_fra), 
       by=ID][, range(V1)]
```
In this case, there seems to be an issue. Please adjust.

## Check stream dist
```{r Bdist, echo=T}
sdistB <- fread(file.path(path_B, '1_Attributes','Stream_dist.csv'))
nrow(sdistB)
summary(sdistB)
```

# Check C_basins_diff_upstrm_lowimp
## Check catchment attributes
```{r Cstructure, echo=T}
path_C <- file.path(path_lamah, 'C_basins_diff_upstrm_lowimp')


cattriC <- fread(file.path(path_C, '1_Attributes', 'Catchment_attributes.csv'))
nrow(cattriC)
str(cattriC)
names(cattriC) %in% names(cattri)
summary(cattriC)
```
###Check that max. sum of LCs doesn't exceed 1 beyond rounding.
```{r Clc, echo=T}
cattriC[, sum(c(agr_fra, bare_fra, forest_fra, glac_fra, lake_fra, urban_fra)), 
       by=ID][, range(V1)]
```

###Check that max. sum of geology doesn't exceed 1 beyond rounding.
```{r Cgeo, echo=T}
range(
  cattriC[,rowSums(.SD, na.rm=T), 
         .SDcols= grep('gc.*', names(cattri), value=T)[-1]]
)
```

###Check that max. sum of texture doesn't exceed 1 -> issue here
```{r Csoil, echo=T}
cattriC[, sum(sand_fra, silt_fra, clay_fra, grav_fra, oc_fra), 
       by=ID][, range(V1)]

cattriC[, sum(sand_fra, silt_fra, clay_fra, grav_fra), 
       by=ID][, range(V1)]
```
In this case, there seems to be an issue. Please adjust.

## Check stream dist
```{r Cdist, echo=T}
sdistC <- fread(file.path(path_C, '1_Attributes','Stream_dist.csv'))
nrow(sdistC)
summary(sdistC)
```
Why are there

## Check gauge hierarchy
```{r Chier, echo=T}
ghierC <- fread(file.path(path_C, '1_Attributes','Gauge_hierarchy.csv'))
nrow(ghierC)
str(ghierC)
summary(ghierC)
```

## Check daily time series
Inspect 2 time series at random
```{r Cdaily}
Cdpath <- file.path(path_C, '2_timeseries', 'daily')


Cdlist <- list.files(Cdpath) %>%
  sample(2) %>%
  file.path(Cdpath, .)


Cdtab1 <- fread(Cdlist[1])
str(Cdtab1)
summary(Cdtab1)

Cdtab2 <- fread(Cdlist[2])
str(Cdtab2)
summary(Cdtab2)
```

## Inspect shapefiles
```{r c shapefile}
UAD <- readOGR(dsn = file.path(path_C, '3_shapefiles'), 
               layer = 'Upstrm_area_diff_lowimp')
names(cattriC) %in% names(UAD) #Check that all attributes are there

#Check CRS
UAD@proj4string

# Plot it
UADformat <- fortify(UAD)

ggplot(UADformat, aes(x = long, y = lat, group=group, fill=group)) + 
  geom_polygon(alpha=0.2, color='black') + 
  theme_bw() +
  theme(legend.position='none')

```

# Check D_gauges
## Check gauges 1_attributes
### Gauge_attributes.csv
```{r Dattri}
path_D <- file.path(path_lamah, 'D_gauges')
Dattridir<- file.path(path_D, '1_attributes')

gattri <- fread(file.path(Dattridir, 'Gauge_attributes.csv'))
str(gattri)
summary(gattri)
which(duplicated(gattri$ID))
which(!gattri$ID %in% cattri$ID)
min(gattri$obsbeg_day)

qplot(gattri$area_ratio)
length(which(abs(1-gattri$area_ratio) >= 0.1))
```

###Hydro_indices_1989_2009.csv
```{r Dindices8909}
gindices_8909 <- fread(file.path(Dattridir, 'Hydro_indices_1989_2009.csv'))
gindices_8909[, .N]
gindices_8909[!is.na(q_mean), .N]
str(gindices_8909)
summary(gindices_8909)
which(duplicated(gindices_8909$ID))
which(!gindices_8909$ID %in% cattri$ID)

```
###Hydro_indices_beg_2017.csv
```{r Dindices17}
gindices_beg17 <- fread(file.path(Dattridir, 'Hydro_indices_beg_2017.csv'))
gindices_beg17[, .N]
str(gindices_beg17)
summary(gindices_beg17)
gindices_beg17[is.na(q_mean), .N]
which(duplicated(gindices_beg17$ID))
which(!gindices_beg17$ID %in% cattri$ID)

gindices_join <- gindices_beg17[gindices_8909, on='ID']
gindices_join[is.na(q_mean) & is.na(i.q_mean), .N]
gindices_join[is.na(q_mean) & !is.na(i.q_mean), ID]

gattri[ID %in% gindices_join[is.na(q_mean), ID], 
       list(obsbeg_day, obsend, gaps)]
```

## Check gauges 2_timeseries
###Daily
```{r Dtimeseries_daily}
Dseriesdir<- file.path(path_D, '2_timeseries')
Ddaily <- file.path(Dseriesdir, 'daily')
gdlist <- list.files(Ddaily, '.', ) %>%
  file.path(Ddaily, .)

plotf <- function(tab) {
  ggplot(tab, aes(x=date, y=qobs+0.001)) + 
  labs(title=tab$ID) +
  #geom_point(alpha=1/10) +
  geom_line() +
  geom_vline(data=tab[qobs==-999,], aes(xintercept=date), color='red') +
  scale_y_log10() + 
  theme_bw()
}

procf <- function(path) {
  fread(path) %>%
  .[, `:=`(date = as.Date(paste(YYYY, MM, DD, sep='-')),
           ID = basename(path))]
}

#Check minimum year of record
min(
  lapply(gdlist, function(path) {
    fread(path, fill=T)[, c(min(YYYY))]
  }) %>%
    unlist
)

#Sample five gauges
gdlist <- sample(gdlist, 5)

  
gdtab1 <- procf(gdlist[1]) 
summary(gdtab1)
plotf(gdtab1)

gdtab2 <- procf(gdlist[2]) 
summary(gdtab2)
plotf(gdtab2)

gdtab3 <- procf(gdlist[3]) 
summary(gdtab3)
plotf(gdtab3)

#Check those that do not have indices
lapply(gindices_join[is.na(q_mean), ID],
       function(x)  plotf(procf(file.path(Ddaily, paste0('ID_',x, '.csv'))))
)

```

###Gaps
```{r Dtimeseries_gaps}
Dgaps <- file.path(Dseriesdir, 'gaps')
gdlist <- list.files(Dgaps, '.', ) %>%
  sample(5) %>%
  file.path(Dgaps, .)

summary(fread(gdlist[1]) )
summary(fread(gdlist[2]) )
summary(fread(gdlist[3]) )
```

## Gauge shapefiles
```{r D shapefiles}
gshp <- readOGR(dsn = file.path(path_D, '3_shapefiles'), 
               layer = 'Gauges')
nrow(gshp)
names(gshp)
names(gattri) %in% names(gshp) #Check that all attributes are there
which(duplicated(gshp$ID))
gshp[duplicated(gshp$lon, gshp$lat) | duplicated(gshp$lon, gshp$lat, fromLast = T),]
str(gshp)

#Check CRS
identical(gshp@proj4string, UAD@proj4string)

# Plot it
ggplotly(
  ggplot(UADformat, aes(x = long, y = lat)) + 
  geom_point(data=as.data.frame(gshp), aes(x=as.numeric(lon), y=as.numeric(lat)),
            color='black') +
  geom_polygon(aes(group=group, fill=group), alpha=0.2, color='black') + 
  theme_bw() +
  theme(legend.position='none')
)
```