FIESTA Overview
The R package, FIESTA (Forest Inventory ESTimation and
Analysis) is a research estimation tool for analysts that work with
sample-based inventory data from the U.S. Department of Agriculture,
Forest Service, Forest Inventory and Analysis (FIA) Program to
accommodate: unique population boundaries, different evaluation time
periods, customized stratification schemes, non-standard variance
equations, integration of multi-scale remotely-sensed data and other
ancillary information, and interaction with other modeling and
estimation tools from CRAN R’s library of packages. FIESTA
contains a collection of functions that can access FIA databases,
summarize and compile plot and spatial data, and generate estimates with
associated sampling errors.
Functions are organized by type or objective and are named with a corresponding prefix:
Core Functions
DB) - functions for querying and
extracting data from FIA’s national database.dat) - functions for summarizing and
exploring FIA data.sp) - functions for manipulating and
summarizing spatial data.Estimation Modules
GB) - functions for FIA’s standard
‘Green-Book’ estimators.PB) - functions for supplementary
photo-based estimators.SA) - functions for integration with
available small area estimators (SAE).MA) - functions for integration with
available Model-Assisted estimators.Analysis Tools
an) - wrapper functions for
stream-lining estimation processes.FIESTA spatial (sp) tools
FIESTA’s spatial tools allow for summarizing and
manipulating spatial data for use in FIESTA’s estimation
modules.
| FUNCTION | DESCRIPTION |
|---|---|
| spImportSpatial() | Imports a spatial layer to a sf object. |
| spExportSpatial() | Exports a sf object to a spatial layer. |
| spMakeSpatialPoints() | Generates an S4 SpatialPoints object from X/Y coordinates. |
| spReprojectSpatial() | Reprojects a sf object. |
| spClipPoint() | Subset points with a Spatial polygon layer. |
| spClipPoly() | Subset Spatial polygons layer with another Spatial polygons layer. |
| spClipRast() | Subset raster layer with a Spatial polygons layer. |
| spExtractPoly() | Extracts point attribute values from Spatial polygons layer(s). |
| spExtractRast() | Extracts point attribute values from raster layer(s). |
| spGetXY() | Wrapper: Extracts XY coordinates and subsets to boundary. |
| spGetPlots() | Wrapper: Extracts plot data and subsets to boundary. |
| spGetAuxiliary() | Wrapper: Extracts point attribute values, area for estimation unit(s), and zonal statistics for strata predictor layers. |
| spGetEstUnit() | Wrapper: Extracts point attribute values and area for estimation unit(s). |
| spGetStrata() | Wrapper: Extracts point attribute values, area for estimation unit(s), and pixel counts for strata spatial layer. |
| spUnionPoly() | Generates one Spatial polygons object from two Spatial polygons layers. |
| spZonalRast() | Extracts summary statistics by polygon (i.e., zone) for a raster. |
The objective of this tutorial is to demonstrate how to use
FIESTA’s spatial tools for manipulating and summarizing
spatial data. The examples use data from three inventory years of field
measurements in the state of Wyoming, from FIADB_1.7.2.00, last updated
June 20, 2018, downloaded on June 25, 2018 and stored as internal data
objects in FIESTA.
| Data Frame | Description |
|---|---|
| WYplt | WY plot-level data |
| WYcond | WY condition-level data |
| WYtree | WY tree-level data |
| External data | Description |
|---|---|
| WYbighorn_adminbnd.shp | Polygon shapefile of WY Bighorn National Forest Administrative boundary* |
| WYbighorn_districtbnd.shp | Polygon shapefile of WY Bighorn National Forest District boundaries** |
| WYbighorn_forest_nonforest_250m.tif | GeoTIFF raster of predicted forest/nonforest (1/0) for stratification*** |
| WYbighorn_dem_250m.img | Erdas Imagine raster of elevation change, in meters**** |
*USDA Forest Service, Automated Lands Program (ALP). 2018. S_USA.AdministrativeForest (). Description: An area encompassing all the National Forest System lands administered by an administrative unit. The area encompasses private lands, other governmental agency lands, and may contain National Forest System lands within the proclaimed boundaries of another administrative unit. All National Forest System lands fall within one and only one Administrative Forest Area.
**USDA Forest Service, Automated Lands Program (ALP). 2018. S_USA.RangerDistrict (http://data.fs.usda.gov/geodata/edw). Description: A depiction of the boundary that encompasses a Ranger District.
***Based on MODIS-based classified map resampled from 250m to 500m resolution and reclassified from 3 to 2 classes: 1:forest; 2:nonforest. Projected in Albers Conical Equal Area, Datum NAD27 (Ruefenacht et al. 2008). Clipped to extent of WYbighorn_adminbnd.shp.
****USGS National Elevation Dataset (NED), resampled from 30m resolution to 250m. Projected in Albers Conical Equal Area, Datum NAD27 (U.S. Geological Survey 2017). Clipped to boundary of WYbighorn_adminbnd.shp.
First, you’ll need to load the FIESTA library:
Next, you’ll have to load some external data from the
FIESTA package and set up objects in your R
global environment:
# 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")
# Other spatial layers used for examples, extracted using the geodata package, gadm function.
# County-level boundaries for USA and subset for Wyoming (Note: must have internet connection)
USAco <- geodata::gadm(country = "USA", level = 2, path=tempdir())
WYco <- USAco[USAco$NAME_1 == "Wyoming",]Finally, you’ll need to set up an “outfolder”. This is just a file
path to a folder where you’d like FIESTA to send your data
output. For this vignette, we have saved our outfolder file path as the
outfolder object.
In the remainder of this vignette, we provide examples of using the
functions in the sp portion of FIESTA.
spImportSpatial()
The spImportSpatial function imports a spatial layer
(e.g., ESRI Shapefile) to a simple feature (sf) object.
## Import external data shapefiles
WYbh <- spImportSpatial(WYbhfn)
WYbhdist <- spImportSpatial(WYbhdistfn)
## Display boundary
plot(sf::st_geometry(WYbhdist), border="blue")
plot(sf::st_geometry(WYbh), add=TRUE)spExportSpatial()
The spExportSpatial function exports an sf
object.
## Export Spatial Polygons layer to a shapefile
spExportSpatial(WYbh,
savedata_opts = list(out_dsn = "WYbh.shp",
outfolder = outfolder,
overwrite_dsn = TRUE)
)spMakeSpatialPoints()
The spMakeSpatialPoints function generates an
sf points object with a defined projection. Note: The
coordinate reference system (crs) is: prj = "longlat" and
datum = "NAD83".
We can use EPSG code with spMakeSpatialPoints.
WYspplt <- spMakeSpatialPoints(xyplt = WYplt,
xy.uniqueid = "CN",
xvar = "LON_PUBLIC",
yvar = "LAT_PUBLIC",
xy.crs = 4269
)To view our output we can run the following code.
## Display output
plot(sf::st_geometry(WYbhdist))
plot(sf::st_geometry(WYspplt), add=TRUE)
## NOTE: To display multiple layers, all layers must be in the same coordinate system.
lapply(list(WYbh, WYbhdist, WYspplt), sf::st_crs)And finally we export the spatial points to an outfolder.
WYspplt <- spMakeSpatialPoints(xyplt = WYplt,
xy.uniqueid = "CN",
xvar = "LON_PUBLIC",
yvar = "LAT_PUBLIC",
xy.crs = 4269,
exportsp = TRUE,
savedata_opts = list(
out_dsn = "spplt",
out_fmt = "shp",
outfolder = outfolder,
out_layer = "WYplots",
overwrite_layer = TRUE)
)spReprojectVector()
The spReprojectVector function reprojects a spatial
vector layer. Note: The layer must have a defined coordinate reference
system (test using sf::st_crs).
sf::st_crs(WYspplt)
prj <- "+proj=utm +zone=12 +ellps=GRS80 +datum=NAD83 +units=m +no_defs"
WYspplt.utm12 <- spReprojectVector(layer = WYspplt,
crs.new = prj)
## Check results
sf::st_crs(WYspplt.utm12)spClipPoint()
The spClipPoint function subsets
SpatialPoints or point file with a Spatial polygon
boundary.
## Get points within Bighorn National Forest boundary (project on the fly)
WYbhptslst <- spClipPoint(xyplt = WYplt,
uniqueid = "CN",
clippolyv = WYbh,
spMakeSpatial_opts=list(xvar = "LON_PUBLIC",
yvar = "LAT_PUBLIC",
xy.crs = 4269)
)
WYbhptslst <- spClipPoint(xyplt = WYplt,
uniqueid = "CN",
clippolyv = WYbh,
savedata = TRUE,
exportsp = TRUE,
spMakeSpatial_opts=list(xvar = "LON_PUBLIC",
yvar = "LAT_PUBLIC",
xy.crs = 4269),
savedata_opts = list(outfolder=outfolder,
out_layer = "WYbh")
)
names(WYbhptslst)
WYbhspplt <- WYbhptslst$clip_xyplt
WYbhprj <- WYbhptslst$clip_polyvNext we check and display the output:
WYbhspplt
plot(sf::st_geometry(WYbhprj), border="red", lwd=2)
plot(sf::st_geometry(WYbhspplt), add=TRUE)Note: If the projection of spplt is not the same as the points layer,
the points layer will be reprojected to the same projection as
clippolyv (See notes in help file for more info).
Now we generate a sf object first, then clip sf points
with the spClipPoint function:
WYspplt <- spMakeSpatialPoints(xyplt = WYplt,
xy.uniqueid = "CN",
xvar = "LON_PUBLIC",
yvar = "LAT_PUBLIC",
xy.crs = 4269)
WYbhptslst <- spClipPoint(xyplt = WYspplt,
uniqueid = "CN",
clippolyv = WYbh)We can also subset other tales with clipped points.
WYbhptslst <- spClipPoint(xyplt = WYspplt,
uniqueid = "CN",
clippolyv = WYbh,
othertabnms = c("WYcond", "WYtree"))
names(WYbhptslst)Next, we can export clipped points.
## Export clipped points
spExportSpatial(WYbhspplt,
savedata_opts = list(out_layer = "WYbhpts",
outfolder = outfolder,
overwrite_layer = TRUE)
)Finally, we can get points within Bighorn National Forest boundary (project on the fly) and save to an outfolder.
WYbhptslst <- spClipPoint(xyplt = WYplt,
uniqueid = "CN",
clippolyv = WYbh,
othertabnms = c("WYcond", "WYtree"),
exportsp = TRUE,
spMakeSpatial_opts=list(xvar = "LON_PUBLIC",
yvar = "LAT_PUBLIC",
xy.crs = 4269),
savedata_opts = list(
outfolder = outfolder,
overwrite_layer = TRUE,
outfn.pre = "clip")
)
names(WYbhptslst)spClipPoly()
spClipPoly function clips (intersects) a polygon vector
layer with another polygon vector layer.
We first subsect the WYco SpatialPolygons layer with
WYbighorn. Note: you must download USAco from the
geodata package and subset to Wyoming. See the “Set Up”
section above.
WYbhco <- spClipPoly(polyv = WYco,
clippolyv = WYbh)We can now check and display WYbighorn_co with
labels.
head(WYbhco)
plot(sf::st_geometry(WYbhco['NAME_2']), col = sf::sf.colors(nrow(WYbhco)))
coords <- sf::st_coordinates(sf::st_centroid(sf::st_geometry(WYbhco)))
text(coords[,"X"], coords[,"Y"], WYbhco[["NAME_2"]])spClipRast()
The spClipRast function subsets a raster to a polygon
extent or boundary.
First, we subset the strata layer with the Medicine Wheel district boundary
WYbhdist
WYbhMW <- WYbhdist[WYbhdist$DISTRICTNA == "Medicine Wheel Ranger District",]
plot(sf::st_geometry(WYbhdist))
plot(sf::st_geometry(WYbhMW), border="red", add=TRUE)Next, we can clip the raster. Note: If the projection of
polyv is not the same as rast,
polyv will by reprojected to the same projection as
rast before clipping (See note in the help file for more
details).
WYbhMW.fornf <- spClipRast(fornffn,
clippolyv = WYbhMW,
outfolder = outfolder)Finally, we display the results:
WYbhMWprj <- FIESTAutils::crsCompare(WYbhMW,
FIESTAutils::rasterInfo(WYbhMW.fornf)$crs)$x
terra::plot(terra::rast(WYbhMW.fornf))
plot(sf::st_geometry(WYbhMWprj),
border = "red",
add = TRUE)spExtractPoly()
The spExtractPoly function subsets a
SpatialPolygons layer with another
SpatialPolygons layer.
First, we extract polygon attributes from WYbighorn to
WYpts, keeping NULL values. Note: If the
projection of spplt is not the same as the
SpatialPoints, the SpatialPoints layer will be
reprojected to the same projection as clippolyv before the
evaluation points (See the note in help file for more details).
extpolylst <- spExtractPoly(WYspplt,
xy.uniqueid = "CN",
polyvlst = WYbhdist)
WYspplt_bh <- extpolylst$spxyext
dim(WYspplt)
dim(WYspplt_bh)
head(WYspplt_bh)
plot(WYspplt_bh["DISTRICTNA"], pch = 8)Next we can extract a subset of polygon attributes from
WYbighorn to WYpts, not keeping
NULL values.
extpolylst <- spExtractPoly(WYspplt,
xy.uniqueid = "CN",
polyvlst = WYbh,
polyvarlst = c("FORESTNUMB", "FORESTNAME"),
keepNA = FALSE)
WYspplt_bh2 <- extpolylst$spxyext
dim(WYspplt)
dim(WYspplt_bh2)
head(WYspplt_bh2)spExtractRast()
The spExtractRast function subsets a
SpatialPolygons layer with another
SpatialPolygons layer.
We can extract raster values from WYdem to
WYpts, not keeping NULL values.
extrastlst <- spExtractRast(WYspplt,
rastlst = c(fornffn, demfn),
xy.uniqueid = "CN",
keepNA = FALSE)
WYspplt_dem <- extrastlst$sppltext
dim(WYspplt)
dim(WYspplt_dem)
head(WYspplt_dem)spGetAuxiliary()
The spGetAuxiliary function extracts data and zonal
statistics for model-assisted or model-based (small area) estimation.
The major steps are as follows:
dunit_layerdomlayer (if
areacalc = TRUE)First we generate an sf object from the WY plot data,
public coordinates (xvar = "LON_PUBLIC",
yvar = "LAT_PUBLIC"). Note: the public coordinate
projection information is: prj = "longlat",
datum = "NAD83".
WYspplt <- spMakeSpatialPoints(xyplt = WYplt,
xy.uniqueid = "CN",
xvar = "LON_PUBLIC",
yvar = "LAT_PUBLIC",
xy.crs = 4269)Next, we derive new 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)Finally, we extract estimation unit and zonal raster statistics (i.e., mean).
rastlst.cont <- c(demfn, slp, asp)
rastlst.cont.name <- c("dem", "slp", "asp")
rastlst.cat <- fornffn
rastlst.cat.name <- "fornf"
modeldat <- spGetAuxiliary(xyplt = WYspplt,
uniqueid = "CN",
unit_layer = WYbhfn,
unitvar = NULL,
rastlst.cont = rastlst.cont,
rastlst.cont.name = rastlst.cont.name,
rastlst.cat = rastlst.cat,
rastlst.cat.name = rastlst.cat.name,
rastlst.cont.stat = "mean",
asptransform = TRUE,
rast.asp = asp,
keepNA = FALSE,
showext = FALSE,
savedata = FALSE)
names(modeldat)
pltassgn <- modeldat$pltassgn
unitzonal <- modeldat$unitzonal
unitvar <- modeldat$dunitvar
inputdf <- modeldat$inputdf
unitarea <- modeldat$unitarea
areavar <- modeldat$areavar
inputdf <- modeldat$inputdf
prednames <- modeldat$prednames
zonalnames <- modeldat$zonalnames
unitvar
areavar
unitzonal
unitarea
head(pltassgn)
prednames
zonalnamesspGetXY()
The spGetXY function extracts XY data within a given
boundary.
We can extract public coordinates within WY Bighorn NF boundary, returning spatial (spxy) and nonspatial plot identifiers (pltids).
WYbhxy <- spGetXY(bnd = WYbhfn,
xy_datsource = "datamart",
eval = "FIA",
eval_opts = eval_options(Cur = TRUE),
returnxy = TRUE)
names(WYbhxy)
pltids <- WYbhxy$pltids
head(pltids)
spxy <- WYbhxy$spxy
plot(sf::st_geometry(spxy))Or, returning only nonspatial plot identifiers (pltids).
WYbhxyids <- spGetXY(bnd = WYbhfn,
xy_datsource = "datamart",
eval = "FIA",
eval_opts = eval_options(Cur = TRUE),
returnxy = FALSE)
names(WYbhxyids)
pltids <- WYbhxyids$pltids
head(pltids)spGetPlots()
The spGetPlots function extracts plot data within a
given boundary.
We can extract public coordinates within WY Bighorn NF boundary.
WYbhdat <- spGetPlots(bnd = WYbhfn,
states = "Wyoming",
datsource = "datamart",
eval = "FIA",
eval_opts = eval_options(Cur = TRUE),
istree = FALSE)
names(WYbhdat)spGetEstUnit()
The spGetEstUnit function extracts estimation unit
values to plots and calculate area by estimation unit.
First, we create a SpatialPoints object from
WYplt
WYspplt <- spMakeSpatialPoints(xyplt = WYplt,
xy.uniqueid = "CN",
xvar = "LON_PUBLIC",
yvar = "LAT_PUBLIC",
xy.crs = 4269)
xyplt <- WYsppltWe now get estimation unit acres for Bighorn National Forest.
unitdat.bh <- spGetEstUnit(xyplt = WYplt,
uniqueid = "CN",
unit_layer = WYbhfn,
spMakeSpatial_opts=list(xvar = "LON_PUBLIC",
yvar = "LAT_PUBLIC",
xy.crs = 4269)
)
names(unitdat.bh)
unitarea.bh <- unitdat.bh$unitarea
unitvar.bh <- unitdat.bh$unitvar
areavar.bh <- unitdat.bh$areavar
unitarea.bh
unitvar.bh
areavar.bhWe now get estimation unit acres and strata information for Bighorn National Forest.
unitdat.bhdist <- spGetEstUnit(xyplt = WYplt,
uniqueid = "CN",
unit_layer = WYbhdistfn,
unitvar = "DISTRICTNA",
spMakeSpatial_opts=list(xvar = "LON_PUBLIC",
yvar = "LAT_PUBLIC",
xy.crs = 4269)
)
names(unitdat.bhdist)
unitarea.bhdist <- unitdat.bhdist$unitarea
unitvar.bhdist <- unitdat.bhdist$unitvar
areavar.bhdist <- unitdat.bhdist$areavar
unitarea.bhdist
unitvar.bhdist
areavar.bhdistspGetStrata()
The spGetStrata function is a wrapper to extract
attribute and area from a polygon or raster estimation unit layer and a
polygon or raster layer with strata pixel categories.
First, we generate an sf object from the WY plot data,
public coordinates (xvar = "LON_PUBLIC",
yvar = "LAT_PUBLIC"). Note: The public coordinate
projection information is: prj = "longlat",
datum = "NAD83".
WYspplt <- spMakeSpatialPoints(xyplt = WYplt,
xy.uniqueid = "CN",
xvar = "LON_PUBLIC",
yvar = "LAT_PUBLIC",
xy.crs = 4269)Next, we extract polygon attributes from WYbighorn to
WYpts, keeping NULL values.
stratlst <- spGetStrata(WYspplt,
uniqueid = "CN",
unit_layer = WYbhfn,
strattype = "RASTER",
strat_layer = fornffn)
names(stratlst)
stratlst$stratalutspUnionPoly()
The spUnionPoly function generates a unioned
sf object with polygons and attributes from two
sf polygon objects.
First, we import spatial data with spImportSpatial, and
then extract polygon attributes from WYbighorn to
WYpts, keeping NULL values.
WYbh <- spImportSpatial(WYbhfn)
polyUnion <- spUnionPoly(polyv1 = USAco[USAco$NAME_1 == "Wyoming",],
polyv2 = WYbh,
areacalc = TRUE)
plot(sf::st_geometry(polyUnion))
head(polyUnion)spZonalRast()
The spZonalRast function extracts summary statistics by
polygon (i.e., zone).
First, we import spatial data with spImportSpatial, and
then extract polygon attributes from WYbighorn to
WYpts, keeping NULL values.
WYbhdist <- spImportSpatial(WYbhdistfn)
zonallst <- spZonalRast(polyv = WYbhdist,
polyv.att = "DISTRICTNA",
rastfn = demfn,
zonalstat = c("mean", "sum"))
names(zonallst)
zonalext <- zonallst$zonalext
outname <- zonallst$outname
rasterfile <- zonallst$rasterfile
head(zonalext)
outname
rasterfile