R/lsat_clean_data.R

In logan-berner/lsatTS: An R package to facilitate retrieval, cleaning, cross-calibration, and phenological modeling of Landsat time-series data

#' Clean Landsat surface reflectance data
#'
#' @description
#' This function enables users to filter out surface reflectance 
#' measurements that exhibit:
#' (1) clouds, cloud shadows, snow, or water flagged by the CFMask algorithm;
#' (2) surface water over the Landsat record;  
#' (3) impossibly high reflectance (>1.0) and abnormally low reflectance (<0.005);
#' (4) scene cloud cover above a user-defined threshold;
#' (5) geometric uncertainty above a user-defined threshold;
#' (6) solar zenith angle above a user-defined threshold.
#' 
#' @param dt Data.table generated by calling lsat_format_data().
#' @param cloud.max Maximum allowable cloud cover in Landsat scene (percentage).
#' @param geom.max Maximum allowable geometric uncertainty (meters).
#' @param sza.max Maximum allowable solar zenith angle (degrees).
#' @param filter.cfmask.snow (TRUE/FALSE) Remove measurements with CFmask flag = snow.
#' @param filter.cfmask.water (TRUE/FALSE) Remove measurements with CFmask flag = water.
#' @param filter.jrc.water (TRUE/FALSE) Remove sample sites that were ever inundated 
#'    based on the maximum surface water extent variable from the 
#'    JRC Global Surface Water Dataset.
#' @return A data.table that includes Landsat measurements that met the quality control criteria.
#' @import data.table
#' @export lsat_clean_data
#'
#' @examples 
#' data(lsat.example.dt)
#' lsat.dt <- lsat_format_data(lsat.example.dt)
#' lsat.dt <- lsat_clean_data(lsat.dt)
#' lsat.dt

lsat_clean_data <- function(dt, 
                            cloud.max=80, 
                            geom.max=30, 
                            sza.max=60, 
                            filter.cfmask.snow = T, 
                            filter.cfmask.water = T, 
                            filter.jrc.water = T){
  
  dt <- data.table::data.table(dt)
  n.orig <- nrow(dt)

  # pixel flags for clear sky
  dt[, clear := mapply(clear_value, qa.pixel)]
  dt <- dt[clear == 1]

  # pixel flags for snow
  if (filter.cfmask.snow == T){
    dt[, snow := mapply(snow_flag, qa.pixel)]
    dt <- dt[snow == 0]
  }

  # pixel flags for water and JRC Max Water Extent
  if (filter.cfmask.water == T){
    dt[, water := mapply(water_flag, qa.pixel)]
    dt <- dt[water == 0]
  }
  
  if (filter.jrc.water == T){
    dt[, jrc.water := as.numeric(jrc.water)]
    dt <- dt[jrc.water == 0]
  }

  # scene flags
  dt <- dt[cloud.cover <= cloud.max]
  dt <- dt[geometric.rmse.model <= geom.max]
  dt <- dt[90-sun.elevation <= sza.max]
  dt <- dt[qa.radsat == 0]

  # filter out unrealistic band values
  dt <- dt[blue > 0.005][green > 0.005][red > 0.005][nir > 0.005]
  dt <- dt[blue < 1][green < 1][red < 1][nir < 1]

  n.final <- nrow(dt)
  n.removed <- n.orig - n.final
  print.msg <- paste0('removed ', n.removed, ' of ', n.orig, ' observations (',
                      round(n.removed / n.orig * 100, 2), '%)')
  print(print.msg)
  dt
}

# LANDSAT CLEAR SKY BIT VALUES (returns 1 if clear and 0 if not clear)
clear_value = function(x) {
  # reverse order of bits so read from left to right
  bit_str = paste(as.integer(intToBits(x)), collapse="")
  # conditions
  filled = substr(bit_str, 1, 1) == '1'
  cloud_shadow = substr(bit_str, 5, 5) == '1'
  not_clear = substr(bit_str, 7, 7) == '0'

  if(filled | cloud_shadow | not_clear){
    return(0)
  }  else{
    return(1)
  }
}

# filter snow
snow_flag = function(x) {
  # reverse order of bits, left to right
  bit_str = paste(as.integer(intToBits(x)), collapse="")
  snow = substr(bit_str, 6, 6) == '1'
  if(snow){return(1)}  else{return(0)}
}

# filter water
water_flag = function(x) {
  # reverse order of bits, left to right
  bit_str = paste(as.integer(intToBits(x)), collapse="")
  water = substr(bit_str, 8, 8) == '1'
  if(water){return(1)}  else{return(0)}
}