---
title: "01_PullMunge"
author: "Simon Topp"
date: "6/29/2020"
output: html_document
editor_options: 
  chunk_output_type: console
---

```{r setup, include=FALSE}
library(tidyverse)
library(lubridate)
library(purrr)
library(furrr)
library(data.table)
library(feather)
library(sf)
library(ggpmisc)
library(ggpubr)
library(Hmisc)

knitr::opts_chunk$set(echo = TRUE)
```

## Read in the data and aggregate some initial stats to see how many observations

```{r eval = F}
## Read in whatever path you saved the raw pull data to
files <- list.files('data/raw', full.names = T)
files <- files[file.size(files) > 1] ## 6 WRS tile returned no values

## Just look how many obs per lake we're getting  
pullStats <- function(file){
  df <- read_csv(file) %>%
    group_by(Hylak_id) %>%
    summarize(count = n())
  return(df)
}

## Do it all in parallel because it takes awhile with ~30 million obs
plan(multiprocess)
stats <- files %>% future_map_dfr(pullStats, .progress = T)
plan(sequential)

stats.agg <- stats %>%
  group_by(Hylak_id) %>%
  summarize(tile.count = n(),
            obs = sum(count))

meanObs <- mean(stats.agg$obs)
medObs <- median(stats.agg$obs)
```

# Looks like we have a mean of `r meanObs` observations and median of `r medObs` obeservations per lake in hydrolakes

```{r}
ggplot(stats.agg, aes(x = obs, fill = factor(tile.count))) + 
  geom_histogram() + 
  scale_fill_viridis_d() +
  theme_bw() +
  labs()
```

## Do a conservative munge on the dataset to reduce the noise

```{r eval = F}
## Munge it conservatively, basically no hillshade, zero clouds/shadow/ice and
## at least 8 pixels
## at the same time extract the date and Landsat ID to join with the metadata
munger <- function(file){
  df <- read_csv(file) %>%
    filter(!is.na(Blue),
           hillShadow == 1,
           cScore_clouds == 0,
           pCount_dswe1 > 8) %>%
    mutate(LandsatID = map_chr(`system:index`, ~str_split(.,'_0000')[[1]][1]),
         date =  map_chr(LandsatID, ~tail(strsplit(., "_")[[1]], n =1)),
         date = ymd(date)) %>%
  select(-`system:index`)
  return(df)
}

## Again, do it all in Parallel
plan(multiprocess)
munged <- files %>% future_map_dfr(munger, .progress = T)
plan(sequential)

## Save the output
write_feather(munged, 'data/sr_us_hydrolakes_dp_20200628.feather')
```

## Look at the post-munging summary stats

```{r}
munged <- read_feather('data/processed/sr_us_hydrolakes_dp_20200628.feather')
## Look at new post-munged values
summary <- setDT(munged)[,.N,.(Hylak_id)]

postMungeMean <- mean(summary$N)
sd(summary$N)
IQR(summary$N)

postMungeMed <- median(summary$N)

ggplot(summary, aes(x = N)) + 
  geom_histogram() +
  geom_vline(aes(xintercept = postMungeMean)) +
  theme_bw() +
  labs(title = '~22.3 Million Obs over 57k Lakes with an Average of 392 Obs per lake')


##Spatial Distribution and count
dp <- st_read('data/in/HydroLakes_DP/HydroLakes_DP.shp')

hl <- st_read('../USLakeClarityTrendr/in/hydroLakes/HydroLAKES_polys_v10_shp/HydroLAKES_polys_v10.shp') %>%
  filter(Hylak_id %in% dp$Hylak_id) %>%
  st_centroid()

summary.sf <- summary %>% inner_join(hl) %>% st_as_sf() %>% st_transform(102003) 

usa <- maps::map('usa', plot = F) %>% st_as_sf() %>% st_transform(102003) 

grid <- st_make_grid(usa, cellsize = c(50000,50000), square = F) %>% st_as_sf() %>% mutate(ID = row_number())
 
counts.sf <- grid %>% st_join(summary.sf, left = F) %>%
  group_by(ID) %>%
  summarise(sumLakes = n(),
            meanObs = median(N))

## Temporal Distribution
tempDist <- as.data.table(munged)[, year := year(date)
               ][, month := month(date)
                 ][,.N, by = .(year, month)
                   ][,date := ymd(paste0(year,'-',month,'-01'))
                     ][, monthLab := lubridate::month(date, label = T)]

p1 <- ggplot() + 
  geom_sf(data = usa) +
  geom_sf(data = counts.sf %>% mutate(sumLakes = ifelse(sumLakes > 100, 100, sumLakes)), aes(fill = sumLakes)) +
  scale_fill_viridis_c('Number of Lakes', trans = 'log10', breaks = c(1,10,100), labels = c(1,10,'>100')) +
  ggthemes::theme_map(base_size = 11) +
  #labs(tag = 'B') +
  theme(legend.position = 'top',
        plot.margin = unit(c(0,0,-.3,0), 'null')) +
  guides(fill = guide_colourbar(title.position="top"))


p2 <- ggplot() + 
  geom_sf(data = usa) +
  geom_sf(data = counts.sf %>%
            mutate(meanObs = ifelse(meanObs > 900, 900, meanObs)), aes(fill = meanObs)) +
  scale_fill_viridis_c('Mean Observations per Lake', option = 'plasma', direction = -1, breaks = c(300,600,900), labels = c(300, 600, '>900')) +
  ggthemes::theme_map(base_size = 11) +
  #labs(tag = 'C') +
  theme(legend.position = 'top',
        plot.margin = unit(c(0,0,-.3,0), 'null')) +
  guides(fill = guide_colourbar(title.position="top"))


p3 <- ggplot(tempDist, aes(x = date, y = N, fill = monthLab)) + 
  geom_col() + 
  annotate('text', x = ymd('1990-01-01'), y = 100000, label = 'Total Observations', fontface = 'bold') +
  scale_fill_viridis_d(option = 'cividis') + 
  scale_y_continuous(breaks = c(0, 50000, 100000), labels = c(0, '50k', '100k')) + #labels=c(10^1, 10^) +
  labs(y = 'Count', fill = 'Month') +
  ggthemes::theme_few() + 
  theme(panel.border = element_blank(),
        axis.title = element_blank(),
        legend.position = c(.46,-.4),
        legend.title = element_blank(),
        plot.margin = unit(c(0,0,.4,0), 'null'),
        legend.text = element_text(size = 10),
        legend.key.size = unit(1, 'lines'),
        legend.background = element_blank()) +
  guides(fill = guide_legend(nrow = 1, override.aes = list(size = 1)))
p3

layout <- rbind(c(2,3),
                c(2,3),
                c(2,3),
                c(1,1))

g <- gridExtra::grid.arrange(p3,p1,p2, layout_matrix = layout)

ggsave('figures/ObsDistribution.png', plot = g, width = 7, units = 'in', dpi = 300)
```

# After munging, we have  a mean of `r postMungeMean` observations and median of `r postMungeMed` observations per lake in hydrolakes.  Now join the data to the metadata and look at potential biases between sensors.

```{r}
## Metadata code available at
##https://code.earthengine.google.com/561b310819705b739ca740f9c6a231e9
meta <- read_csv('data/SceneMetadata.csv')
meta <- meta %>% rename(LandsatID = `system:index`)

refCompPre <- munged %>%
  left_join(meta %>% select(LandsatID, sat = SATELLITE, Clouds_Scene = CLOUD_COVER)) %>%
  select(LandsatID, Hylak_id, date, Aerosol, Blue, Red, Green, Nir, Swir1, Swir2, TIR1, TIR2, sat, pCount_dswe1, pCount_dswe3) 

refCompPre <- as.data.table(refCompPre)[, year := year(date)
               ][,sat := factor(sat, levels = c('LANDSAT_5','LANDSAT_7','LANDSAT_8'),
                      labels = c('l5','l7','l8'))]


## Create polynomial functions based on the 1-99th percentile of each sensor
## for overlapping periods

lm8 <- function(band){
  y <- refCompPre %>% 
    filter(year > 2012, sat == 'l7') %>% 
    .[,band] %>% 
    quantile(., seq(.01,.99, .01))
  
  x = refCompPre %>% 
    filter(year > 2012, sat == 'l8') %>% 
    .[,band] %>% 
    quantile(., seq(.01,.99, .01))
  
  lm <- lm(y~poly(x, 2, raw = T))
  
  df <- tibble(band = band, intercept = lm$coefficients[[1]], B1 = lm$coefficients[[2]], B2 = lm$coefficients[[3]])

  return(df)
}

lm5 <- function(band){
  y <- refCompPre %>% 
    filter(year > 1999, year < 2012, sat == 'l7') %>% 
    .[,band] %>% 
    quantile(., seq(.01,.99, .01))
  
  x = refCompPre %>% 
    filter(year > 1999, year < 2012, sat == 'l5') %>% 
    .[,band] %>% 
    quantile(., seq(.01,.99, .01))
  
  lm <- lm(y~poly(x, 2, raw = T))
  
  df <- tibble(band = band, intercept = lm$coefficients[[1]], B1 = lm$coefficients[[2]], B2 = lm$coefficients[[3]])

  return(df)
}


## Create the functions and look at the resulting corrections
bands <-  c('Blue', 'Green', 'Nir', 'Red', 'Swir1','Swir2', 'TIR1')
funcs.8 <- bands %>% map_dfr(lm8)
funcs.5 <- bands %>% map_dfr(lm5)


funcs.5 %>% mutate(SatCor = 'l5') %>%
  bind_rows(funcs.8 %>% mutate(SatCor = 'l8')) %>%
  write_feather(.,'data/out/landsat_correction_coefficients.feather')

write_feather(refCompPre, 'data/processed/srMunged_Original_us_hydrolakes_dp_20200628.feather')
```

## We know from previous work that even minor differences between sensors can translate to meaningful differences in dominant wavelength, so we'll apply the above corrections to the entire dataset

```{r}
l5corr <- refCompPre[sat == 'l5'
                     ][,Blue := funcs.5[1,2][[1]] +
                                      funcs.5[1,3][[1]]*Blue +
                                      funcs.5[1,4][[1]]*Blue^2
                       ][,Green := funcs.5[2,2][[1]] +
                           funcs.5[2,3][[1]]*Green + 
                           funcs.5[2,4][[1]]*Green^2
                         ][,Nir := funcs.5[3,2][[1]] + 
                             funcs.5[3,3][[1]]*Nir + 
                             funcs.5[3,4][[1]]*Nir^2
                          ][,Red := funcs.5[4,2][[1]] + 
                             funcs.5[4,3][[1]]*Red + 
                             funcs.5[4,4][[1]]*Red^2
                          ][,Swir1 := funcs.5[5,2][[1]] +
                              funcs.5[5,3][[1]]*Swir1 +
                              funcs.5[5,4][[1]]*Swir1^2
                            ][,Swir2 := funcs.5[6,2][[1]] +
                              funcs.5[6,3][[1]]*Swir2 +
                              funcs.5[6,4][[1]]*Swir2^2
                            ][,TIR1 := funcs.5[5,2][[1]] + 
                             funcs.5[5,3][[1]]*TIR1 + 
                             funcs.5[5,4][[1]]*TIR1^2]

l8corr <- refCompPre[sat == 'l8'
                     ][,Blue := funcs.8[1,2][[1]] +
                                      funcs.8[1,3][[1]]*Blue +
                                      funcs.8[1,4][[1]]*Blue^2
                       ][,Green := funcs.8[2,2][[1]] + 
                           funcs.8[2,3][[1]]*Green + 
                           funcs.8[2,4][[1]]*Green^2
                         ][,Nir := funcs.8[3,2][[1]] + 
                             funcs.8[3,3][[1]]*Nir + 
                             funcs.8[3,4][[1]]*Nir^2
                          ][,Red := funcs.8[4,2][[1]] + 
                             funcs.8[4,3][[1]]*Red + 
                             funcs.8[4,4][[1]]*Red^2
                          ][,Swir1 := funcs.8[5,2][[1]] +
                              funcs.8[5,3][[1]]*Swir1 +
                              funcs.8[5,4][[1]]*Swir1^2
                            ][,Swir2 := funcs.8[6,2][[1]] +
                              funcs.8[6,3][[1]]*Swir2 +
                              funcs.8[6,4][[1]]*Swir2^2
                              ][,TIR1 := funcs.8[5,2][[1]] + 
                             funcs.8[5,3][[1]]*TIR1 + 
                             funcs.8[5,4][[1]]*TIR1^2]

srCor <- refCompPre %>%
  filter(sat == 'l7') %>%
  bind_rows(l5corr) %>%
  bind_rows(l8corr)

rm(l5corr, l8corr)

srCor <- as.data.table(srCor)[,dWL := fui.hue(Red, Green, Blue)]
write_feather(srCor, 'data/processed/srCorrected_us_hydrolakes_dp_20200628.feather')
write_csv(srCor, 'data/processed/srCorrected_us_hydrolakes_dp_20200628.csv')
```


## Make a function for comparing corrected vs uncorrected values

```{r}
correctionPlot <- function(band, sat, dataPre, dataPost){
  if(sat == 'l8'){

    df <- tibble(l7 = quantile(dataPost[sat == 'l7' & year > 2012, ..band
                           ][[1]], seq(.01,.99,.01)),
    Original = quantile(dataPre[sat == sat & year > 2012, ..band
                                 ][[1]], seq(.01,.99,.01)),
    PostCorrection = quantile(dataPost[sat == sat & year > 2012, ..band
                                       ][[1]], seq(.01,.99,.01)))
  }else if(sat == 'l5'){

    df <- tibble(l7 = quantile(dataPost[sat == 'l7' & year < 2012 & year > 1999, ..band
                           ][[1]], seq(.01,.99,.01)),
    Original = quantile(dataPre[sat == sat & year < 2012 & year > 1999, ..band
                                ][[1]], seq(.01,.99,.01)),
    PostCorrection = quantile(dataPost[sat == sat & year < 2012 & year > 1999, ..band
                                       ][[1]], seq(.01,.99,.01)))
  }
  
  ogBias <- round(Metrics::bias(df$l7, df$Original), 3)
  CBias <- round(Metrics::bias(df$l7, df$PostCorrection), 3)
  
  df <- df %>% gather(Original, PostCorrection, key = "Correction", value = 'Reflectance')
  
  ggplot(df, aes(x = l7, y = Reflectance, color = Correction)) + geom_point(alpha = .8) + 
    geom_abline(color = 'red') + 
    scale_color_viridis_d(end = .7, labels = c('Original', 'Post\nCorrection')) +
    #stat_regline_equation(aes(label =  paste(..adj.rr.label..))) +
    annotate('text', x= Inf, y = Inf, vjust = 3.8, hjust = 1, 
             label = paste0('Original Bias: ', ogBias,'\nCorrected Bias: ', CBias)) +
    theme_bw() +
    theme(axis.title = element_blank()) +
    scale_y_continuous(trans = 'log10') +
    scale_x_continuous(trans = 'log10') +
    labs(title = capitalize(band))
}
```


## Plot up some examples

```{r}
#srCor <- read_feather('data/processed/srCorrected_us_hydrolakes_dp_20200628.feather')
#refCompPre <- read_feather('data/processed/srMunged_Original_us_hydrolakes_dp_20200628.feather')

srCor <- as.data.table(srCor)
refCompPre <- as.data.table(refCompPre)
p1 <- correctionPlot('Blue', 'l5', refCompPre, srCor)
p2 <- correctionPlot('Green', 'l5', refCompPre, srCor)
p3 <- correctionPlot('Red', 'l5', refCompPre, srCor)
p4 <- correctionPlot('Nir', 'l5', refCompPre, srCor)
#p5 <- correctionPlot('TIR1', 'l5', refCompPre, srCor)
p6 <- correctionPlot('Blue', 'l8', refCompPre, srCor)
p7 <- correctionPlot('Green', 'l8', refCompPre, srCor)
p8 <- correctionPlot('Red', 'l8', refCompPre, srCor)
p9 <- correctionPlot('Nir', 'l8', refCompPre, srCor)
#p10 <- correctionPlot('TIR1', 'l8', refCompPre, srCor)

g1 <- ggarrange(p1, p2, p3, p4, nrow = 2, ncol = 2, common.legend = T) 
g2 <- ggarrange(p6,p7,p8,p9, nrow = 2, ncol = 2, common.legend = T)

g11 <- gridExtra::grid.arrange(g1, bottom = 'Landsat 7 Surface Reflectance', left = 'Landsat 5 Surface Reflectance')
g22 <- gridExtra::grid.arrange(g2, bottom = 'Landsat 7 Surface Reflectance', left = 'Landsat 8 Surface Reflectance')

g <- gridExtra::grid.arrange(g11,g22, nrow = 2)               

ggsave('figures/srCorrections.png', plot = g, width = 5.5, height = 9, units = 'in')

rm(g1, g2, g11, g22, g)
```

## Look at some general patterns in dWL and Reflectance.

```{r}
## By period

srCorSummary.sf <- merge(srCor, hl) %>% st_as_sf()

usa <- maps::map('usa', plot = F) %>% st_as_sf() %>% st_transform(st_crs(hl)) 

grid <- st_make_grid(usa, cellsize = c(100000,1000000), square = F) %>% st_as_sf() %>% mutate(ID = row_number())
 
gridSummary <- grid %>% st_join(srCorSummary.sf, left = F) %>%
  set_geometry(NULL) %>%
  as.data.table(.)[,ref := (Blue + Green + Red + Nir)/4
                   ][,.(medianRef = median(ref, na.rm = T), mediandWL = median(dWL, na.rm = T)), by = ID]


```

#Visualize the Deepest point for a figure

```{r}
hl <- st_read('../USLakeClarityTrendr/in/hydroLakes/HydroLAKES_polys_v10_shp/HydroLAKES_polys_v10.shp') %>%
  filter(Hylak_id %in% dp$Hylak_id) 

hl_sample <- dp %>% distinct(Hylak_id) %>% sample_n(1000)

dpSamp <- dp %>% filter(Hylak_id %in% hl_sample$Hylak_id) %>%
  mutate(Measurement = factor(type, labels = c('Centroid', 'Deepest Point'))) %>%
  st_transform(102003)


addBuff <- function(i){
  point <- dpSamp[i,]
  st_buffer(point, dist = point$distance)
}

dpBuff <- addBuff(dpSamp[1,])

for(i in c(2:nrow(dpSamp))){
  dpBuff <- rbind(dpBuff, addBuff(dpSamp[i,]))
}

dpBuff <- c(1:nrow(dpSamp)) %>% map_dfr(addBuff)

mapview(dpBuff)

dpBuff <- dpSamp %>%
  mutate(geometry = split(., 1:5) %>% purrr::map(st_bbox))
  rowwise() %>%
  st_as_sf() %>%
  st_buffer(., dist = distance) %>%
  ungroup()


hl_poly <- hl %>% filter(Hylak_id %in% hl_sample$Hylak_id) %>% st_transform(102003)
  
mapview(dpBuff, zcol = 'Measurement', color = 'grey60', col.regions = c('cyan','magenta'), layer.name = 'Sample Point and Associated\nDistance to Nearest Shore') + 
  mapview(dpSamp, zcol = 'Measurement',layer.name = 'Sample Point', color = 'black', col.regions = c('cyan','magenta')) + mapview(hl_poly, col.regions = 'black', layer.name = 'Lake Polygon')
```