0_surfaceReflectanceCorrection.Rmd

---
title: "SensorCorrection"
output: html_document
date: "2023-03-13"
editor_options: 
  chunk_output_type: console
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

# Purpose of file
Read in lake and channel reflectance data and adjust Landsat 5 and Landsat 8 data to Landsat 7 reflectances

# Import Libraries
```{r}
library(data.table)
library(tidyverse)
library(dtplyr)
library(feather)
library(Metrics)
library(ggpubr)
```

# Functions
```{r}
#takes RGB to calculate dominant wavelength
chroma <- function(R, G, B) {  
  require(colorscience)

# Convert R,G, and B spectral reflectance to dominant wavelength #based
# on CIE chromaticity color space

# see Wang et al 2015. MODIS-Based Radiometric Color Extraction and
# Classification of Inland Water With the Forel-Ule
# Scale: A Case Study of Lake Taihu

# chromaticity.diagram.color.fill()
Xi <- 2.7689*R + 1.7517*G + 1.1302*B
Yi <- 1.0000*R + 4.5907*G + 0.0601*B
Zi <- 0.0565*G + 5.5943*B

# calculate coordinates on chromaticity diagram
x <-  Xi / (Xi + Yi +  Zi)
y <-  Yi / (Xi + Yi +  Zi)
z <-  Zi / (Xi + Yi +  Zi)

# calculate hue angle
alpha <- atan2( (x - (1/3)), (y - (1/3))) * 180/pi

# make look up table for hue angle to wavelength conversion
cie <- cccie31 %>%
  dplyr::mutate(a = atan2( (x - (1/3)), (y - (1/3))) * 180/pi) %>%
  dplyr::filter(wlnm <= 700) %>%
  dplyr::filter(wlnm >=380) 

# find nearest dominant wavelength to hue angle
wl <- cie[as.vector(sapply(alpha,function(x) which.min(abs(x - cie$a)))) , 'wlnm']

return(wl)
}

## Make a function for comparing and plotting corrected vs uncorrected values
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))
}

```

# code adapted from https://github.com/GlobalHydrologyLab/LakeReflectanceRepo/blob/master/01_LakeExtractor.Rmd 
```{r} 
# 1. import lake & channel reflectance (only keep mean reflectance, we don't need the 10th & 90th% reflectance--using those didn't improve results)
setwd("C:/Users/whyana/OneDrive/DocumentsLaptop/001_GraduateSchool/Research/Connectivity/Mackenzie/Data/GEE Downloads")

all.lakes = fread('MackenzieLakeExport_20230323_2.csv')  %>% 
  select(!matches("p10") & !matches("p90"))
all.channels = fread('MackenzieChannelExport_20230313_allYears.csv')  %>% 
  select(!matches("p10") & !matches("p90"))

# 2. import image metadata & prep for analysis
metadata.df = fread("SceneMetadata.csv") %>% lazy_dt() %>% 
  rename(Landsat = SPACECRAFT_ID) %>% 
  mutate(Landsat = parse_number(Landsat),
         dateTime = as_datetime(`system:time_start`*0.001),
         date = as_date(dateTime)) %>% 
  select(-`system:index`, -.geo)

# 3. Filter lake data to only low cloud observations
lakes.maxPix = all.lakes %>% group_by(OBJECTID) %>% summarise(max.pix = max(Red_count))# count max number of pix--> likely the clear sky total pixel count of the lake
lakes.filter = all.lakes %>% lazy_dt() %>% 
  left_join(lakes.maxPix, by="OBJECTID") %>% 
  mutate(dateTime = as_datetime(`system:time_start_mean`*0.001),
         date = as_date(dateTime),
         year = year(date),
         month = month(date),
         doy = yday(date)) %>% 
  dplyr::filter(Red_count/max.pix >= 0.5) %>% 
  rename(Landsat = constant_mean, WRS_PATH = WRS_PATH_mean, WRS_ROW=WRS_ROW_mean,
         `system:time_start`=`system:time_start_mean`) %>% 
  filter(Landsat!=9) %>% as.data.table() # remove landsat 9 bcs not enough data


# 4. Filter channel data to only low cloud observations
channels.maxPix = all.channels %>% group_by(fid) %>% summarise(max.pix = max(Red_count))# count max number of pix--> likely the clear sky total pixel count of the lake
channels.filter = all.channels %>% lazy_dt() %>%
  left_join(channels.maxPix, by="fid") %>%
  mutate(dateTime = as_datetime(`system:time_start_mean`*0.001),
         date = as_date(dateTime),
         fid=fid*(-1)) %>% # make it so we can tell fids from objectids
  dplyr::filter(Red_count/max.pix >= 0.5) %>%
  rename(Landsat = constant_mean, WRS_PATH = WRS_PATH_mean, WRS_ROW=WRS_ROW_mean,
         `system:time_start`=`system:time_start_mean`,
         OBJECTID = fid) %>%
  filter(Landsat!=9) # remove landsat 9 bcs not enough data


# 4. Join lake and channel data to related image metadata
## Lake
data.join = lakes.filter%>% lazy_dt() %>% 
  left_join(metadata.df, 
            by=c("system:time_start", "dateTime", "date", "WRS_ROW", "WRS_PATH", "Landsat")) %>% 
  mutate(year=year(date))

data.join <- as.data.table(data.join)[, year := year(date)
               ][,sat := factor(Landsat, levels = c(5,7,8),
                      labels = c('l5','l7','l8'))]
## Channel
data.join.chan = channels.filter%>% lazy_dt() %>% 
  left_join(metadata.df, 
            by=c("system:time_start", "dateTime", "date", "WRS_ROW", "WRS_PATH", "Landsat")) %>% 
  mutate(year=year(date))
data.join.chan <- as.data.table(data.join.chan)[, year := year(date)
               ][,sat := factor(Landsat, levels = c(5,7,8),
                      labels = c('l5','l7','l8'))]

# 5.  Export joined data as an intermediary step if you don't want to run the above code every time
#setwd("C:/Users/whyana/OneDrive/DocumentsLaptop/001_GraduateSchool/Research/Connectivity/Mackenzie/Data/intermediaryDownloads")
#data.join.save = data.join[,dom_wv_mean := chroma(Red_mean, Green_mean, Blue_mean)]
#write_feather(data.join.save, 'srMunged_Original_mackLakes_20230324.feather')
#data.join.save = read_feather('srMunged_Original_mackLakes_20230324.feather')

# 6. Define functions to compare all observations to landsat 7
## Lake functions
lm5 <- function(band){
  y <- data.join %>% 
    filter(year > 1999, year < 2012, sat == 'l7') %>% 
    .[,..band]
   
  y = quantile(y[[1]], seq(.01,.99, .01))
  x = data.join %>% 
    filter(year > 1999, year < 2012, sat == 'l5') %>% 
    .[,..band]
  x = quantile(x[[1]], seq(.01,.99, .01))
  
  lm <- lm(y~stats::poly(x, 2, raw = T))
  
  df <- tibble(band = band, intercept = lm$coefficients[[1]], B1 = lm$coefficients[[2]], B2 = lm$coefficients[[3]])
  return(df)
}
lm8 <- function(band){
  y <- data.join %>% 
    filter(year > 2012, sat == 'l7') %>% 
    .[,..band] 
    y = quantile(y[[1]], seq(.01,.99, .01))
  
  x = data.join %>% 
    filter(year > 2012, sat == 'l8') %>% 
    .[,..band] 
  x = quantile(x[[1]], seq(.01,.99, .01))
  
  lm <- lm(y~stats::poly(x, 2, raw = T))
  
  df <- tibble(band = band, intercept = lm$coefficients[[1]], B1 = lm$coefficients[[2]], B2 = lm$coefficients[[3]])
  return(df)
}

## Channel functions
lm5_chan <- function(band){
  y <- data.join.chan %>% 
    filter(year > 1999, year < 2012, sat == 'l7') %>% 
    .[,..band]
   
  y = quantile(y[[1]], seq(.01,.99, .01))
  x = data.join.chan %>% 
    filter(year > 1999, year < 2012, sat == 'l5') %>% 
    .[,..band]
  x = quantile(x[[1]], seq(.01,.99, .01))
  
  lm <- lm(y~stats::poly(x, 2, raw = T))
  
  df <- tibble(band = band, intercept = lm$coefficients[[1]], B1 = lm$coefficients[[2]], B2 = lm$coefficients[[3]])
  return(df)
}
lm8_chan <- function(band){
  y <- data.join.chan %>% 
    filter(year > 2012, sat == 'l7') %>% 
    .[,..band] 
    y = quantile(y[[1]], seq(.01,.99, .01))
  
  x = data.join.chan %>% 
    filter(year > 2012, sat == 'l8') %>% 
    .[,..band] 
  x = quantile(x[[1]], seq(.01,.99, .01))
  
  lm <- lm(y~stats::poly(x, 2, raw = T))
  
  df <- tibble(band = band, intercept = lm$coefficients[[1]], B1 = lm$coefficients[[2]], B2 = lm$coefficients[[3]])
  return(df)
}



# 7. Apply correction functions to get coefficients for the corrections
band <-  c('Blue_mean', 'Green_mean',  'Red_mean', 'Nir_mean', 'Swir1_mean', 'Swir2_mean')
funcs.5 <- band %>% map_dfr(lm5)
funcs.8 = band %>% map_dfr(lm8)
band_chan <-  c('Blue_mean', 'Green_mean',  'Red_mean', 'Nir_mean')
funcs.5.chan <- band_chan %>% map_dfr(lm5_chan)
funcs.8.chan = band_chan %>% map_dfr(lm8_chan)

# 8. Use the generated coefficients to correct the reflectance data
## Lakes
### Landsat 5 data
l5corr = data.join[sat == 'l5'][,Blue_mean := funcs.5[1,2][[1]] +
                                      funcs.5[1,3][[1]]*Blue_mean +
                                      funcs.5[1,4][[1]]*Blue_mean^2
                       ][,Green_mean := funcs.5[2,2][[1]] +
                           funcs.5[2,3][[1]]*Green_mean + 
                           funcs.5[2,4][[1]]*Green_mean^2
                         ][,Red_mean := funcs.5[3,2][[1]] + 
                             funcs.5[3,3][[1]]*Red_mean + 
                             funcs.5[3,4][[1]]*Red_mean^2
                          ][,Nir_mean := funcs.5[4,2][[1]] + 
                             funcs.5[4,3][[1]]*Nir_mean + 
                             funcs.5[4,4][[1]]*Nir_mean^2
                          ][,Swir1_mean := funcs.5[5,2][[1]] + 
                             funcs.5[5,3][[1]]*Swir1_mean + 
                             funcs.5[5,4][[1]]*Swir1_mean^2
                          ][,Swir2_mean := funcs.5[6,2][[1]] + 
                             funcs.5[6,3][[1]]*Swir2_mean + 
                             funcs.5[6,4][[1]]*Swir2_mean^2
                          ]

### Landsat 8 data
l8corr = data.join[sat == 'l8'][,Blue_mean := funcs.8[1,2][[1]] +
                                      funcs.8[1,3][[1]]*Blue_mean +
                                      funcs.8[1,4][[1]]*Blue_mean^2
                       ][,Green_mean := funcs.8[2,2][[1]] +
                           funcs.8[2,3][[1]]*Green_mean + 
                           funcs.8[2,4][[1]]*Green_mean^2
                         ][,Red_mean := funcs.8[3,2][[1]] + 
                             funcs.8[3,3][[1]]*Red_mean + 
                             funcs.8[3,4][[1]]*Red_mean^2
                          ][,Nir_mean := funcs.8[4,2][[1]] + 
                             funcs.8[4,3][[1]]*Nir_mean + 
                             funcs.8[4,4][[1]]*Nir_mean^2
                          ][,Swir1_mean := funcs.8[5,2][[1]] + 
                             funcs.8[5,3][[1]]*Swir1_mean + 
                             funcs.8[5,4][[1]]*Swir1_mean^2
                          ][,Swir2_mean := funcs.8[6,2][[1]] + 
                             funcs.8[6,3][[1]]*Swir2_mean + 
                             funcs.8[6,4][[1]]*Swir2_mean^2
                          ]

srCor <- data.join %>%
  filter(sat == 'l7') %>%
  bind_rows(l5corr) %>%
  bind_rows(l8corr) 
rm(l5corr, l8corr)
srCor <- as.data.table(srCor)[,dom_wv_mean := chroma(Red_mean, Green_mean, Blue_mean)]

setwd("C:/Users/whyana/OneDrive/DocumentsLaptop/001_GraduateSchool/Research/Connectivity/Mackenzie/Data/intermediaryDownloads")
write_feather(srCor, 'srCorrected_mackLakes_20230324.feather')
rm(all.lakes) 
gc()

## Channel data
### Landsat 5
l5corr.chan = data.join.chan[sat == 'l5'][,Blue_mean := funcs.5.chan[1,2][[1]] +
                                      funcs.5.chan[1,3][[1]]*Blue_mean +
                                      funcs.5.chan[1,4][[1]]*Blue_mean^2
                       ][,Green_mean := funcs.5.chan[2,2][[1]] +
                           funcs.5.chan[2,3][[1]]*Green_mean + 
                           funcs.5.chan[2,4][[1]]*Green_mean^2
                         ][,Red_mean := funcs.5.chan[3,2][[1]] + 
                             funcs.5.chan[3,3][[1]]*Red_mean + 
                             funcs.5.chan[3,4][[1]]*Red_mean^2
                          ][,Nir_mean := funcs.5.chan[4,2][[1]] + 
                             funcs.5.chan[4,3][[1]]*Nir_mean + 
                             funcs.5.chan[4,4][[1]]*Nir_mean^2
                          ]

### Landsat 8
l8corr.chan = data.join.chan[sat == 'l8'][,Blue_mean := funcs.8.chan[1,2][[1]] +
                                      funcs.8.chan[1,3][[1]]*Blue_mean +
                                      funcs.8.chan[1,4][[1]]*Blue_mean^2
                       ][,Green_mean := funcs.8.chan[2,2][[1]] +
                           funcs.8.chan[2,3][[1]]*Green_mean + 
                           funcs.8.chan[2,4][[1]]*Green_mean^2
                         ][,Red_mean := funcs.8.chan[3,2][[1]] + 
                             funcs.8.chan[3,3][[1]]*Red_mean + 
                             funcs.8.chan[3,4][[1]]*Red_mean^2
                          ][,Nir_mean := funcs.8.chan[4,2][[1]] + 
                             funcs.8.chan[4,3][[1]]*Nir_mean + 
                             funcs.8.chan[4,4][[1]]*Nir_mean^2
                          ]
srCor.chan <- data.join.chan %>%
  filter(sat == 'l7') %>%
  bind_rows(l5corr.chan) %>%
  bind_rows(l8corr.chan) 
rm(l5corr.chan, l8corr.chan)
setwd("C:/Users/whyana/OneDrive/DocumentsLaptop/001_GraduateSchool/Research/Connectivity/Mackenzie/Data/intermediaryDownloads")
write_feather(srCor.chan, 'srCorrected_mackChans_20230424.feather')


```

# correction plotted results for lakes  - Supplemental Figure 3
```{r}
setwd("C:/Users/whyana/OneDrive/DocumentsLaptop/001_GraduateSchool/Research/Connectivity/Mackenzie/Data/intermediaryDownloads")
srCor=read_feather('srCorrected_mackLakes_20230324.feather')
srCor <- as.data.table(srCor)
refCompPre <- as.data.table(data.join)
#ls5
p1 <- correctionPlot('Blue_mean', 'l5', refCompPre, srCor)
p2 <- correctionPlot('Green_mean', 'l5', refCompPre, srCor)
p3 <- correctionPlot('Red_mean', 'l5', refCompPre, srCor)
p4 <- correctionPlot('Nir_mean', 'l5', refCompPre, srCor)
p5 <- correctionPlot('Swir1_mean', 'l5', refCompPre, srCor)
p6 <- correctionPlot('Swir2_mean', 'l5', refCompPre, srCor)
#ls8
p7 <- correctionPlot('Blue_mean', 'l8', refCompPre, srCor)
p8 <- correctionPlot('Green_mean', 'l8', refCompPre, srCor)
p9 <- correctionPlot('Red_mean', 'l8', refCompPre, srCor)
p10 <- correctionPlot('Nir_mean', 'l8', refCompPre, srCor)
p11 <- correctionPlot('Swir1_mean', 'l8', refCompPre, srCor)
p12 <- correctionPlot('Swir2_mean', 'l8', refCompPre, srCor)
# combine plots
g1 <- ggarrange(p1, p2, p3, p4,p5,p6, common.legend = T) 
g2 <- ggarrange(p7,p8,p9,p10,p11, p12,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)       

setwd("C:/Users/whyana/OneDrive/DocumentsLaptop/001_GraduateSchool/Research/Connectivity/Mackenzie/images")
ggsave('srCorrections_20230423.pdf', plot = g, width = 8, height = 10, units = 'in')
rm(g1, g2, g11, g22, g)
```