1_nhd_join_and_munge.Rmd

---
title: "1_nhd_join_and_munge"
author: "Simon Topp"
date: "11/15/2018"
output: html_document
editor_options: 
  chunk_output_type: console
---

# This script takes the downloaded reflectance data from Google Earth Engine and joins it back up with the original water quality information as well as lake, catchement, and watershed data from NHDPlusV2 and LakeCat.

## Load in the downloaded data from AquaSat

```{r}
# This initial dataset is archived on figshare (10.6084/m9.figshare.12227273)
sr <- read.csv('../Aquasat/3_rswqjoin/data/out/sr_wq_rs_join.csv', stringsAsFactors = F) %>%
  filter(type == 'Lake') %>%
  select(-c(p_sand, tis, doc)) %>%
  gather(secchi, tss, chl_a, key = parameter, value = value) %>%
  filter(!is.na(value)) %>%
  mutate(UniqueID = row_number())

# Create Ecoregions shapefile the same as those used in the NLA lake assessments to join to our training data following Herlihy  et al. 2008 https://bioone.org/journals/Freshwater-Science/volume-27/issue-4/08-081.1/Striving-for-consistency-in-a-national-assessment--the-challenges/10.1899/08-081.1.full

# ecoRegs <- st_read('in/NLA/us_eco_l3/us_eco_l3.shp') %>%
#   mutate(region = ifelse(US_L3CODE %in% c(33,35,65, 74, 63, 73, 75, 84,34,76), 'Coastal Plain',
#                           ifelse(US_L3CODE %in% c(58, 62, 82, 59, 60, 61, 83), "Northern Appalachians",
#                                  ifelse(US_L3CODE %in% c(45, 64, 71, 36, 37, 38, 39, 66, 67, 68, 69, 70), 'Southern Appalachians',
#                                         ifelse(US_L3CODE %in% c(49,50, 51,52,56), 'Upper Midwest',
#                                                ifelse(US_L3CODE %in% c(53, 54, 55, 57, 72, 40, 46, 47, 48, 28), 'Temperate Plains', ifelse(US_L3CODE %in% c(42,43), 'Northern Plains', ifelse(US_L3CODE %in% c(44, 25, 26, 27, 29, 30, 32, 31), 'Southern Plains', ifelse(US_L3CODE %in% c(10, 12, 13, 18, 20, 22, 80, 14, 81, 24, 6, 7, 79, 85), 'Xeric West', ifelse(US_L3CODE %in% c(11,15,16,17,19,21,4,41,11,15,16,17,29,21,4,41,11,5,77,78,9,1,2,3,23,8), 'Western Mountains', NA)))))))))) %>%
#   group_by(region) %>%
#   summarize(inputs = glue::collapse(unique(US_L3CODE), sep = ', '))
# 
#st_write(ecoRegs, 'in/NLA/NLA_Ecoregions/EcoRegsMerged.shp')

region <- st_read('in/NLA/NLA_Ecoregions/EcoRegsMerged.shp')

sr.sf <- st_as_sf(sr,coords=c('long','lat'), crs=4326) %>%
  st_transform(st_crs(region))

```

## Download lakeCat data
```{r, eval = F}
# library(RCurl)
# 
# url = "ftp://newftp.epa.gov/EPADataCommons/ORD/NHDPlusLandscapeAttributes/LakeCat/FinalTables/"
# filenames = getURL(url, ftp.use.epsv = FALSE, dirlistonly = TRUE) %>%
#   strsplit(., "\n") %>%
#   unlist() %>%
#   grep(.,pattern = '.zip', value = T)
# 
# for (file in filenames) {
#   download.file(paste0(url, file), paste0('lakeCat/', file))
# }
# 
# for (file in filenames){
#   unzip(zipfile = paste0('lakeCat/',file), exdir = 'lakeCat/')
# }
```

## Join the data to NHD to extract additional landscape and lake specific variables

```{r packages_sites, eval=TRUE, echo=TRUE, warning=FALSE, include=TRUE, message=FALSE}
#Make an SF object of individual sites
sites.sf <- sr %>%
  distinct(SiteID, lat, long) %>%
  mutate(x = long, y = lat) %>% #Preserve lat and long in ouput object
  st_as_sf(.,coords=c('x','y'), crs=4326) %>%
  st_transform(st_crs(region))

##Join sites to the regions
plan(multiprocess)
reg.Join <- sites.sf %>%
  split(., c(1:10)) %>%
  future_map(~st_join(., region %>% st_simplify(dTolerance = 500), join = st_nearest_feature), .progress = T)
plan(sequential)

## For some reason map_dfr and row_bind don't play well with spatial data so join them the clunky way.
reg.Join <- rbind(reg.Join[[1]], reg.Join[[2]], reg.Join[[3]], reg.Join[[4]], 
                  reg.Join[[5]], reg.Join[[6]],reg.Join[[7]], reg.Join[[8]], 
                  reg.Join[[9]], reg.Join[[10]])

## Load in nhd lake data downlowded from, http://www.horizon-systems.com/NHDPlus/V2NationalData.php
nhdLakes <- st_read('D:/GIS_Data/NHD/nhdPlusV2_WaterBodiesFull.shp') %>%
  filter(FTYPE == 'LakePond' | FTYPE == 'Reservoir') %>%
  st_transform(st_crs(region))

plan(multiprocess)
centerJoin <- reg.Join %>%
  st_buffer(90) %>% ##Add buffer to capture sample sites at lake edges
  split(., c(1:10)) %>%
  future_map(~st_join(.,nhdLakes), .progress = T)
plan(sequential)

centerJoin <- rbind(centerJoin[[1]], centerJoin[[2]], centerJoin[[3]], centerJoin[[4]], centerJoin[[5]],
                    centerJoin[[6]], centerJoin[[7]], centerJoin[[8]], centerJoin[[9]], centerJoin[[10]])

##2680 sites aren't within 90 meters of a waterbody and ~178 sites are within 90 meters of 2 or more bodies
check <- centerJoin %>%
  group_by(SiteID) %>%
  summarize(count = n()) %>%
  filter(count > 1)
check <- centerJoin %>% filter(SiteID %in% check$SiteID)

##Do join with no buffer for sites with multiple returns
dup.join <- reg.Join %>% filter(SiteID %in% check$SiteID) %>% st_join(.,nhdLakes)

## Remove duplicates and join up with the unbuffered join.
centerJoin <- centerJoin %>% 
  st_centroid() %>%
  filter(!SiteID %in% check$SiteID) %>% rbind(dup.join)

## Real quick look at how many obs we loose from no COMID 
missing <- centerJoin %>% filter(is.na(COMID))
check <- sr %>% filter(SiteID %in% missing$SiteID)
## Its about 24k obs, pretty close to what we get with HydroLakes Join, of those ~ 12k are secchi.

##Finally, munge the landscape variables from LakeCat and connect it to our lake data

## Connect it all to the lakeCat data previously downloaded using the COMID identifier
lakeCat <- tibble(COMID = unique(centerJoin$COMID)) 
lc.files <- list.files('in/lakeCat/unzip', full.names = T)

for(i in lc.files){
  if(i == first(lc.files)){
    lc <- read.csv(i)
    lakeCat.full <- lakeCat %>%
      left_join(lc, by ='COMID')}
  else{
    lc <- read.csv(i) %>%
      select(-c(CatAreaSqKm, WsAreaSqKm, CatPctFull,WsPctFull,inStreamCat))
    lakeCat.full <- lakeCat.full %>%
      left_join(lc, by = 'COMID')
  }
}


#write_feather(lakeCat.full,'out/lakeCatFull.feather')
lakeCat.full <- read_feather('out/lakeCatFull.feather')


#############  Munge the lakeCat data a little:
# Look at the names of the variables to select appropriate inputs. See 
# ftp://newftp.epa.gov/EPADataCommons/ORD/NHDPlusLandscapeAttributes/LakeCat/Documentation/VariableList-QuickReference.html 
# for variable descriptions.

names.in <- names(lakeCat.full)

# Check how many Ws vs Cat variables there are.
length(grep(names.in, pattern = '*Ws'))
length(grep(names.in, pattern = '*Cat'))

#Double check to make sure watershed and catchment variables are the same.
names.ws <- names.in[grep(names.in, pattern = '*Ws')] %>%
  gsub(., pattern = 'Ws', replacement = '')

names.cat <- names.in[grep(names.in, pattern = '*Cat')] %>%
  gsub(., pattern = 'Cat', replacement = '')

# Look at what we end up with if we remove all the watershed ones
names.cat[!(names.cat %in% names.ws) == T]
grep(names.in, pattern = '*Ws', value = T, invert = T)


## Remove unwanted parameters and munge a couple other ones.  Unfortunately a lot of the date related ones are at too coarse a scale to really use in the modelling pipeline.  As a result, we'll keep the 2006 ones (in case we want to use them as 'summary' type variables), but drop the rest.  We'll also assume that catchment variables are better predictors than coarser watershed ones.

lakeCat.munged <- lakeCat.full %>%
  filter(!is.na(CatAreaSqKm)) %>%  ## ~125 sites from SR weren't within 90 meters of an NHD waterbody.
  select(-c(grep(names.in, pattern = '*Ws', value = T), 
            grep(names.in, pattern = '2011Cat', value = T),
            ### 2011 values aren't representative of the study period
            grep(names.in, pattern = 'PctFire', value = T), 
            grep(names.in, pattern = 'PctFrstLoss', value = T),
            ### While potentially important, the forest metrics lack 
            ### the temporal resolution we need.
            grep(names.in, pattern = 'NABD', value = T), 
            ### Dam info acts as identifiers and are redundant with lake_type/area
            grep(names.in, pattern = 'Dam', value = T), 
            ### Dam info acts as identifiers and are redundant with lake_type/area
            'Precip08Cat', 'Precip09Cat', 'Tmean08Cat', 'Tmean09Cat')) %>%  
            ### Use long term averages not these
  mutate(inStreamCat = factor(inStreamCat))

## Look at variables and round everything to a reasonable number of digits to avoid ID like variables.
summary(lakeCat.munged)

#Round the majority of values to the nearest 1 to avoid variables become 'unique ID's)
round1 <- names(lakeCat.munged %>% select(c(CatAreaSqKm:CatPctFull,PctImp2006Cat, PctCarbResidCat:WetIndexCat)))

## Round a few to 10 because they're variation is real big and larger values with 
#Round a subset to the nearest 10th because there values more or less range from 0 to 1.
round.1 <- names(lakeCat.munged %>% select(c(AgKffactCat, KffactCat, MineDensCat)))

lakeCat.munged <- lakeCat.munged %>%
  mutate_at(round1, round, digits = 0) %>%
  mutate_at(round.1, round, digits = 1)

## Double check
summary(lakeCat.munged)

site.data.full <- centerJoin %>% left_join(lakeCat.munged)

#write_feather(site.data.full %>% as.data.frame() %>% select(-geometry), 'out/siteInfoFull.feather')
site.data.full <- read_feather('out/siteInfoFull.feather')
```

# Check for duplicates between lagos and wqp now that we have lakes linked to SiteIDs

```{r}
# Now that we have everything linked up to lakes 
# we can check for observations that may be getting 
# double counted due to WQP and Lagos inclusion.

obs.check <- sr %>% 
  left_join(site.data.full %>% as_data_frame(.) %>% select(-geometry)) %>%
  mutate(doubleID = paste0(COMID, parameter, value, landsat_id, date, lat, long)) %>%
  group_by(doubleID) %>%
  summarize(count = n()) %>%
  filter(count > 1)

sum(obs.check$count)

# Look at the duplicates to make sure they're all in LAGOs area
check <- sr %>% 
  left_join(site.data.full %>% as_data_frame(.) %>% select(-geometry)) %>%
  mutate(doubleID = paste0(COMID, parameter, value, landsat_id, date, lat, long)) %>%
  filter(doubleID %in% obs.check$doubleID) %>%
  st_as_sf(coords=c('long','lat'),crs=4326) %>%
  st_transform(.,2163)

# They aren't :(.  This just means there are duplicates in WQP which isn't surprising.
mapview(check)

# Look at the duplicates a little
duplicates <- check %>%
  group_by(doubleID) %>%
  summarise(count = n(),
            red_sd = sd(red),
            blue_sd = sd(blue),
            green_sd = sd(green),
            nir_sd = sd(nir)) %>%
  mutate(SameRef = ifelse(red_sd + blue_sd + green_sd + nir_sd == 0, 'Same', 'Different'))

ggplot(duplicates) + geom_bar(aes(x = SameRef, y = ..count..))
ggplot(check) + geom_bar(aes(x = source))

## They all have the exact same reflectance values which means they're truly duplicates (they were entered into WQP twice with different SiteIDs.  We'll just get rid of them by calling distinct further on.
```


## Workflow for creating munged parameter dataframe from fully joined reflectance, NHD, and hydrolakes data.
```{r}
# Function for munging full data based on reasonable values.  
munger <- function(df, param, dswe.lim = 1, dswe.sd.lim = 0, minPix = 8, maxClouds = 50, maxTimeDiff = 129600, maxRef = 2000, minRef = 1, h.shadow = 1, h.shadow_sd = 0){
  
  minValue = 0.001
  maxValue = quantile(df$value, .99)
  
  if(param == 'doc'){
    minValue = 0.1  #from sobek 2007
    maxValue = 350}
  if(param == 'chl_a'){
    minValue = 0.001 #from NLA 2012 minimum detection limit
    maxValue = 764}  # NLA max reported
  if(param == 'tss'){
    minValue = 0.001 # Using the same as Chl-a cause I'm not sure
    maxValue = 5000} # Just leaving um all in there because I don't really know.
  if(param == 'secchi'){
    minValue = 0.1
    maxValue = 45} #deepest recorded lake sdd by Larson et al 2007

  data <- df %>%
    filter(pixelCount > minPix, # We want at least 9 good water pixels
         hillshadow == h.shadow, # No shadow effects
         hillshadow_sd <= h.shadow_sd, 
         clouds < maxClouds,
         abs(timediff) < maxTimeDiff,
         dswe == dswe.lim, ##Only high confidence water
         dswe_sd <= dswe.sd.lim,
         value > minValue,
         value < maxValue) %>%
    filter_at(vars(blue,green,red,nir,swir1,swir2),all_vars(.>minRef & .<maxRef)) %>%
    mutate(date = ymd(date),
           date_unity = ymd_hms(date_unity),
           month = month(date),
           year = year(date),
           sat = factor(sat)) %>%
    select(-c(system.index, pixelCount, qa, qa_sd, date_only)) 
return(data)
}

files <- list.files('in/MERRA2_data2', full.names = T)

munged <- sr %>%
  group_by(parameter) %>%
  nest() %>%
  mutate(munge = purrr::map2(.x = data, .y = parameter, ~munger(.x, param = .y, dswe.sd = .3))) %>%
  select(-data) %>%
  unnest(munge) %>%
  left_join(site.data.full) %>%
  filter(!is.na(COMID)) %>%
         ## We loose ~20k obs which aren't within 90 meters of a NHD Lake 
  mutate(Reservoir = factor(ifelse(FTYPE == "Reservoir", 1, 0))) %>%
  distinct(COMID, value, parameter, landsat_id, date, lat, long, .keep_all = T) %>% #remove duplicate values (see workflow at line ~200)
## Finally, tag on Aerosol Optical Depth from the MERRA-2 Reanalysis
  # rowwise() %>%
  # mutate(AOD = pullMerra(date = date, lat = lat, long =long)) %>% ##This awhile
  # ungroup() %>%
  select(-.geo)
  ##Join to NLA EcoRegions


munged <- munged %>% 
  #select(-geometry) %>%
  rename_at(vars(FDATE:LakeArea), tolower) #Because of previous framework.
write_feather(munged, 'out/srMunged.feather')
#munged <- read_feather('out/srMunged.feather')
```

##Take a look at the dist and count of the final

```{r}
hucs <- st_read('in/hucs/hucs.shp')

sr.sf <- read_feather('out/srMunged.feather') %>%
  filter(parameter == 'secchi',
         value <= 10,
         !is.na(CatAreaSqKm)) %>%
  st_as_sf(coords = c('long', 'lat'), crs = 4326) %>%
  st_transform(st_crs(hucs))


usa = st_as_sf(maps::map('usa',plot=F,fill=T)) %>%
    st_transform(.,2163) %>%
    st_buffer(.,0) 

huc.sr.sf <- st_join(hucs,sr.sf) %>%
  as_tibble() %>%
  filter(!is.na(value)) %>%
  group_by(HUC_8) %>%
  summarise(Overpasses=n()) %>%
  filter(Overpasses > 0) %>%
  left_join(hucs, by = 'HUC_8')

ggplot() +
  geom_sf(data = usa, aes(geometry = geometry), fill = 'grey20') + 
  geom_sf(data = huc.sr.sf,aes(fill=Overpasses, geometry = geometry), color=NA)  + 
  scale_fill_gradient2(low='#a50026',mid='#abd9e9', high='#313695',na.value='black',midpoint = 2, trans='log10', name = 'No. of Coincident\nCloud Free Observations') +
  theme_bw() +
  #ggtitle('Cloud free coincident field and Landsat Observations') +
  theme(legend.position = 'bottom', 
        legend.direction = 'horizontal',
        panel.grid.major = element_line(color = 'white'),
        panel.border = element_blank(),
        axis.text = element_blank(),
        axis.ticks = element_blank(),
        plot.title = element_text(hjust = .5))

ggsave('figures/MungedOverpasses.png', width = 6, height = 5, units = 'in', dpi = 300)
```

```{r}
## Finally, we need to standardize our reflectance values.  To do so we'll use the entire stack of images from the NLA and standardize reflectance values between the percentiles of distributions for each satellite during their coincident overpass period.
files <- list.files('lake_data/NLA2012_cntr/', full.names = T)

refCompPre <- purrr::map_dfr(files, read.csv, stringsAsFactors = F) %>%
  filter(!is.na(blue),
         dswe == 1,
         dswe_sd < .4,
         hillshadow == 1,
         cScore == 0,
         Cloud_Cover < 50,
         pixelCount > 8) %>%
  filter_at(vars(red, blue, green, nir, swir1, swir2), all_vars(.>0 & .<2000)) %>%
  mutate(dWL = fui.hue(red, green, blue),
         date = ymd_hms(date),
         year = year(date),
         sat = factor(sat, levels = c(5,7,8), labels = c('l5','l7','l8')),
         source = 'NLA') %>%
  select(year, blue, red, green, nir, dWL, sat, source, COMID)

## Create functions for calculating the polynomial regression equation
lm8 <- function(band, x){
  x <- x %>% filter(year > 2012, source == 'NLA')
  df <- tibble(l7 = quantile(x %>% filter(sat == 'l7') %>% .[,band], 
                             seq(.01,.99, .01)), 
               l8 = quantile(x %>% filter(sat == 'l8') %>% .[,band], 
                             seq(.01,.99, .01)))
  lm <- lm(l7~poly(l8, 2), data = df)
  #lut <- tibble(quantile = names(lm$fitted.values), value = lm$fitted.values)
  return(lm)
}

lm5 <- function(band, x){
  x <- x %>% filter(year > 1999, year < 2012, source == 'NLA') 
  df <- tibble(l7 = quantile(x %>% filter(sat == 'l7') %>% .[,band], 
                             seq(.01,.99, .01)), 
               l5 = quantile(x %>% filter(sat == 'l5') %>% .[,band], 
                             seq(.01,.99, .01)))
  lm <- lm(l7~poly(l5, 2), data = df)
  return(lm)
}

##Check out the resulting corrections a little
## Function for plotting up the resutls
correctionPlot <- function(band, sat){
  if(sat == 'l8'){
    x <- refCompPre %>% filter(source == 'NLA', year > 2012)
    }else if(sat == 'l5'){
      x <- refCompPre %>% filter(source == 'NLA', year > 1999, year < 2012)
    }
  max.x <- max(x[,band])
  df <- tibble(y = quantile(x %>% filter(sat == 'l7') %>% .[,band], 
                               seq(.01,.99, .01)), 
               x = quantile(x %>% filter(sat == sat) %>% .[,band], 
                               seq(.01,.99, .01)))
  lm <- lm(y~poly(x, 2), data = df)
  r2 <- summary(lm)$adj.r.squared %>% round(4)
  
  df <- df %>% bind_cols(Corrected = lm$fitted.values) %>%
    gather(x, Corrected, key = 'Reflectance', value = x) %>%
    mutate(Reflectance = ifelse(Reflectance == 'x', 'Original', Reflectance))
  
  
  ggplot(df, aes(x = x, y = y, color = Reflectance)) + geom_point(alpha = .5) + 
    geom_abline(color = 'red') + 
    annotate('text', label = paste0("italic(R)^2 ==", r2), parse = T, x = -Inf, y = Inf, vjust = 1.2, hjust = -.2) +
    scale_color_viridis_d(begin = .2, end = .8) +
    theme_bw() +
    #scale_y_continuous(sec.axis = sec_axis(~.*(100/max.x), name = 'Quantile')) + 
    #labs(x = paste0(capitalize(sat), " SR"), y = 'L7 SR', title = paste0(capitalize(sat),' ', capitalize(band), " Correction")) +
    labs(title = paste0(capitalize(band))) +
    theme(legend.position = 'none')#,
          #axis.title = element_blank())
}

p1 <- correctionPlot('blue', 'l5')
p2 <- correctionPlot('green', 'l5')
p3 <- correctionPlot('red', 'l5')
p4 <- correctionPlot('nir', 'l5')
p5 <- correctionPlot('blue', 'l8')
p6 <- correctionPlot('green', 'l8')
p7 <- correctionPlot('red', 'l8')
p8 <- correctionPlot('nir', 'l8')

legend <-  cowplot::get_legend(p1 + theme(legend.position = 'bottom')) 
g1 <- arrangeGrob(p1,p2,p3,p4, nrow = 2, bottom = 'L5 Surface Reflectance')
g2 <- arrangeGrob(p5,p6,p7,p8, nrow = 2, bottom = 'L8 Surface Reflectance')

g <- grid.arrange(g1,g2, nrow = 2, left = "L7 Surface Reflectance", top = legend)
ggsave('figures/CorrectionPlots.png', plot = g, height = 8, width = 4, units = 'in')

## Example for AGU talk
correctionPlot('nir', 'l8') + labs(x = 'Landsat 8 NIR Reflectance', y = 'Landsat 7 NIR Reflectance', title = 'Reflectance Correction Procedure') + theme(legend.position = c(.8, .2))
ggsave('figures/NirCorrExample.png', height = 4, width = 4, units = 'in')


## Generate functions to actually apply to the data in the analysis
bands <- tibble(band = c('blue', 'green', 'nir', 'red')) 
funcs.8 <- bands %>%
  mutate(lm8 = purrr::map(band, lm8, x = refCompPre))
funcs.5 <- bands %>%
  mutate(lm5 = purrr::map(band, lm5, x = refCompPre))

save(funcs.8, funcs.5, file = 'models/landsat_poly_corrs.Rdata')

##Apply the correction to the entire dataset.  Done to SR data in 00_ParameterConfigurations
refCompPost <- refCompPre %>%
  mutate(uniqueID = row_number()) %>%
  select(uniqueID, blue:nir, sat) %>%
  gather(blue:nir, key = 'band' , value = 'value') %>%
  spread(sat, value) %>%
  group_by(band) %>%
  nest() %>%
  left_join(funcs.8) %>%
  left_join(funcs.5) %>%
  mutate(pred8 = map2(lm8, data, predict),
         pred5 = map2(lm5, data, predict)) %>%
  select(-c(lm8,lm5)) %>%
  unnest(data, pred8, pred5) %>%
  select(-c(l8,l5)) %>%
  rename(l8 = pred8, l5 = pred5) %>% gather(l5,l7,l8, key = 'sat', value = 'value') %>%
  spread(band, value) %>%
  na.omit() %>%
  left_join(refCompPre[,c('year','source','COMID','uniqueID')]) %>%
  mutate(correction = 'Post.Correction',
         NR = nir/red,
         BG = blue/green,
         dWL = fui.hue(red, green, blue))
```