FIESTA: Forest Inventory Estimation and Analysis

Documented in modMApop

#' Model-Assisted module - Generate population data for MA module.
#' 
#' Generates population data for generating model-assisted estimation.  Plots
#' that are totally nonsampled are excluded from estimation dataset. Next, an
#' adjustment factor is calculated by strata (if MAmethod="PS") or by
#' estimation unit to adjust for nonsampled (nonresponse) conditions that have
#' proportion less than 1. Attributes adjusted to a per-acre value are summed
#' by plot, divided by the adjustment factor, and averaged by stratum and/or
#' estimation unit. Note: population data must be generated by MA method.
#' 
#' If variables are NULL, then it will prompt user to input variables.
#' 
#' Necessary variables:\cr 
#' \tabular{llll}{ 
#' \tab \bold{Data} \tab \bold{Variable} \tab \bold{Description}\cr 
#' \tab tree \tab tuniqueid 
#' \tab Unique identifier for each plot, to link to pltassgn (e.g. PLT_CN).\cr 
#' \tab \tab CONDID \tab Unique identifier of each condition on plot, to link to
#' cond.  Set CONDID=1, if only 1 condition per plot.\cr 
#' \tab \tab TPA_UNADJ \tab Number of trees per acre each sample tree represents 
#' (e.g. DESIGNCD=1: TPA_UNADJ=6.018046 for trees on subplot; 74.965282 for 
#' trees on microplot).\cr 
#' \tab cond \tab cuniqueid \tab Unique identifier for each plot, to link to 
#' pltassgn (e.g. PLT_CN).\cr 
#' \tab \tab CONDID \tab Unique identifier of each condition on plot. Set 
#' CONDID=1, if only 1 condition per plot.\cr 
#' \tab \tab CONDPROP_UNADJ \tab Unadjusted proportion of condition on
#' each plot.  Set CONDPROP_UNADJ=1, if only 1 condition per plot.\cr 
#' \tab \tab COND_STATUS_CD \tab Status of each forested condition on plot 
#' (i.e. accessible forest, nonforest, water, etc.)\cr 
#' \tab \tab NF_COND_STATUS_CD \tab If ACI=TRUE. Status of each nonforest 
#' condition on plot (i.e. accessible nonforest, nonsampled nonforest)\cr 
#' \tab \tab SITECLCD \tab If landarea=TIMBERLAND. Measure of site 
#' productivity.\cr 
#' \tab \tab RESERVCD \tab If landarea=TIMBERLAND. Reserved status.\cr 
#' \tab \tab SUBPROP_UNADJ \tab Unadjusted proportion of subplot conditions 
#' on each plot. Set SUBPROP_UNADJ=1, if only 1 condition per subplot.\cr 
#' \tab \tab MICRPROP_UNADJ \tab If microplot tree attributes. Unadjusted 
#' proportion of microplot conditions on each plot. Set MICRPROP_UNADJ=1, 
#' if only 1 condition per microplot.\cr 
#' \tab \tab MACRPROP_UNADJ \tab If macroplot tree attributes. Unadjusted 
#' proportion of macroplot conditions on each plot. Set MACRPROP_UNADJ=1, 
#' if only 1 condition per macroplot.\cr \tab pltassgn \tab puniqueid 
#' \tab Unique identifier for each plot, to link to cond (e.g. CN).\cr
#' \tab \tab STATECD \tab Identifies state each plot is located in.\cr 
#' \tab \tab INVYR \tab Identifies inventory year of each plot.\cr 
#' \tab \tab PLOT_STATUS_CD \tab Status of each plot (i.e. sampled, 
#' nonsampled). If not included, all plots are assumed as sampled.\cr }
#' 
#' For available reference tables: sort(unique(FIESTAutils::ref_codes$VARIABLE)) \cr
#' 
#' @param popType String. Type of evaluation(s) to include in population data.
#' Note: currently only c('CURR', 'VOL', 'LULC') are available. See details
#' below for descriptions of each.
#' @param popTabs List of population tables the user would like returned.
#'  See help(popTables) for a list of options.
#' @param popTabIDs List of unique IDs corresponding to the population tables
#' that the user has requested. See help(popTableIDs) for a list of
#' options.
#' @param popFilter List of population filters. See help(popFilters) for a 
#' list of options. 
#' @param pltassgn DF/DT, Optional. R object, sf R object, comma-delimited
#' file(.csv), layer or spatial layer in dsn, or shapefile(.shp). Plot-level
#' assignment of estimation unit and/or strata, with one record for each plot.
#' @param pltassgnid String.
#' @param datsource String. Name of data source ('obj', 'sqlite', 'postgres').
#' @param dsn String. Name of database where tree, cond, and plot-level tables
#' reside.  The dsn varies by driver. See gdal OGR vector formats
#' (https://www.gdal.org/ogr_formats.html).
#' @param dbconn Open database connection.
#' @param pjoinid String. Join variable in plot to match pltassgnid. Does not
#' need to be uniqueid. If using most current XY coordinates for plot
#' assignments, use identifier for plot (e.g., PLOT_ID).
#' @param areawt String. Name of variable for summarizing area weights (e.g.,
#' CONDPROP_UNADJ).
#' @param adj String. How to calculate adjustment factors for nonsampled
#' (nonresponse) conditions based on summed proportions for by plot ('samp',
#' 'plot', 'none').  'samp' - adjustments are calculated at strata/estimation unit
#' level; 'plot' - adjustments are calculated at plot-level. Adjustments are
#' only calculated for annual inventory plots (DESIGNCD=1).
#' @param defaultVars Logical. If TRUE, a set of default variables are selected.
#' @param unitvar String. Name of the estimation unit variable in unitarea and
#' cond or pltassgn data frame with estimation unit assignment for each plot
#' (e.g., 'ESTN_UNIT'). Optional if only one estimation unit.
#' @param unitarea Numeric or DF. Total area by estimation unit. If only 1
#' estimation unit, include number of total acreage for the area of interest or
#' a data frame with area and estimation unit. If more than one estimation
#' unit, provide a data frame of total area by estimation unit, including
#' unitvar and areavar.
#' @param areavar String. Name of area variable in unitarea. Default="ACRES".
#' @param unitzonal DF/DT. Table with zonal auxiliary information by estimation
#' unit. For continuous data, means by estimation unit; for categorical data,
#' proportion of class by estimation unit.
#' @param prednames String vector. Name(s) of predictor variables to include in
#' model.
#' @param predfac String vector. Name(s) of prednames that are factors (i.e.,
#' categorical). Names will change in output depending on number of categories.
#' @param standardize Logical. If TRUE, predictors are standardized. 
#' @param returndata Logical. If TRUE, returns data objects.
#' @param savedata Logical. If TRUE, saves table(s) to outfolder.
#' @param saveobj Logical. If TRUE, saves returned list object to outfolder.
#' @param objnm String. Name of *.rds object.
#' @param unit_opts List. See help(unit_options()) for a list of options.
#' @param savedata_opts List. See help(savedata_options()) for a list
#' of options. Only used when savedata = TRUE.  
#' @param database_opts List. See help(database_options()) for a list
#' of options. Only used when datsource = 'postgres'.  
#' @param MAdata List. Data output from FIESTA::MAdata().
#' @param pltdat R List object. Output data list components from
#' FIESTA::spGetPlots().
#' @param auxdat List. Auxiliary data output from FIESTA::spGetAuxiliary().
#' @param ... For extendibility.
#' @return A list with population data for Green-Book estimates.
#' 
#' \item{condx}{ Data frame. Condition-level data including plot-level
#' assignment of estimation unit and stratum (if strata=TRUE) and adjusted
#' condition proportion. } 
#' \item{pltcondx}{ Data frame. Condition-level data, merged with plot data. } 
#' \item{cuniqueid}{ String. Unique identifier of plot in condx and pltcondx. } 
#' \item{condid}{ String. Unique identifier of condition in condx and pltcondx. } 
#' \item{treex}{ Data frame. If esttype='TREE', tree-level data, including 
#' sample adjustment factor. }
#' \item{tuniqueid}{ String. If esttype='TREE', unique identifier of plot in
#' treex. } 
#' \item{ACI.filter}{ String. If ACI=FALSE, ACI.filter="COND_STATUS_CD
#' == 1" . } 
#' \item{unitarea}{ String. Returned table of area by estimation unit. } 
#' \item{unitvar}{ String. Variable name for estimation unit. }
#' \item{expcondtab}{ String. If ACI=FALSE, ACI.filter="COND_STATUS_CD == 1". } 
#' \item{plotsampcnt}{ Data frame. Number of plots by PLOT_STATUS_CD. } 
#' \item{condsampcnt}{ Data frame. Number of conditions by COND_STATUS_CD. } 
#' \item{states}{ String. State names in dataset. } 
#' \item{invyrs}{ String. Range of inventory years in dataset. }
#' 
#' \tabular{lll}{ \tab \bold{Variable} \tab \bold{Description}\cr 
#' \tab unitvar \tab estimation unit \cr 
#' \tab n.total \tab number of plots for estimation unit \cr 
#' \tab CONDPROP_UNADJ_SUM \tab summed condition proportion by strata and 
#' estimation unit \cr 
#' \tab CONDPROP_ADJFAC \tab adjusted condition proportion by strata after 
#' nonsampled plots removed \cr 
#' \tab AREA_USED \tab total area of estimation unit \cr 
#' \tab expfac \tab strata-level expansion factor after nonsampled plots and 
#' conditions removed (AREA_USED/n.strata) \cr 
#' \tab EXPNS \tab strata-level area expansions (expfac * strwt)\cr }
#' 
#' Table(s) are also written to outfolder.
#' @note
#' 
#' ADJUSTMENT FACTOR:\cr The adjustment factor is necessary to account for
#' nonsampled conditions. It is calculated for each estimation unit by strata.
#' by summing the unadjusted proportions of the subplot, microplot, and
#' macroplot (i.e. *PROP_UNADJ) and dividing by the number of plots in the
#' strata/estimation unit).
#' 
#' An adjustment factor is determined for each tree based on the size of the
#' plot it was measured on. This is identified using TPA_UNADJ as follows:
#' 
#' \tabular{llr}{ \tab \bold{PLOT SIZE} \tab \bold{TPA_UNADJ} \cr 
#' \tab SUBPLOT \tab 6.018046 \cr 
#' \tab MICROPLOT \tab 74.965282 \cr 
#' \tab MACROPLOT \tab 0.999188 \cr }
#' 
#' If ACI=FALSE, only nonsampled forest conditions are accounted for in the
#' adjustment factor. \cr 
#' If ACI=TRUE, the nonsampled nonforest conditions are
#' removed as well and accounted for in adjustment factor.  This is if you are
#' interested in estimates for all lands or nonforest lands in the
#' All-Condition-Inventory.
#' 
#' unitcombine:\cr If TRUE and less than 2 plots in any one estimation unit,
#' all estimation units with 10 or less plots are combined. The current method
#' for combining is to group the estimation unit with less than 10 plots with
#' the estimation unit following in consecutive order (numeric or
#' alphabetical), restrained by survey unit (UNITCD) if included in dataset,
#' and continuing until the number of plots equals 10. If there are no
#' estimation units following in order, it is combined with the estimation unit
#' previous in order.
#' 
#' stratcombine:\cr If TRUE and less than 2 plots in any one strata class
#' within an estimation unit, all strata classes with 2 or less plots are
#' combined. The current method for combining is to group the strata with less
#' than 2 plots with the strata class following in consecutive order (numeric
#' or alphabetical), restrained by estimation unit (if unitcombine=FALSE), and
#' continuing until the number of plots equals 10. If there are no strata
#' classes following in order, it is combined with the estimation unit previous
#' in order.
#' @author Tracey S. Frescino, Paul L. Patterson
#' @references Scott, Charles T.; Bechtold, William A.; Reams, Gregory A.;
#' Smith, William D.; Westfall, James A.; Hansen, Mark H.; Moisen, Gretchen G.
#' 2005. Sample-based estimators used by the Forest Inventory and Analysis
#' national information management system. Gen. Tech. Rep. SRS-80. Asheville,
#' NC: U.S. Department of Agriculture, Forest Service, Southern Research
#' Station, p.53-77.
#' @keywords data
#' @examples 
#' \donttest{
#' # NOTE: FIA data objects used in these examples are stored in `FIESTA`, but
#' # can be generated for populations of interest by the user with functions in
#' # `FIESTA` such as `spGetPlots()`, `spGetAuxiliary()`, etc. For more
#' # information, see `FIESTA`'s extensive vignettes.
#'  
#' # Population data for counties in Wyoming
#' modMApop(popTabs = list(tree = FIESTA::WYtree,
#'                         cond = FIESTA::WYcond),
#'          pltassgn = FIESTA::WYpltassgn,
#'          pltassgnid = "CN",
#'          unitarea = FIESTA::WYunitarea,
#'          unitvar = "ESTN_UNIT",
#'          unitzonal = FIESTA::WYunitzonal,
#'          prednames = c("dem", "tcc", "tpi", "tnt"),
#'          predfac = "tnt")
#'          
#' # Adding seedling data as well
#' modMApop(popTabs = list(tree = FIESTA::WYtree,
#'                         cond = FIESTA::WYcond,
#'                         seed = FIESTA::WYseed),
#'          pltassgn = FIESTA::WYpltassgn,
#'          pltassgnid = "CN",
#'          unitarea = FIESTA::WYunitarea,
#'          unitvar = "ESTN_UNIT",
#'          unitzonal = FIESTA::WYunitzonal,
#'          prednames = c("dem", "tcc", "tpi", "tnt"),
#'          predfac = "tnt")
#' }
#' @export modMApop
modMApop <- function(popType="VOL",
                     popTabs = popTables(),
                     popTabIDs = popTableIDs(), 
                     popFilter = popFilters(),
                     pltassgn = NULL,
                     pltassgnid = "PLT_CN", 
                     datsource = "sqlite",
                     dsn = NULL, 
                     dbconn = NULL,
                     pjoinid = "CN", 
                     areawt = "CONDPROP_UNADJ", 
                     adj = "plot", 
                     defaultVars = TRUE,
                     unitvar = NULL, 
                     unitarea = NULL, 
                     areavar = "ACRES", 
                     unitzonal = NULL, 
                     prednames = NULL, 
                     predfac = NULL, 
                     standardize = TRUE,
                     returndata = TRUE,
                     savedata = FALSE, 
                     saveobj = FALSE, 
                     objnm = "MApopdat", 
                     unit_opts = NULL, 
                     savedata_opts = NULL, 
                     database_opts = NULL,
                     MAdata = NULL, 
                     pltdat = NULL, 
                     auxdat = NULL,
                     ...){
  
  ##################################################################################
  ## DESCRIPTION:
  ## Generates population data 'on-the-fly', including strata weights, number
  ## of plots by strata and estimation unit, strata-level expansion factors,
  ## and sample-based area adjustment factors.
  ## - checks input parameters and data tables, including removing nonsampled
  ##   plots and conditions (see check.popdata for details).
  ## - checks auxiliary data (i.e., stratification data).
  ## - calculates adjustment factors for nonresponse and appends an adjustment
  ##   variable to condition and tree data.
  ##################################################################################
  
  ## CHECK GUI - IF NO ARGUMENTS SPECIFIED, ASSUME GUI=TRUE
  gui <- FALSE
  
  ## If gui.. set variables to NULL
  if (gui) {
    areavar=cuniqueid=ACI=tuniqueid=savedata=unitvar=projectid <- NULL
  }
  
  ## Set parameters
  adjtree <- FALSE
  nonsamp.pfilter=nonsamp.cfilter=schema=vcondstrx=vcondsppx=outlst <- NULL
  returnlst <- list(module = "MA")
  
  ## Set global variables
  ONEUNIT=n.total=expcondtab=bndx <- NULL
  strata <- FALSE
  condid <- "CONDID"
  areawt2 <- NULL
  pvars2keep <- NULL
  pltidsadjindb=savepltids=dsnreadonly <- FALSE
  
  
  ##################################################################
  ## CHECK PARAMETER NAMES
  ##################################################################
  
  ## Check input parameters
  input.params <- names(as.list(match.call()))[-1]
  formallst <- names(formals(modMApop)) 
  if (!all(input.params %in% formallst)) {
    miss <- input.params[!input.params %in% formallst]
    stop("invalid parameter: ", toString(miss))
  }
  
  ## Check parameter lists
  pcheck.params(input.params = input.params,
                unit_opts = unit_opts, 
                savedata_opts = savedata_opts, database_opts = database_opts)
 
  ## Check parameter option lists
  optslst <- pcheck.opts(optionlst = list(
                         popFilter = popFilter,
                         unit_opts = unit_opts, 
                         savedata_opts = savedata_opts,
                         database_opts = database_opts))
  savedata_opts <- optslst$savedata_opts  
  unit_opts <- optslst$unit_opts  
  database_opts <- optslst$database_opts  
  popFilter <- optslst$popFilter
  
  for (i in 1:length(unit_opts)) {
    assign(names(unit_opts)[[i]], unit_opts[[i]])
  }
  for (i in 1:length(savedata_opts)) {
    assign(names(savedata_opts)[[i]], savedata_opts[[i]])
  }
 
  
  ##################################################################
  ## CHECK PARAMETER INPUTS
  ##################################################################
  
  ## Check returndata 
  returndata <- FIESTAutils::pcheck.logical(returndata, varnm="returndata", 
        title="Return data as objectsd?", first="YES", gui=gui, stopifnull=TRUE)
  
  ## Check savedata 
  savedata <- pcheck.logical(savedata, varnm="savedata", 
        title="Save data tables?", first="YES", gui=gui, stopifnull=TRUE)
  ## Check saveobj 
  saveobj <- pcheck.logical(saveobj, varnm="saveobj", 
        title="Save SApopdat object?", first="YES", gui=gui, stopifnull=TRUE)
  
  ## Check output
  ########################################################
  if (savedata || saveobj) {
    outlst <- pcheck.output(savedata_opts = savedata_opts)
    if (savedata) {
      if (outlst$out_fmt == "sqlite" && is.null(outlst$out_dsn)) {
        outlst$out_dsn <- "MApopdat.db"
      }
      outlst$add_layer <- TRUE
    }
  }
  
  if (saveobj) {
    outobj_fmtlst <- c('rds', 'rda')
    outobj_fmt <- pcheck.varchar(var2check=outobj_fmt, varnm="outobj_fmt", gui=gui,
                                 checklst=outobj_fmtlst, caption="outobj_fmt",
                                 multiple=FALSE, stopifnull=TRUE)
    
    if (is.null(objnm)) {
      objnm <- "MApopdat"
    }
    #if (append_layer) overwrite_layer <- FALSE
    if (append_layer) message("currently cannot append to object lists")
    objfn <- getoutfn(outfn = objnm,
                      ext = outobj_fmt,
                      outfolder = outfolder, 
                      overwrite = overwrite_layer,
                      outfn.pre = outfn.pre,
                      outfn.date = outfn.date)
  }
  
 
  ## Check popType
  ########################################################
  #evalTyplst <- c("ALL", "CURR", "VOL", "LULC", "P2VEG", "INV", "GRM", "DWM")
  DWM_types <- c("CWD", "FWD_SM", "FWD_LG", "DUFF")
  evalTyplst <- c("ALL", "CURR", "VOL", "LULC", "P2VEG", "INV", "DWM",
                  "CHNG", "GRM", "GROW", "MORT", "REMV")
  popType <- pcheck.varchar(var2check=popType, varnm="popType", gui=gui,
                            checklst=evalTyplst, caption="popType", multiple=FALSE, 
                            stopifinvalid=FALSE)
  if (is.null(popType)) {
    message("popType is invalid... must be from following list:\n", toString(evalTyplst))
  }
  popevalid <- popFilter$evalid
  if (!is.null(popevalid)) {
    popevalid <- as.character(popevalid)
    substr(popevalid, nchar(popevalid)-1, nchar(popevalid)) <- 
      formatC(FIESTAutils::ref_popType[FIESTAutils::ref_popType$popType %in% popType, "EVAL_TYP_CD"], 
              width=2, flag="0")
    #evalid <- as.character(evalid)
    #substr(evalid, nchar(evalid)-1, nchar(evalid)) <- "01"
  } 
  if (popType %in% c("GROW", "MORT", "REMV")) {
    popType <- "GRM"
  }
  
  ###################################################################################
  ## Load data
  ###################################################################################
  if (!is.null(MAdata)) {
    list.items <- c("tabs", "unitarea", "unitvar", "unitzonal")
    MAdata <- pcheck.object(MAdata, "MAdata", list.items=list.items)
    #bnd <- MAdata$bnd
    popTabs <- MAdata$tabs
    popTabIDs <- MAdata$tabIDs
    pltassgn <- MAdata$pltassgn
    pltassgnid <- MAdata$pltassgnid
    unitarea <- MAdata$unitarea
    areavar <- MAdata$areavar
    unitzonal <- MAdata$unitzonal
    puniqueid <- MAdata$puniqueid
    pjoinid <- MAdata$pjoinid
    
    if (is.null(unitvar)) {
      unitvar <- MAdata$unitvar
      unitvar2 <- MAdata$unitvar2
    } 
    if (is.null(npixelvar)) {
      npixelvar <- MAdata$npixelvar
    }
    if (is.null(prednames)) {
      prednames <- MAdata$prednames
    } else {
      if (!all(prednames %in% MAdata$prednames)) {
        stop("invalid prednames: ", toString(prednames[!prednames %in% MAdata$prednames]))
      }
    }
    if (is.null(predfac)) {
      predfac <- MAdata$predfac
    }
    predfac <- predfac[predfac %in% prednames]
    
  } else {
    if (!is.null(pltdat)) {
      datsource <- "obj"
      tabnames <- if (sum(names(pltdat$tabs) %in% names(popTables())) == 0) {
        stop("no tables exist in pltdat")
      }
      popTabs <- pltdat$tabs
      popTabIDs <- pltdat$tabIDs
      pjoinid <- pltdat$pjoinid
      spxy <- pltdat$spxy
      xy.uniqueid <- pltdat$xy.uniqueid
    }
    if (!is.null(auxdat)) {
      list.items <- c("pltassgn", "unitzonal", "unitvar", "prednames", "unitarea")
      auxdat <- pcheck.object(auxdat, "auxdat", list.items=list.items)
      pltassgn <- auxdat$pltassgn
      pltassgnid <- auxdat$pltassgnid
      unitzonal <- auxdat$unitzonal
      unitvar <- auxdat$unitvar
      unitvar2 <- auxdat$unitvar2
      unitarea <- auxdat$unitarea
      areavar <- auxdat$areavar
      
      if (is.null(npixelvar)) {
        npixelvar <- auxdat$npixelvar
      }
      if (is.null(prednames)) {
        prednames <- auxdat$prednames
      } else {
        if (!all(prednames %in% auxdat$prednames))
          stop("invalid prednames: ", toString(prednames[!prednames %in% auxdat$prednames]))
      }
      if (is.null(predfac)) {
        predfac <- auxdat$predfac
      }
      predfac <- predfac[predfac %in% prednames]
    }
  } 
  
  ## Set user-supplied popTable values 
  popTables_defaults_list <- formals(popTables)[-length(formals(popTables))]
  if (length(popTabs) > 0) {
    for (i in 1:length(popTabs)) {
      if (names(popTabs)[[i]] %in% names(popTables_defaults_list)) {
        assign(names(popTabs)[[i]], popTabs[[i]])
      } else {
        stop(paste("Invalid parameter: ", names(popTabs)[[i]]))
      }
    }
  } else {
    stop("need to include popTabs")
  }
  
  list.items <- c("cond")
  if (popType == "VOL") {
    list.items <- c(list.items, "tree")
  }
  if (popType == "P2VEG") {
    list.items <- c(list.items, "vsubpstr", "subplot", "subp_cond")
  }
  if (popType == "DWM") {
    list.items <- c(list.items, "cond_dwm_calc")
  }
  if (popType == "CHNG") {
    list.items <- c(list.items, "sccm")
  }
  popTabs <- pcheck.object(popTabs, "popTabs", list.items=list.items)
  
  
  ## Set user-supplied popTabIDs values
  ### Check for invalid parameters first
  popTableIDs_defaults_list <- formals(popTableIDs)[-length(formals(popTableIDs))]
  for (i in 1:length(popTabIDs)) {
    if (!(names(popTabIDs)[[i]] %in% names(popTableIDs_defaults_list))) {
      stop(paste("Invalid parameter: ", names(popTabIDs)[[i]]))
    }
  }
  ### Then actually set the values
  for (nm in names(popTabs)) {
    if (!any(names(popTabIDs) == nm)) {
      popTabIDs[[nm]] <- popTableIDs_defaults_list[[nm]]
    }
  }
  
  ###################################################################################
  ## CHECK PLOT PARAMETERS AND DATA
  ## Generate table of sampled/nonsampled plots and conditions
  ## Remove nonsampled plots (if nonsamp.pfilter != "NONE")
  ## Applies plot filters
  ###################################################################################
  pltcheck <- 
    check.popdataPLT(dsn = dsn, dbconn = dbconn, schema = schema,
                     datsource = datsource, 
                     tabs = popTabs, tabIDs = popTabIDs, 
                     pltassgn = pltassgn,
                     pltassgnid = pltassgnid, pjoinid = pjoinid, 
                     module = "MA", popType = popType,
                     popevalid = popevalid, adj = adj, 
                     popFilter = popFilter,
                     nonsamp.pfilter = nonsamp.pfilter, 
                     unitarea = unitarea, areavar = areavar,
                     unitvar = unitvar, unitvar2 = unitvar2,
                     areaunits = areaunits, 
                     unit.action = unit.action, 
                     auxlut = unitzonal, 
                     defaultVars = defaultVars,
                     prednames = prednames, predfac = predfac,
                     pvars2keep = pvars2keep, 
                     dsnreadonly = dsnreadonly)
  if (is.null(pltcheck)) return(0)
  pltassgnx <- pltcheck$pltassgnx
  pltassgnid <- pltcheck$pltassgnid
  pltassgn. <- pltcheck$pltassgn.
  plotlst <- pltcheck$plotlst
  pltidsWITHqry <- pltcheck$pltidsWITHqry
  pltidsid <- pltcheck$pltidsid
  pltidvars <- pltcheck$pltidvars
  projidvars <- pltcheck$projidvars
  pdoms2keep <- pltcheck$pdoms2keep
  ACI <- pltcheck$ACI
  unitvar <- pltcheck$unitvar
  unitvar2 <- pltcheck$unitvar2
  unitarea <- pltcheck$unitarea
  areavar <- pltcheck$areavar
  areaunits <- pltcheck$areaunits
  unit.action <- pltcheck$unit.action
  P2POINTCNT <- pltcheck$P2POINTCNT 
  plotsampcnt <- pltcheck$plotsampcnt
  states <- pltcheck$states
  invyrs <- pltcheck$invyrs
  dbconn <- pltcheck$dbconn
  SCHEMA. <- pltcheck$SCHEMA.
  pltaindb <- pltcheck$pltaindb
  datindb <- pltcheck$datindb
  POP_PLOT_STRATUM_ASSGN <- pltcheck$POP_PLOT_STRATUM_ASSGN
  getdataWITHqry <- pltcheck$getdataWITHqry
  getdataCNs <- pltcheck$getdataCNs
  plotunitcnt <- pltcheck$plotunitcnt
  prednames <- pltcheck$prednames
  predfac <- pltcheck$predfac
  auxlut <- unitzonal
  getdataWITHqry <- pltcheck$getdataWITHqry
  getdataCNs <- pltcheck$getdataCNs

  if (ACI) {
    nfplotsampcnt <- pltcheck$nfplotsampcnt
  }
  
  ###################################################################################
  ## Check auxiliary data
  ###################################################################################
  makedummy <- TRUE
  auxcheck <- 
    check.auxiliary(module = "MA", 
                    pltx = pltassgnx,
                    puniqueid = pltassgnid, 
                    unitvar = unitvar,
                    unitvar2 = unitvar2, 
                    unitarea = unitarea, 
                    areavar = areavar, 
                    minplotnum.unit = minplotnum.unit, 
                    unit.action = unit.action, 
                    auxlut = auxlut, 
                    prednames = prednames, 
                    predfac = predfac, 
                    makedummy = makedummy, 
                    npixelvar = npixelvar, 
                    standardize = standardize,
                    auxtext = "unitlut", 
                    removetext = "unitarea", 
                    AOI = popFilter$AOIonly)
  pltassgnx <- setDT(auxcheck$pltx)
  unitarea <- auxcheck$unitarea
  unitvar <- auxcheck$unitvar
  unitvars <- auxcheck$unitvars
  unitlut <- auxcheck$auxlut
  prednames <- auxcheck$prednames
  predfac <- auxcheck$predfac
  npixels <- auxcheck$npixels
  unitNA <- auxcheck$unitNA
  if (is.null(key(pltassgnx))) setkeyv(pltassgnx, pltassgnid)
  

  if (popType %in% c("ALL", "CURR", "VOL")) {
    
    ###################################################################################
    ## Check parameters and data for popType AREA/VOL
    ###################################################################################
    areawt <- "CONDPROP_UNADJ"
    popcheck <- 
      check.popdataVOL(tabs = popTabs, tabIDs = popTabIDs,
                       popType = popType,
                       datindb = datindb, pltaindb = pltaindb,
                       pltidsWITHqry = pltidsWITHqry, 
                       pltidsid = pltidsid,
                       pltidvars = pltidvars, projidvars = projidvars,
                       pdoms2keep = pdoms2keep,
                       defaultVars = defaultVars,
                       pltidsadjindb = FALSE, 
                       pltassgnid = pltassgnid, 
                       pltassgnx = pltassgnx,
                       POP_PLOT_STRATUM_ASSGN = POP_PLOT_STRATUM_ASSGN,
                       adj = adj, ACI = ACI, 
                       plotlst = plotlst,  
                       condid = condid, 
                       areawt = areawt, areawt2 = areawt2,
                       unitvars = unitvars,
                       nonsamp.cfilter = nonsamp.cfilter,
                       dbconn = dbconn, SCHEMA. = SCHEMA.,
                       getdataWITHqry = getdataWITHqry,
                       getdataCNs = getdataCNs,
                       returndata = returndata,
                       savedata = savedata,
                       outlst = outlst)
    if (is.null(popcheck)) return(NULL)
    pltidsadj <- popcheck$pltidsadj
    pltcondx <- popcheck$pltcondx
    pltcondflds <- popcheck$pltcondflds
    cuniqueid <- popcheck$cuniqueid
    condid <- popcheck$condid
    adjfactors <- popcheck$adjfactors
    adjvarlst <- popcheck$adjvarlst
    condsampcnt <- popcheck$condsampcnt
    dbqueries <- popcheck$dbqueries
    dbqueriesWITH <- popcheck$dbqueriesWITH
    estfromqry <- popcheck$estfromqry
    adjcase <- popcheck$adjcase
    pjoinid <- popcheck$pjoinid
    
    if(popType == "VOL") {
      treex <- popcheck$treex
      seedx <- popcheck$seedx
      tuniqueid <- popcheck$tuniqueid
      if (is.null(treex) && is.null(seedx)) {
        stop("must include tree data")
      }
    }    
  
  } else {
    stop("invalid popType")
  }
 
  ###################################################################################
  ## Return population data objects
  ###################################################################################
  estvar.area <- ifelse(adj == "none", "CONDPROP_UNADJ", "CONDPROP_ADJ")
  if (is.null(key(unitarea))) {
    setkeyv(unitarea, unitvar)
  }
  
  ###################################################################################
  ## Add new variables to pltcondx for estimation
  ###################################################################################
  if (returndata || savedata) {
    
    ## Get order of pltcondx columns
    pltcondxcols <- names(pltcondx)
    pltcondxkey <- key(pltcondx)
    newcols <- {}
    
    ## Add LANDSTATUSCD based on the following lookup table
    landstatuscdnm <- findnm("LANDSTATUSCD", pltcondxcols, returnNULL=TRUE)
    if (is.null(landstatuscdnm)) {
      condstatusnm <- findnm("COND_STATUS_CD", pltcondxcols, returnNULL=TRUE)
      reservcdnm <- findnm("RESERVCD", pltcondxcols, returnNULL=TRUE)
      siteclcdnm <- findnm("SITECLCD", pltcondxcols, returnNULL=TRUE)
      
      if (all(!sapply(c(condstatusnm, reservcdnm, siteclcdnm), is.null))) {
        lower <- ifelse (condstatusnm == "COND_STATUS_CD", FALSE, TRUE)
        landstatusnm <- ifelse(lower, "landstatus", "LANDSTATUS")
        
        LANDSTATUSlut <- data.frame(LANDSTATUS = c(101:108, 111:117),
                                    LANDSTATUSCD = c(rep(1, 6), rep(2, 2), rep(3, 6), 4),
                                    LANDSTATUSNM = c(rep("Timberland", 6), 
                                                     rep("Other forestland", 2), 
                                                     rep("Reserved productive forestland", 6),
                                                     "Reserved other forestland"))
        if (lower) names(LANDSTATUSlut) <- tolower(names(LANDSTATUSlut))
        
        pltcondx[[landstatusnm]] <- 
          with(pltcondx, get(condstatusnm) * 100 + get(reservcdnm) * 10 + get(siteclcdnm))
        pltcondx <- merge(pltcondx, LANDSTATUSlut, by=landstatusnm, all.x=TRUE)
        pltcondx[[landstatusnm]] <- NULL
        newcols <- c("LANDSTATUSCD", "LANDSTATUSNM")
        if (lower) newcols <- tolower(newcols)
        
        if (popType %in% c("CHNG", "GRM")) {
          prevnm <- ifelse(lower, "prev_", "PREV_")
          names(LANDSTATUSlut) <- paste0(prevnm, names(LANDSTATUSlut))
          
          pltcondx[[paste0(prevnm, landstatusnm)]] <- 
            with(pltcondx, get(paste0(prevnm, condstatusnm)) * 100 + 
                   get(paste0(prevnm, reservcdnm)) * 10 + get(paste0(prevnm, siteclcdnm)))
          pltcondx <- merge(pltcondx, LANDSTATUSlut, by=paste0(prevnm, landstatusnm), all.x=TRUE)
          pltcondx[[paste0(prevnm, landstatusnm)]] <- NULL
          newcols <- c(newcols, paste0(prevnm, newcols))
        }
      }
    }
    
    ## Add FORTYPGRPCD to pltcondx if not already in dataset
    fortypgrpnm <- findnm("FORTYPGRPCD", pltcondxcols, returnNULL=TRUE)

    if (is.null(fortypgrpnm)) {
      fortypnm <- findnm("FORTYPCD", pltcondxcols, returnNULL=TRUE)
      
      if (!is.null(fortypnm)) {
        lower <- ifelse (fortypnm == "FORTYPCD", FALSE, TRUE)
        
        ref_fortyp <- ref_codes[ref_codes$VARIABLE == "FORTYPCD", c("VALUE", "GROUPCD")]
        names(ref_fortyp) <- c("FORTYPCD", "FORTYPGRPCD")
        if (lower) names(ref_fortyp) <- tolower(names(ref_fortyp))
        
        pltcondx <- merge(pltcondx, ref_fortyp, by=fortypnm, all.x=TRUE)
        newcols <- c(newcols, ifelse(lower, "fortypgrpcd", "FORTYPGRPCD"))
        
        if (popType %in% c("CHNG", "GRM")) {
          prevnm <- ifelse(lower, "prev_", "PREV_")
          names(ref_fortyp) <- paste0(prevnm, names(ref_fortyp))
          
          pltcondx <- merge(pltcondx, ref_fortyp, by=paste0(prevnm, fortypnm), all.x=TRUE)
          newcols <- c(newcols, ifelse(lower, "prev_fortypgrpcd", "PREV_FORTYPGRPCD"))
        }  
      }
    }
    
    ## Add DSTRBGRP to pltcondx if not already in dataset
    dstrgrpnm <- findnm("DSTRBGRP", pltcondxcols, returnNULL=TRUE)
    
    if (is.null(dstrgrpnm)) {
      dstrbcd1nm <- findnm("DSTRBCD1", pltcondxcols, returnNULL=TRUE)
      
      ref_dstrbcd <- ref_codes[ref_codes$VARIABLE == "DSTRBCD", c("VALUE", "GROUPCD")]
      names(ref_dstrbcd) <- c("DSTRBCD1", "DSTRBGRP")
      if (lower) names(ref_dstrbcd) <- tolower(names(ref_dstrbcd))
      
      pltcondx <- merge(pltcondx, ref_dstrbcd, by=dstrbcd1nm, all.x=TRUE)
      newcols <- c(newcols, ifelse(lower, "dstrbgrp", "DSTRBGRP"))
      
      if (popType %in% c("CHNG", "GRM")) {
        prevnm <- ifelse(lower, "prev_", "PREV_")
        names(ref_dstrbcd) <- paste0(prevnm, names(ref_dstrbcd))
        
        pltcondx <- merge(pltcondx, ref_dstrbcd, by=paste0(prevnm, dstrbcd1nm), all.x=TRUE)
        newcols <- c(newcols, ifelse(lower, "prev_dstrbgrp", "PREV_DSTRBGRP"))
      }  
    }
    
    
    ## Move new columns to end of table
    setcolorder(pltcondx, c(pltcondxcols, newcols))
    pltcondflds <- c(pltcondflds, newcols)
    setkeyv(pltcondx, pltcondxkey)
  }

  ## Save pltids, including adjustment factors
  if (savepltids) {
    ## Add PROJECTID to pltassgnx
    if (!is.null(projectid)) {
      pltidsadj$PROJECTID <- projectid
    }
    
    message("saving pltids...")
    outlst$out_layer <- "pltids"
    if (!append_layer) index.unique.pltids <- c(projectid, puniqueid)
    datExportData(pltidsadj, 
                  savedata_opts = outlst)
  }
  
  
  
  ## Build list of data to return
  ###################################################################################
  returnlst$popType <- popType
  if(!is.null(bndx)) {
    returnlst$bndx <- bndx
  }
  returnlst <- append(returnlst, list(
    pltidsadj = pltidsadj, pltcondx=pltcondx, 
    pltcondflds = pltcondflds, pjoinid = pjoinid,
    cuniqueid = cuniqueid, pltassgnid = pltassgnid,
    condid = condid, ACI = ACI,
    areawt = areawt, areawt2 = areawt2, adjcase = adjcase,
    dbqueries = dbqueries, dbqueriesWITH = dbqueriesWITH,
    pltassgnx = pltassgnx, unitlut = data.table(unitlut),
    unitarea = unitarea, npixels = npixels,
    npixelvar = npixelvar, estvar.area = estvar.area,
    areavar = areavar, areaunits = areaunits, 
    unitvar = unitvar, unitvars = unitvars,
    plotsampcnt = plotsampcnt, condsampcnt = condsampcnt,
    states = states, invyrs = invyrs, adj = adj,
    P2POINTCNT = P2POINTCNT, plotunitcnt = plotunitcnt))
  
  if (popType == "VOL") {
    if (!is.null(treex)) {
      returnlst$treex <- treex
      returnlst$tuniqueid <- tuniqueid
      returnlst$adjtree <- adjtree
    }
    if (!is.null(seedx)) {
      returnlst$seedx <- seedx
    }
  }
  
  if (!is.null(popevalid)) {
    returnlst$evalid <- popevalid
  }
  if (adj != "none") {
    returnlst$adjfactors <- adjfactors
    returnlst$adjvarlst <- adjvarlst
  }

  returnlst$prednames <- prednames
  returnlst$predfac <- predfac
  
  
  ## Save data frames
  ##################################################################
  if (returndata) {
    returnlst$popdatindb <- FALSE
  } else {
    returnlst$popdatindb <- TRUE
    
    if (savedata) {
      if (outlst$out_fmt == "sqlite") {
        returnlst$pop_fmt <- "sqlite"
        returnlst$pop_dsn <- file.path(outlst$outfolder, outlst$out_dsn)
        returnlst$pop_schema <- NULL
      }
      
      message("saving pltassgnx...")
      outlst$out_layer <- "pltassgn"
      datExportData(pltassgnx, 
                    savedata_opts = outlst)
      
      message("saving unitarea...")
      outlst$out_layer <- "unitarea"
      datExportData(unitarea, 
                    savedata_opts = outlst)
      
      rm(pltassgnx)
      rm(unitarea)

      # if (popType %in% c("TREE", "GRM")) {
      #   message("saving REF_SPECIES...")
      #   outlst$out_layer <- "REF_SPECIES"
      #   datExportData(REF_SPECIES,
      #                 savedata_opts = outlst)
      # }
      
      
      if (!is.null(vcondsppx)) {
        message("saving vcondsppx...")
        outlst$out_layer <- "vcondsppx"
        datExportData(vcondsppx, 
                      savedata_opts = outlst)
        rm(vcondsppx)
        # gc()
      }
      if (!is.null(vcondstrx)) {
        message("saving vcondstrx...")
        outlst$out_layer <- "vcondstrx"
        datExportData(vcondstrx, 
                      savedata_opts = outlst)
        rm(vcondstrx)
        # gc()
      }
      
    } else if (datindb) {
      
      returnlst$pop_fmt <- datsource
      returnlst$pop_dsn <- dsn
      returnlst$pop_schema <- schema
      returnlst$popconn <- dbconn
    }
  }
  
  ## Save list object
  ##################################################################
  if (saveobj) {
    if (getext(objfn) == "rds") {
      message("saving list object to: ", objfn)
      saveRDS(returnlst, objfn)
    } else if (getext(objfn) == "rda") {
      message("saving list object to: ", objfn)
      save(returnlst, objfn)
    } else {
      message("invalid object name... must end in: ", toString(c("rds", "rda")))
    } 
  } 
  
  

  return(returnlst)
  
  
}