High-resolution mapping of aboveground biomass for forest carbon monitoring system in the Tri-State region of Maryland, Pennsylvania and Delaware, USA

The study area, Tri-State region of Pennsylvania (PA; 119 280 km²), Maryland (MD; 32 133 km²), and Delaware (DE; 6452 km²) in the USA (hereafter Tri-State region) covers a land area of approximately 157 865 km². The Tri-State region contains four major ecological regions that are generally uniform with respect to species-composition and environmental gradients (figure 2). These regions are the Outer Coastal Plain Mixed, Eastern Broadleaf, Central Appalachian Broadleaf-Coniferous, and Northeastern Mixed Forest, adapted from Bailey's ecoregions for the Conterminous United States (Bailey 1995). Forests, including woody wetlands, comprise 56% of the total area based on National Land Cover Dataset (NLCD) 2011 (Homer et al 2015). The topography varies substantially with elevations ranging from about sea level near the Delaware River to 1027 m in Western Maryland.

Zoom In Zoom Out Reset image size

The field data collected over the domain is built from two sections: USDA forest inventory analysis (FIA) plot data for model development, and independent field measurements collected by the University of Maryland (UMD) for product validation. For this study, 2243 forest plots that most closely matched the lidar year of acquisition over the domain were extracted from the FIA database (EVALIDator Version 1.6.0.03a). These plots, consist of measurements for individual tree height, diameter and species, were collected from year 2004 to 2013 in MD and DE, and from year 2005 to 2009 in PA, 92% of them had a less than 2 years difference (supplement S2 and table S2 is available online at stacks.iop.org/ERL/14/095002/mmedia).

Additionally, a total of 1163 variable-radius plots (VRP) were collected by UMD plots) between 2010 and 2013 in MD (Dubayah 2012, Huang et al 2015), in the summers of 2014 and 2015 for the northeastern forest in PA and DE, respectively. VRP locations were determined by a modified random stratified sampling scheme broken down by ecoregions (Bailey 1995), NLCD 2011 land cover classes (Home et al 2015), and lidar canopy height (figure 2). Individual tree species, dbh and maximum height were measured in the field.

Plot-level aboveground live biomass was estimated using allometric models. For the sake of comparison with different map products, two types of allometric relationships were prepared: (1) regional models from FIA's Component Ratio Method (CRM) (Heath et al 2008, Woodall et al 2011) and (2) national generalized models (JKS) from Jenkins et al (2003). County and state scale total biomass was calculated based on sample probability for each allometric equation (FIAcrm and FIAjks).

Lastly, both the FIA and UMD plot data were screened for outliers. The screening omitted plots located in no lidar data areas (i.e. missing tiles in Maryland) and excluded mismatched plots between biomass and lidar heights with three filtering rules: had zero biomass but lidar heights (>10 m); had large biomass (>500 Mg ha⁻¹) but low lidar heights (<10 m); and small biomass (<50 Mg ha⁻¹) but high lidar heights (>30 m). These filters reduced the number of FIA plots from 2243 to 2215, and the number of UMD plots from 1163 to 1046. Summary statistics of the FIA and UMD plots are given in table 1.

Discrete return lidar data were collected through different county and state government agencies (summarized in supplement S2 and table S1) over the Tri-State region (figure 2(b)). These airborne lidar collections were originally dedicated to surveying digital elevation model (DEM). Therefore leaf-off season (December–January) were picked for flight. Lidar first returns were interpolated into digital surface model (DSM); last returns were utilized to derive a bare earth DEM using quick terrain modeler and other GIS-based software (O'Neil-Dunne et al 2013, 2014). Then DSM and DEM were differenced to obtain a normalized difference surface model (nDSM).

High-resolution tree canopy cover map at 1 m resolution for circa 2011 over the domain was created using an object-based approach that integrated the lidar data and high-resolution multi-spectral images from the National Agricultural Imagery Program or NAIP (O'Neil-Dunne et al 2013, 2014, USDA 2018, supplement S3). A canopy height model (CHM) was then extracted by applying the canopy cover mask to the nDSM, thereby creating a height model that excluded non-canopy height.

The definition of tree canopy used for high-resolution [1 m] product was a minimum height of 2.5 m and a minimum area of 4 m². This definition is wider than the FIA definition of forest, which only refers to forested areas that are not developed for nonforest land uses, at least an acre in size and 120' in width, and contain a live plus missing canopy cover of at least 10% (USFS 2018a).

A set of percentage canopy cover and mean canopy height maps were generated from aggregating the 1 m tree canopy cover and tree height maps. The percentage tree canopy cover was aggregated to 30 m. From the CHM, the max tree height within each 10 m cell was first calculated and then averaged to 30 m resolution (hereafter referred to as average-maximum height).

A set of forestry metrics were prepared, including percentile heights, density metrics, mean heights, canopy covers, and topographical indices (table 2). Specifically, metrics were calculated at plot-scale and pixel-scale for modeling and mapping purpose. At plot-scale, we located the extent of a given plot using the plot centroid and calculated metrics within a given size of plot (e.g. 30 m). These metrics have been found high correlated with biomass (Goetz and Dubayah 2011, Swatantran et al 2012, Huang et al 2013, Whitehurst et al 2013, Zhao et al 2018, Hurtt et al 2019), shown in figure S2. At pixel-scale, we calculated metrics at 30 m resolution for each pixel within the domain. Countywide maps were created first and merged into a seamless statewide map aligned to a common extent.

Four national and a global scale biomass maps were prepared for map comparisons, summarized in table 3. The national biomass and carbon dataset (NBCD) 2000 was the first 30 m national product developed using medium resolution satellite optical imagery from Landsat enhanced thematic mapper plus (ETM+) and InSAR data from shuttle radar topography mission (SRTM) (Kellndorfer et al 2012). Two versions of the NBCD biomass maps were extracted for consistency with other national maps, including the FIA database version using CRM for allometry (NBCDcrm) and the nationally consistent allometric models (NCE) version using JKS allometry models (NBCDjks). The Saatchi map is a 90 m national-scale biomass map generated using optical, radar and lidar data and forest inventory data (Saatchi et al 2012). Both the Blackard and Wilson were 250 m national biomass datasets generated by USDA forest service two generated by USFS (Blackard et al 2008a, 2008b, Wilson et al 2013b). The GlobBiomass map was a 100 m global biomass product converted from growing stock volume (Santoro 2018). Detailed descriptions for these maps are given in supplement S6.

For forest mask, the National Land Cover Dataset or NLCD 2011 (30 m) (Homer et al 2015) were utilized to differentiate between forested and non-forested areas for map comparisons. Specifically, woody forestlands (deciduous, conifers and mixed) and woody wetlands were included as forested area. From the NLCD forest mask, we applied a 20% threshold to aggregate the mask from 30 to 250 m. For consistency and statistical measurement purposes, all biomass and land cover maps were projected into Albers Equal-area Conic projection. All pixels were aligned to NLCD 2011 30 m cell.

Random Forest models were developed using 80% of FIA plots to predict field-estimated biomass as a function of forestry and topographical metrics. Quantile Regression Forests (QRF), which predicts a particular quantile (10% and 90%), was employed to provide model error bounds (Meinshausen 2006). The R packages 'randomForest' and 'quantregForest' were used for modeling analyses (Liaw and Wiener 2002, R Development Core Team 2016). We employed several statistical measurements such as coefficient of determination (R²), root-mean-squared error (RMSE), RMSE CV (%), and mean estimated bias error (MBE) to validate results at plot-scale.

$$\begin{eqnarray}&&{\rm{R}}{\rm{M}}{\rm{S}}{\rm{E}}=\sqrt{\displaystyle \sum _{i=1}^{n}\displaystyle \frac{{\left({y}_{i}-{x}_{i}\right)}^{2}}{n}}\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&{\rm{R}}{\rm{M}}{\rm{S}}{\rm{E}}\,{\rm{C}}{\rm{V}}=100\times \displaystyle \frac{{\rm{R}}{\rm{M}}{\rm{S}}{\rm{E}}}{\bar{x}}\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&{\rm{M}}{\rm{B}}{\rm{E}}=\displaystyle \frac{\displaystyle {\sum }_{i=1}^{n}\left({y}_{i}-{x}_{i}\right)}{n}\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&{R}^{2}=1-\displaystyle \frac{\displaystyle {\sum }_{i=1}^{n}{\left({y}_{i}-{x}_{i}\right)}^{2}}{\displaystyle {\sum }_{i=1}^{n}{\left({x}_{i}-\bar{x}\right)}^{2}},\end{eqnarray} \tag{ 4 }$$

where x_i is the value from field measurements, y_i is the predicted value from RF model; i is the sample index; and $\bar{x}$ and $\bar{y}$ are the means of field measurements and predicted values respectively; and n is the sample size.

Models developed for each ecoregion were applied to make countywide predictions for a total of 94 counties. Specifically, Outer Coastal model was applied to 16 counties, Eastern Broadleaf model was applied to 31 counties, Central Appalachian model was applied to 21 counties, and Northeastern Mixed model was applied to 26 counties. Details about random forest models were given in supplement S2.

Lastly, county-scale mean biomass density (Mg ha⁻¹) and biomass stock (Tg) over DE, MD, and PA were summarized using average and sum methods. The Maryland canopy height map had gaps comprising around 1.5% of the state because of data collection restrictions and quality issues. To overcome this, these data gaps were filled with values from the NBCD product. Specifically, the height values were from weighted basal area heights (NBCD_BAW) and biomass values were from the national consistent allometry equation version of biomass (NBCDjks) (Huang et al 2015).

At the plot-scale, we compared estimates of biomass density from this study (RFjks) to with-held FIA plots and UMD plots. Explicitly, we conducted plot-scale validation using the 20% with-held FIA plots. Then, we used all UMD plots as independent plot-scale validations. In addition, we test the scale effect by comparing the RF predicts from 30 m and 90 m metrics to estimates from FIA plots (supplement S2). Maps of uncertainty were generated for every 30 m pixel. RF predicts the mean response for any given location. Prediction intervals were developed using QRF. QRF is similar to RF but it retains all responses and allows for constructing prediction intervals. We generated the 10th and 90th prediction percentiles and analyzed the spatial distribution of errors across the state. The uncertainty index is calculated as:

$$\begin{eqnarray}&&U=\displaystyle \frac{{\rm{P}}{\rm{e}}{\rm{r}}{\rm{c}}{\rm{t}}{90}_{i}-{\rm{P}}{\rm{e}}{\rm{r}}{\rm{c}}{\rm{t}}{10}_{i}}{{y}_{i}},\end{eqnarray} \tag{ 5 }$$

where y_i is the predicted value from RF model; Perct10_i and Perct90_i are the predicted values from 10th and 90th prediction percentiles, respectively.

Estimates of biomass from this study (RFjks) were compared to four national and one global products (NBCDjks, Saatchi, Blackard, Wilson, and GlobBiomass) at the pixel, county and state scales following the framework we setup in previous studies in Maryland and California (Huang et al 2015, 2017).

At the pixel scale, we compared biomass maps at 250 m spatial resolution, the coarsest resolution among all maps in this study. We choose to use this resolution because previous analysis found that the map comparisons at coarser resolution are more reliable (Huang et al 2017). The aggregation process reduces the difference caused by geolocation and misalignment between different maps. Statistical measurements such as R², root-mean-squared difference (RMSD, calculated similar to RMSE) and mean estimated bias difference (MBD, calculated similar to MBE), were selected to evaluate results at pixel-scale.

At the county to the state scale, we compared estimates from five map products to FIA estimates. Especially, we compared the map-based county and state-scale estimates of aboveground biomass from this study, four national products (NBCDjks, Saatchi, Wilson, and Blackard), a global product (GlobBiomass) and sample-based estimates from FIA plots using two different allometric equations (FIAjks and FIAcrm). Additionally, at the state- to Tri-State region scales, we compared estimates of mean carbon density and stock from map-based results to sample-based estimates from FIA. A coefficient of 0.5 was used to convert from biomass to carbon (Chave et al 2005). We used similar statistical indicators (i.e. R², RMSE, and MBE) as described in section 2.3 to evaluate the map-based results.

Visually, the spatial patterns of canopy cover, mean canopy height and biomass matched land cover and topography (figure 3). In general, counties in western Maryland and northern Pennsylvania had the highest mean canopy heights, largest forest cover, and total and mean biomass. The Blue Ridge region had large biomass while the coastal plains had smaller biomass in general. Biomass distribution was patchy and heterogeneous as expected from the patchy nature of the landscape. The map of uncertainty index derived from QRF models shows similar spatial patterns matched canopy cover and height (figure S2).

Zoom In Zoom Out Reset image size

High-resolution maps reveal more details (figure 4). The canopy cover map at 1 m resolution detailed delineating out important urban forests that traditionally defined as 'non-forest', such as open, low and medium high developed land cover types. The canopy height map revealed not only the extent of forest but also vertical information, which is critical for estimating carbon in forests. Larger biomass values were observed in urban areas that have tall trees and residential neighborhoods that have tree plantation programs (figure S6(a)). Biomass values were larger near water bodies, as canopy heights were higher in protected wetlands (figure S6(b)).

Zoom In Zoom Out Reset image size

Plot-scale validation of biomass density showed different levels of agreement when using FIA plots (figure 4) and UMD plots (figure S4). The Random Forest models explained moderate variability in biomass from FIA plot data (R² = 0.46–0.54, RMSE = 51.4–54.7 Mg ha⁻¹, RMSE CV = 34.1%–39.9%) with small estimated biases (MBE = 4.2–12.7 Mg ha⁻¹) (a). These values are slightly better than RMSE values reported by NBCDjks, ranged between 52.0 and 58.0 Mg ha⁻¹ (i.e. 5.2–5.8 kg m⁻²) in their bootstrap validation for five mapping zones (zone 60 to zone 64) that covers the study area. Unexpectedly, we noticed saturation effects in all ecoregion models. Validation using 20% with-held FIA plots indicated that four ecoregion models performed differently (b). Independent validation using UMD plots (figure S3) showed larger errors (R² = 0.33–0.61, RMSE = 85.3–100.9 Mg ha⁻¹, RMSE CV = 43.2%–85.0%) and different scales of under estimations (MBE = −8.9 ∼ −45.5 Mg ha⁻¹), but followed similar patterns as the results using FIA plots.

Similar patterns in biomass distribution were found when comparing the RF products from 30 and 90 m metrics to estimates from FIA plots (figure S5). Two of the four ecoregion models (Eastern Broadleaf and Northeastern Mixed) have smaller estimated biases in biomass predictions (MBE = 0.6–1.6 Mg ha⁻¹) using metrics generated at 30 m or 90 m scales. While the other two ecoregion biomass models had either underestimation (Outer Coaster, MBE = −17.4 Mg ha⁻¹) or overestimation (Central Appalachian, MBE = 11.6 Mg ha⁻¹) when using metrics at 30 m scales.

The random forest derived biomass map (RFjks) was compared to four national products (NBCD, Saatchi, Blackard, Wilson) and one global map product (GlobBiomass) over the Tri-State region. Although all map products followed general patterns in topography and landscape, fine-scale details varied especially over urban fringes (figure 5). Over forested area, less disagreements were observed among the three finer scaled map products (NBCDjks 30 m, Saatchi 90 m and GlobBiomass 100 m) and RFjks, compared to larger disagreements found among the two coarser map products (Wilson and Blackard, 250 m) and RFjks, as illustrated in the scatter plots and associated errors (figure 6). Biomass estimates had stronger correlations (R² > 0.65) among the three InSAR-optically fused national maps (NBCDjks, Saatchi and GlobBiomass) and the RFjks map over all regions. The two coarser resolution maps (Wilson and Blackard) had the poorest agreement with the RFjks.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The relationships at county- and state-scale varied between biomass estimates among five map products and those produced by the USFS FIA for the Tri-State region (table 4). All map-based estimates agreed well with the FIA estimates at county-scale (R² > 0.92, RMSD ranged between 3.0 and 9.0 Tg; figure S6). Details in county-scale of biomass density and stock were illustrated in figures S7 and S8. At the state-scale, the estimated aboveground biomass among different maps and USFS FIA showed large differences over forest and non-forest area. In Pennsylvania and Tri-State region (PA dominated), the four map-based estimates (RFjks, NBCDjks, Saatchi and GlobBiomass) that used either lidar or InSAR derived variables as model predictors, achieved closer values to each other as well as to the FIAjks over forest regions, but were more diverse over non-forest regions. In Maryland and Delaware, the estimates from this study (RFjks) matched best to FIA estimates over forest regions, followed by the estimates from the other three finer scaled maps (NBCDjks, GlobBiomass, and Saatchi). The estimated aboveground carbon density and stock among different maps and USFS FIA estimates showed large differences for each state and the Tri-State region. The map-based estimates in aboveground biomass ranged between 527.5 and 680.2 Tg across the Tri-State region. Whereas the USFS FIA estimates ranged between 638.4 (CRM) and 695 (Jenkins) Tg. Our estimates in aboveground carbon (563.8, 103.4 and 13 Tg) matched best to USFS for Pennsylvania, Maryland and Delaware (562.4, 116.7 and 15.9 Tg). The carbon estimates from the two coarser maps (Wilson & Blackard) were consistently ∼20% less than FIAcrm estimates in forested areas.