RSFIA: FIA Tree Mortality Analysis

#' @include generics.R
#' @include nc_helper.R
#'
#' @title nc_helper extension that stores FIA variables
#'
#' @details
#'
#' For `remeasure`:
#'     `0`: return all plots
#'     `1`: return plots which are an initial census
#'     `2`: initial census plots that HAVE NOT BEEN remeasured at least once
#'     `3`: only plot remeasurements, no initial census
#'     `4`: initial census of remeasurements + remeasurements
#'     `5`: initial census plots that have been remeasured at least once
#'
#' @author Brandon McNellis
#'
#' @name FIA
#' @rdname FIA
NULL
#'
#' An S4 class for FIA data
#'
#' @rdname FIA
FIA <- setClass(
  'FIA',
  slots = list(
    plot_df = 'data.frame',
    cond_df = 'data.frame',
    tree_df = 'data.frame',
    FIA_dir = 'character',
    db_ver = 'numeric',
    provinces = 'character',
    remeasure = 'integer',
    plot_status = 'integer'
  ),
  contains = 'nc_helper'
)
#' @export
setValidity('FIA', function(object) {
  errors <- character()

  if (object@db_ver != 7.2) {
    msg <- paste0('Unsupported database version')
    errors <- c(errors, msg)
  }

  if (!(all(object@plot_status %in% c(1L, 2L, 3L)))) {
    msg <- paste0('Bad plot_status input, must be in 1:3')
    errors <- c(errors, msg)
  }

  if (!(object@remeasure %in% c(0L, 1L, 2L, 3L, 4L, 5L))) {
    msg <- paste0('Bad remeasure input, must be in 1:5')
    errors <- c(errors, msg)
  }

  if (length(object@provinces) > 0) {
    full_provs <- unique(ClelandEcoregions::Cleland_meta_df$province_code)
    if (!all(object@provinces %in% full_provs)) {
      msg <- paste0('Provinces are not in Cleland metadata')
      errors <- c(errors, msg)
    }
  }

  # returns
  if (length(errors) == 0) {
    TRUE
  } else {
    errors
  }
})
#' @rdname FIA
#' @export
setMethod('initialize',
          signature(.Object = 'FIA'),
          function (.Object, ...) {
            params <- list(...)

            if ('db_ver' %in% names(params)) {
              .Object@db_ver <- params$db_ver
            } else {
              .Object@db_ver <- 7.2
            }

            if ('plot_status' %in% names(params)) {
              .Object@plot_status <- params$plot_status
            } else {
              .Object@plot_status <- c(1L, 2L)
            }

            if ('remeasure' %in% names(params)) {
              .Object@remeasure <- params$remeasure
            } else {
              .Object@remeasure <- 0L
            }

            .Object@time_units <- 'years'

            # returns
            .Object <- callNextMethod()
            mt <- validObject(.Object)
            if (isTRUE(mt)) {
              return(.Object)
            } else {
              return(mt)
            }
          }
)
#' @rdname FIA
setMethod('SetupDataFile',
          signature(object = 'FIA'),
          function(object, overwrite = F) {
            vars <- c('PLOT_', 'COND_', 'TREE_')
            fn <- object@file_name
            for (i in seq_along(vars)) {
              ii <- vars[i]
              object@file_name <- paste0(ii, fn)
              nc_path <- paste0(object@nc_dir, '/', object@file_name)
              cat(nc_path, '\n')

              if (file.exists(nc_path)) {
                if (!overwrite) {
                  stop('file already exists, set overwrite = T to update')
                }
              }

              callNextMethod()
            }
            object@file_name <- fn
            invisible(object)
          })
#' @rdname FIA
#' @export
setMethod('GetMeta',
          signature(class = 'FIA'),
          function (class) {
            NULL
          }
)
#' @rdname FIA
#' @export
ImportPlots <- function(object, overwrite = F) {

  # input checks:
  validObject(object)
  path <- object@FIA_dir
  nc_path <- paste0(object@nc_dir, '/PLOT_', object@file_name)
  provs <- object@provinces
  rm0 <- object@remeasure

  # these can maybe be changed for other db versions?
  ec_col <- 'ECOSUBCD'
  yr_col <- 'INVYR'
  st_col <- 'PLOT_STATUS_CD'

  # work loop
  fl_list <- list.files(path = path, pattern = '_PLOT.csv')
  df0 <- data.frame(stringsAsFactors = F)
  okfl <- paste0(state.abb, '_PLOT.csv')

  message('ImportPlots:')
  for (i in seq_along(fl_list)) {
    ii <- fl_list[i]
    ii_path <- paste0(path, '/', ii)
    if (!(ii %in% okfl)) {
      next
    }

    ifl <- read.csv(ii_path, stringsAsFactors = F)
    if (nrow(ifl) == 0) {
      next
    }
    ifl <- data.frame(lapply(ifl, function(x) {
      if (is.character(x)) {
        y <- trimws(x)
      } else {
        y <- x
      }
      return(y) # this only works in long form for some reason?
    }), stringsAsFactors = F)

    if (length(provs) > 0) {
      subsecs <- unique(ClelandEcoregions::ScaleDownClelandCode(provs)$subsection_code)
      ifl <- ifl[which(ifl[, ec_col] %in% subsecs), ]
    }

    # other options:
    ifl <- ifl[which(ifl$DESIGNCD == 1), ]
    ifl <- ifl[which(ifl[[yr_col]] %in% object@time), ]
    ifl <- ifl[which(ifl[[st_col]] %in% object@plot_status), ]

    df0 <- rbind.fill(df0, ifl)

    cat('\r', format(i / length(fl_list) * 100, digits = 2, nsmall = 2), '%  ')
  } # end i

  if (nrow(df0) == 0) {
    stop('didnt find data')
  }

  df0 <- MakeUniquePlot(df0)

  df0 <- FilterRemeasured(df0, remeasure = object@remeasure)

  df0_nodup <- df0[, c('LON', 'LAT', 'unique_plot_key')]
  df0_nodup <- df0_nodup[!duplicated(df0_nodup), ]
  object <- AddCoords(object, lon = df0_nodup$LON, lat = df0_nodup$LAT, sample = df0_nodup$unique_plot_key)
  object@variables <- colnames(df0)[which(colnames(df0) %in% c('LON', 'LAT', 'INVYR', 'unique_plot_key') == F)]

  #SetupDataFile(object, overwrite = T)

  colnames(df0)[which(colnames(df0) == 'unique_plot_key')] <- 'sample'
  colnames(df0)[which(colnames(df0) == 'INVYR')] <- 'time'
  #load('/media/bem/F_STORAGE/scratch/df0.rda')

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

  df0$ECOSUBCD <- ClelandEcoregions::SubAsIntger(df0$ECOSUBCD)
  df0 <- data.frame(lapply(df0, function(x) {
    if (is.character(x)) {
      y <- suppressWarnings(as.numeric(x))
    } else {
      y <- x
    }
    return(y)
  }), stringsAsFactors = F)

  #df0 <- as.data.frame(tidyr::complete(df0, sample, time))
  #for (i in seq_along(object@time)) {
  #  ii <- object@time[i]
  #  dfi <- df0[which(df0$time == ii), ]
  #  FillArray(object, df = dfi, nc = nc0)
  #}
  if (TRUE) {
    # delete this when you get the nc to work
    cat('\nWriting .csv to', object@nc_dir)
    write.csv(df0, file = paste0(object@nc_dir, '/PLOT_DF.csv'))
  }

  #ncdf4::nc_close(nc0)
  object@plot_df <- df0
  object <- UpdateTimeStamp(object)
  return(object)

}
#' @rdname FIA
#' @export
FilterRemeasured <- function(FIA_df, remeasure = 0L) {
  stopifnot(
    'PREV_PLT_CN' %in% colnames(FIA_df),
    'CN' %in% colnames(FIA_df)
  )
  if (remeasure == 0L) {
    cat('No remeasurement filtering, returning all plots.')
    return(FIA_df)
  }

  # m1: initial establishment plot measurements
  # m2: remeasurement visits of previously established plots
  # m3: initial establishment plots that have been remeasured at least once in current dataset
  # m4: initial establishment plots that have yet to be remeasured in current dataset

  m1 <- FIA_df[which(is.na(FIA_df$PREV_PLT_CN)), ] # old var: first_meas
  m2 <- FIA_df[which(!is.na(FIA_df$PREV_PLT_CN)), ] # old var: was_remeas
  if (any(!(m2$PREV_PLT_CN %in% FIA_df$CN))) {
    w0 <- which(!(m2$PREV_PLT_CN %in% FIA_df$CN))
    m2 <- m2[-w0, ]
    cat('\nDropped', length(w0), 'plots with no initial visit plot in input dataframe\n')
  }

  m3 <- m1[which(m1$CN %in% m2$PREV_PLT_CN), ] # old var: remeas
  m4 <- m1[which(!(m1$CN %in% m2$PREV_PLT_CN)), ] # old var: not_remeas

  df0 <- data.frame(stringsAsFactors = F)
  switch(remeasure,
         '1' = {
           message('Returning initial-census plots (remeasured + not remeasured).')
           df0 <- m1
           nr <- table(is.na(df0[['PREV_PLT_CN']]))['TRUE']
           stopifnot(nr == nrow(df0))
         },
         '2' = {
           message('Returning not-remeasured plots.')
           df0 <- m4
         },
         '3' = {
           message('\nReturning remeasured plots, no initial census.')
           df0 <- m2
         },
         '4' = {
           message('Returning remeasured plots, including initial census.')
           df0 <- rbind(m2, m3, stringsAsFactors = F)
         },
         '5' = {
           message('Returning initial census plots that have been remeasured at least once.')
           df0 <- m3
         })

  #stopifnot(nrow(df0) > 0, nrow(df0) <= nrow(FIA_df))
  return(df0)
}
#' @rdname FIA
UpdateSubsections <- function(object, overwrite = F) {
  validObject(object)
  if (length(object@subsections) > 0) {
    if (overwrite == F) {
      stop('FIA object has subsections already, use overwrite = T?')
    }
  }
  provs <- object@provinces
  subs <- unique(ClelandEcoregions::ScaleDownClelandCode(provs)$subsection_code)
  object@subsections <- subs
  return(object)
}
#' @rdname FIA
#' @export
MakeUniquePlot <- function(FIA_df) {

  cat('\n')
  message('MakeUniquePlots:')
  check_vars <- c('LAT', 'LON', 'STATECD', 'COUNTYCD', 'PLOT')
  stopifnot(all(check_vars %in% colnames(FIA_df)))
  if (any(is.na(FIA_df[, check_vars]))) {
    stop('broke plot maker')
  }
  CNs <- unique(FIA_df$CN)

  unique_plot_key <- rep(1, nrow(FIA_df))
  for (i in seq_along(check_vars)) {
    ii <- check_vars[i]
    unique_plot_key <- unique_plot_key * FIA_df[, ii]
  }
  unique_plot_key <- abs(round(unique_plot_key))
  unique0 <- unique(unique_plot_key)

  check_out <- numeric()
  for (i in seq_along(unique0)) {
    ii <- unique0[i]

    cn1 <- na.omit(FIA_df$CN[which(unique_plot_key == ii)]) # ith plot, CNs
    cn2 <- na.omit(FIA_df$PREV_PLT_CN[which(unique_plot_key == ii)]) # ith plot, previous CNs
    cn3 <- na.omit(FIA_df$CN[-which(unique_plot_key == ii)]) # all other plots, CNs
    # all cn2 must be in cn1, and none of cn1 or cn2 can be in cn3

    if (any(
      !all(cn2 %in% cn1),
      any(c(cn1, cn2) %in% cn3)
      )) {
      #test_df <- FIA_df[which(FIA_df$CN %in% c(cn1, cn2)), ]
      #browser()
      check_out <- c(check_out, cn1)
    }
    cat('\r', format(i / length(unique0) * 100, nsmall = 2, digits = 2), '%  ')
  }

  df_out <- data.frame(unique_plot_key, FIA_df, stringsAsFactors = F)

  if (length(check_out) > 0) {
    cat('\nBad plots:\n')
    cat(check_out, '\n')
    df_out <- df_out[-which(df_out$CN %in% check_out), ]
  }

  return(df_out)
}
#' @rdname FIA
#' @export
ImportConditions <- function(object, overwrite = F) {

  # input checks:
  validObject(object)
  path <- object@FIA_dir
  c0 <- CoordVecsToList(object@coords)
  sample <- object@sample
  nc_cond_path <- paste0(object@nc_dir, '/', 'COND_', object@file_name)
  nc_plot_path <- paste0(object@nc_dir, '/', 'PLOT_', object@file_name)

  filt_PLT_CN <- object@plot_df$CN
  #if (file.exists(nc_cond_path)) {
  #  if (!overwrite) {
  #    stop('cond file exists, use overwrite = T to update')
  #  }
  #} else {
  #  stop('need to run nc setup using ImportPlots first')
  #}

  #ncc <- ncdf4::nc_open(nc_cond_path, write = T)
  #ncp <- ncdf4::nc_open(nc_plot_path, write = F)

  #varCN <- ncdf4::ncvar_get(ncp, varid = 'CN', start = c(1, 1), count = c(-1, -1))
  #varCN <- ncdf4::ncvar_get(ncp, varid = 'CN')
  #varCN <- data.frame(object@sample, varCN)
  #colnames(varCN) <- c('sample', as.character(2000:2015))
  #cn0 <- tidyr::gather(varCN, time, value = CN, as.character(2000:2015))
  #filt_PLT_CN <- unique(cn0$CN)

  fls <- list.files(path = path, pattern = '_COND.csv')
  fls <- fls[-grep(fls, pattern = 'SUBP')]
  if (length(fls) < 1) {
    stop('No COND files found.')
  }
  df0 <- data.frame(stringsAsFactors = F)

  # work loop
  message('COND file read:')
  for (i in seq_along(fls)) {
    ii <- fls[i]
    rp <- paste0(path, '/', ii)
    print(ii)
    i_fl <- read.csv(rp, stringsAsFactors = F)

    if (length(filt_PLT_CN) > 0) {
      rel_con <- which(i_fl[['PLT_CN']] %in% filt_PLT_CN)
    } else {
      rel_con <- c(1:nrow(i_fl))
    }

    if (length(rel_con) < 1) {
      next
    } else {
      df0 <- rbind.fill(df0, i_fl[rel_con, ])
    }
  }
  cat('\n')

  filt_PLT_CN <- object@plot_df$CN
  df0 <- df0[which(df0$PLT_CN %in% filt_PLT_CN), ]
  if (TRUE) {
    write.csv(df0, file = paste0(object@nc_dir, '/COND_DF.csv'))
    object@cond_df <- df0
  }

  return(object)

  #df0 <- as.data.frame(tidyr::complete(df0, sample, time))
  #for (i in seq_along(object@time)) {
  #  ii <- object@time[i]
  #  dfi <- df0[which(df0$time == ii), ]
  #  FillArray(object, df = dfi, nc = nc0)
  #}

}
#' @rdname FIA
#' @export
ImportTrees <- function(object) {
  validObject(object)
  path <- object@FIA_dir
  filt_PLT_CN <- object@plot_df$CN

  fls <- list.files(path = path, pattern = '_TREE.csv')
  #fls <- fls[-grep(fls, pattern = 'SUBP')]
  if (length(fls) < 1) {
    stop('No TREE files found.')
  }

  df0 <- data.frame(stringsAsFactors = F)

  # work loop
  message('TREE file read:')
  for (i in seq_along(fls)) {
    ii <- fls[i]
    rp <- paste0(path, '/', ii)
    print(ii)
    i_fl <- read.csv(rp, stringsAsFactors = F)

    if (length(filt_PLT_CN) > 0) {
      rel_con <- which(i_fl[['PLT_CN']] %in% filt_PLT_CN)
    } else {
      rel_con <- c(1:nrow(i_fl))
    }

    if (length(rel_con) < 1) {
      next
    } else {
      df0 <- rbind.fill(df0, i_fl[rel_con, ])
    }
  }
  cat('\n')

  filt_PLT_CN <- object@plot_df$CN
  df0 <- df0[which(df0$PLT_CN %in% filt_PLT_CN), ]
  write.csv(df0, file = paste0(object@nc_dir, '/TREE_DF.csv'))
  object@tree_df <- df0

  return(object)
}
#' @rdname FIA
#' @export
Smallify <- function(object) {
  cond_fname <- paste0(object@nc_dir, '/COND_DF.csv')
  plot_fname <- paste0(object@nc_dir, '/PLOT_DF.csv')
  tree_fname <- paste0(object@nc_dir, '/TREE_DF.csv')
  stopifnot(
    file.exists(cond_fname),
    file.exists(plot_fname),
    file.exists(tree_fname)
  )
  object@cond_df <- data.frame(cond_fname, stringsAsFactors = F)
  object@plot_df <- data.frame(plot_fname, stringsAsFactors = F)
  object@tree_df <- data.frame(tree_fname, stringsAsFactors = F)
  return(object)
}
#' @rdname FIA
#' @export
Biggify <- function(object) {
  cond_fname <- as.character(object@cond_df)
  plot_fname <- as.character(object@plot_df)
  tree_fname <- as.character(object@tree_df)
  stopifnot(
    file.exists(cond_fname),
    file.exists(plot_fname),
    file.exists(tree_fname)
  )
  object@cond_df <- read.csv(cond_fname)
  object@plot_df <- read.csv(plot_fname)
  object@tree_df <- read.csv(tree_fname)
}
#' @rdname FIA
#' @export
PullDF <- function(object, type) {
  if (type == 'plots') {
    return(read.csv(as.character(object@plot_df)))
  } else if (type == 'conds') {
    return(read.csv(as.character(object@cond_df)))
  } else if (type == 'trees') {
    return(read.csv(as.character(object@trees_df)))
  } else {
    stop('bad type')
  }
}