Skip to contents

Spatial Tools

Source: vignettes/FIESTA_tutorial_sp.Rmd
FIESTA_tutorial_sp.Rmd

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

  • Database tools (DB) - functions for querying and extracting data from FIA’s national database.
  • Data tools (dat) - functions for summarizing and exploring FIA data.
  • Spatial tools (sp) - functions for manipulating and summarizing spatial data.

Estimation Modules

  • Green-Book (GB) - functions for FIA’s standard ‘Green-Book’ estimators.
  • Photo-Based (PB) - functions for supplementary photo-based estimators.
  • Small Area (SA) - functions for integration with available small area estimators (SAE).
  • Model-Assisted (MA) - functions for integration with available Model-Assisted estimators.

Analysis Tools

  • Analysis tools (an) - wrapper functions for stream-lining estimation processes.

Overview of 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.

Objective of tutorial

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.

Example data - Wyoming (WY), Inventory Years 2011-2012

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.

Set up

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")

# County-level boundaries for USA and subset for Wyoming
WYco <- stunitco[stunitco$STATENM == "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 set our outfolder file path as a temporary directory. .

Examples

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.

View Example 
# 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.

View Example 
# Export Spatial Polygons layer to a shapefile
spExportSpatial(sfobj = WYbh, 
                savedata_opts = savedata_options(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".

View Example

We can use EPSG code with spMakeSpatialPoints.

WYspplt <- spMakeSpatialPoints(xyplt = WYplt, 
                               xy.uniqueid = "CN", 
                               xvar = "LON_PUBLIC", 
                               yvar = "LAT_PUBLIC", 
                               xy.crs = 4269)

We can look at the results to verify that WYspplt is now an sf points object with the desired projection:

WYspplt[ ,c("CN", "geometry")]

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 = savedata_options(out_dsn = "spplt",
                                                                out_fmt = "shp",
                                                                outfolder = outfolder,
                                                                out_layer = "WYplots",
                                                                overwrite_layer = TRUE))

spReprojectVector()

The spReprojectVector function re-projects a spatial vector layer. Note: The layer must have a defined coordinate reference system.

View Example 
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.

View Example 
# 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_polyv

Next we check and display the output:

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 tables 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 = savedata_options(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()

The spClipPoly function clips (intersects) a polygon vector layer with another polygon vector layer.
View Example

We first subset 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.

plot(sf::st_geometry(WYbhco['COUNTYCD']), col = sf::sf.colors(nrow(WYbhco)))
coords <- sf::st_coordinates(sf::st_centroid(sf::st_geometry(WYbhco)))
text(coords[,"X"], coords[,"Y"], WYbhco[["COUNTYCD"]])

spClipRast()

The spClipRast function subsets a raster to a polygon extent or boundary.

View Example

First, we subset the strata layer with the Medicine Wheel district boundary

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(rast = 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.

View Example

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

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

head(WYspplt_bh2)

spExtractRast()

The spExtractRast function extracts values from one or more raster layers and appends to a SpatialPoints layer or data frame.

View Example

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

We can verify that the values from the raster layers were indeed added to our SpatialPoints layer.

head(WYspplt_dem[ , c("WYbighorn_forest_nonforest_250m", "WYbighorn_dem_250m")])

spGetAuxiliary()

The spGetAuxiliary function extracts data and zonal statistics for model-assisted or model-based (small area) estimation. The major steps are as follows:

  1. Check parameters
  2. Extract point values from dunit_layer
  3. Set up output data structures
  4. Extract point values and get zonal statistics from continuous raster layers
  5. Extract point values and get zonal statistics from categorical raster layers
  6. Get total acres from domlayer (if areacalc = TRUE)
View Example

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 call spGetAuxiliary with out points object, boundary layer, and raster files.

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)

We start by getting a sense for all that spGetAuxiliary outputs by looking at the names of the modeldat object:

names(modeldat)

The pltassgn object returned by this function contains all of the necessary plot identifiers as well as the extracted values from each of the auxiliary rasters.

head(modeldat$pltassgn)

unitzonal contains zonal statistics for each domain in the input layer. In this case we are treating the entire layer as a single domain so we just get one row of output.

head(modeldat$unitzonal)

unitarea contains area values by domain in the input layer.

head(modeldat$unitarea)

Finally, we can look at the resulting predictor names (prednames).

prednames <- modeldat$prednames

spGetXY()

The spGetXY function extracts XY data within a given boundary.

View Example

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(WYbhprj))
plot(sf::st_geometry(spxy), add=TRUE)

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.

View Example

We can extract plots using the public coordinates within WY Bighorn NF boundary. In this example we utilize the FIA datamart to pull all of the necessary tables in and draw plots from the most recent FIA evaluation.

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.

View Example

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 <- WYspplt

By Bighorn National Forest District

We now get acres and strata information for Bighorn National Forest using the districts as our estimation units.

unitdat.bhdist <- spGetEstUnit(xyplt = WYplt, 
                               uniqueid = "CN", 
                               unit_layer = WYbhdistfn, 
                               unitvar = "DISTRICTNA", 
                               spMakeSpatial_opts = list(xvar = "LON_PUBLIC", 
                                                         yvar = "LAT_PUBLIC", 
                                                         xy.crs = 4269))

head(unitdat.bhdist$unitarea)
head(unitdat.bhdist$pltassgn[ ,c("CN", "DISTRICTNA")])

spGetStrata()

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.

View Example

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$stratalut

spUnionPoly()

The spUnionPoly function generates a unioned sf object with polygons and attributes from two sf polygon objects.

View Example

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 = stunitco[stunitco$STATENM == "Wyoming", ], 
                         polyv2 = WYbh, 
                         areacalc = TRUE)

plot(sf::st_geometry(polyUnion))
head(polyUnion)

spZonalRast()

The spZonalRast function extracts summary statistics by polygon (i.e., zone).

View Example

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
head(zonalext)