21st Century Planning Techniques for Creating Fire-Resilient Forests in the American West

: Data-driven decision making is the key to providing effective and efﬁcient wildﬁre protection and sustainable use of natural resources. Due to the complexity of natural systems, management decision(s) require clear justiﬁcation based on substantial amounts of information that are both accurate and precise at various spatial scales. To build information and incorporate it into decision making, new analytical frameworks are required that incorporate innovative computational, spatial, statistical, and machine-learning concepts with ﬁeld data and expert knowledge in a manner that is easily digestible by natural resource managers and practitioners. We prototyped such an approach using function modeling and batch processing to describe wildﬁre risk and the condition and costs associated with implementing multiple prescriptions for risk mitigation in the Blue Mountains of Oregon, USA. Three key aspects of our approach included: (1) spatially quantifying existing fuel conditions using ﬁeld plots and Sentinel 2 remotely sensed imagery; (2) spatially deﬁning the desired future conditions with regards to fuel objectives; and (3) developing a cost/revenue assessment (CRA). Each of these components resulted in spatially explicit surfaces describing fuels, treatments, wildﬁre risk, costs of implementation, projected revenues associated with the removal of tree volume and biomass, and associated estimates of model error. From those spatially explicit surfaces, practitioners gain unique insights into tradeoffs among various described prescriptions and can further weigh those tradeoffs against ﬁnancial and logistical constraints. These types of datasets, procedures, and comparisons provide managers with the information needed to identify, optimize, and justify prescriptions across the landscape.


Introduction
Wildfires are increasing in frequency [1], cost [2,3], and impact [4]. While this increase can be attributed to many factors such as a warming climate [5], recreation [6], increased number of severe events [7], and various anthropogenic-related activities [8], the focus, with regards to the impact of fire, tends to revolve around anthropogenic impacts. Most recently, we have witnessed those consequences in terms of life, property, and social unrest [9], which in turn has helped to focus the wildfire discussions around efficient and effective wildfire protection [10].
Novel approaches to framing wildfire protection (e.g., potential operational delineations) [11] and advancements in physically based modeling tools [12] have been critical to improving our understanding of fire and its potential impacts. Additionally, the successful use of those tools requires data that are accurate at fine resolution, spatially explicit, and current. However, such data often do not exist or are extremely time intensive and costly to develop [13,14]. Moreover, much of the information generated from various existing fire models can be difficult to translate to other metrics related to managing natural resources. In large part, these two limitations, fine resolution, accurate data, and the ability to translate Simulation Software (FSim) exist, these resources alone do not provide the wall-to-wall information needed, at the appropriate scales, to address the finely grained spatially explicit questions or the justification for related prescriptions designed to restore the landscape to a fire-resilient condition.
However, coupled with finely grained remotely sensed information, such as Sentinel 2 imagery [34] and the National Elevation Dataset (NED) [35], data such as multiparty monitoring plots and TIGER road networks and statistical and machine-learning relationships can be leveraged to produce finely grained surfaces depicting the existing forest condition [23] and the costs to move biomass from the forest to a given sawmill [24]. In this study, we demonstrate how existing forest monitoring plots [32], TIGER line files [33], Sentinel 2 imagery [34], and PODs derived from fire managers and common fire modeling tools [11] were used to quantify various forest metrics and treatment costs which in turn were used to inform, justify, and improve restoration decisions. Due to the variety of data sources used, the methodological steps outlined, and the specialized language of the fire, forest, geospatial, and statistical and machine-learning modeling communities, many acronyms are presented in this study. To aid in reading, we provide a list of all acronyms used in this study in Appendix A. Topologically Integrated Geographic Encoding and Referencing (TIGER) road inventories [33], and derived products such as those produced by the Landfire project and Wildfire Risk Simulation Software (FSim) exist, these resources alone do not provide the wall-to-wall information needed, at the appropriate scales, to address the finely grained spatially explicit questions or the justification for related prescriptions designed to restore the landscape to a fire-resilient condition. However, coupled with finely grained remotely sensed information, such as Sentinel 2 imagery [34] and the National Elevation Dataset (NED) [35], data such as multiparty monitoring plots and TIGER road networks and statistical and machine-learning relationships can be leveraged to produce finely grained surfaces depicting the existing forest condition [23] and the costs to move biomass from the forest to a given sawmill [24]. In this study, we demonstrate how existing forest monitoring plots [32], TIGER line files [33], Sentinel 2 imagery [34], and PODs derived from fire managers and common fire modeling tools [11] were used to quantify various forest metrics and treatment costs which in turn were used to inform, justify, and improve restoration decisions. Due to the variety of data sources used, the methodological steps outlined, and the specialized language of the fire, forest, geospatial, and statistical and machine-learning modeling communities, many acronyms are presented in this study. To aid in reading, we provide a list of all acronyms used in this study in Appendix A.

Study Site & Primary Data
Our study site was located in eastern Oregon, USA and consists of both private and public lands, covering approximately 1.2 million ha ( Figure 2). The study area encompasses 1 of 23 high-priority landscapes managed by the Forest Service that receives augmented funding from Congress to accelerate the pace and scale of restoration treatments under the Collaborative Forest Landscape Restoration Program (CFLRP) [36]. Vegetation types across this landscape include dry grasslands, shrubs, and pine and mixed conifer forest types [37]. This research utilized data from a long-term multiparty monitoring program, the Forest Vegetation and Fuels (FVF) program, led by Forest Ser-

Study Site & Primary Data
Our study site was located in eastern Oregon, USA and consists of both private and public lands, covering approximately 1.2 million ha ( Figure 2). The study area encompasses 1 of 23 high-priority landscapes managed by the Forest Service that receives augmented funding from Congress to accelerate the pace and scale of restoration treatments under the Collaborative Forest Landscape Restoration Program (CFLRP) [36]. Vegetation types across this landscape include dry grasslands, shrubs, and pine and mixed conifer forest types [37]. This research utilized data from a long-term multiparty monitoring program, the Forest Vegetation and Fuels (FVF) program, led by Forest Service managers, Oregon State University researchers, and the Blue Mountains Forest Partners, a stakeholder group based in John Day, Oregon that convenes diverse stakeholders to help plan and implement restoration treatments on Forest Service lands [32]. The FVF program provides data that informs accelerated restoration in an adaptive management framework [38].
Historically, chronic surface fire across the study area removed fine surface fuels and promoted low-severity fire effects [39]. Fire exclusion policies have resulted in extremely high wildfire fuel loads [40] that, along with a warming climate, drive large fast-moving fires that threaten human communities and wildlife habitat [41][42][43]. Reducing forest density to facilitate reintroduction of low-severity fire is a major goal of accelerated restoration funded by the CFLRP [44].  [32].
The FVF program provides data that informs accelerated restoration in an adaptive management framework [38]. Historically, chronic surface fire across the study area removed fine surface fuels and promoted low-severity fire effects [39]. Fire exclusion policies have resulted in extremely high wildfire fuel loads [40] that, along with a warming climate, drive large fast-moving fires that threaten human communities and wildlife habitat [41][42][43]. Reducing forest density to facilitate reintroduction of low-severity fire is a major goal of accelerated restoration funded by the CFLRP [44]. Primary raster datasets used in the case study included NED 30 m digital elevation model [35] and Sentinel 2 satellite level 2a processed imagery [45] (Table 1). Sentinel 2 satellite images were acquired for a spring and fall season during 2019 and consist of a total of eight tiles for spring and fall seasons ( Figure 2, Table 1). Acquisition dates for spring and fall tiles fell on the same day within each season, had minimal cloud cover (<1%), and visually appeared to be at a common relative radiometric scale. As such, level 2a bottom of atmosphere products were used without further relative or absolute normalization. To address various grain sizes of the imagery and elevation datasets, all raster surface metrics used as predictor variables for modeling were resampled to a 30 m spatial resolution based on the environmental research institute's (ESRI's) nearest neighbor algorithm [46]. Primary vector datasets used included 2019 TIGER/Line files Primary raster datasets used in the case study included NED 30 m digital elevation model [35] and Sentinel 2 satellite level 2a processed imagery [45] (Table 1). Sentinel 2 satellite images were acquired for a spring and fall season during 2019 and consist of a total of eight tiles for spring and fall seasons ( Figure 2, Table 1). Acquisition dates for spring and fall tiles fell on the same day within each season, had minimal cloud cover (<1%), and visually appeared to be at a common relative radiometric scale. As such, level 2a bottom of atmosphere products were used without further relative or absolute normalization. To address various grain sizes of the imagery and elevation datasets, all raster surface metrics used as predictor variables for modeling were resampled to a 30 m spatial resolution based on the environmental research institute's (ESRI's) nearest neighbor algorithm [46]. Primary vector datasets used included 2019 TIGER/Line files acquired from the US census bureau [33], Flowline and Waterbody line and polygon features from the National Hydrography Dataset (NHD) [47], FVF monitoring plots [32], POD polygons [11], and the location of the Malheur Lumber Company (Table 1).

Existing Condition
To quantify the existing condition of the landscape within a classical forest management perspective, at a fine spatial resolution, we employed an enhanced general additive modeling approach (EGAM) [23] that related summarized plot data to texture metrics derived from multitemporal Sentinel 2 level 2a processed satellite imagery [45] and topographic metrics derived from the NED digital elevation model (DEM) dataset [35]. The EGAM procedure extends generalized additive modeling (GAM) [48] by incorporating a Monte Carlo sampling scheme to create and independently test multiple GAM models as an ensemble. The ensemble of GAMs was then aggregated to produce estimates of response mean and standard error. As part of the EGAM procedure, predictor variables are iteratively selected from all potential predictor variables using all sample observations [49]. After selecting predictor variables, datasets are internally partitioned into training and test subsets to both build and independently test each model, respectively. For each response metric, a total of 50 GAMs were created using random subsets of 80% of the data to build the model and 20% of the data to independently test the model. To guard against extremely poor predictive models becoming part of the EGAM, we assigned a k-factor of 10 for regression-based models and an improvement factor of 0.02 for multinomial models. K and improvement factors were conservatively assigned to identify extremely poor performing models due to extreme samples. Once built, EGAMs were applied back to predictor surface to create continuous estimates of response variable mean and standard error. An implementation of the EGAM procedure can be found in [50].
Key response metrics derived from the monitoring plot data and used to describe the landscape included: cover type (Pine, Fir, or Nonforest), tree basal area measured in m 2 ha −1 (BAH), the number of trees ha −1 (TPH), bone dry tonnes (i.e., a metric ton or megagram) of stem wood biomass ha −1 (SWBH), and bone dry tonnes of total above ground biomass ha −1 (AGBH). BAH, TPH, SWBH, and AGBH were estimated for trees greater than 12.7 cm in diameter at breast height based on commonly used equations ( Table 2). Field plot data collection protocol is described in [32], and expansion factors (EF) for plot estimates were used to bring BAH, TPH, SWBH, and AGBH to the scale of a ha. BAH = basal area ha −1 (m 2 ), TPH = tree counts ha −1 , SWBH = stem wood biomass ha −1 , AGBH = above ground biomass ha −1 , EF = expansion factor for a given fixed plot radius, and b1 and b2 are species-specific coefficients in [52].
Likewise, predictor metrics were created at approximately the same spatial resolution as the extent of each monitoring plot and included mean and standard deviation derived from a 3 by 3 cell neighborhood window passed through each cell, within each image band of the Sentinel 2 imagery used in the analyses (Table 3). Sentinel 2 image bands acquired at a spatial resolution of 20 m (i.e., 5-7, 8a, 11, and 12) were resampled to a spatial resolution of 10 m, using ESRI's nearest neighbor algorithm, prior to performing neighborhood analyses. After performing neighborhood analysis, all raster surfaces were resampled to a 30 m spatial resolution to match the spatial resolution of topographic features and approximate the spatial resolution of the monitoring plots. Topographic features created for the EGAM included slope, northing, and easting and were created from the NED DEM and combined with Sentinel 2 images as potential predictor variables (Table 3). Table 3. List of raster metrics and calculation used to create predictor variable values collected at the location of each plot [25].

Mean
Band cell mean value within a 3 by 3 moving window Standard Deviation Band cell standard deviation within a 3 by 3 moving window Slope Percent slope calculated from a digital elevation model Northing Northing calculation from azimuth Easting Easting calculation from azimuth FVF monitoring plots represent a collection of various systematically located field plots in randomly selected treatments units and nearby untreated controls designed to track changes in vegetation given different implemented treatments [32]. As such, when aggregated together they represent a sample biased to the treatments monitored across the landscape. To address this sampling bias and attempt to center our estimates to the population of pixels that make up our study site while also spreading our dataset across predictor variable space, we implemented a Generalized Random Tessellation Stratified (GRTS) [53] strategy to develop a comparison sample (CS). Our GRTS strategy uses the first two component scores of a correlation-based principal component analysis (PCA) [25,[54][55][56] resulting from predictor variable values found at 100,000 randomly selected locations across the study area to create a sample of observations spread and balanced across predictor variable space [57][58][59]. Using the CS and the predictor variable values found at those locations, we then performed an unsupervised K-means clustering [25,54,60] of the observations, splitting the sample into 100 different K-mean classes. Next, we applied the same K-means algorithm [25,54] to our monitoring sample and compared the proportion of K-mean class counts found within the monitoring plot sample to the CS. K-mean classes that were over-or under-represented in the monitoring plots were then identified and labeled. Observations falling within K-mean classes labeled as over-represented in the monitoring plot sample were then randomly selected until class count proportions matched those of the CS.
K-mean classes identified as being under-represented were then supplemented with additional monitoring plot locations randomly chosen across the study area, which met K-mean class criterion, until class count proportions matched. Newly assigned plots were spatially identified and scheduled to be collected during future field data collection campaigns.

Desired Future Condition (DFC) & Biomass Removal
The ideal or desired condition of the landscape can vary based on management objectives. For this study we narrowed the management objectives to bring Pine dominant forest cover types (DomType, Table 2) to BAH densities that generally resulted in reduced flame lengths (fire-resilient condition) and spread across POD containment zones. POD containment zones were previously delineated for the study area [11] and were prioritized into strategic response zones with the following labels: Protect, Restore, Maintain, Exclude, and High complexity (Appendix B). PODs serve as a natural mechanism to spatially compartmentalize the landscape and build prescriptions designed to generally reduce tree BAH within POD boundaries and further reduce tree densities at and around those boundaries. Using the results of [61] that describe the impacts to flame length based on basal area variable density thinning prescription across similar landscapes, we generated a series of rules that spatially identify desired future BAH, for each 30 m cell within the study area (Table 4).
Primarily Pine dominant forest cover types (DomType, Table 2) were selected for potential treatment. However, next to POD boundaries, cells classified as Fir dominant cover types were also treated. Generally, variable density BAH thresholds for north facing slopes were 27 m 2 ha −1 , while south facing slope thresholds were 20 m 2 ha −1 . For areas within 0.4 km of a POD boundary, BAH thresholds were 15 m 2 ha −1 and prescriptions were designed to harden POD boundaries and reduce the ability of fire to move across those boundaries. All perennial stream side management zones (SMZs, 30 m from a stream) were excluded from mechanical treatment (deferred). Likewise, areas with a slope greater than 50% were deferred due to operational constraints. To be relatively sure that modeled estimates of BAH were accurate, we excluded all cells that had BAH mean estimate that were not statistically different than zero at the 95% confidence level.
Using these rules, Batch Processing [25,62], elevation, slope, the results from the EGAM surfaces, and raster map algebra [63], we quantified how much biomass would need to be removed from each 30 m cell to meet the desired future condition (DFC, Figure 3). Removal values were calculated on a cell basis by subtracting the existing BAH from the DFC BAH (∆BAH). The proportion of DFC BAH removed from the existing condition was then used to estimate the amount of SWBH (∆SWBH) and AGBH (∆AGBH) removed from the landscape as follows: ∆SWBH and ∆AGBH were converted into sawlog (ST) and biomass (BT) tonne products removed from the landscape by incorporating leakage and cell area as follows: In Equations (3) and (4), % leakage represents the amount of woody biomass left on the landscape due to felling and processing and the constant 0.09 was used to convert ha −1 estimates to cell area estimates. For ST, leakage was held constant at 5%. For BT, leakage was held constant at 50%. Cells with negative ∆SWBH and ∆AGBH that did not meet topographic or SMZ thresholds, or that had BAH estimates that were not statistically different than zero, were assigned ∆SWBH and ∆AGBH values of 0. In Equations (3) and (4), % leakage represents the amount of woody biomass left on the landscape due to felling and processing and the constant 0.09 was used to convert ha −1 estimates to cell area estimates. For ST, leakage was held constant at 5%. For BT, leakage was held constant at 50%. Cells with negative ΔSWBH and ΔAGBH that did not meet topographic or SMZ thresholds, or that had BAH estimates that were not statistically different than zero, were assigned ΔSWBH and ΔAGBH values of 0. Figure 3. The process flow used to spatially identify the desired future condition (DFC) expressed as basal area ha −1 (BAH). A digital elevation model (DEM), potential operational delineations (PODs), stream side management zones, and the dominant forest classification were used as inputs within a batch processing spatial modeling framework [25,62] to identify cells meeting Table 4 criterion. Raster surface cells meeting specified criterion were populated with the corresponding desired BAH. Orange collaborative forest landscape restoration program polygons are displayed for reference in the DFC graphic. To address estimation variation in the decision-making process for existing condition continuous EGAMs (Section 3.1), raster surface cell standard error values were used to create a binary mask of estimates that were statistically different than zero. Using cell mean estimates, standard error estimates, and a z-score for two standard deviations from the mean (1.96), we created a binary mask identifying raster cells in which the 95% confidence interval of the mean included an estimate of zero. Cell values with estimated mean values including zero were remapped to a value of zero while cells with mean estimates statistically greater than zero (95% confidence interval) were given a value of one. Continuous response variable masks were then multiplied by EGAM estimates in Sections 3.2 and 3.3 to zero out mean estimates and increase the confidence that the existing condition did not include a value of zero or below. A digital elevation model (DEM), potential operational delineations (PODs), stream side management zones, and the dominant forest classification were used as inputs within a batch processing spatial modeling framework [25,62] to identify cells meeting Table 4 criterion. Raster surface cells meeting specified criterion were populated with the corresponding desired BAH. Orange collaborative forest landscape restoration program polygons are displayed for reference in the DFC graphic. To address estimation variation in the decision-making process for existing condition continuous EGAMs (Section 3.1), raster surface cell standard error values were used to create a binary mask of estimates that were statistically different than zero. Using cell mean estimates, standard error estimates, and a z-score for two standard deviations from the mean (1.96), we created a binary mask identifying raster cells in which the 95% confidence interval of the mean included an estimate of zero. Cell values with estimated mean values including zero were remapped to a value of zero while cells with mean estimates statistically greater than zero (95% confidence interval) were given a value of one. Continuous response variable masks were then multiplied by EGAM estimates in Sections 3.2 and 3.3 to zero out mean estimates and increase the confidence that the existing condition did not include a value of zero or below.

Cost Revenue Assessment (CRA)
Delivered cost was calculated using the model described in [24]. Cost in this context referred to all costs associated with felling, processing, and moving woody biomass from the forest to the Malheur Lumber Company conversion facility ( Figure 2). The delivered cost model in [24] uses a least-cost path algorithm [64] to transform surface distances, both on-and off-road, into travel times based on specified harvesting systems and associated machine and operation rates (Table 5). For our study, rubber-tire-skidder-based systems were only considered for the implementation of mechanically based prescriptions. Of note, the delivered cost model can include barriers to motion for off-road travel. We used NHD perennial streams and waterbodies and TIGER\line U.S. interstate and state highways as barriers to motion for off-road travel.
Travel time was then converted to dollars tonne −1 (bone dry) using cost hr −1 operation rates (Table 5). Outputs from the delivered cost model included medium-to fine-grained raster surfaces that estimate total round-trip travel costs of off-road skidding and onroad transportation for each cell across the landscape. Inputs for this model included a TIGER\linear road network with rates of speed defined for each road segment (Table 6), the NED DEM, and cost components depicting the rates of speed For moving biomass from the site of harvest to the roadside, hourly rates for machines use an average payload (tonnes trip −1 ) for different types of equipment, fixed operations, administration cost (as a dollar tonne −1 constant), and a vector layer of off-road barriers. Table 6. Queries used to allocate the average speed to each road segment and its derived raster cells, considering the Master Address File (MAF) TIGER\Line Feature Class Code (MTFCC) 1 .

Query
Speed km h −1 MTFCC = "S1400": Local Road, Rural Road, City Street 88 MTFCC = "S1200": Secondary road 56 MTFCC = "S1100": Primary road 40 NOT (MTFCC = "S1400" OR MTFCC = "S1200" OR MTFCC = "S1100") 40 Operation costs of USD 65 tonne −1 included felling (USD 15) and processing (USD 50), but not skidding, which is included in off-road travel cost accounting. The administrative costs were assumed to be zero to provide a normalized comparison across all forest ownerships. Together, on-road, off-road, and operation costs constituted the optimal potential total cost (OPTC) of moving biomass, measured in dollars tonne −1 , to Malheur Lumber Company in a spatially explicit manner, at the spatial resolution of 30 m 2 across the study site. The estimated actual cost of prescription implementation (ACPI) was then calculated by multiplying OPTC by the amount of biomass removed at each cell, where 0.09 converts ∆AGBH estimates to cell area estimates as follows: Like ACPI, estimated revenue is dependent on ∆AGBH. However, estimated revenue (ER) also depends on gate prices for various products removed from forests as follows: ER biomass = BT * Gate Price biomass (7) We fixed gate prices for sawlogs at USD 55 green tonne −1 and USD 28 green tonne −1 for biomass products [65]. Assuming an average green moisture content of 94% [66], dry tonne gate price equivalents were estimated at USD 106.70 tonne −1 for sawlogs and USD 54.32 tonne −1 for biomass products. Using map algebra, function modeling, and Equations (6) and (7), we calculated saw and biomass ER for each 30 m cell and created spatially explicit raster surfaces. Adding saw and biomass ERs together, we calculated total ER. Finally, we subtracted ACPI from total ER for overlapping raster cells to create a raster surface of estimated anticipated income (EAI) of implementing a prescription as follows:

Quantifying Existing Condition
Our GRTS selection strategy selected 429 plots balanced and spread across the first two principal components of transformed predictor variable space. The first two principal components of our GRTS selection strategy accounted for 45% of the correlation within data. Of the 429 selected plots, 301 plots had been visited and were used to calibrate EGAMs. In total, 128 plots were identified for future field data collection campaigns. Ideally, EGAMs would be created after collecting the additional 128 plots. However, given the spread and approximate balance of this sample (Figure 4), we chose to use these plots to calibrate our EGAMs. From this sample, we expect global model bias, due to sampling, to be minimal and that interpolated model estimates across predictor variable space will approximate the relationship occurring within unsampled K-mean classes. Of importance with regards to estimation error in under-sampled regions of predictor variable space is the variability across GAMs within our EGAMs. To spatially identify those regions of predictor variable space we projected K-mean class into coordinate space and provided that binary map as a band within each of the raster surfaces.  Using our sample of 301 observations, we created 50 GAMs from random subsets of the data and combined them together as an ensemble for each of our response variables. Statistically selected predictor variables identified in the selection process are presented in Table 7. On average the selection process reduced the potential number of predictor variables from 44 to 8.4 for each EGAM. The mean and standard deviation of estimates from all GAMs within each EGAM were used to populate raster surface band cell mean and standard error values (Figures 5 and 6). EGAM fit statistics are presented in Table 7 and Figure 7 and on average explain 55.99% of the variation in the data for continuous responses and had an average map accuracy of 83.38% for the dominant forest classification. Of particular interest in Figure 7 is the overall strength of the relationship between observed and predicted values and the nonconstant variation in estimated values. Using our sample of 301 observations, we created 50 GAMs from random subsets of the data and combined them together as an ensemble for each of our response variables. Statistically selected predictor variables identified in the selection process are presented in Table 7. On average the selection process reduced the potential number of predictor variables from 44 to 8.4 for each EGAM. The mean and standard deviation of estimates from all GAMs within each EGAM were used to populate raster surface band cell mean and standard error values (Figures 5 and 6). EGAM fit statistics are presented in Table 7 and Figure 7 and on average explain 55.99% of the variation in the data for continuous responses and had an average map accuracy of 83.38% for the dominant forest classification. Of particular interest in Figure 7 is the overall strength of the relationship between observed and predicted values and the nonconstant variation in estimated values.
Using our sample of 301 observations, we created 50 GAMs from random subsets of the data and combined them together as an ensemble for each of our response variables. Statistically selected predictor variables identified in the selection process are presented in Table 7. On average the selection process reduced the potential number of predictor variables from 44 to 8.4 for each EGAM. The mean and standard deviation of estimates from all GAMs within each EGAM were used to populate raster surface band cell mean and standard error values (Figures 5 and 6). EGAM fit statistics are presented in Table 7 and Figure 7 and on average explain 55.99% of the variation in the data for continuous responses and had an average map accuracy of 83.38% for the dominant forest classification. Of particular interest in Figure 7 is the overall strength of the relationship between observed and predicted values and the nonconstant variation in estimated values.    Table 7. Ensemble of generalized additive models (EGAM) predictor variables and fit statistics using all observations. Root mean squared error and overall map accuracy are reported for continuous (a) and categorical (b) response variables.

Desired Future Condition & Biomass Removals
Using the Pine dominant forest type class, the DEM, stream networks, POD boundaries, batch processing, function modeling, and the rules described in Table 4, we were able to create a raster surface that spatially defined the DFC of a landscape in terms of BAH (Figure 3). Subtracting DFC from the existing BAH raster surfaces, we spatially . Accuracy assessment with user, producer, and overall accuracies, chi-squared (X 2 ) statistics, and Cohen's Kappa statistics for the most likely DomType class (left) and general trends between observed and predicted values for basal area ha −1 (BAH), trees ha −1 (TPA), stem woody biomass ha −1 (SWBH), and above ground biomass ha −1 (AGBH). For BAH, TPH, SWBH, and AGBH graphics, a loess trend line (black dashed) and corresponding 99% confidence region (grey region) around that line were overlaid with a linear 1 to 1 trend line (red line) to visually aid in evaluating model fit. Additionally, the mean (x) and standard deviation (σ) of observed values and Pearson's correlation coefficient (r) comparing observed and predicted values are reported in the bottom right corner of BAH, TPH, SWBH, and AGBH graphics.

Desired Future Condition & Biomass Removals
Using the Pine dominant forest type class, the DEM, stream networks, POD boundaries, batch processing, function modeling, and the rules described in Table 4, we were able to create a raster surface that spatially defined the DFC of a landscape in terms of BAH (Figure 3). Subtracting DFC from the existing BAH raster surfaces, we spatially quantified the reduction in BAH needed to meet DFC. Applying the proportion change in BAH to existing SWBH and AGBH surfaces and incorporating leakage (Equations (3) and (4)), we created two raster surfaces depicting the ST and BT removed from the landscape to meet DFC (Figure 8). Approximately 527,000 ha across the study area had negative ΔBAH. Those cells represented areas in which the existing BAH was less than the DFC and as such were remapped to zero tonnes removed. Conversely, there were approximately 236,000 ha identified across the study area as having more BAH than desired. In many locations (approximately 86,000 ha) ΔBAH was relatively small (<4 m 2 ha −1 ), suggesting minimal deviation from the desired future condition that may be of less importance to treat than cells with larger ΔBAH.
To meet DFC, approximately 6.5 and 0.9 million tonnes of sawlog and biomass products, respectively, will need to be removed from the landscape. Interestingly, while POD hardening zones account for only 22.66% of the total study area, the majority of tonnes removed to meet DFC came from those zones (approximately 54%). Fortunately, many of the POD boundaries were drawn along road networks, potentially reducing implementation costs associated with meeting wildfire objectives.

Cost Revenue Assessment
We spatially quantified OPTC and travel time for on-and off-road costs (Figure 9). Using batch processing, function modeling, machine rates, and payloads (Table 5) Approximately 527,000 ha across the study area had negative ∆BAH. Those cells represented areas in which the existing BAH was less than the DFC and as such were remapped to zero tonnes removed. Conversely, there were approximately 236,000 ha identified across the study area as having more BAH than desired. In many locations (approximately 86,000 ha) ∆BAH was relatively small (<4 m 2 ha −1 ), suggesting minimal deviation from the desired future condition that may be of less importance to treat than cells with larger ∆BAH.
To meet DFC, approximately 6.5 and 0.9 million tonnes of sawlog and biomass products, respectively, will need to be removed from the landscape. Interestingly, while POD hardening zones account for only 22.66% of the total study area, the majority of tonnes removed to meet DFC came from those zones (approximately 54%). Fortunately, many of the POD boundaries were drawn along road networks, potentially reducing implementation costs associated with meeting wildfire objectives.

Cost Revenue Assessment
We spatially quantified OPTC and travel time for on-and off-road costs (Figure 9). Using batch processing, function modeling, machine rates, and payloads (Table 5), and the outputs from the delivered cost tool, we transformed travel times into dollars tonne −1 raster surfaces that can be easily manipulated within batch processing and function modeling. . Geospatial processes and outputs from the delivered optimum potential total cost (OPTC) model [24]. Raster cell values depicting travel times and delivered costs spatially quantify the potential time and corresponding cost using machine rates and payloads ( Approximately, 77% of the study area had a delivered OPTC between USD 60 and USD 70 tonne −1 . Only 6% of the landscape had a delivered OPTC less than USD 60 tonne −1 and 17% of the study area had a delivered OPTC more than USD 70 tonne −1 . Moreover, there were no locations across our study area that had an estimated delivered OPTC less than USD 54.32, suggesting biomass products alone will produce negative profits. While the delivered OPTC identified costs on a tonne basis, the ACPI was dependent on the amount of total biomass tonne removed. On a cell basis we used map algebra and Equation (6) to spatially quantify ACPI (Figure 10). To meet DFC across the landscape would cost approximately USD 648 million. However, only 52% of that total cost was attributed to hardening POD boundaries. Moreover, 94% of the total ACPI area occurred within cells having an OPTC less than USD 70 tonne −1 , suggesting that restoring the landscape to DFCs would result in positive profits if a component of the prescription produced sawlog products (i.e., sawlog gate prices will offset delivered costs).
Like ACPI, ER for sawlog and biomass products was dependent on the amount of total biomass tonne removed from the landscape. Unlike ACPI, though, leakage and product type were factored into the estimated amount of ST and BT delivered to Malheur Lumber Company (Equations (3) and (4)). On a cell basis, we accurately estimated ER and ER using Equations (6) and (7) (Figure 10). Profit margins were then calculated and spatially depicted by subtracting ACPI from total ER ER + ER on a cell basis (Figure 10). . Geospatial processes and outputs from the delivered optimum potential total cost (OPTC) model [24]. Raster cell values depicting travel times and delivered costs spatially quantify the potential time and corresponding cost using machine rates and payloads (Table 5) to move woody biomass from a given location on the landscape to the Malheur Lumber Company (yellow point within each graphic). Linear road network (black lines) is displayed for reference in the various graphics.
Approximately, 77% of the study area had a delivered OPTC between USD 60 and USD 70 tonne −1 . Only 6% of the landscape had a delivered OPTC less than USD 60 tonne −1 and 17% of the study area had a delivered OPTC more than USD 70 tonne −1 . Moreover, there were no locations across our study area that had an estimated delivered OPTC less than USD 54.32, suggesting biomass products alone will produce negative profits.
While the delivered OPTC identified costs on a tonne basis, the ACPI was dependent on the amount of total biomass tonne removed. On a cell basis we used map algebra and Equation (6) to spatially quantify ACPI ( Figure 10). To meet DFC across the landscape would cost approximately USD 648 million. However, only 52% of that total cost was attributed to hardening POD boundaries. Moreover, 94% of the total ACPI area occurred within cells having an OPTC less than USD 70 tonne −1 , suggesting that restoring the landscape to DFCs would result in positive profits if a component of the prescription produced sawlog products (i.e., sawlog gate prices will offset delivered costs).
Like ACPI, ER for sawlog and biomass products was dependent on the amount of total biomass tonne removed from the landscape. Unlike ACPI, though, leakage and product type were factored into the estimated amount of ST and BT delivered to Malheur Lumber Company (Equations (3) and (4)). On a cell basis, we accurately estimated ER saw and ER biomass using Equations (6) and (7) (Figure 10). Profit margins were then calculated and spatially depicted by subtracting ACPI from total ER (ER saw + ER biomass ) on a cell basis ( Figure 10). Approximately 99% of the area identified as needing treatment had a positive profit margin, suggesting that most of the landscape that departed from DFC had enough sawlog products to compensate for delivered costs. If all prescriptions were implemented across the landscape, we estimated a net positive profit of USD 121 million with an average positive return of USD 513 ha −1 . However, more than 80% of all Profit cell values had dollar ha −1 estimates less than USD 513, suggesting marginal returns for most operational units. Relatively speaking, hardening POD boundaries returned larger net profits than prescriptions designed to increase fire resilience (Figure 11). Approximately 99% of the area identified as needing treatment had a positive profit margin, suggesting that most of the landscape that departed from DFC had enough sawlog products to compensate for delivered costs. If all prescriptions were implemented across the landscape, we estimated a net positive profit of USD 121 million with an average positive return of USD 513 ha −1 . However, more than 80% of all Profit cell values had dollar ha −1 estimates less than USD 513, suggesting marginal returns for most operational units. Relatively speaking, hardening POD boundaries returned larger net profits than prescriptions designed to increase fire resilience (Figure 11). Figure 11. Supply chain by product type, implementation cost zones (x-axis), and prescription (left graphic) and area treated (y-axis 1) and profit margin (y-axis 2) by implementation cost zones (x-axis) and prescription (right graphic) graphs.

Discussion
Optimal decision making requires detailed and accurate information. We present a methodology that uses function modeling [25], batch processing [62], field plot data, and readily available information to create spatially explicit metrics used within the field of forestry to manage forests and fuels. Our approach builds on the outputs of common fire risk and simulation models by converting field plots and passively acquired remotely sensed data into raster surfaces that contain values of commonly used and well-known forest management metrics, spatially tied to fire risk and spread through PODs. Additionally, using a spatially explicit modeling approach we quantified the cost of delivering wood products to a specified processing location [24]. Moreover, using batch processing and map algebra, we transformed prescriptions designed to increase fire resilience and harden POD permitters into spatially explicit representations of DFC and further used those DFCs to quantify tree biomass removal. Finally, we combined each of those products at the cell level (spatial resolution of 30 m) to spatially identify and quantify various aspects of implementing treatments to meet DFCs.
Our analyses provide unique insights into both scale and application. Of note, the process of linking field plot data with Sentinel 2 multitemporal imagery using EGAMs allowed us to quickly and easily create 30 m raster surfaces depicting various mean metrics and estimates of standard error that can be used to manage forests at various spatial scales. These mean metrics and estimates of mean error provide practitioners with accurate, spatially explicit information pertaining to the forested condition at a spatial resolution that can be easily evaluated based on known forest practices and can be used to inform decision making. In our study, we demonstrated how these estimates (Section 3.1) coupled with DFCs (Section 3.2) can be used to quantify the amount of biomass removed from the landscape to both harden POD boundaries and restore forests to a more fire-resilient condition. Moreover, we demonstrated how estimates of mean error can be used to focus our restoration efforts on locations that have mean BAH estimates different than zero (95% confidence limits).
While our EGAMs accounted for much of the variation in the data (Section 3.1), improvements in model estimation error can be made by collecting additional field plots within under-represented portions of predictor variable space and portions of predictor variable space with large standard error estimates. Moreover, mean and standard error estimates can further be used with additional sampling and design-based estimators (e.g., ratio estimation [67]) to improve mean and total basal area, tree count, and biomass estimates for subregions (stands) within the study area. However, given the spread and Figure 11. Supply chain by product type, implementation cost zones (x-axis), and prescription (left graphic) and area treated (y-axis 1) and profit margin (y-axis 2) by implementation cost zones (x-axis) and prescription (right graphic) graphs.

Discussion
Optimal decision making requires detailed and accurate information. We present a methodology that uses function modeling [25], batch processing [62], field plot data, and readily available information to create spatially explicit metrics used within the field of forestry to manage forests and fuels. Our approach builds on the outputs of common fire risk and simulation models by converting field plots and passively acquired remotely sensed data into raster surfaces that contain values of commonly used and well-known forest management metrics, spatially tied to fire risk and spread through PODs. Additionally, using a spatially explicit modeling approach we quantified the cost of delivering wood products to a specified processing location [24]. Moreover, using batch processing and map algebra, we transformed prescriptions designed to increase fire resilience and harden POD permitters into spatially explicit representations of DFC and further used those DFCs to quantify tree biomass removal. Finally, we combined each of those products at the cell level (spatial resolution of 30 m) to spatially identify and quantify various aspects of implementing treatments to meet DFCs.
Our analyses provide unique insights into both scale and application. Of note, the process of linking field plot data with Sentinel 2 multitemporal imagery using EGAMs allowed us to quickly and easily create 30 m raster surfaces depicting various mean metrics and estimates of standard error that can be used to manage forests at various spatial scales. These mean metrics and estimates of mean error provide practitioners with accurate, spatially explicit information pertaining to the forested condition at a spatial resolution that can be easily evaluated based on known forest practices and can be used to inform decision making. In our study, we demonstrated how these estimates (Section 3.1) coupled with DFCs (Section 3.2) can be used to quantify the amount of biomass removed from the landscape to both harden POD boundaries and restore forests to a more fire-resilient condition. Moreover, we demonstrated how estimates of mean error can be used to focus our restoration efforts on locations that have mean BAH estimates different than zero (95% confidence limits).
While our EGAMs accounted for much of the variation in the data (Section 3.1), improvements in model estimation error can be made by collecting additional field plots within under-represented portions of predictor variable space and portions of predictor variable space with large standard error estimates. Moreover, mean and standard error estimates can further be used with additional sampling and design-based estimators (e.g., ratio estimation [67]) to improve mean and total basal area, tree count, and biomass estimates for subregions (stands) within the study area. However, given the spread and balance of our sample and the fit of our existing condition models, we feel confident that our estimates can be used to help inform and prioritize prescriptions across the study area to reduce wildfire risk and increase resilience to wildfire.
Coupled with known information related to variable density treatments [61], fire severity [68], wildfire risk [29,30], and the conceptual idea behind PODs [11], these types of spatially explicit raster surfaces provide information that can be easily manipulated within a geographic information system (GIS) to prioritize multiple objectives and compare and contrast alternatives [69]. For example, with these types of datasets we can prioritize and summarize biomass removal and retention, cost and revenue of implementation, and profit using POD boundaries. We can further use these datasets to develop and summarize quality risk assessment metrics and a series of spatial filters that take into consideration contiguous areas with enough biomass to justify operations ( Figure 12). Moreover, we can compare coarse wildfire restoration strategies such as Maintain, Protect, and Restore, across strategies and among PODs at fine spatial scales. Additionally, using spatial adjacency rules and estimates of standard errors, we can filter and define operational treatment units that are at least 5 ha in size and that have at least three tonne of biomass on a cell basis (95% confidence interval threshold). Finally, treatment unit geometry can be used to summarize removals, costs, revenues, and profits at the treatment unit level. balance of our sample and the fit of our existing condition models, we feel confident that our estimates can be used to help inform and prioritize prescriptions across the study area to reduce wildfire risk and increase resilience to wildfire. Coupled with known information related to variable density treatments [61], fire severity [68], wildfire risk [29,30], and the conceptual idea behind PODs [11], these types of spatially explicit raster surfaces provide information that can be easily manipulated within a geographic information system (GIS) to prioritize multiple objectives and compare and contrast alternatives [69]. For example, with these types of datasets we can prioritize and summarize biomass removal and retention, cost and revenue of implementation, and profit using POD boundaries. We can further use these datasets to develop and summarize quality risk assessment metrics and a series of spatial filters that take into consideration contiguous areas with enough biomass to justify operations ( Figure 12). Moreover, we can compare coarse wildfire restoration strategies such as Maintain, Protect, and Restore, across strategies and among PODs at fine spatial scales. Additionally, using spatial adjacency rules and estimates of standard errors, we can filter and define operational treatment units that are at least 5 ha in size and that have at least three tonne of biomass on a cell basis (95% confidence interval threshold). Finally, treatment unit geometry can be used to summarize removals, costs, revenues, and profits at the treatment unit level. Figure 12. Example of using potential operational delineations (PODs) to perform a coarse scale prioritization (Maintain, Protect, Restore) combined with fine scale thresholding (Criterion) to identify contiguous regions that have cell removal values statistically greater than three tonne (95% confidence level threshold) and that are at least five ha in size. These regions, within and along a given POD, represent operational treatment units for which tonne, cost, revenue, and profit Figure 12. Example of using potential operational delineations (PODs) to perform a coarse scale prioritization (Maintain, Protect, Restore) combined with fine scale thresholding (Criterion) to identify contiguous regions that have cell removal values statistically greater than three tonne (95% confidence level threshold) and that are at least five ha in size. These regions, within and along a given POD, represent operational treatment units for which tonne, cost, revenue, and profit can be quickly summarized and reported. Restoration PODs A and B in the right graphic identify two units with restoration foci that emphasize within-unit concerns (A) versus controlling spread into neighboring protected regions (B).
Although Section 3.2 provides estimates of woody biomass removed to meet DFC, implementation costs often constrain where and the types of treatments that can occur. Our CRA estimates spatially quantify those costs and revenues for a snapshot in time given described machine rates and gate prices. As such, these estimates are dependent on the systems specified and can vary by operation. However, using batch processing, function modeling, and the delivered cost model, machine rates and gate prices can be adjusted quickly and easily. Additionally, these rates can be spatially allocated as raster surfaces within the delivered cost model to specify unique cost structures given various conditions [69]. Moreover, because CRA surfaces are spatially defined, additional costs not defined in our analyses can easily be incorporated into the CRA analysis on a ha −1 , tonne −1 , and spatially explicit basis. For example, both the cost of prescribed burning and road maintenance were not defined within our CRA. These costs can easily be added to the overall cost associated with implementation in a spatially explicit manner using batch processing and map algebra. Likewise, gate prices can vary from one sawmill to another and can be modified easily within the batch processing modeling framework. Because estimates of removed biomass are spatially identified, gate prices can also be allocated in a spatially explicit manner, providing flexibility in quantifying various prescriptions.
While changes to machine rates and gate prices and additional costs will change the overall cost, revenue, and profit surfaces, the spatial nature of these estimates provides flexibility in determining actual treatment costs and revenues. Moreover, these surfaces can be used to identify areas with limited access and financially justify treatments.
Generally, we recommend using the CRA surfaces to quantify the costs and revenue of prescriptions designed to meet forest objectives. However, cost and revenue surfaces could also be used for treatment optimization. Regardless of how they are used, it is important to recognize that machine rates and gate prices vary based on markets, and as such, the variability of markets, machine rates, and assumed prices should be factored into the planning and decision-making process.
In our CRA, the cost of felling and processing non-sawlog biomass was greater than gate prices. This suggests that to have positive net profits, sawlog biomass must be a component of the treatment. All prescriptions in our study included sawlog biomass products, so net positive profits were realized for approximately 99% of the landscape identified with more BAH than DFC thresholds. However, this does not mean that all cells with a net positive profit are practical to treat. Further thresholding, such as is depicted in Figure 12, may need to be employed to help identify operational treatment units that are practical.
Likewise, delivered cost for our CRA represents an optimal cost. That is, landings to which biomass are skidded are identified based on the least cost path to each individual 30 m road cell. Thus, the true cost associated with skidding material to a given landing will be more for all cells that are not defined as the least cost path. If landing locations are known prior to performing the delivered cost model, they can be incorporated into the analyses.
Finally, our data-driven approach to identifying, describing, and quantifying existing conditions, desired future conditions, and costs and revenues in a finely grained, spatially explicit manner across broad landscapes necessitates the use of spatial modeling, remote sensing, machine learning, and statistics. In many instances, this dependency will require translating GIS analyses into actionable information that practitioners can use to inform decision making [70]. While we recognize this to be a substantial hurdle to adoption, we argue that a larger emphasis must be placed on building those skillsets and integrating them within natural resource management.

Conclusions
We present an accurate, finely grained, spatially explicit approach to quantifying forest characteristics, desired future conditions, and the expected cost and revenues of implementing prescriptions designed to reduce wildfire risk and increase forest resistance to wildfire. Our approach leverages field plots, spatial, statistical, and machine-learning tools, and readily available datasets such as remotely sensed imagery, existing road and water networks, and elevation models to accurately quantify various aspects of restoring the forest landscape to a more fire-resilient condition. These spatially explicit datasets provide the types of information needed to meaningfully inform data-driven decision making.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.

Acknowledgments:
We would like to thank the reviewers for their helpful comments and suggestions and would like to state that the views, thoughts, and opinions expressed within the article belong solely to the authors and do not necessarily represent the views of the United States Agricultural Department or the United States government.

Conflicts of Interest:
The authors declare no conflict of interest.