Methods

Data inputs

The LANDFIRE Program provides 25+ national geo-spatial layers (e.g. vegetation, fuel, disturbance, etc.), databases, and ecological models that are available to the public for the U.S. and insular areas.

We relied on Three LANDFIRE Products as inputs for this initial analysis:

  1. Biophysical Settings (BpS) Models, Documents and Reference Condition tables: For estimating “reference” (historic) amounts of open- and closed late-succession classes.
  2. Biophysical Settings Spatial Data: Were used to map ‘ecosystem footprints’.
  3. Succession Class Spatial Data: To create a map of current succession classes. These are linked to the BpS products.

*NOTE: LANDFIRE data is delivered as geotifs with 30M pixel size. All maps were created in ArcGIS Pro. All charts were plotted in R/R-Studio.

Area of interest:

Show the code 
# load packages for page
library(sf)
library(terra)
library(tidyverse)
library(tmap)
library(janitor)
library(scales)
library(kableExtra)

#  read shape
MapZone2<- st_read("data/mz2.shp", quiet = TRUE) %>% 
  st_transform(crs = 5070) %>%
  st_union() %>%
  st_sf()

# toggle tmap mode to interactive viewing
tmap_mode("view")

# create a quick interactive map
quickmap <- qtm(MapZone2, 
                borders = "darkgreen", 
                fill = NULL, 
                check.and.fix = TRUE, 
                basemaps = 'https://server.arcgisonline.com/ArcGIS/rest/services/World_Topo_Map/MapServer/tile/{z}/{y}/{x}',
                title = 'LANDFIRE Map Zone 2',)

quickmap

Figure 2. Study area is the western coast of Oregon along the western side of the coastal range of the Cascade Mountains. This area is also known in LANDFIRE as Map Zone 2. See map of all LANDFIRE Map Zones here.

Below we describe the data inputs and GIS operations.

GIS operations

Clipping ecosystem footprint data

We relied on LANDFIRE’s Biophysical Settings (BpS) to obtain “potential ecosystem footprints” to represent what would have been present prior to European colonization of North America. LANDFIRE’s BpS suite is based on both the current biophysical environment and an approximation of the historical disturbance regime.

For example, each BpS, (i.e. “western hemlock douglas fir forest”) is accompanied by a description that includes information on its succession classes and a model which captures historical disturbances and estimated amounts of succession classes.

We used the Clip Raster tool in ArcGIS pro to extract the area of interest from a broader extent of the raster data.

Obtaining historical percentages of late-successional forests

For each BpS, we needed to get estimates how much late-open and late-closed forest would have historically been present. This information was developed in BpS state-and-transition models, and available in a couple formats:

  1. The BpS Reference Condition Table obtainable by clicking “Reference Conditions” here.
  2. BpS description documents. We downloaded them by filtering for Map Zone 2 at https://landfirereview.org/search.php.

Example - Below is a screenshot of a BpS model: 1. Boxes represent the succession classes, or structural stages of an ecosystem 2. Lines represent succession and disturbance pathways 3. What’s missing from the diagram? The quantitative, disturbance inputs from the model, or the outputs which include percent of each succession class that would have occurred historically (visit the LANDFIRE Vegetation Modeling site to learn more about the modeling process).

Assessing current succession classes

We used the LANDFIRE succession class dataset to estimate the percentage of each succession class that occurs (per BpS) on today’s landscape. Each pixel within a BpS is labeled “A through E” for the “reference” states, and includes labels for other categories such as agriculture and urban.

Step 1: Use the “clip raster” tool to extract the area of interest from a broader extent of the raster data.

Step 2: Do a ‘combine’ using the Combine Tool of the BpS and succession class data to get a dataframe with pixel counts of each unique combination of BpS and succession class to allow for estimating of the amounts of current succession classes per BpS.

Step 3: Use the Join Field Tool to join the ‘BPS_NAME’, ‘BPS_MODEL’ attributes from the BpS dataset, and the ‘Label’ attribute from the Succession Class dataset (see output data example below).

The rules for defining the structural characteristics of each succession class are in the BpS descriptions, and are essentially combinations of height, cover and type.

Wrangling data exported from ArcGIS Pro

Once we have clipped, combined and attributed spatial data in ArcGIS pro, we need to do a little data wrangling in R to further prep the spatial data for mapping and to make some charts. In summary we:

  1. Read the combined data attribute table and supporting data tables from LANDFIRE into R (all located in ‘data’ directory in Git repo)
    • BpS-Succession Class Combined attribute table (named “bps_scl_cmbn.csv”)
    • Succession Class Descriptions (named “scls_descriptions.csv”)
    • Reference condition percentages (named “ref_con_long.csv”)
    • Map Zone 2 Biophysical Settings (named “bps_mz2.csv”)
  2. Build and join in succession class ‘descriptions’ (e.g., ‘Late1-Closed’)
  3. Do a little cleaning, and make a new ‘join_field’
  4. Calculate reference and current acres per succession class per BpS

Click the “Code” button to see code.

Show the code 
## read in raw data and sclass descriptions

raw_bps_scls <- read.csv("data/bps_scls_cmbn.csv")

sclass_descriptions <- read.csv("data/scls_descriptions.csv")

reference_percents <- read.csv("data/ref_con_long.csv")

bps_mz2 <- read.csv("data/bps_mz2.csv")

##  clean and prep raw combined data

clean_bps_scls_cmbn <- raw_bps_scls %>%
  clean_names() %>%
  dplyr::select(-c(4, 5)) %>%
  unite("join_field", bps_model,label, sep = "_", remove = FALSE )

## clean and prep sclass descriptions

sclass_descriptions_clean <- sclass_descriptions %>%
  dplyr::select(-c(4)) %>% # remove column
  rename("model_code" = "StratumID",
         "scls_label" = "ClassLabelID",
         "state_class_id" = "StateClassID" ) %>% # rename columns
  unite("join_field", model_code:scls_label, sep = "_", remove = FALSE ) %>%
  separate(state_class_id, into = c("age_category", "canopy_category"), sep = ":", remove = FALSE) 


## clean and prep reference percents

unique_sclass_labels_ref <- unique(reference_percents$refLabel)
#print(unique_sclass_labels_ref)

unique_sclass_lables_cmbn <- unique(clean_bps_scls_cmbn$label)
#print(unique_sclass_lables_cmbn)
# does not have barren/sparse and there are differences, e.g., Urban-Developed between this and sclass label
# will assume Barren/Sparse, NoData and Snow/Ice is minimal; will change "Developed" to "Urban" in reference df cleaning code 

clean_ref_percents <- reference_percents %>%
  clean_names() %>%
  mutate(across('ref_label', str_replace, 'Developed', 'Urban')) %>%
  mutate(across('model_label', str_replace, 'Developed', 'Urban')) %>%
  rename("join_field" = "model_label",
         "bps_name" = "bp_s_name" )

## need to winnow this df to only the bps model codes in MZ 2

clean_ref_percents_mz2 <- clean_ref_percents %>%
  filter(model_code %in% bps_mz2$BPS_MODEL)


## create 'final' dataframe with reference and current sclass percents, acres and labels

## first ref con and sclass descriptions
final_df <- left_join(clean_ref_percents_mz2, sclass_descriptions_clean) 


# looks OK, now full join to add reference percents then clean a bit

final_df2 <- full_join(final_df, clean_bps_scls_cmbn) %>%
  dplyr::select(-c(6, 10, 13, 14)) %>%
  rename( "cur_scls_count" = "count")

# now for the math: need count/acres per bps, cur sclass percents and differences

final_df3 <- final_df2 %>%
  group_by(bps_name) %>%
  mutate(bps_count = sum(cur_scls_count, na.rm = TRUE)) %>%
  ungroup() %>%
  mutate(bps_acres = bps_count*0.2223945,
         ref_scls_acres = bps_acres*(ref_percent/100),
         cur_scls_acres = cur_scls_count*0.2223945,
         cur_percent = (cur_scls_acres/bps_acres)*100) %>%
  mutate(across(12:15, round, 0))

# save to csv to explore in Excel, and reorder columns

#write.csv(final_df3, file = "data/final_df.csv", row.names = FALSE)

The resulting final data table was then joined back to the combined BpS-Succession Class data in ArcGIS Pro for mapping.