The LANDFIRE map legend for existing vegetation and land cover encompasses the continuum from natural to ruderal and cultural vegetation types. The hierarchical structure of the USNVC includes broad units at upper levels defined by vegetation physiognomy, followed by progressively narrow units at lower levels defined by vegetation floristic composition (Federal Geographic Data Committee [FGDC] 2008). The full spectrum, from “natural” to “cultural” vegetation types are encompassed by the USNVC, but here we will reference use of “natural” to “ruderal” units. Table 1 provides an example of the USNVC hierarchy, with defining characteristics. Here, tallgrass prairie types have been well described at all levels of the hierarchy down to the association level, where multiple dominant and diagnostic species are used to define a given type. Over 6,000 associations describe natural vegetation types within the conterminous United States (CONUS). While this level of thematic detail is not currently feasible to map beyond relatively local scales, Group and Alliance levels are increasingly feasible to target in regional and national map legends. Within the CONUS, LANDFIRE mapped nearly 300 natural USNVC Group concepts. Descriptions of this classification hierarchy and these units may be found at the USNVC website (http://USNVC.org/) and NatureServe Explorer (https://explorer.natureserve.org/).

While the USNVC Group level provided a useful classification of natural vegetation units for the map legend, additional map legend categories were used to provide robust map product. First, the revised USNVC includes “ruderal” units that are defined as including plant assemblages with no natural analog. These commonly result from prior land conversion and subsequent abandonment, so they encompass what are often referred to as “old fields” and secondary forests where exotic species and/or native species are present in abundances not found where prior human influence is less discernable. Several ruderal vegetation units, approximating the Group level, were documented for use in LANDFIRE map products.

Second, a series of map class modifiers were used to facilitate mapping structural variants within each USNVC Group-based map class informed by the using the National Land Cover Database (NLCD; Homer et al. 2015). For example, where “evergreen” “deciduous” and “mixed” variants of a given forest type were discernable, they were mapped separately. Additionally, where “forest” vs. “shrubland” or “herbaceous” structural stages in forest succession occurred in discernable pattern, they were also differentiated in the map legend.

Modern mapping methods include use of georeferenced sample locations – each labeled to the intended map units they represent – to train models that will combine predictor layers to generate a vegetation map. Due the the very large number of georeferenced samples needed for national land cover mapping at thematic levels like the USNVC Group, LANDFIRE produced algorithmic tools called “autokeys” for processing vegetation sample plot data for subsequent modeling and mapping. The autokey algorithm scans the content of each sample plot to detect species presence and abundance as well as structural categories to determine to which map legend class the sample belongs. It then applies the appropriate label for use in subsequent modeling steps. Autokeys were designed and implemented within regions determined by clustering ecologically similar ecoregions modifed from US Forest Service (Cleland et al. 2007) and US Environmental Protection Agency (EPA) (Omernik and Griffith 2014) sources. In the CONUS project area 278 USNVC Groups were processed for use in the LANDFIRE legend. Each autokey pertains to one of 17 regions, including the southern tip of Florida which was treated along with adjacent Caribbean islands.

Expert ecologists reviewed and hand labeled nearly 18,000 samples to assess autokey performance. The total number of samples labeled by autokeys varied by region, from a high of 80,148 in the Rocky Mountains to 3,517 in the North Coast region. For most regions, the proportion of plots used in assessment that were reviewed by experts was 4% to 8%. Validation statistics for each map legend category were used throughout the development and final evaluation of each autokey. The overall validation statistic is a useful measure of how well each autokey performed across all types. It is the number of matches between expert and autokey labels divided by the total number of expert plots × 100, and the overall validation statistic was calculated for each autokey. Overall agreement for the USNVC Group keys ranged from a high of 90% (Texas-Oklahoma Hill Prairie) to a low of 39.9% (Coastal Plain). In most cases, lower performance occurred where substantial proportions of the landscape are dominated by ruderal vegetation, and distinguishing among very similar vegetation types using sample plots becomes more challenging.

Although over 500,000 vegetation samples were labeled through autokeys, there were still hundreds of thousands of samples with insufficient quantitative information to run through the autokeys. These often included documented locations from local natural resource inventories where an existing classification was used to label the location without including vegetation composition and structure. A series of classification crosswalks were used to reconcile these differences and label samples to the intended unit on the LANDFIRE map legend. In the 2016 LANDFIRE map, over one million samples were processed either by autokeys or through expert labeling across the CONUS.

A key factor in ecosystem modeling is determining the boundaries within which to build and apply models for existing vegetation types. Initial prototyping found that Omernik Level III Ecoregions (Omernik and Griffith 2014) provided a more ecologically meaningful framework for modeling existing vegetation types than map zones used in previous LANDFIRE mapping efforts (Picotte et al. 2019). The ecoregions were grouped to create 13 vegetation production units (Figure 1) across the CONUS that were of a manageable size for efficiently preparing satellite imagery and georeferenced sample inputs. These were similar to, but not identical to the ecoregion clusters used for design and implementation of autokeys.

The modeling process began by removing sample plots in recently disturbed areas collected prior to the disturbance using LANDFIRE annual disturbance products. Spectral outliers were identified by summing Landsat bands one through six for each class and sample plots, those plots greater than two standard deviations from the mean were removed. The resultant filtered plots were used to model lifeform and vegetation structure. The spectral test was performed separately for vegetation types prior to witholding samples for map validation.

Vegetation structure is important for fire behavior fuel models. Therefore, existing vegetation products were designed to nest by lifeform. For example, pixels identified as tree in the lifeform mask will be assigned a tree cover, height, and vegetation type. The lifeform modeling process began with an initial output using the filtered sample plots. The initial lifeform model output was improved through an iterative process by adding expert-labeled training samples based on desktop review of aerial photos to correct obvious mapping errors.

Using USNVC Group concepts as a guide, sample plots were separated into three lifeforms: tree, shrub, and herbaceous vegetation types, as well as barren or sparsely vegetated types (<10% total cover). Plots were further separated into wetland vs non-wetland categories, and alpine vs non-alpine categories where they existed. Classification tree models were generated with the See5 algorithm using raster predictor variables (Table 2). Models were generated within individual ecoregions to produce categorical outputs for each lifeform layer. The layers were then combined using the lifeform mask, wetland mask, barren/sparse mask, and alpine mask created with a separate modeling process to restrict the mapping of certain vegetation in appropriate locations where applicable. After modeling, vegetation types that comprised a mix of evergreen and deciduous dominant/co-dominant species to varying degrees were separated using the NLCD (Homer et al. 2015) categories to further refine the map and aid in assessing fire behavior fuel models. An overview of the modeling process is shown in Figure 2.

The draft thematic map was edited using rulesets based on geography or topography, or manual pixel reclassification with hand-drawn polygons based on expert opinion and review. Draft maps were also revised by removing problematic plots identified during the modeling process, reclassifying plots to a better fit, or adding sample plots based on expert opinion to correct modeling errors and improve mapping of problematic classes. Draft maps were reviewed by regional experts with NatureServe and staff from state agencies. Wildland Urban Interface maps produced by the Forest Service (Radeloff et al. 2017) were used in combination with rulesets based on existing vegetation to identify ruderal vegetation in proximity to developed areas. Recently disturbed areas (within previous 10-years) were identified using LANDFIRE disturbance products to appropriately label transitional vegetation. For example, it was more appropriate to label regenerating clearcuts in the Pacific Northwest as recently disturbed herbaceous cover rather than native montane grassland.

Vegetation percent cover and height training samples were derived from field estimates of vegetation cover and height. In addition, tree canopy percent cover estimates were calculated from the percentage of Light Detection and Ranging (LiDAR) point cloud above 3 m and tree height estimates were derived from the 90^th percentile of LiDAR returns. Sample plots were separated into three lifeforms: tree, shrub, and herbaceous cover and height. Regression tree models predicting percent cover and height of dominant vegetation were generated with the Cubist algorithm using the following predictor layers identified in Table 2: Seasonal Landsat imagery, tasseled cap, and NDVI 5-year statistics. Models were generated within a vegetation production unit to predict continuous outputs for each lifeform layer. The layers were then combined using the lifeform mask. Recently disturbed areas (within previous 3 years) did not model well because the satellite image composites spanned multiple years and comprised a mixture of pre and post-disturbance pixels. Disturbance severity and timing rulesets based on LANDFIRE annual disturbance products were used to assign lifeform and estimated cover.

Several masks were developed to identify open water, barren land, sparse vegetation, developed, and agricultural lands. Open water was identified using custom modeling methods based on the Landsat Level 3 Dynamic Surface Water Extent (DSWE) Science Product from the U.S. Geological Survey. Fragmented segments along streams and rivers were connected using National Hydrography Dataset (USGS 2016a) flowlines to form a continuous network. Barren land (0% vegetated cover) and sparse vegetation (<9% total cover) were identified using NDVI 5-year median thresholds and calibrated by location based on known barren and sparsely vegetated areas. Developed land and snow/ice were identified using NLCD. Ruderal vegetation classes were modeled within NLCD Class 21 (Developed – Open Space) in order to assign the appropriate fire behavior fuel models. Agricultural land was identified using the Cropland Data Layer (USDA 2015) and summarized by Common Land Unit polygons (USDA 2006) using a zonal majority.

The mapping process for existing vegetation based on the USNVC Group concepts generated a 30 m pixel resolution map raster with 499 natural, ruderal, and cultural map classes (Figure 3). The majority of these map classes are reflected in the more than 300 USNVC Groups. In addition to these, map class modifiers distinguish structural variants within each natural USNVC Group. Then a series of map classes reflect USNVC Groups for ruderal vegetation plus cultural land cover derived from other sources as described above.

The LANDFIRE Program implemented a map assessment approach that utilized existing information because no resources were available to collect additional data. The goal was to provide assessment results in tandem with data delivery. The assessment sampling strategy for the LANDFIRE 2016 Remap randomly withdrew 10% of the available plots for each of the terrestrial ecological systems classification developed by NatureServe (Comer et al. 2003) within a vegetation production unit, with two caveats. First, no more than 300 plots were withdrawn for any ecological system type. Second, no plots were withdrawn for an ecological system type within an individual production unit if less than 30 plots would be available for the assessment in an effort to maximize training data for modeling categories with few plots. This same set of withdrawn plots was used for the USNVC assessment, although the plot labels were assigned using independent methods and therefore did not always represent conceptually equivalent types. Withdrawn plots were never used in the modeling process so they represent an independent assessment sample.

Confusion tables were created for each of the six LANDFIRE GeoArea delivery packages across CONUS by cross-tabulating the autokey USNVC Group assignment for each assessment plot against the LANDFIRE USNVC Group assignment for map pixels at the plot location. Category agreement focused tables were then generated from each GeoArea contingency table. No stratification, spatial buffering, or category weighting was used (Table 3). A summary report was also generated for each GeoArea to provide an initial indication of map performance (available from https://www.landfire.gov/remapevt_assessment.php).

The assessment sample was based on plots previously available to the program so the sample size and distribution reflected the overall plot numbers and categorical distribution present in any GeoArea. Across the continent the assessment sample was not sufficient for most of the mapped categories. While the agreement results were not high, there was variation in the results across GeoAreas and across categories within each GeoArea. No consistent error patterns were identified, although there is some indication that forest types tend to have lower error rates than shrub and herbaceous types. The opportunities for comparing category error rates across GeoAreas are limited by the sample sizes. For example, Southern Rocky Mountain Ponderosa Pine was mapped in both the Southwest and Northwest GeoAreas but the limited extent in the Northwest GeoArea resulted in too small an assessment set for this GeoArea. Results were not specifically linked to the number of categories assessed. For example, while the Southeast GeoArea had the lowest number of assessed categories and the lowest percentage of assessed categories with more than 70% agreement, the North Central had the lowest percentage of assessed categories with more than 50% agreement. Readers should note that these are absolute errors. If the plot assignment did not match the map assignment exactly it was designated as an error, so errors between floristically similar groups are counted the same as errors between floristically dissimilar groups. Users can review results for USNVC Groups of specific interest to fully understand the results of the assessment analysis.

To understand the results and ramifications for the USNVC Group-based map, a small portion of the Category Agreement Table for the Northwest GeoArea is presented in Table 4. One row represents higher agreement, one row represents moderate agreement and one row represents low agreement. The key information presented in this table is what type of errors were made, not just overall agreement.

For example, the most prevalent misclassification for Intermountain Basins Dry Tall Sagebrush Shrubland & Steppe was Intermountain Low & Black Sagebrush Shrubland & Steppe, followed by Mesic Tall Sagebrush Shrubland & Steppe. These types can occur immediately adjacent to each other across the western landscapes where they are found, and share substantial floristic composition, while the Dry-Mesic Spruce - Fir Forest & Woodland is much less similar so those errors may be more substantial depending on the application. This type of variation in agreement results was common across the GeoAreas, so users can review results for USNVC Groups of specific interest to fully understand the results of the assessment analysis.

Map users should also note that, in addition to the issues with sample size and distribution, these results do not indicate the scope of misclassifications, e.g., how much area within a GeoArea had agreement greater than or less than 50% or 70%.

Table 5 provides a summary of map inputs, including existing map sources for potential distribution modeling. Here we emphasize project components outside the USA. Existing national and regional maps, along with georeferenced field sample data for vegetation types, were all reconciled thematically to the IVC and NatureServe ecological systems classifications (Comer et al. 2003, Josse et al. 2003). Again, given the intent to map potential distribution of “natural” vegetation, only these types were sampled from existing sources. Given limitations of available field samples, randomized samples were also gathered from existing local maps in order to provide a robust and spatially balanced representation of each target map class where there was an acceptable level of confidence in the map source. Here we define “acceptable” as being judged sufficiently reliable by project ecologists experienced in the region and familiar with each map source.

Next, we screened the patch sizes of a given type, patches > 10 km² in area provided the pool of source areas for sample selection. Selection of the 10 km² is again an expert judgment, having evaluated existing maps and concluded that sampling from types depicted in smaller areas risked introducing substantial error. We acknowledge that this risks exclusion of naturally rare ecosystem types, but we judged this risk was warranted given the quality of existing map information for this purpose. This pool of map polygons encompasses 95% of natural landscapes. Stratified random sample selection was weighted by continent-wide area of each type using the log₁₀ (area)*100, providing a sample total weighted towards types of lesser area. A total of 595,951 georeferenced samples were generated for the Americas, with an additional 70,380 held aside for map validation.

Explanatory variables, represented as map surfaces, included a series of biophysical factors, such as bioclimate, landform, slope, and aspect, as well as surface flow accumulation (Table 6).

Bioclimates, as modeled by Metzger et al. (2013), reflect the categorization of temperature and precipitation regime to globally-available remotely sensed data, resulting in a total of 125 unique bioclimates at 1 km² spatial resolution. Geophysical map surfaces were developed using 90 m × 90 m digital elevation data (Jarvis et al. 2008). Slope and aspect were measured in terms of degrees. The methodology for the landform class derivation used a variable moving window to assess relative relief and followed other regional scale approaches to model macro-landforms (Dikau et al. 1991, True et al. 2000). Landforms as discrete units were derived from Weiss (2001), who used the 90 m continental digital elevation to assign pixels into one of the following regional physiographic types: “canyons, deeply incised streams,” “midslope drainages, shallow valleys,” “upland drainages, headwaters,” “u-shaped valleys plains,” “open slopes,” “upper slopes, mesas,” “local ridges, hills in valleys,” “midslope ridges, small hills in plains,” and “mountain tops, high ridges.” Land surface flow accumulation was derived from existing continental data (Lehner et al. 2006), based on a 90 m × 90 m hydrologically conditioned digital elevation. This data set specifically aims to use topographic surfaces to indicate stream flow direction and flow accumulation for use in stream, wetland, and riparian ecosystem analysis. EarthSat NaturalVue (1998–2002) multi-year composites of a red-green-blue translation of 6 bands from Landsat 5–7 at 150 m pixel resolution served as the only spectral inputs to the model. These data served to differentiate the locations of natural types tending to occur in proximity, and therefore, similar to geophysical settings. An overview of the hierarchical modeling process to map potential distibrutions of IVC Macrogroups is shown in Figure 4.

We used a sequential mapping process where maps derived for multiple broader levels of the IVC classification hierarchy were then used as input to modeling distributions of types defined at lower hierarchical levels. In this application, the first thematic level for inductive modeling was the IVC Division (Level 4 from Table 1). For example, in South America the 80 types may be viewed as continental expressions of vegetation formations; with vegetation responding most directly to global climate patterns. On average, 1,032 (Min = 53, Max = 5,724) samples per map class were used to generate the South America portion of this map. No satellite imagery was used in the development of the IVC division-level map output, the EarthSat NaturalVue (1998–2002) imagery was used for both MacroGroup and Ecological Systems level of classfication. Macrogroups were subsequently modeled using an average of 1,234 (Min = 33, Max = 3,433) samples per map class. Both the IVC division map output and EarthSat NaturalVue (1998–2002) imagery were used as map inputs for the macrogroup model. Terrestrial ecological systems, being most numerous and most finely differentiated among the classification units used in this effort, were modeled using the macrogroup map as an additional model input. Once completed, the modeled terrestrial ecological systems layer is the finest thematic scale achievable using this technique. Because these units could be conceptually nested into IVC macrogroup concepts, the “bottom-up” aggregation of maps depicting these units is expected to provide the most reliable map of macrogroups. This aggregation was then reviewed and edited to finalize the distribution of each IVC macrogroup (Figure 5). Numbers of map classes by region and classification level are listed in Table 7.

Over-prediction of more common (over rare, or low sample size) land cover types is a common source of error in CART-based inductive modeling of land cover (Weiss 1995, Lowry et al. 2007). That is, more common land cover classes can be over-mapped at the expense of less common classes. This could be anticipated in this particular application where there is high similarity in predictor variable combinations for vegetation types that are naturally adjacent. In these instances, predicted distributions may be skewed in favor of some over other types in portions of their range of co-occurrence. Because of this, some form of expert-based review and map refinement is unavoidable. We assumed that over-prediction would be concentrated in regional landscapes with extensive land use history and only fragmentary remnants of natural vegetation types. However, because the generalized distribution of each terrestrial ecological system type had been previously documented by country and World Wildlife Fund ecoregion (Josse et al. 2003), this knowledge was used in expert type-by-type review and refinement. Draft model outputs were attributed as extent measures per WWF ecoregion. These distributions were compared against known ecoregion distributions to identify likely error. Types found to be in error had their pixel distributions recoded to most-likely correct types for each WWF ecoregion. In turn, a second phase expert review and refinement followed the procedure used in ecoregions but was applied to each type using a common grid of 100 km² hexagons (Sahr 2013). Again, with each type attributed to the hexagon grid, type-by-type review led to final recoding of pixels to most-likely correct types. Final map products were produced at 90 m and 270 m pixel resolutions with resampling to the unified pixel size accomplished using bilinear interpolation technique suitable for continuous data.

As noted above, during initial sample data collection from map sources, georeferenced samples of each vegetation type were gathered and set aside for use in map validation. These samples were gathered for types that had existing polygons in regional/local source maps > 10 km² in size in South America and > 5 hectares for temperate and tropical North America. Of the 315 Macrogroup map classes in North and South America, 284 had sufficient samples to be quantitatively assessed.

Once map edits were finalized for the 90 m products, validation samples were used to score the degree of agreement between samples and map classes for each map class at three spatial scales. Circular buffers around each sample encompassed 1-km² (within 6 pixels of center) and 5-km² (within 28 pixels of center). A point sample was defined from the centroid of each pixel of the 6 × 6 neighborhood of the 90 m product and is equivalent to the 270 m version of each map. Overlay of these samples on the final map product generated tabular summaries to determine whether or not the mapped class present matched the type labeled to each sample; i.e., the same types co-occur within the buffered area. While truly independent samples could not be acquired to evaluate a spatial model depicting “potential/historical” extent of these vegetation types, this technique provides one initial measure of map quality, and serves as a primary input to decisions regarding use of the map for type-by-type assessment. Thus, the percentage of agreement between validation samples and maps can indicate the degree of map reliability for use with a practical minimum map unit of 270 m vs. 1 km² vs. 5 km².

Table 8 provides a summary of validation statistics for potential distribution maps of macrogroups. Additional detail is found in Suppl. material 1. For the 315 mapped macrogroup types, per class numbers of samples for the 1-km² validation sample area varied from a high of 972 to a low of 10. Summary of validation statistics indicate high (> 80%) to moderate (> 60%) map accuracy overall, and on a per map class basis at 270 m vs. 1 km² vs. 5 km² map resolutions. Using the most demanding “point” (or 270 m) validation sample area, 8 types scored at 90–100% agreement, 11 types scored 80–90% agreement, 16 types scored 70–80% agreement, 39 types scored 60–70% agreement and 52 types scored 50–60% agreement. A total of 158 types (56% of all assessed map classes) scored below 50% agreement.

The inclusion of the 1 km² was limited to the North American portion of the map product for two reasons. First, the sample sizes available for CONUS in North America was substantially higher than that for adjacent countries. Secondly, the inclusion of the 1 km² allowed the examination of the gradient of model performance over a spatial gradient of neighborhoods. Using the 1-km² validation sample area in North America only, 44 types scored at 90–100% agreement, 32 types scored 80–90% agreement, 14 types scored 70–80% agreement, 8 types scored 60–70% agreement, and 9 types scored 50–60% agreement. A total of 12 types (10% of all assessed map classes) scored < 50% agreement. For 1 km² samples, the total sample agreement was 85% and the median level of map class agreement for the types assessed was 88%.

Using the 5 km² validation sample area, 160 types scored at 90–100% agreement, 50 types scored 80–90% agreement, 23 types scored 70–80% agreement, 20 types scored 60–70% agreement and 15 types scored 50–60% agreement. A total of 16 types (6% of all assessed map classes) scored below 50% agreement. For 5-km² samples, the total sample agreement was 85% and the median level of map class agreement was 92%.

These results indicate that map reliability is limited on a per pixel basis (at 270 m pixels), but within relatively small clusters of adjacent pixels, the reliability of the map increases for most map classes.