RSFIA: FIA Tree Mortality Analysis

#' @include generics.R
#' @include nc_helper.R
#'
#' @title Class, methods, and functions for remote sensing objects
#'
#' @description
#'
#' Wrapper/processing functions for remotely-sensed data.
#'
#' @param object Object of class 'remote'
#' @param bands Integer vector of bands to include (in Landsat 8 form)
#'
#'
#' @details
#'
#' @author Brandon McNellis
#'
#' @name remote
#' @rdname remote
NULL
#'
#' An S4 class for analysis results
#'
#' @slot ft_dir Directory location of fusion table/GEE output files
#' @slot bands Bands from BandKey() to include, must be in GEE filename as band*
#'
#' @rdname remote
remote <- setClass(
  'remote',
  slots = list(
    bands = 'integer',
    ft_dir = 'character'
  ),
  contains = c('nc_helper')
)
#' @export
setValidity('remote', function(object) {
  errors <- character()

  # ft_dir
  if (!dir.exists(object@ft_dir)) {
    msg <- paste0('ft_dir does not exist')
    errors <- c(errors, msg)
  }

  # bands
  ba <- BandKey()[, 1]
  if (all(object@bands %in% ba)) {
    # do nothing
  } else {
    msg <- paste0('bands out of range')
    errors <- c(errors, msg)
  }

  # Returns
  if (length(errors) == 0) {
    TRUE
  } else {
    errors
  }
})
#' @rdname remote
#' @export
setMethod('initialize',
          signature(.Object = 'remote'),
          function (.Object, ...) {

            params <- list(...)

            # ft_dir
            if ('ft_dir' %in% names(params)) {
              .Object@ft_dir <- params$ft_dir
            } else {
              .Object@ft_dir <- getwd()
            }

            # bands
            if ('bands' %in% names(params)) {
              .Object@bands <- params$bands
            } else {
              .Object@bands <- seq(11)
            }

            # returns
            .Object <- callNextMethod()
            mt <- validObject(.Object)
            if (isTRUE(mt)) {
              return(.Object)
            } else {
              return(mt)
            }
          }
)
#' @rdname remote
#' @export
setMethod('SetupDataFile',
          signature(object = 'remote'),
          function(object, overwrite = F, backup = F) {

            # ft section:
            fname <- paste0(object@ft_dir, '/coordinates.csv')

            lon <- sapply(object@coords, function(x) x[1])
            lat <- sapply(object@coords, function(x) x[2])
            df <- data.frame(object@sample, lon, lat, stringsAsFactors = F)
            colnames(df) <- c('id', 'longitude', 'latitude')

            write.csv(df, fname, row.names = F)
            message('Wrote fusion table coordinates to:')
            cat(fname, '\n')

            callNextMethod()
          }
)
#' @rdname remote
#' @export
ReadEarthEngineOutputs <- function(object, satellite) {
  validObject(object)
  key <- object@sample
  time <- object@time
  bd <- object@bands
  bk <- BandKey()
  if (satellite %in% colnames(bk)) {
    #bd0 <- data.frame(bk$band, bk[[satellite]][match(bd, bk[, 1])], stringsAsFactors = F)
    bdind <- match(bd, bk[, 1])
    bd0 <- data.frame(bk$band[bdind], bk[[satellite]][bdind])
    bindx <- which(!(is.na(bd0[, 2])))
    bd1 <- paste0('band', bd0[bindx, 2])
    bd0[bindx, 'ftag'] <- bd1
    colnames(bd0) <- c('band', satellite, 'ftag')
  } else {
    stop('satellite name must be a column in BandKey()')
  }
  bands <- bd0$ftag

  fl_list <- list.files(path = object@ft_dir)
  nc_path <- paste0(object@nc_dir, '/', object@file_name)
  lon <- sapply(object@coords, function(x) x[1])
  lat <- sapply(object@coords, function(x) x[2])

  # Magic numbers for column #s that contain coordinates as well as measurement
  indM <- 2 # measurement
  indC <- 5 # coordinates

  if (!file.exists(nc_path)) {
    SetupDataFile(object)
  }

  nc <- ncdf4::nc_open(nc_path, write = T)

  yrs <- object@time
  sc <- 0L
  i0 <- 0L
  f0 <- length(yrs) * length(bands)
  for (i in seq_along(yrs)) {
    ii <- yrs[i]
    i_fls <- fl_list[grep(ii, fl_list)]

    for (j in seq_along(bands)) {
      jj <- bands[j]
      j_fls <- i_fls[grep(jj, i_fls)]

      if (length(j_fls) > 1) {
        stop('too many files matched')
      } else {
        if (length(j_fls) < 1) {
          sc <- sc + 1L
          if (sc == length(yrs)) {
            stop('didnt match any files')
          }
          next
        }
        jfls0 <- paste0(object@ft_dir, '/', j_fls)
        jfl <- read.csv(jfls0, stringsAsFactors = F)
        jfl <- jfl[, c(indC, indM)]

        for (k in seq(nrow(jfl))) {
          st <- regexpr('\\[', jfl[k, 1])[1]
          en <- regexpr('\\]', jfl[k, 1])[1]
          jfl[k, 1] <- substr(jfl[k, 1], st + 1, en - 1)
        }

        c0 <- data.frame(CoordVecsToList(strsplit(jfl[, 1], ',')), jfl[, 2])
        c1 <- data.frame(CoordVecsToList(object@coords), object@sample, stringsAsFactors = F)
        fdf <- dplyr::left_join(c1, c0, by = c('lon', 'lat'))
        bn <- sapply(strsplit(jj, 'band'), function(x) x[2])
        colnames(fdf) <- c('lon', 'lat', 'sample', bn)
        fdf <- fdf[order(fdf$sample), ]
        fdf$time <- rep(ii, nrow(fdf))
        if (all(fdf$sample == object@sample)) {
          # do nothing
        } else {
          stop('object sample doesnt match input sample')
        }
        fdf <- fdf[, c('time', 'sample', bn)]

        cat('\r', format(i0 / f0 * 100, digits = 2, nsmall = 2), '%')
        cat(', year =', ii, ', band =', jj, '     ')
        FillArray(object, df = fdf, nc = nc)
        i0 <- i0 + 1L

      } # end ifelse (=file length match check)
    } # end j (=band)
  } # end i (=time)

  return(object)
}
#' @rdname remote
#' @export
BandKey <- function() {
  # ref doi:10.3390/rs61010232
  key <- data.frame(matrix(nrow = 11, ncol = 0))
  key$band <- c(seq(11))
  key$`Landsat 8 Name` <- c('coastal aerosol', 'blue', 'green', 'red', 'NIR', 'SWIR1', 'SWIR2',
                            'panchromatic', 'cirrus', 'TIRS1', 'TIRS2')
  key$`Landsat 8` <- seq(11)
  key$`Landsat 7` <- c(2, 3, 4, 5, 6, 10, 7, 8, NA, NA, NA)

  key$`Landsat 5_TM` <- c(2, 3, 4, 5, 6, 10, 7, NA, NA, NA, NA)
  key$`Landsat 5_MSS` <- c(3, 4, NA, 3, 5, NA, NA, NA, NA, NA, NA)

  key$`Landsat 4_TM` <- c(2, 3, 4, 5, 6, 10, 7, NA, NA, NA, NA)
  key$`Landsat 4_MSS` <- c(3, 4, NA, 3, 5, NA, NA, NA, NA, NA, NA)

  key$`Landsat 3_MSS` <- c(NA, NA, NA, 3, 4, NA, 5, NA, NA, NA, NA)
  key$`Landsat 2_MSS` <- c(NA, NA, NA, 3, 4, NA, 5, NA, NA, NA, NA)
  key$`Landsat 1_MSS` <- c(NA, NA, NA, 3, 4, NA, 5, NA, NA, NA, NA)

  # Other indicies:
  key[nrow(key) + 1, ] <- c(101, 'NDVI', rep(101, ncol(key) - 2))
  key[nrow(key) + 1, ] <- c(102, 'NDMI', rep(102, ncol(key) - 2))
  key[nrow(key) + 1, ] <- c(103, 'RGI', rep(103, ncol(key) - 2))
  # MIR is SWIR #2, or 'longer' SWIR, from Meddens et al. (2012)
  key[nrow(key) + 1, ] <- c(104, 'MIR', 7, 7, rep(NA, ncol(key) - 4))

  key[, -2] <- data.frame(lapply(key[, -2], as.integer))

  return(key)
}
#' @rdname remote
#' @export
DistillRasterStack <- function(coords, geo_dir = getwd()) {
  # Input validity checks
  stopifnot(
    class(geo_dir) == 'character',
    length(geo_dir) == 1,
    class(coords) == 'data.frame',
    ncol(coords) == 2
  )
  # Check lat/long range
  if (any(
    range(coords[, 1])[1] < -90,
    range(coords[, 1])[2] > 90,
    range(coords[, 2])[1] < -180,
    range(coords[, 2])[2] > 180
  )) {
    stop('Coordinates out of range.')
  }
  # Find the files and check if they exist
  file_ls <- list.files(path = geo_dir, pattern = '.tif$')
  if (length(file_ls) < 1) stop('No .tif files found for import.')
  # Preallocation
  return_df <- data.frame(matrix(nrow = nrow(coords), ncol = 0))
  for (i in 1:length(file_ls)) {
    i_tif <- raster::raster(file_ls[i])
    i_pts <- raster::extract(i_tif, coords)
    xpix <- res(i_tif)[1]
    if (!is.na(res(i_tif)[2])) ypix <- res(i_tif)[2] else ypix <- xpix
    for (j in 1:nrow(coords)) {
      if (!is.na(i_pts[j])) next
      xinc <- xpix * 0.1
      yinc <- ypix * 0.1
      j_pt <- coords[j, ]
      if (any(
        j_pt[, 1] > i_tif@extent@ymax,
        j_pt[, 1] < i_tif@extent@ymin,
        j_pt[, 2] > i_tif@extent@xmax,
        j_pt[, 2] < i_tif@extent@xmin
      )) {
        cat('Coordinate', j, 'outside raster extent\n')
        next
      }
      wc <- 0
      while(is.na(i_pts[j, ])) {
        j_ext <- raster::extent(c(
          j_pt[, 2] - xinc,
          j_pt[, 2] + xinc,
          j_pt[, 1] - yinc,
          j_pt[, 1] + yinc
        ))
        j_pt_new <- raster::extract(i_tif, j_ext)
        if (length(j_pt_new) == 1) {
          i_pts[j, ] <- j_pt_new
        } else {
          xinc <- xinc + xpix * 0.1
          yinc <- yinc + ypix * 0.1
        }
        if (wc > 100) stop('Broke a while loop')
      }
    } # end j loop
    return_df <- cbind(return_df, i_pts)
    cat(round(i / length(file_ls) * 100), '% complete.\n')
  } # end i loop
  # Clean-up & return
  return_df <- data.frame(coords, return_df, stringsAsFactors = F)
  colnames(return_df) <- c('lat', 'lon', file_ls)
  return(return_df)
}
#' @rdname remote
#' @export
DropFirePerimeters <- function(df, fire_dir,
                               yrs = c(2000, 2015), as_logical = F) {
  # Sanity checks:
  if (length(yrs) > 2) stop('yrs must be length-2')

  warning('DropFirePerimeters isnt generalized... be careful')
  cat('\n')
  col_loc <- which(colnames(df) %in% c('LAT', 'LON'))
  pts <- df[, col_loc]
  colnames(pts) <- c('lat', 'lon')
  a0 <- yrs[1]
  a1 <- yrs[2]
  in_perims <- FilterPointsByGeoMac(coords = pts, years = c(a0:a1), dir = fire_dir)
  if (ncol(in_perims) == 2) {
    warning('FilterPointsByGeoMac failed, returning coordinates')
    return(in_perims)
  }
  if (as_logical) {
    in_fire_perim <- as.logical(rowSums(in_perims[, -c(1, 2)]))
    in_perims <- data.frame(in_perims[, c(1:2)], in_fire_perim,
                            stringsAsFactors = F)
    colnames(in_perims)[1:2] <- c('LAT', 'LON')
    df <- dplyr::left_join(df, in_perims, by = c('LAT', 'LON'))
    return(df)
  } else {
    stop('i think the more complex method is broken?')
  }
  yr <- 2000
  col <- 3
  drop_out <- numeric()
  while (a0 < (a1 + 1)) {
    drop_plot_keys <- df$unique_plot_key[in_perims[, col]]
    drop_out <- c(drop_out, drop_plot_keys)
    a0 <- a0 + 1
    col <- col + 1
  }
  if (length(which(df$unique_plot_key %in% drop_out)) > 0) {
    df <- df[-which(df$unique_plot_key %in% drop_out), ]
  }
  return(df)
}
#' @rdname remote
#' @export
FilterPointsByGeoMac <- function(coords, years, dir = getwd(), tag = 'perimeters') {
  # C:/Users/Brandon/Documents/docs/PHD/RSFIA_project/GeoMAC

  # Setup:
  require(sp)
  if (!(all(colnames(coords) %in% c('lat', 'lon')))) {
    if (all(c('LAT', 'LON') %in% colnames(coords))) {
      coords <- coords[, c('LAT', 'LON')]
      colnames(coords) <- c('lat', 'lon')
    } else {
      stop('coords need proper colnames')
    }
  }
  if (as.numeric(length(years)) < 1) stop('Bad years input')
  years <- as.character(years)
  exit_dir <- getwd()
  on.exit(setwd(exit_dir))
  setwd(dir)
  out_df <- data.frame(matrix(ncol = 2, nrow = nrow(coords)))
  out_df[, c(1, 2)] <- coords
  colnames(out_df)[1:2] <- c('lat', 'lon')
  coords <- data.frame(coords[, 2], coords[, 1])
  colnames(coords) <- c('lon', 'lat')
  coords_sp <- sp::SpatialPoints(coords)

  # Find folders:
  relv_dir <- list.files(pattern = tag)
  if (length(relv_dir) == 0) stop('didnt find folders')
  for (i in relv_dir) {
    if (!(unlist(strsplit(i, '_'))[1] %in% years)) next
    in_rast <- rgdal::readOGR(dsn = i)
    proj4string(coords_sp) <- sp::proj4string(in_rast)
    overlap_coords <- as.data.frame(coords_sp[in_rast]@coords)
    if (nrow(overlap_coords) == 0) next
    colnames(overlap_coords) <- c('lon', 'lat')
    overlap_coords[, 3] <- TRUE
    colnames(overlap_coords)[3] <- paste0('pt_in_fire_perim_', i)
    out_df <- dplyr::left_join(out_df, overlap_coords, by = c('lat', 'lon'))
  }
  if (ncol(out_df) == 2) {
    warning('didnt filter any points')
    return(out_df)
  }
  out_df[, 3:ncol(out_df)] <- data.frame(lapply(out_df[, 3:ncol(out_df)], function(x) {
    y <- ifelse(is.na(x), 0, x)
    y <- as.logical(y)
    return(y)
  }), stringsAsFactors = F)

  return(out_df)
}