## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(warning = FALSE, message = FALSE)

## ----include=FALSE------------------------------------------------------------
# Sets up output folding
hooks = knitr::knit_hooks$get()
hook_foldable = function(type) {
  force(type)
  function(x, options) {
    res = hooks[[type]](x, options)
    
    if (isFALSE(options[[paste0("fold.", type)]])) return(res)
    
    paste0(
      "<details><summary>", type, "</summary>\n\n",
      res,
      "\n\n</details>"
    )
  }
}
knitr::knit_hooks$set(
  output = hook_foldable("output"),
  plot = hook_foldable("plot")
)

## ----echo=-1------------------------------------------------------------------
data.table::setDTthreads(2)

## ----warning = F, message = F-------------------------------------------------
library(FIESTA)

## -----------------------------------------------------------------------------
outfolder <- tempdir()

## -----------------------------------------------------------------------------
# File names for external spatial data
WYbhfn <- system.file("extdata", "sp_data/WYbighorn_adminbnd.shp", package="FIESTA")
WYbhdistfn <- system.file("extdata", "sp_data/WYbighorn_districtbnd.shp", package="FIESTA")
WYbhdist.att <- "DISTRICTNA"

fornffn <- system.file("extdata", "sp_data/WYbighorn_forest_nonforest_250m.tif", package="FIESTA")
demfn <- system.file("extdata", "sp_data/WYbighorn_dem_250m.img", package="FIESTA")

# Derive new predictor layers from dem
library(terra)

dem <- rast(demfn)

slpfn <- paste0(outfolder, "/WYbh_slp.img")
slp <- terra::terrain(dem,
                      v = "slope",
                      unit = "degrees",
                      filename = slpfn, 
                      overwrite = TRUE, 
                      NAflag = -99999.0)

aspfn <- paste0(outfolder, "/WYbh_asp.img")
asp <- terra::terrain(dem,
                      v = "aspect",
                      unit = "degrees", 
                      filename = aspfn,
                      overwrite = TRUE, 
                      NAflag = -99999.0)

## -----------------------------------------------------------------------------
smallbnd <- WYbhdistfn
smallbnd.domain <- "DISTRICTNA"

## -----------------------------------------------------------------------------
SApltdat <- spGetPlots(bnd = WYbhdistfn,
                       xy_datsource = "obj",
                       xy = WYplt,
                       xy_opts = xy_options(xy.uniqueid = "CN",
                                            xvar = "LON_PUBLIC", 
                                            yvar = "LAT_PUBLIC",
                                            xy.crs = 4269),
                       datsource = "obj",
                       dbTabs = dbTables(plot_layer = WYplt,
                                         cond_layer = WYcond,
                                         tree_layer = WYtree, 
                                         seed_layer = WYseed),
                       eval = "custom",
                       eval_opts = eval_options(invyrs = 2011:2013),
                       showsteps = TRUE,
                       returnxy = TRUE,
                       savedata_opts = savedata_options(outfolder = outfolder))

## -----------------------------------------------------------------------------
str(SApltdat, max.level = 1)

## ----results='hide'-----------------------------------------------------------
rastlst.cont <- c(demfn, slpfn, aspfn)
rastlst.cont.name <- c("dem", "slp", "asp")
rastlst.cat <- fornffn
rastlst.cat.name <- "fornf"

unit_layer <- WYbhdistfn
unitvar <- "DISTRICTNA"

auxdat <- spGetAuxiliary(xyplt = SApltdat$spxy,
                         uniqueid = "PLT_CN",
                         unit_layer = unit_layer,
                         unitvar = "DISTRICTNA",
                         rastlst.cont = rastlst.cont,
                         rastlst.cont.name = rastlst.cont.name,
                         rastlst.cont.stat = "mean",
                         rastlst.cont.NODATA = 0,
                         rastlst.cat = rastlst.cat,
                         rastlst.cat.name = rastlst.cat.name,
                         asptransform = TRUE,
                         rast.asp = aspfn,
                         keepNA = FALSE,
                         showext = FALSE,
                         savedata = FALSE)

## -----------------------------------------------------------------------------
str(auxdat, max.level = 1)

## -----------------------------------------------------------------------------
SApopdat <- modSApop(pltdat = SApltdat, 
                     auxdat = auxdat,
                     smallbnd = WYbhdistfn,
		                 smallbnd.domain = smallbnd.domain)

## -----------------------------------------------------------------------------
str(SApopdat, max.level = 1)

## -----------------------------------------------------------------------------
all_preds <- c("slp", "dem", "asp_cos", "asp_sin", "fornf")

## -----------------------------------------------------------------------------
area1 <- modSAarea(SApopdatlst = SApopdat,        # pop - population calculations for WY, post-stratification
                   prednames = all_preds,         # est - character vector of predictors to be used in the model
                   SApackage = "JoSAE",           # est - character string of the R package to do the estimation
                   SAmethod = "unit",             # est - method of small area estimation. Either "unit" or "area"
                   multest = FALSE)               # est - whether to also run all other available small area estimators

## -----------------------------------------------------------------------------
area1$est

## -----------------------------------------------------------------------------
str(area1$raw, max.level = 1)

## -----------------------------------------------------------------------------
area2 <- modSAarea(SApopdatlst = SApopdat,   # pop - population calculations for WY, post-stratification
                   prednames = "dem",        # est - character vector of predictors to be used in the model
                   SApackage = "JoSAE",      # est - character string of the R package to do the estimation
                   SAmethod = "area",        # est - method of small area estimation. Either "unit" or "area"
                   multest = TRUE)           # est - whether to also run all other available small area estimators

## -----------------------------------------------------------------------------
area2$est

## -----------------------------------------------------------------------------
area2$multest

## -----------------------------------------------------------------------------
area3 <- modSAarea(
  SApopdatlst = SApopdat,   # pop - population calculations for WY, post-stratification
  prednames = all_preds,    # est - character vector of predictors to be used in the model
  SApackage = "hbsae",      # est - character string of the R package to do the estimation
  SAmethod = "unit",        # est - method of small area estimation. Either "unit" or "area"
  multest = TRUE
  )


## -----------------------------------------------------------------------------
area3$est
area3$raw$SAmethod
area3$raw$SApackage

## -----------------------------------------------------------------------------
area4 <- modSAarea(
  SApopdatlst = SApopdat,     # pop - population calculations for WY, post-stratification
  prednames = all_preds,      # est - character vector of predictors to be used in the model
  SApackage = "hbsae",        # est - character string of the R package to do the estimation
  SAmethod = "unit",          # est - method of small area estimation. Either "unit" or "area"
  na.fill = "DIR",
  prior = function(x) 1       # est - prior on ratio of between and within area variation
  )


## -----------------------------------------------------------------------------
area3$est
area4$est

## -----------------------------------------------------------------------------


## -----------------------------------------------------------------------------
area5 <- modSAarea(
  SApopdatlst = SApopdat,      # pop - population calculations for WY, post-stratification
  prednames = all_preds,       # est - character vector of predictors to be used in the model
  SApackage = "JoSAE",         # est - character string of the R package to do the estimation
  SAmethod = "unit",           # est - method of small area estimation. Either "unit" or "area"
  modelselect = TRUE           # est - elastic net variable selection
  )

## -----------------------------------------------------------------------------
area5$est
area5$raw$predselect.unit

## -----------------------------------------------------------------------------
estvar <- "VOLCFNET"
live_trees <- "STATUSCD = 1"

## -----------------------------------------------------------------------------
tree1 <- modSAtree(
    SApopdatlst = SApopdat,      # pop - population calculations for WY, post-stratification
    prednames = all_preds,       # est - character vector of predictors to be used in the model
    SApackage = "JoSAE",         # est - character string of the R package to do the estimation
    SAmethod = "unit",           # est - method of small area estimation. Either "unit" or "area" 
    landarea = "FOREST",         # est - forest land filter
    estvar = estvar,             # est - net cubic-foot volume
    estvar.filter = live_trees   # est - live trees only
    )


## -----------------------------------------------------------------------------
tree1$est
tree1$multest

## -----------------------------------------------------------------------------
tree2 <- modSAtree(
    SApopdatlst = SApopdat,      # pop - population calculations for WY, post-stratification
    prednames = all_preds,       # est - character vector of predictors to be used in the model
    SApackage = "JoSAE",         # est - character string of the R package to do the estimation
    SAmethod = "unit",           # est - method of small area estimation. Either "unit" or "area"  
    landarea = "FOREST",         # est - forest land filter
    estvar = estvar,             # est - net cubic-foot volume
    estvar.filter = live_trees,   # est - live trees only
    modelselect = TRUE
    )

## -----------------------------------------------------------------------------
tree2$raw$predselect.unit
tree2$est

## -----------------------------------------------------------------------------
tree3 <- modSAtree(
    SApopdatlst = SApopdat,      # pop - population calculations for WY, post-stratification
    prednames = all_preds,       # est - character vector of predictors to be used in the model
    SApackage = "JoSAE",         # est - character string of the R package to do the estimation
    SAmethod = "unit",           # est - method of small area estimation. Either "unit" or "area"  
    landarea = "FOREST",         # est - forest land filter
    estvar = "DRYBIO_AG",               # est - net cubic-foot volume
    estvar.filter = live_trees,  # est - live trees only
    returntitle = TRUE
    )

## -----------------------------------------------------------------------------
tree3$est

## -----------------------------------------------------------------------------
tree3$titlelst

## -----------------------------------------------------------------------------
tree4 <- modSAtree(
    SApopdatlst = SApopdat,      # pop - population calculations for WY, post-stratification
    prednames = "dem",       # est - character vector of predictors to be used in the model
    SApackage = "sae",         # est - character string of the R package to do the estimation
    SAmethod = "area",           # est - method of small area estimation. Either "unit" or "area"  
    landarea = "FOREST",         # est - forest land filter
    estvar = "DRYBIO_AG",               # est - net cubic-foot volume
    estvar.filter = live_trees,  # est - live trees only
    returntitle = TRUE
    )

## -----------------------------------------------------------------------------
tree4$est