Next Article in Journal
Tensile Properties of Natural and Synthetic Rattan Strips Used as Furniture Woven Materials
Previous Article in Journal
Effect of Seasonal Rains and Floods on Seedling Recruitment and Compositional Similarity in Two Lowland Tropical Forests
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Estimation of Forest Disturbance from Retrospective Observations in a Broad-Scale Inventory

1
USDA Forest Service, Southern Research Station, Center for Integrated Forest Science, Blacksburg, VA 24060, USA
2
Department of Forest Resources, University of Minnesota, St. Paul, MN 55108, USA
3
USDA Forest Service, Northern Research Station, Forest Inventory and Analysis, York, PA 17402, USA
4
USDA Forest Service, Southern Research Station, Forest Inventory and Analysis, Blacksburg, VA 24060, USA
*
Author to whom correspondence should be addressed.
Submission received: 12 November 2020 / Revised: 30 November 2020 / Accepted: 1 December 2020 / Published: 3 December 2020
(This article belongs to the Section Forest Inventory, Modeling and Remote Sensing)

Abstract

:
Understanding the extent and timing of forest disturbances and their impacts is critical to formulating effective management and policy responses. Broad-scale inventory programs provide key estimates of forest parameters that indicate the extent and severity of disturbance impacts. Here, we review the use of a post-stratified estimator in a panelized design, in the context of disturbance observations that are collected retrospectively. We further develop a sample weight adjustment that is requisite for proper estimation of the extent and timing of disturbances. Using populations from areas of Arkansas, California, and Maine in the US, the weight adjustment technique was tested in a Monte Carlo simulation. We found that the estimated area of disturbance using the weight adjustment technique had satisfactory agreement with the true population values and performed considerably better than the conventional post-stratified estimation approach. The proliferation of panelized forest inventory designs globally suggests that accurate estimates of areal extent and timing of disturbances will often require that weighting adjustment techniques be employed in the estimation process.

1. Introduction

Disturbance is a universal component of forests and quantifying the areal extent and timing of forest disturbances is key to understanding and managing forested lands. Quality data on forest impacts from disturbance events aids in identifying management opportunities, restoration, and more accurately estimating effects on carbon stocks, sequestration rates, and other forest structure parameters. The primary causes of disturbance in the United States are fire, insects and diseases, drought, and harvesting [1], all of which have transformative effects on forest function and total biomass even at low levels of mortality [2]. Abiotic and biotic natural disturbance agents and interactions between them have resulted in amplified levels of forest disturbance [3,4,5], providing a pressing need for information on the extent, timing, and effects of forest disturbance. Likewise, there have been shifts in harvest frequencies over time. While some disturbance events are stand replacing, others result in selective mortality to a species or cohort. This research aims to develop improved techniques for disturbance estimation under rotating panel, repeated measure sampling designs for forest inventory, such as those used in the United States, France, and China [6].
Broad-scale forest disturbance is monitored through two main sources: the field-based component of national forest inventories (NFI) and remotely sensed imagery. Remotely sensed imagery has been used in forestry since the 1950s. More recent remote sensing applications include producing annual disturbance maps based on the Landsat time series stacks [7,8,9]. Since 1982, Landsat has collected spectral reflectance data at a 30 m spatial scale and 16-day temporal density, which is particularly applicable to forestry given the ability to detect changes across relatively short time intervals. With the application of robust models, these data allow for evaluating the timing and geographic extent of major forest disturbances. However, with the exceptions of Moisen and colleagues [10], Schleeweis and colleagues [11], and the Monitoring Trends in Burn Severity (MTBS) project [12,13], there has been little focus on developing models that identify the type of disturbance (e.g., drought, fire, harvest, etc.). Additionally, the impacts on forest structure and biomass are often not quantified in existing remote sensing approaches to disturbance assessment.
Forest inventory data from NFIs provide an alternative to remote sensing approaches or may be combined with remotely sensed data using model-assisted or model-based estimators. Through NFIs, in situ data on disturbance occurrence, disturbance type, and forest structure pre- and post-disturbance are collected, providing a valuable base for disturbance related analytics. However, it is common with longitudinal designs to record mortality-causing disturbances that have occurred since the previous plot measurement. This retrospective observation of disturbance, when coupled with rotating panel repeated measure designs, creates issues for estimating annual disturbance totals. To our knowledge, these issues are not addressed in the classic sampling and estimation literature offered by Cochran [14], Kish [15], and Särndal and colleagues [16], or the NFI specific literature (e.g., [17,18]). However, these estimation issues (described in detail in the following section) must be rectified, not only to construct correct annual estimates based on a direct estimator or model-assisted estimator, but also to construct correct annual estimates when employing small area estimation techniques where remotely sensed information serves as ancillary data (e.g., the Fay–Herriot area model [19]).
The objective of this research is to develop an appropriate approach for estimating the annual disturbance parameters such as areal extent and mortality volume when multiple panels of a field-based inventory in a rotating panel sample design are combined for estimation. To accomplish this objective, we first develop a technique to modify the sample weights of the rotating panel design to account for the retrospective observation and then test the technique in a Monte Carlo framework using three hypothetical populations with different disturbance levels and temporal trajectories. We further demonstrate the technique with a case study.

2. Materials and Methods

For this research we rely on the statistical design implemented by the USDA Forest Service, Forest Inventory and Analysis (FIA) program. Specifics about the sampling design and estimators are provided by Bechtold and Patterson [17]. Details regarding data collection protocols are provided by the USDA [20].

2.1. Plot-Level Disturbance

The FIA program employs a 5 or 7 panel design by state in the eastern United States and a 10-panel design in the western United States. One panel is measured each year, so the remeasurement schedule is approximately equal to the number of panels (P). Up to three observable disturbances (e.g., fire) and the year of each disturbance are recorded. Further, up to three forest management treatments (e.g., timber harvesting) are also recorded, along with the year of each treatment [20]. For simplicity, we refer to both disturbances and treatments simply as disturbances.
Plot-level disturbances are identified in a retrospective fashion, in that a field crew visiting a permanent sample plot determines whether the plot is disturbed relative to the previous visit. Disturbance type and the year of disturbance are the basic data collected, but we must think of the data in a different way for estimation purposes. In fact, under a 5-panel design, there are five different observations collected, one observation for each year since the last measurement. This scenario is expressed in Table 1, which represents a hypothetical example based on a 5-year remeasurement period, with “1” indicating a disturbance, “0” indicating no disturbance, and “--” indicating unobservable (i.e., a year in the future from the field visit or a year preceding the previous visit). In this example, plot 1 was visited in 2008, and the data recorded were binary indicators for whether the plot was disturbed for each year since the last measurement. As such, there are five retrospective disturbance observations recorded; one for each year between field visits. Plots from other measurement years, or panels in the five-panel rotation, have the same pattern. In many cases, no disturbances are observed, as illustrated by plot 4 in Table 1, having zeros recorded for each of the possible disturbance years.

2.2. Post-Stratified Estimator

Bechtold and Patterson [17] (2005) describe the post-stratified estimator used by the FIA program. The technique is post-stratified in that permanent plot locations are used and strata are determined based on ancillary data after the sample was drawn. The ancillary data are typically derived from classified satellite imagery [21] such as the National Land Cover Database (NLCD) land cover map or the NLCD tree canopy cover map [22,23]. The estimate of the population mean for a domain (d) of interest for parameter y is:
y d ¯ = h H W h y h d ¯
where,
y h d ¯ = n h 1 i = 1 n h y h i d .
Wh is the post-stratification weight, nh is the sample size in stratum h = 1…H, and yi is the observation from the ith plot of variable y. The estimated variance of y d ¯ is:
v ( y d ¯ ) = n 1 ( h H W h n h v ( y h d ¯ ) + h H ( 1 W h ) n h n 1 v ( y h d ¯ ) )
where,
v ( y h d ¯ ) = n h 1 ( i = 1 n h ( y h i d y h d ¯ ) 2 n h 1 ) .
The estimate of the total is then Y d ^ = A T y d ¯ with estimated variance v ( Y d ^ ) = A T 2 v ( y d ¯ ) where AT is the total area of the population. The standard error of the estimate is S E ( Y d ^ ) = v ( Y d ^ ) 0.5 .
In practice, yhid is either a proportion (e.g., proportion of the ith plot in stratum h and domain d that was forested) or on an areal unit basis (e.g., tree volume m3 ha−1). The sampling frame is area-based and the sample is equal probability, so the sampling weights are AT/n and the estimation weights (assuming Wh is known) are ATWh/nh. For estimation, all panels are combined, and the same stratification is applied throughout and, thus, weights are constant across panels and disturbance years for a given stratum.

2.3. Sampling Weight Adjustment

The rotating panel design is assumed to produce an equal probability sample of the population for each panel and across panels. However, we contend that the retrospective observation of disturbance combined with the panel design results in an unequal probability sample when multiple panels are used for estimating the areal extent of disturbance for a single year (Table 2). Table 2 shows an example with field measurements (indicated by M) collected from 2006 through 2010 and the year that a disturbance could be reported for each p panel (indicated by x). Note that disturbances that occurred in 2002 can only be observed during the measurement of panel 1 in 2006 and, likewise, disturbances that occurred in 2010 can only be observed in panel 5. However, disturbances that occurred in 2006 can be observed in all panels in the example.
The overall sampling weight Ω = (AT/n) arises from the sample design and is the inverse of the inclusion probability. Given that AT is fixed for any particular population, it is inappropriate to alter the weight without justifying a change in n. However, from Table 2, we note that the realized sample is not balanced across disturbance years. For example, if the domain of interest is disturbance in 2003 then nd = 2003 = np = 1 + np = 2 because 2003 disturbances could only be observed in panels 1 and 2 of 5 panels. The recalculation of n for each disturbance year alters the sampling weights by disturbance year. Based on the effective sample size (nd) the sampling weights differ by disturbance year. The effective Ωd are approximately 5 times larger than Ω for disturbance years d = 2002 and d = 2010, 2.5 times larger for d = 2003 and d = 2009, 1.67 times larger for d = 2004 and d = 2008, 1.25 times larger for d = 2005 and d = 2007, and the same for d = 2006 ( A T p = 1 5 n p = A T n ) , the only year in which all panels in the example have a field observation.
The effective sample size can be calculated as nd = nD, where n is a 1 × P matrix of the number of field plots in each panel and D is a P × D design matrix, where each row corresponds to the disturbance years that can be observed for each panel. D can take on numerous forms, depending on how many years of retrospective observations are considered. If one assumes that a disturbance and the disturbance year can be collected since the last plot measurement then
D = 1 1 1 1 1 0 0 0 0 0 1 1 1 1 1 0 0 0 0 0 1 1 1 1 1 0 0 0 0 0 1 1 1 1 1 0 0 0 0 0 1 1 1 1 1
This same approach can be used to determine the sample size for post-stratified estimation of disturbance. In this case nhd = nhD where nh is an H × P matrix of the number of plots by panel for each h stratum and the resulting nhd matrix is H × D:
n h d = n h = 1 , d = 1 n h = 1 , d = D n h = H , d = 1 n h = H , d = D
The total sample size for domain d is nd = 1nhd where 1 is a 1 × D matrix of ones. The elements from nhd and nd corresponding to h and d in Equations (2)–(4) are used:
y h d ¯ = n h d 1 i = 1 n h d y h i d ,
v ( y d ¯ ) = n d 1 ( h H W h n h d v ( y h d ¯ ) + h H ( 1 W h ) n h d n d 1 v ( y h d ¯ ) ) ,
v ( y h d ¯ ) = n h d 1 ( i = 1 n h d ( y h i d y h d ¯ ) 2 n h d 1 ) .
The post-stratified estimator as presented in Equations (1)–(4) has not changed; rather, the correct within stratum and across strata sample size is used for the time domain of interest. It is important to note that the original design-based sampling weights (AT/n) and the post-stratified weights (Ah/nh) are incorrect and the correct post-stratified weights are Ah/nhd.

2.4. Monte Carlo Simulation

We considered two estimation approaches. The PS (post-stratified) approach employed the post-stratified estimator as presented in equations 1 through 4. The WAPS (weight adjusted post-stratified) method employed the weight adjustment as presented in equations 1, and 5 through 7. To compare these approaches, we identified three populations from the forest disturbance time-series maps [11,24], which were considered the true population totals for area of disturbance. These maps identified forest disturbance, causal agent, and disturbance year from 1994–2010 at a 30 m spatial resolution. We used these populations for a Monte Carlo simulation. We selected populations in southwest Arkansas, the northern interior portion of California, and the Aroostook region in Northern Maine (Figure 1). The disturbances of interest were fire in the California population, and harvest in the Arkansas and Maine populations. In the Arkansas population (AT = 3.6 million ha), there were 2.87 million ha of forest and harvest occurrence steadily decreased from 2002–2010. The California population (AT = 5.8 million ha) had 3.44 million ha of forest and experienced episodic fire in 2006 and 2008. The Maine population (AT = 1.8 million ha) had 1.58 million ha of forest and experienced relatively stable harvest from 2002–2010. The NLCD Percent Tree Canopy Cover layer (Homer et al. 2016) was used to create three strata for each population. For the Arkansas population, the strata were 0–9% canopy cover (h = 1), 10–79% canopy cover (h = 2), and 80–100% canopy cover (h = 3). The strata for the California population were 0–9% canopy cover (h = 1), 10–49% canopy cover (h = 2), and 50–100% canopy cover (h = 3). For the Maine population the strata were 0–9% canopy cover (h = 1), 10–69% canopy cover (h = 2), and 70–100% canopy cover (h = 3). The strata boundaries for each population were identified following the Westfall et al. (2011) recommendations on minimum within strata and total sample sizes. Weights (Wh) were recorded for each stratum in each population.
Each population was sampled based on the FIA hexagonal sampling frame (Bechtold and Patterson 2005), where each hexagon was approximately 2400 ha and contained one sample plot. Each hexagon was assigned to 1 panel in a 5 panel design, and, for our purposes, we assumed panel 1 was measured in 2006, panel 2 was measured in 2007, etc., and panel 5 was measured in 2010. Step 1 was to generate a random plot center location within each hexagon. Step 2 was to assign the plot to a stratum. Step 3 was to generate the plot design based on each plot center. The plot design was a four-point cluster with point 1 being subplot 1 and plot center. Three additional subplots were located 36.6 m away from plot center at azimuths of 0 degrees, 120 degrees, and 240 degrees. Each subplot was 0.017 ha in size. In Step 4, the four subplots were then intersected with the populations, and the proportion of area disturbed across subplots by disturbance year was recorded as the observation for the ith plot. In calculating the proportion disturbed we only allowed for disturbances that could have been observed given the measurement year of the panel. For example, if the measurement year was 2010 (panel 5) then only disturbances from 2006–2010 were recorded. In step 5, the PS and the WAPS estimates (total ( Y ^ ) and the standard error of the total ( S E ( Y ^ ) ) ) of disturbed area by disturbance year were constructed for each population. Steps 1–5 were repeated B = 20,000 times. The root mean square error (RMSE) and bias were calculated based on the 20,000 replicates for each population. Specifically, R M S E = b B ( Y ^ b Y ) 2 B and b i a s = b B ( Y ^ b Y ) B where Y was the true population value for each disturbance year and population.
To understand the characteristics of the PS and WAPS approaches we examined the results in the following ways. First, for each approach, we compared the distribution of Y ^ from the Monte Carlo replicates to the true population total (Y). Second, we examined the bias of the PS and WAPS approaches. Third, for each approach, we compared the distribution of SE( Y ^ ) from the Monte Carlo replicates to the RMSE across Monte Carlo replicates.

3. Results

There were distinct differences between the PS and WAPS estimate of the total. Overall, the PS method tended to underestimate the true population total where the inter-quartile range (IQR) of PS estimates generally failed to include the true population value (Figure 2). There were two exceptions. Exception 1 was for the 2006 disturbance year. This exception was because the sample was properly weighted for that year. Exception 2 was for the CA population, where estimates in the beginning of the period (2002–2007) exhibited little underestimation. This was because there was little fire disturbance at the beginning of the time-period, and proper weighting of the sample was, hence, less consequential. The underestimation of the PS approach generally increased towards the beginning of the time-period (2002) and the end of the time period (2010) considered. The width of the IQR was smaller for the PS estimates of the total as compared to the WAPS estimates. The WAPS estimates of the total generally followed the true population total where in all cases the true total fell within the IQR of the WAPS estimates.
Neither the WAPS nor PS method was unbiased (Figure 3). The PS method yielded biased estimates across years and populations, however, estimates for 2002 and 2010 generally showed the largest bias. The exception for the CA population was due to little disturbance in the beginning and end of the time period. Estimates arising from the PS method typically were negatively biased (i.e., underestimates). While the WAPS method produced estimates with less bias than the PS method, the WAPS method did not result in a mean bias of zero for any particular estimate. However, the WAPS method did not tend to either underestimate or overestimate disturbance area.
The standard errors (SE) of the PS estimates were generally smaller than the SE of the WAPS estimates (Figure 4). This was because under the WAPS approach the sample size was adjusted downward for each disturbance year, except 2006, when the sample was the same for the two methods. One would generally expect the distribution of the SE of the estimate to be similar to the RMSE from a simulation approach. The IQR of the SE of the WAPS estimate generally contained the RMSE. However, there were large differences between the IQR of the SE of the PS estimate, except for the 2006 disturbance years, and when there was little disturbance in the CA population (Figure 4).

Case Study

To demonstrate the WAPS method, we constructed estimates of the extent and timing of disturbances in east Texas as well as tree mortality volume. We selected east Texas as the case study area, owing to the occurrence of several prominent disturbance events that occurred during the period of annual data collection 2001–2018 [25]. East Texas is an area encompassing approximately 9 million hectares, of which about 4.9 million hectares is forest land [26]. The temperate forest biome that covers much of the southeastern US reaches its southwestern extent in east Texas [27,28]. In 2005, Hurricane Rita passed through east Texas with sufficient force to impact 6 percent of the total timber volume in east Texas [29]. In 2008 a second hurricane, Ike, impacted 4 percent of the total timber volume [30]. In 2011, an exceptional drought occurred that resulted in the mortality of approximately 4 percent of all trees in the region [31,32,33]. Fire was also prominent in 2011, during the drought event and in subsequent years. The expectation is properly constructed estimates will show spikes of wind disturbance in 2005 and 2008 and drought and fire in 2011.
Based on the design and results of the simulation study one might consider just using the panels needed to construct estimates for the year(s) of interest, however, the execution of the rotating panel design is rarely as simple as described in the simulation. The actual measurement of a panel may span more than one year and, hence, the remeasurement period rarely equals the number of panels. Further, some plots may be measured out of panel because of temporary hazardous conditions and other practical problems that field crews and inventory implementers face.
Trends are similar for both disturbance area and mortality volume (Figure 5). For consistency with area estimates, volumes reported are mortality (any cause) on the conditions of noted disturbance. Estimates show large spikes of wind disturbance occurring in 2005 and 2008, which coincide with Hurricanes Rita and Ike, respectively. Hurricane Rita is generally believed to have had greater impact on timber resources than Hurricane Ike [29,30], but estimates here suggest similar impacts, at least in terms of area disturbed and mortality volume. The first sign of drought impacts is in 2011, the year in which the worst 1-year drought in Texas history occurred. Drought impacts were more severe than either hurricanes, a finding consistent with another study that looked at the impacts of three events on standing dead trees [25]. The drought affected at least 48% more area and 62% more mortality volume than either hurricane. Substantial drought impacts are estimated for 2012, but much reduced in 2013–2015. Fire is also most prominent in 2011, as expected. Interestingly, fire has continued to affect sizeable area of east Texas 2012–2018, but impacts on mortality volume have not been very large. Fire estimates reported here are combined estimates of three different fire categories: crown and ground fire, ground fire, and crown fire. Crown fires are not that common, with 2011 being a notable exception where approximately 50% of the fires were crown fires. More typical is ground fire, which likely explains the smaller mortality volume even when sizeable areas have been impacted by fire as is the case in recent years (e.g., approximately 93% of fires were ground fires in 2018).

4. Discussion

The goal of this research was to develop an appropriate approach for estimating the annual areal extent of disturbance under a rotating panel sample design when multiple panels are combined for estimation. Our literature review revealed no guidance for this problem, a finding that we contend is because the issue is not with the estimator, but rather the sampling weights. Most of the available literature focuses on estimation assuming that sampling weights are defined during the initial survey design. The design-based weights of the FIA sample are equal within a population; however, considering the nature of retrospective observation, the correct sampling weights are unequal when multiple panels are considered. The change in the sampling weights for estimating annual disturbance area are not immediately clear to the practitioner. Further, standard estimation tools do not take the correct sampling weights into account when constructing estimates of annual disturbance extent. The exception is when the disturbance year of interest is the mid-point (i.e., 2006 in the example). We have demonstrated that the post-stratified estimator is suitable to estimate annual disturbance totals, but the effective sample size and, hence, sampling weights, must be calculated. It turns out to be a rather simple task to correctly calculate the effective sample size and sampling weights for annual disturbance estimates. These results are relevant to other NFIs that use rotating panel design and collect disturbance information in a retrospective fashion.
While we considered the WAPS methods suitable to estimate annual disturbance totals, it is important to note that the estimates were biased based on the simulation results. The percent bias of the WAPS method was 5.8%, across years and populations. The bias arose because the disturbances were rare events or small domains and, under traditional NFI sampling intensities, only a small number of disturbed plots may be observed in a given year. We note that bias was typically minimized for 2006 estimates (Figure 3). This was because the effective sample size included all plots from all panels (Table 2), which further suggests that increased sample sizes may be required to achieve unbiased estimates of disturbance extent by year. As suggested by Kish [15] and Cochran [14], there are sampling techniques that are appropriate for rare events. However, most NFIs were not designed as rare event surveys and, further, their sample designs and intensities are relatively fixed at the national level. Alternative estimators such as, model-assisted, model based, and small area techniques may prove valuable for reducing the bias and increasing precision of annual disturbance estimates.
The analysis presented examined the performance of the WAPS approach for estimating disturbance extent and timing. The reason for this emphasis was because information was available to serve as a true population in the Monte Carlo simulation. However, our results are applicable to a suite of domain estimates that can be constructed given the variables collected as part of broad-scale forest inventories. For example, one may use the WAPS technique to estimate the volume of trees that died due to fires in a given year, because tree mortality year and cause are also collected, retrospectively. Likewise, annual estimates of removal volumes can also be constructed. There are numerous other estimates relating to disturbance impacts that can be constructed using the proposed technique.
There are other methods to account for the retrospective observations when all panels are combined for estimation. One straightforward approach is to annualize the observations. Observations can be annualized by simply dividing the observed value from the plot by the length of the measurement interval for the plot. This technique is currently used by FIA when constructing estimates of growth, removals, and mortality [34]. If the retrospective observations of disturbance are annualized then no adjustment to the sampling weights is needed. However, the average annual estimate of the disturbed area for any population would cover an approximately 2P year span. Actual annual estimates are more informative for addressing current resource issues than an estimate of average annual disturbance area over a 10+ year time span. Further, as Roesch [35] points out, these average annual estimates, calculated over different measurement intervals, are different population parameters.
Our focus here was on developing an estimation technique for retrospective observations under a rotating panel design, and we did not address measurement error issues related to plot-level retrospective observation. Most NFIs have published quality objectives, and those of the FIA program suggest that the year of disturbance should be within one-year of the actual disturbance year 99% of the time [20]. Pollard and colleagues [36] found that this quality objective is generally met. However, the study by Pollard and colleagues [36] should be updated given the expanding FIA data volume available. In some cases, a practitioner may wish to only consider disturbances (and year of disturbance) recorded within the three-years of plot measurement (for example). In this case, it is relatively straightforward to modify the D matrix to take into account that only three years of retrospective observations are being used to construct the estimates.

5. Conclusions

Disturbance at varying frequencies, intensities, and spatial scales is an important component of forest systems. Understanding the extent of disturbances and their impacts is critical to formulating effective management and policy responses. A key information source is large-area forest inventory data, such as those produced across the US by the FIA program. Due to the retrospective nature of the observations in a panelized design setting, the use of standard estimators that generate typical forest resource statistics can lead to inaccurate representations of the areal extent and timing of disturbance. In this paper, the post-stratified estimator was extended to account for the specific manner in which disturbance data are recorded. In addition to the theoretical development, the estimates exhibited satisfactory agreement with the true population values when assessed in a simulation environment. Broadly, the increasing use of panelized forest inventory designs globally suggests it will often be critical for researchers to account for the necessary weighting adjustments by observation year of the panel when estimates of areal extent and timing of disturbances are desired. Thus, we recommend the methods presented in this paper be adopted or the concepts generally incorporated into other types of forest inventory estimation frameworks.

Author Contributions

J.W.C. conceived the research, developed the methodology, performed the formal analysis, wrote the original draft and the wrote final draft of the manuscript. C.B.E. performed the formal analysis, wrote the original draft, and wrote the final draft of the manuscript. J.A.W. developed methodology, wrote the original draft, and wrote the final draft of the manuscript. M.E.T. wrote the original draft and wrote the final draft of the manuscript. All authors have read and agreed to the published version of the manuscript.

Funding

Funding to support CBE in this effort was provided by the USDA Forest Service Project 15 FIA–Forest Biometrics Program Support (RJVA 15-JV-11242305-100). J.W.C., J.A.W., and M.E.T. received no external funding in support of this research.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Masek, J.G.; Goward, S.N.; Kennedy, R.E.; Cohen, W.B.; Moisen, G.G.; Schleeweis, K.; Huang, C. United States forest disturbance trends observed using Landsat time series. Ecosystems 2013, 16, 1087–1104. [Google Scholar] [CrossRef] [Green Version]
  2. Franklin, J.F.; Spies, T.A.; Van Pelt, R.; Carey, A.B.; Thornburgh, D.A.; Berg, D.R.; Lindenmayer, D.B.; Harmon, M.E.; Keeton, W.S.; Shaw, D.C.; et al. Disturbances and structural development of natural forest ecosystems with silvicultural implications, using Douglas-fir forests as an example. For. Ecol. Manag. 2002, 155, 399–423. [Google Scholar] [CrossRef]
  3. Flannigan, M.D.; Stocks, B.J.; Wotton, B.M. Climate change and forest fires. Sci. Total Environ. 2000, 262, 221–229. [Google Scholar] [CrossRef]
  4. Seidl, R.; Thom, D.; Kautz, M.; Martin-Benito, D.; Peltoniemi, M.; Vacchiano, G.; Wild, J.; Ascoli, D.; Petr, M.; Honkaniemi, J.; et al. Forest disturbances under climate change. Nat. Clim. Chang. 2017, 7, 395–402. [Google Scholar] [CrossRef] [Green Version]
  5. Karel, T.H.; Man, G. (Eds.) Forest Service. Major Forest Insect and Disease Conditions in the United States: 2015; Department of Agriculture Forest Service, Forest Health Protection: Washington, DC, USA, 2017; p. 45.
  6. Vidal, C.; Alberdi, I.; Hernández, L.; Redmond, J. (Eds.) National Forest Inventories: Assessment of Wood Availability and Use; Springer: Cham, Switzerland, 2016. [Google Scholar]
  7. Brooks, E.B.; Thomas, V.A.; Wynne, R.H.; Coulston, J.W. Fitting the multitemporal curve: A Fourier series approach to the missing data problem in remote sensing analysis. IEEE Trans. Geosci. Remote Sens. 2012, 50, 3340–3353. [Google Scholar] [CrossRef] [Green Version]
  8. Goward, S.N.; Waring, R.H.; Dye, D.G.; Yang, J. Ecological remote sensing at OTTER: Satellite macroscale observations. Ecol. Appl. 1994, 4, 322–343. [Google Scholar] [CrossRef]
  9. Huang, C.; Goward, S.N.; Masek, J.G.; Thomas, N.; Zhu, Z.; Vogelmann, J. An automated approach for reconstructing recent forest disturbance history using dense Landsat time series stacks. Remote Sens. Environ. 2010, 114, 183–198. [Google Scholar] [CrossRef]
  10. Moisen, G.G.; Meyer, M.C.; Schroeder, T.A.; Liao, X.; Schleeweis, K.G.; Freeman, E.A.; Toney, C. Shape selection in Landsat time series: A tool for monitoring forest dynamics. Glob. Chang. Biol. 2016, 22, 3518–3528. [Google Scholar] [CrossRef]
  11. Schleeweis, K.G.; Moisen, G.G.; Schroeder, T.A.; Toney, C.; Freeman, E.A.; Goward, S.N.; Huang, C.; Dungan, J.L. US National Maps Attributing Forest Change: 1986–2010. Forests 2020, 11, 653. [Google Scholar] [CrossRef]
  12. Eidenshenk, J.; Schwind, B.; Brewer, K.; Zhu, Z.; Quayle, B.; Howard, S. A project for monitoring trends in burn severity. Fire Ecol. 2007, 3, 3–21. [Google Scholar] [CrossRef]
  13. Kolden, C.A.; Smith, A.M.; Abatzoglou, J.T. Limitations and utilization of Monitoring Trends in Burn Severity products for assessing wildfire severity in the USA. Int. J. Wildland Fire 2015, 24, 1023–1028. [Google Scholar] [CrossRef]
  14. Cochran, W.G. Sampling Techniques; John Wiley & Sons: New York, NY, USA, 1977; 428p. [Google Scholar]
  15. Kish, L. Survey Sampling; John Wiley & Sons: New York, NY, USA, 1995; 643p. [Google Scholar]
  16. Särndal, C.E.; Swensson, B.; Wretman, J. Model Assisted Survey Sampling; Springer: New York, NY, USA, 1992; 694p. [Google Scholar]
  17. Bechtold, W.A.; Patterson, P.L. (Eds.) The Enhanced Forest Inventory and Analysis Program–National Sample Design and Estimation Procedures; Gen. Tech. Rep. SRS-80; U.S. Department of Agriculture, Forest Service, Southern Research Station: Asheville, NC, USA, 2005; 85p.
  18. Diaz-Ramos, S.; Stevens, D.L.; Olsen, A.R. EMAP Statistics Methods Manual; U.S. Environmental Protection Agency, Office of Research and Development, National Health and Environmental Effects Research Laboratory: Corvallis, OR, USA, 1996; p. 13.
  19. Fay, R.E., III; Herriot, R.A. Estimates of income for small places: An application of James-stein procedures to census data. J. Am. Stat. Assoc. 1979, 74, 269–277. [Google Scholar] [CrossRef]
  20. Northern Research Station. Forest Inventory and Analysis Nation Core Field Guide Volume 1: Field Data Collection Procedures for Phase 2 Plots; United States Department of Agriculture, Forest Service: Madison, WI, USA, 2018.
  21. Westfall, J.A.; Patterson, P.L.; Coulston, J.W. Post-stratified estimation: Within-strata and total sample size recommendations. Can. J. For. Res. 2011, 41, 1130–1139. [Google Scholar] [CrossRef] [Green Version]
  22. Homer, C.; Dewitz, J.; Yang, L.; Jin, S.; Danielson, P.; Xian, G.; Coulston, J.; Herold, N.; Wickham, J.; Megown, K. Completion of the 2011 National Land Cover Database for the conterminous United States–representing a decade of land cover change information. Photogramm. Eng. Remote Sens. 2015, 81, 345–354. [Google Scholar]
  23. Coulston, J.W.; Moisen, G.G.; Wilson, B.T.; Finco, M.V.; Cohen, W.B.; Brewer, C.K. Modeling percent tree canopy cover: A pilot study. Photogramm. Eng. Remote Sens. 2012, 78, 715–727. [Google Scholar] [CrossRef]
  24. Zhao, F.; Huang, C.; Goward, S.N.; Schleeweis, K.; Rishmawi, K.; Lindsey, M.A.; Denning, E.; Keddell, L.; Cohen, W.B.; Yang, Z.; et al. Development of Landsat-based annual US forest disturbance history maps (1986–2010) in support of the North American Carbon Program (NACP). Remote Sens. Environ. 2018, 209, 312–326. [Google Scholar] [CrossRef]
  25. Edgar, C.B.; Westfall, J.A.; Klockow, P.A.; Vogel, J.G.; Moore, G.W. Interpreting effects of multiple, large-scale disturbances using national forest inventory data: A case study of standing dead trees in east Texas, USA. Forest Ecol. Manag. 2019, 437, 27–40. [Google Scholar] [CrossRef]
  26. Cooper, J.A.; Bentley, J.W. East Texas, 2011—Forest Inventory and Analysis Factsheet; E-Science Update SRS-052; U.S. Department of Agriculture Forest Service: Asheville, NC, USA, 2012; p. 5.
  27. Burns, R.M.; Honkala, R.H. Silvics of North America: Volume 1 Conifers; Agriculture Handbook 654; U.S. Department of Agriculture, Forest Service: Washington, DC, USA, 1990; p. 675.
  28. Burns, R.M.; Honkala, R.H. Silvics of North America: Volume 2 Hardwoods; Agriculture Handbook 654; U.S. Department of Agriculture, Forest Service: Washington, DC, USA, 1990; p. 877.
  29. Texas Forest Service, Hurricane Rita Timber Damage Assessment. Available online: http://tfsweb.tamu.edu/uploadedFiles/TFSMain/Data_and_Analysis/Forest_Economics_and_Resource_Analysis/Resource_Analysis/Resource_Analysis_publications/HurricaneRitaTimberDamageAssessment2005.pdf (accessed on 13 November 2017).
  30. Texas Forest Service, Hurricane Ike Timber Damage Assessment. Available online: http://tfsweb.tamu.edu/uploadedFiles/TFSMain/Data_and_Analysis/Forest_Economics_and_Resource_Analysis/Resource_Analysis/Resource_Analysis_publications/HurricaneIke_Timber-Damage-Assessment.pdf (accessed on 13 November 2017).
  31. Nielsen-Gammon, J.W. The 2011 Texas drought. Tex. Water J. 2012, 3, 59–95. [Google Scholar]
  32. Moore, G.W.; Edgar, C.B.; Vogel, J.G.; Washington-Allen, R.A.; March, R.G.; Zehnder, R. Tree mortality from an exceptional drought spanning mesic to semiarid ecoregions. Ecol. Appl. 2016, 26, 602–611. [Google Scholar] [CrossRef]
  33. Klockow, P.A.; Vogel, J.G.; Edgar, C.B.; Moore, G.W. Differences in lagged mortality among tree species four years after an exceptional drought in east Texas. Ecosphere 2018, 9, e02455. [Google Scholar] [CrossRef] [Green Version]
  34. Scott, C.T.; Bechtold, W.A.; Reams, G.A.; Smith, W.D.; Westfall, J.A.; Hansen, M.H.; Moisen, G.G. Sampled-based estimators used by the Forest Inventory and Analysis national information management system. In The Enhanced Forest Inventory and Analysis Program–National Sampling Design and Estimation Procedures; General Technical Report SRS-80; Bechtold, W.A., Patterson, P.L., Eds.; U.S. Department of Agriculture, Forest Service, Southern Research Station: Asheville, NC, USA, 2005; pp. 43–67. [Google Scholar]
  35. Roesch, F.A. Continuous inventories and the components of change. In Proceedings of the 2nd International Conference on Forest Measurements and Quantitative Methods and Management, Hot Springs, AK, USA, 15–18 June 2004. [Google Scholar]
  36. Pollard, J.E.; Westfall, J.A.; Patterson, P.L.; Gartner, D.L.; Hansen, M.; Kuegler, O. Forest Inventory and Analysis National Data Quality Assessment Report 2000–2003; General Technical Report RMRS-181; U.S. Department of Agriculture, Forest Service, Rocky Mountain Research Station: Fort Collins, CO, USA, 2006; p. 43.
Figure 1. 2002 to 2010 disturbance trajectories for the Arkansas, California, and Maine populations.
Figure 1. 2002 to 2010 disturbance trajectories for the Arkansas, California, and Maine populations.
Forests 11 01298 g001
Figure 2. Distribution of the estimated total ( Y ^ ) for the post-stratified (PS) method (black box and whisker plot) and for the weight adjusted post-stratified (WAPS) method (gray box and whisker plots) for each disturbance year and population. The solid grey line is the true population total (Y).
Figure 2. Distribution of the estimated total ( Y ^ ) for the post-stratified (PS) method (black box and whisker plot) and for the weight adjusted post-stratified (WAPS) method (gray box and whisker plots) for each disturbance year and population. The solid grey line is the true population total (Y).
Forests 11 01298 g002
Figure 3. Bias for each estimate of the population total using the PS and WAPS method.
Figure 3. Bias for each estimate of the population total using the PS and WAPS method.
Forests 11 01298 g003
Figure 4. Distribution of the standard error of estimated total (SE( Y ^ )) for the PS method (black box and whisker plot) and for the WAPS method (gray box and whisker plots) for each disturbance year and population. The dotted black line is the root mean square error (RMSE) from the PS method and the solid grey line is the RMSE from the WAPS method.
Figure 4. Distribution of the standard error of estimated total (SE( Y ^ )) for the PS method (black box and whisker plot) and for the WAPS method (gray box and whisker plots) for each disturbance year and population. The dotted black line is the root mean square error (RMSE) from the PS method and the solid grey line is the RMSE from the WAPS method.
Forests 11 01298 g004
Figure 5. Estimates of extent for drought, wind, and fire disturbances (top row) and corresponding estimates of mortality volume for each disturbance (bottom row). Error bars are one standard error.
Figure 5. Estimates of extent for drought, wind, and fire disturbances (top row) and corresponding estimates of mortality volume for each disturbance (bottom row). Error bars are one standard error.
Forests 11 01298 g005
Table 1. Hypothetical example showing retrospective disturbance data collected for a five-year rotating panel forest inventory. “1” indicates a disturbance was recorded, “0” indicates a disturbance wasn’t recorded, and “--“ denotes an unobservable.
Table 1. Hypothetical example showing retrospective disturbance data collected for a five-year rotating panel forest inventory. “1” indicates a disturbance was recorded, “0” indicates a disturbance wasn’t recorded, and “--“ denotes an unobservable.
Disturbance Year
PlotInventory YearPanel2004200520062007200820092010
12008110100----
220092--10001--
320103----01000
420103----00000
Table 2. Schematic depicting the Forest Inventory and Analysis (FIA) panel design where “M” identifies the year each panel (p) is measured and “x” identifies the corresponding years where disturbance can be retrospectively recorded by a field measurement. The total number of plots (n) from panels p = 1 to 5 are decomposed by panel (np) and disturbance year (nd). The overall sampling weight Ω = AT/n is also decomposed into panel sampling weight (Ωp) and disturbance year sampling weights (Ωd).
Table 2. Schematic depicting the Forest Inventory and Analysis (FIA) panel design where “M” identifies the year each panel (p) is measured and “x” identifies the corresponding years where disturbance can be retrospectively recorded by a field measurement. The total number of plots (n) from panels p = 1 to 5 are decomposed by panel (np) and disturbance year (nd). The overall sampling weight Ω = AT/n is also decomposed into panel sampling weight (Ωp) and disturbance year sampling weights (Ωd).
Disturbance Year
PanelInventory Year200220032004200520062007200820092010npΩp
12006xxxxM,x n 1 n P A T n 1
22007 xxxxM,x n 2 n P A T n 2
32008 xxxxM,x n 3 n P A T n 3
42009 xxxxM,x n 4 n P A T n 4
52010 xxxxM,x n 5 n P A T n 5
nd n d = p = 1 1 n p n d = p = 1 2 n p n d = p = 1 3 n p n d = p = 1 4 n p n d = p = 1 5 n p n d = p = 2 5 n p n d = p = 3 5 n p n d = p = 4 5 n p n d = p = 5 5 n p
Ωd A T n d A T n p A T n d A T 2 n p A T n d A T 3 n p A T n d A T 4 n p A T n d A T 5 n p A T n d A T 4 n p A T n d A T 3 n p A T n d A T 2 n p A T n d A T n p
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Coulston, J.W.; Edgar, C.B.; Westfall, J.A.; Taylor, M.E. Estimation of Forest Disturbance from Retrospective Observations in a Broad-Scale Inventory. Forests 2020, 11, 1298. https://0-doi-org.brum.beds.ac.uk/10.3390/f11121298

AMA Style

Coulston JW, Edgar CB, Westfall JA, Taylor ME. Estimation of Forest Disturbance from Retrospective Observations in a Broad-Scale Inventory. Forests. 2020; 11(12):1298. https://0-doi-org.brum.beds.ac.uk/10.3390/f11121298

Chicago/Turabian Style

Coulston, John W., Christopher B. Edgar, James A. Westfall, and Marcus E. Taylor. 2020. "Estimation of Forest Disturbance from Retrospective Observations in a Broad-Scale Inventory" Forests 11, no. 12: 1298. https://0-doi-org.brum.beds.ac.uk/10.3390/f11121298

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop