Methodology

###############################################################################################
# Daily temporal combination 
###############################################################################################

# Function to replace all cloud pixels, missing pixels, saturated, fill and no decision pixels to NA in order to cover them 

NA_fill <- function(images){
  for(i in 1:length(images)){
    images[[i]][images[[i]]@data@values == 250] <- NA
    images[[i]][images[[i]]@data@values == 200] <- NA
    images[[i]][images[[i]]@data@values == 201] <- NA
    images[[i]][images[[i]]@data@values == 254] <- NA
    images[[i]][images[[i]]@data@values == 255] <- NA
  }
  return(images)
}

# Function to go through all sets of data and cover all missing pixels possible 
# This function will replace Terra NA pixels with NA from VIIRS first, and then from Aqua

Replace <- function(TERRA, VIIRS, AQUA){
  month_combined <- list()
  for(i in 1:length(TERRA)){
    month_combined[[i]] <- cover(TERRA[[i]], VIIRS[[i]], AQUA[[i]])
  }
  return(month_combined)
}

# This function is the same as the previous one
# However, it is used for the years prior to the availability of VIIRS (prior to 2012)
# Only uses Terra and Aqua

Replace1 <- function(TERRA, AQUA){
  month_combined <- list()
  for(i in 1:length(TERRA)){
    month_combined[[i]] <- cover(TERRA[[i]], AQUA[[i]])
  }
  return(month_combined)
}

# function to add back cloud pixels, note that after this function you must crop function again
# This function is used in order to calculate the remaining cloud pixels after each step in the methodology is implemeneted

cloud_add <- function(na_images){
  for(i in 1:length(na_images)){
    na_images[[i]][is.na(na_images[[i]]@data@values)] <- 2
  }
  return(na_images)
}

###############################################################################################
# Adjecent temporal deduction 
###############################################################################################

# This function takes the combined images and creates an index for i+1, i+2, i-1, i-2
# Note there are 2 functions here, one for a +/-1 day ATD and one for a +/-2 days
# We need both of these functions in order to calcualte the cloud pixels removed from an additional day used in the ATD
# then using the cover function, will replace the NA pixels for one day before and one day after first
# then they take the day two days ahead and two days after, respectively 

# +/- 1 Day
temporal_filter <- function(combined_images){
  temporal_combine <- list()
  a2 <- list()
  a3 <- list()
  a4 <- list()
  a5 <- list()
  for(i in 3:93){
    a2[[i]] <- combined_images[[i+1]]
    a3[[i]] <- combined_images[[i-1]]
    temporal_combine[[i]] <- cover(combined_images[[i]], a3[[i]], a2[[i]], tolerance = 0.05, overlap = T)
  }
  return(temporal_combine)
}

# +/= 2 Days
temporal_filter_2 <- function(combined_images){
  temporal_combine <- list()
  a2 <- list()
  a3 <- list()
  a4 <- list()
  a5 <- list()
  for(i in 3:93){
    a2[[i]] <- combined_images[[i+1]]
    a3[[i]] <- combined_images[[i-1]]
    a4[[i]] <- combined_images[[i+2]]
    a5[[i]] <- combined_images[[i-2]]
    temporal_combine[[i]] <- cover(combined_images[[i]], a3[[i]], a2[[i]], a5[[i]], a4[[i]], tolerance = 0.05, overlap = T)
  }
  return(temporal_combine)
}

###############################################################################################
# Going through all years now with parts of algorithm applicable
###############################################################################################

# From 2003-2011 only Terra and Aqua are used, so we use Replace1 function on NA filled images

Year_2003 <- Replace1(TERRA = NA_fill(Terra_2003), AQUA = NA_fill(Aqua_2003))
Year_2004 <- Replace1(TERRA = NA_fill(Terra_2004), AQUA = NA_fill(Aqua_2004))
Year_2005 <- Replace1(TERRA = NA_fill(Terra_2005), AQUA = NA_fill(Aqua_2005))
Year_2006 <- Replace1(TERRA = NA_fill(Terra_2006), AQUA = NA_fill(Aqua_2006))
Year_2007 <- Replace1(TERRA = NA_fill(Terra_2007), AQUA = NA_fill(Aqua_2007))
Year_2008 <- Replace1(TERRA = NA_fill(Terra_2008), AQUA = NA_fill(Aqua_2008))
Year_2009 <- Replace1(TERRA = NA_fill(Terra_2009), AQUA = NA_fill(Aqua_2009))
Year_2010 <- Replace1(TERRA = NA_fill(Terra_2010), AQUA = NA_fill(Aqua_2010))
Year_2011 <- Replace1(TERRA = NA_fill(Terra_2011), AQUA = NA_fill(Aqua_2011))

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

# For the remainder of the images, since we want to know how much additional cloud is removed with VIIRS
# TAC will be done from 2012-2019 as well 

# TAC 
Year_2012_TAC <- Replace1(TERRA = NA_fill(Terra_2012), AQUA = NA_fill(Aqua_2012))
Year_2013_TAC <- Replace1(TERRA = NA_fill(Terra_2013), AQUA = NA_fill(Aqua_2013))
Year_2014_TAC <- Replace1(TERRA = NA_fill(Terra_2014), AQUA = NA_fill(Aqua_2014))
Year_2015_TAC <- Replace1(TERRA = NA_fill(Terra_2015), AQUA = NA_fill(Aqua_2015))
Year_2016_TAC <- Replace1(TERRA = NA_fill(Terra_2016), AQUA = NA_fill(Aqua_2016))
Year_2017_TAC <- Replace1(TERRA = NA_fill(Terra_2017), AQUA = NA_fill(Aqua_2017))
Year_2018_TAC <- Replace1(TERRA = NA_fill(Terra_2018), AQUA = NA_fill(Aqua_2018))
Year_2019_TAC <- Replace1(TERRA = NA_fill(Terra_2019), AQUA = NA_fill(Aqua_2019))

# Combined 3
# Now we combine terra, aqua and viirs into one image composite using the Replace function 

Year_2012 <- Replace(TERRA = NA_fill(Terra_2012), VIIRS = NA_fill(VIIRS_2012), AQUA = NA_fill(Aqua_2012))
Year_2013 <- Replace(TERRA = NA_fill(Terra_2013), VIIRS = NA_fill(VIIRS_2013), AQUA = NA_fill(Aqua_2013))
Year_2014 <- Replace(TERRA = NA_fill(Terra_2014), VIIRS = NA_fill(VIIRS_2014), AQUA = NA_fill(Aqua_2014))
Year_2015 <- Replace(TERRA = NA_fill(Terra_2015), VIIRS = NA_fill(VIIRS_2015), AQUA = NA_fill(Aqua_2015))
Year_2016 <- Replace(TERRA = NA_fill(Terra_2016), VIIRS = NA_fill(VIIRS_2016), AQUA = NA_fill(Aqua_2016))
Year_2017 <- Replace(TERRA = NA_fill(Terra_2017), VIIRS = NA_fill(VIIRS_2017), AQUA = NA_fill(Aqua_2017))
Year_2018 <- Replace(TERRA = NA_fill(Terra_2018), VIIRS = NA_fill(VIIRS_2018), AQUA = NA_fill(Aqua_2018))
Year_2019 <- Replace(TERRA = NA_fill(Terra_2019), VIIRS = NA_fill(VIIRS_2019), AQUA = NA_fill(Aqua_2019))

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

# Applying the +/- 1 day ATD frmom 2003-2019

ATD_2003 <- temporal_filter(Year_2003)
ATD_2004 <- temporal_filter(Year_2004)
ATD_2005 <- temporal_filter(Year_2005)
ATD_2006 <- temporal_filter(Year_2006)
ATD_2007 <- temporal_filter(Year_2007)
ATD_2008 <- temporal_filter(Year_2008)
ATD_2009 <- temporal_filter(Year_2009)
ATD_2010 <- temporal_filter(Year_2010)
ATD_2011 <- temporal_filter(Year_2011)

ATD_2012 <- temporal_filter(Year_2012)
ATD_2013 <- temporal_filter(Year_2013)
ATD_2014 <- temporal_filter(Year_2014)
ATD_2015 <- temporal_filter(Year_2015)
ATD_2016 <- temporal_filter(Year_2016)
ATD_2017 <- temporal_filter(Year_2017)
ATD_2018 <- temporal_filter(Year_2018)
ATD_2019 <- temporal_filter(Year_2019)

# Applying ATD +/- 2 days from 2003-2019

ATD_2003_2 <- temporal_filter_2(Year_2003)
ATD_2004_2 <- temporal_filter_2(Year_2004)
ATD_2005_2 <- temporal_filter_2(Year_2005)
ATD_2006_2 <- temporal_filter_2(Year_2006)
ATD_2007_2 <- temporal_filter_2(Year_2007)
ATD_2008_2 <- temporal_filter_2(Year_2008)
ATD_2009_2 <- temporal_filter_2(Year_2009)
ATD_2010_2 <- temporal_filter_2(Year_2010)
ATD_2011_2 <- temporal_filter_2(Year_2011)

ATD_2012_2 <- temporal_filter_2(Year_2012)
ATD_2013_2 <- temporal_filter_2(Year_2013)
ATD_2014_2 <- temporal_filter_2(Year_2014)
ATD_2015_2 <- temporal_filter_2(Year_2015)
ATD_2016_2 <- temporal_filter_2(Year_2016)
ATD_2017_2 <- temporal_filter_2(Year_2017)
ATD_2018_2 <- temporal_filter_2(Year_2018)
ATD_2019_2 <- temporal_filter_2(Year_2019)

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

# this function can be used to convert raster values to a tibble for each coordinate set 

daily_tibbles <- function(data_to_tibble){
  data <- lapply(data_to_tibble, as.data.frame,  xy = T)
  daily_tibbles[[i]]<- na.omit(as_tibble(data))
  return(daily_tibbles)
}

sorting_files <- function(MODIS_list){
  for(i in 1:187){
    MODIS_sorted <- sort(MODIS_list[[i]]@data@names, decreasing = F)
  }
  return(MODIS_sorted)
}

# Since there are days before and after the study period used for the ATD
# these NULL images must be removed as they no longer  

ATD_2003 <- ATD_2003[-c(1,2)]
ATD_2004 <- ATD_2004[-c(1,2)]
ATD_2005 <- ATD_2005[-c(1,2)]
ATD_2006 <- ATD_2006[-c(1,2)]
ATD_2007 <- ATD_2007[-c(1,2)]
ATD_2008 <- ATD_2008[-c(1,2)]
ATD_2009 <- ATD_2009[-c(1,2)]
ATD_2010 <- ATD_2010[-c(1,2)]
ATD_2011 <- ATD_2011[-c(1,2)]

ATD_2012 <- ATD_2012[-c(1,2)]
ATD_2013 <- ATD_2013[-c(1,2)]
ATD_2014 <- ATD_2014[-c(1,2)]
ATD_2015 <- ATD_2015[-c(1,2)]
ATD_2016 <- ATD_2016[-c(1,2)]
ATD_2017 <- ATD_2017[-c(1,2)]
ATD_2018 <- ATD_2018[-c(1,2)]
ATD_2019 <- ATD_2019[-c(1,2)]

ATD_2003_2 <- ATD_2003_2[-c(1,2)]
ATD_2004_2 <- ATD_2004_2[-c(1,2)]
ATD_2005_2 <- ATD_2005_2[-c(1,2)]
ATD_2006_2 <- ATD_2006_2[-c(1,2)]
ATD_2007_2 <- ATD_2007_2[-c(1,2)]
ATD_2008_2 <- ATD_2008_2[-c(1,2)]
ATD_2009_2 <- ATD_2009_2[-c(1,2)]
ATD_2010_2 <- ATD_2010_2[-c(1,2)]
ATD_2011_2 <- ATD_2011_2[-c(1,2)]

ATD_2012_2 <- ATD_2012_2[-c(1,2)]
ATD_2013_2 <- ATD_2013_2[-c(1,2)]
ATD_2014_2 <- ATD_2014_2[-c(1,2)]
ATD_2015_2 <- ATD_2015_2[-c(1,2)]
ATD_2016_2 <- ATD_2016_2[-c(1,2)]
ATD_2017_2 <- ATD_2017_2[-c(1,2)]
ATD_2018_2 <- ATD_2018_2[-c(1,2)]
ATD_2019_2 <- ATD_2019_2[-c(1,2)]


###############################################################################################
# 3x3 kernal focal filter 
###############################################################################################

# This function will reclassify rasters to make all cloud pixels == 1 
# this is how we reclassify rasters to create cloud duration maps

reclass_binary_c <- function(raster_list){
  reclassify(raster_list, c(0,240,0, 230, 255, 1), na.omit = T)
}

# This function will reclassify all pixels that are 50% NDSI and above == snow 
# Creating a binary snow map

reclass_binary <- function(raster_list){
  reclassify(raster_list, c(0, 49, 0, 49,100,1, 100, 255, 0), na.omit = T)
}

# this function is used for the focal filter
# creating a neighborhood of 8 pixels surround the pixel needing interpolation 
# it also specifies that the means for interpolation will be the mean of the surrounding points

fill.na <- function(x, i=8) {
  if( is.na(x)[i] ) {
    return( round(mean(x, na.rm=TRUE),0) )
  } 
}  

# The focal filter, specifies a matrix of 3x3 kernel for interpolating NA pixels only 

filter_focal <- function(ras){
  for(i in 1:length(ras)){
    ras[[i]] <- focal(ras[[i]], w = matrix(1, nrow = 3, ncol = 3), na.rm = F, pad = T, fun = fill.na, NAonly = TRUE)
  }
  return(ras)
}

# Applying the spatial filter for each annual dataset
# Using the ATD +/- 2 days as the least cloud containing dataset

SF_2003 <- filter_focal(ATD_2003_2)
SF_2004 <- filter_focal(ATD_2004_2)
SF_2005 <- filter_focal(ATD_2005_2)
SF_2006 <- filter_focal(ATD_2006_2)
SF_2007 <- filter_focal(ATD_2007_2)
SF_2008 <- filter_focal(ATD_2008_2)
SF_2009 <- filter_focal(ATD_2009_2)
SF_2010 <- filter_focal(ATD_2010_2)
SF_2011 <- filter_focal(ATD_2011_2)
SF_2012 <- filter_focal(ATD_2012_2)
SF_2013 <- filter_focal(ATD_2013_2)
SF_2014 <- filter_focal(ATD_2014_2)
SF_2015 <- filter_focal(ATD_2015_2)
SF_2016 <- filter_focal(ATD_2016_2)
SF_2017 <- filter_focal(ATD_2017_2)
SF_2018 <- filter_focal(ATD_2018_2)
SF_2019 <- filter_focal(ATD_2019_2)

###############################################################################################
# Calculating regional snow line and regional landline and filling in NA/ cloud pixels that are
# above and below this threshold accordingly
###############################################################################################
# This code block is taken from Parajka 2010 for regional snowline calculation

SNOWL <- function(DEM, SF_images, shp){
  demmin<-minValue(DEM)
  demmax<-maxValue(DEM)
  
  rsle<--999
  rsl<-10000000000
  
  nc<-freq(mask(SF_images,shp),value=NA) # frequency of cloud pixels
  ns<-freq(mask(SF_images,shp),value=1) #frequency of snow pixels
  nl<-freq(mask(SF_images,shp),value=0) #frequency of land pixels
  
  for (elev in seq(demmin, demmax, by=50)){
    fun <- function(x) { x[x>elev] <- NA; return(x) }
    tmp <- reclassify(calc(DEM, fun),c(0,demmax,1))*SF_images
    
    nss<-freq(tmp,value=1) #frequency of snow pixels
    nll<-nl-freq(tmp,value=0) #frequency of land pixels
    
    if(nll+nss <= rsl) {
      rsl<-nll+nss
      rsle<-elev
    }
  }
  
  
  rlle<--999
  rll<-10000000000
  
  
  for (elev in seq(demmin, demmax, by=50)){
    fun <- function(x) { x[x>elev] <- NA; return(x) }
    tmp <- reclassify(calc(DEM, fun),c(0,demmax,1))*SF_images
    
    nss<-freq(tmp,value=1) #frequency of snow pixels
    nll<-nl-freq(tmp,value=0) #frequency of land pixels
    
    if(nll+nss <= rll) {
      rll<-nll-nss
      rlle<-elev
    }
    
  }
  
  image_1 <- raster::overlay(DEM, SF_images, fun= function(e, s){s[e>= rsle & is.na(s)] <- 1; s})
  image_2 <- raster::overlay(DEM, image_1, fun= function(e, s){s[e<= rlle & is.na(s)] <- 0; s})
  return(image_2)
}

# This SNOWL method is a modification to Parajka 2010 code to calcualtion the regional snowline 
# the vectors return the mean elevation values for snow (1) and no snow (0)
# Then the snow and no snow at the bottom will change the values based on this mean 
# first for snow, then no snow on the same image and return the raster

SNOWL <- function(image, DEM){
  vect <- function(e,s){
    mean(e[s==1], na.rm = T)
  }
  vect1 <- function(e,s){
    mean(e[s==0], na.rm = T)
  }
  t <- raster::overlay(DEM, image, fun = Vectorize(vect))
  t2 <- raster::overlay(DEM, image, fun = Vectorize(vect1))
  
  snow <- raster::overlay(DEM, image, fun= function(e, s){s[e>= mean(t@data@values, na.rm = T) & is.na(s)] <- 1; s})
  no_snow <- raster::overlay(DEM, snow, fun= function(e, s){s[e>= mean(t2@data@values, na.rm = T) & is.na(s)] <- 0; s})
  
  return(no_snow)
}

# this function will output the mean values of elevation for both the landlines and the snowlines daily 

SNOWL_lines <- function(image, DEM){
  vect <- function(e,s){
    mean(e[s==1], na.rm = T)
  }
  vect1 <- function(e,s){
    mean(e[s==0], na.rm = T)
  }
  t <- raster::overlay(DEM, image, fun = Vectorize(vect))
  t2 <- raster::overlay(DEM, image, fun = Vectorize(vect1))
  
  snowline <- mean(t@data@values, na.rm = T)
  landline <- mean(t2@data@values, na.rm = T)
  
  lines <- c(snowline, landline)%>%
    as.table()
  return(lines)
}


# Then we apply this SNOWL method to each of the spatially filtered images
# Note they are reclassified this is the reason for the _c beside each set 

SL_2003 <- lapply(SF_2003_c, SNOWL, DEM=DEM)
SL_2004 <- lapply(SF_2004_c, SNOWL, DEM=DEM)
SL_2005 <- lapply(SF_2005_c, SNOWL, DEM=DEM)
SL_2006 <- lapply(SF_2006_c, SNOWL, DEM=DEM)
SL_2007 <- lapply(SF_2007_c, SNOWL, DEM=DEM)
SL_2008 <- lapply(SF_2008_c, SNOWL, DEM=DEM)
SL_2009 <- lapply(SF_2009_c, SNOWL, DEM=DEM)
SL_2010 <- lapply(SF_2010_c, SNOWL, DEM=DEM)
SL_2011 <- lapply(SF_2011_c, SNOWL, DEM=DEM)
SL_2012 <- lapply(SF_2012_c, SNOWL, DEM=DEM)
SL_2013 <- lapply(SF_2013_c, SNOWL, DEM=DEM)
SL_2014 <- lapply(SF_2014_c, SNOWL, DEM=DEM)
SL_2015 <- lapply(SF_2015_c, SNOWL, DEM=DEM)
SL_2016 <- lapply(SF_2016_c, SNOWL, DEM=DEM)
SL_2017 <- lapply(SF_2017_c, SNOWL, DEM=DEM)
SL_2018 <- lapply(SF_2018_c, SNOWL, DEM=DEM)
SL_2019 <- lapply(SF_2019_c, SNOWL, DEM=DEM)

# The function to output each pair of mean elevation values for land and snow lines 
# Note this is for south thompson, which is why the values of the DEM and crop function are different 

line_2003 <- lapply(SF_2003_c, crop_fun2, thomfile=south_thom)%>%lapply(SNOWL_lines, DEM = DEM_ST)
line_2004 <- lapply(SF_2004_c, crop_fun2, thomfile=south_thom)%>%lapply(SNOWL_lines, DEM = DEM_ST)
line_2005 <- lapply(SF_2005_c, crop_fun2, thomfile=south_thom)%>%lapply(SNOWL_lines, DEM = DEM_ST)
line_2006 <- lapply(SF_2006_c, crop_fun2, thomfile=south_thom)%>%lapply(SNOWL_lines, DEM = DEM_ST)
line_2007 <- lapply(SF_2007_c, crop_fun2, thomfile=south_thom)%>%lapply(SNOWL_lines, DEM = DEM_ST)
line_2008 <- lapply(SF_2008_c, crop_fun2, thomfile=south_thom)%>%lapply(SNOWL_lines, DEM = DEM_ST)
line_2009 <- lapply(SF_2009_c, crop_fun2, thomfile=south_thom)%>%lapply(SNOWL_lines, DEM = DEM_ST)
line_2010 <- lapply(SF_2010_c, crop_fun2, thomfile=south_thom)%>%lapply(SNOWL_lines, DEM = DEM_ST)
line_2011 <- lapply(SF_2011_c, crop_fun2, thomfile=south_thom)%>%lapply(SNOWL_lines, DEM = DEM_ST)
line_2012 <- lapply(SF_2012_c, crop_fun2, thomfile=south_thom)%>%lapply(SNOWL_lines, DEM = DEM_ST)
line_2013 <- lapply(SF_2013_c, crop_fun2, thomfile=south_thom)%>%lapply(SNOWL_lines, DEM = DEM_ST)
line_2014 <- lapply(SF_2014_c, crop_fun2, thomfile=south_thom)%>%lapply(SNOWL_lines, DEM = DEM_ST)
line_2015 <- lapply(SF_2015_c, crop_fun2, thomfile=south_thom)%>%lapply(SNOWL_lines, DEM = DEM_ST)
line_2016 <- lapply(SF_2016_c, crop_fun2, thomfile=south_thom)%>%lapply(SNOWL_lines, DEM = DEM_ST)
line_2017 <- lapply(SF_2017_c, crop_fun2, thomfile=south_thom)%>%lapply(SNOWL_lines, DEM = DEM_ST)
line_2018 <- lapply(SF_2018_c, crop_fun2, thomfile=south_thom)%>%lapply(SNOWL_lines, DEM = DEM_ST)
line_2019 <- lapply(SF_2019_c, crop_fun2, thomfile=south_thom)%>%lapply(SNOWL_lines, DEM = DEM_ST)

# This string of code will melt the lines to include both the snow nad land line and then seperate them via split
# Then we can write this dataframe to an excel file 

L <- melt(line_2019)
L_snow <- split(L, L$Var1 == "A")
df <- data.frame(snow = L_snow$'TRUE', land = L_snow$'FALSE')

write.csv(df, "lines_2003.csv")
###############################################################################################