Use of MODIS NDVI Products to Map Tree Mortality Levels in Forests A ﬀ ected by Mountain Pine Beetle Outbreaks

: Extensive bark beetle outbreaks have recently occurred in western North American forests, resulting in overstory tree mortality across millions of hectares. Annual aerial surveys are currently used to operationally monitor bark beetle induced tree mortality, though this method is subjective and can exclude some forest areas. Daily Moderate Resolution Imaging Spectroradiometer (MODIS) satellite data o ﬀ er a potential alternative means to develop regional tree mortality maps. Accurate methods using such data could aid natural resource managers in surveys of forests with frequent overstory mortality, helping to prioritize forest treatment and restoration activities. This paper discusses a study to test the potential of using MODIS data to detect tree mortality. We developed and tested an approach to use 250-m resolution MODIS Normalized Di ﬀ erence Vegetation Index (NDVI) data products collected during a mountain pine beetle (MPB) outbreak and related tree mortality event in the northern Rocky Mountains of Colorado, USA. The 94 km 2 study area is predominantly lodgepole pine forest with most of the MPB-caused mortality occurring between 2003 and 2008. We used a 2.4-m forest conditions map from 2008 aerial multispectral imagery to calculate percentage of mortality within 240-m pixels for use as reference data. Using either daily or 16-day products, MODIS NDVI change products were calculated for 2008 versus either 2000 or 2003 baselines. MODIS change products were used as predictors in linear regression analysis to assess correlation between MODIS data and the aerial percent forest mortality map. Depending on the MODIS product, linear regression analyses yielded r 2 values ranging from 0.362 to 0.544 without outliers removed and from 0.406 to 0.570 with extreme outliers removed. Daily MODIS NDVI products from 2003 and 2008 were used with exponential regression to improve the r 2 to 0.593. The project showed some MODIS NDVI data potential for mapping percent tree mortality in forests subjected to regional bark beetle outbreaks and severe drought. power for this NDVI change map overall. The scatter plot shows that a linear relationship occurs between a given MODIS NDVI change map and aerial forest mortality map data. However, a range of NDVI change values can occur for a given percent mortality level from the aerial data. Based on the linear trend line, the percent mortality appears to peak at ~80%, though the maximum percent mortality on the aerial map is 85%.


Mapping US Forest Mortality with Aerial Insect and Disease Survey Data
For at least the last 20 years, the primary source of information on the location, severity, and extent of insect damage to forests in the US has been obtained through what is now known as the annual USDA Forest Service aerial Insect and Disease Survey (IDS) mapping products [9]. Previously known as Aerial Detection Survey (ADS) data, IDS products can include estimates of dead trees per acres for bark beetle-related forest mortality. Aerial surveys of disturbed forests have been used to depict areas with forest mortality caused by bark beetles [1,10] and used to assess bark beetle damage impacts on carbon stocks at broad regional scales [11,12]. IDS products have also been used as a reference data set in assessing other geospatial forest disturbance mapping products [13][14][15][16][17][18].
Because the information in the aerial surveys are subjectively gathered and recorded by human observers, these data have some uncertainty. Like electro-optical remote sensing data, IDS data collections can also be affected by impaired visibility from cloud cover, smoke, and haze [19]. Several studies assessed accuracy of IDS in estimating tree mortality, and found errors in both omission and commission, as well as locational errors compared to ground reference data [10,20,21]. The authors of [10] concluded that IDS forest mortality maps contain classification errors, though are still useful. When compared to fine-resolution airborne and satellite imagery, tree mortality from US aerial surveys was underestimated by factors of 3-20 [3]. A final issue that because of budget constraints, not all forested areas may be observed with an aerial IDS every year. While imperfect, IDS does provide basic annual information on the location, extent, severity, and causal agent for many forest disturbances in the United States (CONUS) that are not available otherwise [9].

Remote Sensing of Forest Mortality
Given the uncertainties and subjective nature of IDS regional maps of overstory tree mortality, remotely sensed imagery has the potential to complement such IDS products. Previous work has included efforts to apply high (i.e., fine), medium, and coarse spatial resolution multispectral remote sensing data sets to derive forest mortality maps. Fine and moderate spatial resolution aerial and satellite (e.g., QuickBird, Landsat) imagery with visible and near infrared multispectral bands has been used for detecting MPB forest mortality for various regional study areas [22][23][24][25][26][27][28]. While useful, such satellite data are not typically always available in cloud free form for broad regions on either an annual or seasonal basis at a time of year when insect-related forest disturbance can be easily detected. Acquisition cost can be another factor when it comes to procuring aerial and satellite image data, especially data acquired by commercial systems at high spatial resolutions (e.g., 5-m pixels or finer).
In contrast to sensors like QuickBird and Landsat, satellite imagery from the Moderate Resolution Imaging Spectroradiometer (MODIS) sensors offers advantages in terms of twice-daily data availability over broad areas (or processed to 8-day composites) and is available for free. MODIS data from two satellites, Terra and Aqua, together typically collect two images per day for most areas in the conterminous US at a 250-m spatial resolution for the red and near-infrared bands and provide time series of data on regional forest conditions from 2000 to present. This combination of spatial and temporal resolution has been used for detecting many regional forest disturbances [21,[29][30][31] but may be more limited in its application to detect smaller, localized forest disturbances [20,32].
Different spectral vegetation indices from daily MODIS reflectance data have been used to detect regional forest defoliation from some insects [20,21,29,31]. The authors of [33] employed various spectral indices from 16-day composited MODIS MOD13 data to map forest mortality within pine plantations in Australia damaged by Ips bark beetles. In doing so, 231.7-m Normalized Difference Vegetation Index (NDVI), Normalized Difference Infrared Index (NDII), and Enhanced Vegetation Index (EVI) change products had relatively low r 2 values with a 15-cm resolution aerial image-based map of tree mortality that was used as reference data. The use of the MOD13 NDVI data resulted in a r 2 of 0.39 [33], although relatively few observations were used (n = 62). The authors of [33] also tested the influence of mixed pixels, which did not appear to affect r 2 unless the pixels were less than 25 percent forested. Such mixed pixels were mitigated using a forest mask technique [33].
Despite its high potential, there have been only a few initial studies on the capability of using MODIS for monitoring regional forest mortality from MPB outbreaks. The authors of [26] generated a thematic map of forest mortality from a single date of MODIS multispectral data that yielded an overall accuracy of 78% compared to an accurate, field-checked 2.4-m forest mortality map derived from aerial imagery. A separate study [15] later generated a 2008 versus 2003 MODIS NDVI difference map depicting MPB mortality in Colorado that demonstrated visual agreement with IDS maps of MPB mortality for the same area.
Thus, there exists the need for more studies to assess additional methods using MODIS imagery to map MPB-caused overstory tree mortality. Consequently, we investigated the potential utility of MODIS for this application in two ways. First, given the coarse spatial resolution, we assessed the application of MODIS NDVI data in estimating percent mortality on a per pixel basis (i.e., depicting damage severity), not just classifying pixels as predominantly having mortality present or absent. Secondly, we assessed the potential of using time series of daily satellite remote sensing data (as opposed to 16-day temporally composted data) to supplement or complement IDS-based forest mortality mapping products.

Study Area
The study area is in a portion of the Arapaho National Forest just west of the Continental Divide in the U.S. north-central Colorado Rocky Mountains that experienced a major outbreak of MPB causing extensive forest mortality ( Figure 2). Located in Grand County, the study area bio-physical characteristics discussed herein are based on [34]. With elevations ranging from 2500 to 3400 m, this area includes extensive subalpine forest dominated by lodgepole pine (Pinus contorta Dougl. ex Loud.). Spruce-fir forest stands composed of Engelmann spruce (Picea engelmannii Parry ex Engelm) and subalpine fir (Abies lasiocarpa (Hook.) Nutt.) trees occur along riparian zones and at higher elevations. Herbaceous meadows occur in riparian areas, along with occasional water bodies. Some trembling aspen (Populus tremuloides Michx.) trees occur at intermediate elevations. The authors of [34] also report that the study area has an average annual precipitation of approximately 52 cm, and with mean minimum and maximum temperatures of −10.2 and 9.8 • C, respectively, based on Fraser, CO weather station data from [35].
The MPB outbreak in our study area began in the mid-1990s and accelerated after a multiyear drought occurred in the early years of the 21st century. This outbreak in conjunction with the drought led to much of the mature lodgepole pine west of the continental divide in northern Colorado being killed by August 2008, with dead red forest crowns frequent across the landscape and dead gray crowns also occurring [36]. The MPB-caused tree mortality was related in part to drought and higher temperatures [37][38][39] with the predominance of mature lodgepole pine forest as an additional important factor [39]. Drought in the region during 2001-2003 was noted as being severe by [39] based on Palmer Drought Severity Index data from NOAA. The 2002 drought for many areas in Colorado was particularly severe according to [40]. subalpine fir (Abies lasiocarpa (Hook.) Nutt.) trees occur along riparian zones and at higher elevations. Herbaceous meadows occur in riparian areas, along with occasional water bodies. Some trembling aspen (Populus tremuloides Michx.) trees occur at intermediate elevations. The authors of [34] also report that the study area has an average annual precipitation of approximately 52 cm, and with mean minimum and maximum temperatures of −10.2 and 9.8 • C, respectively, based on Fraser, CO weather station data from [35].

MODIS Data Preparation
For this study, we used Collection 5 MODIS satellite data products acquired from the USGS Land Processes Distributed Active Archive Center (LP DAAC), including 250-m 16-day NDVI MOD13 and MYD13 products, plus 250-m daily surface reflectance data from MOD09 and MYD09 products ( Table 1). The acquired MODIS data sets were temporally processed for the mid to late summer period when MPB-caused mortality was highly distinct [19,41] and when cloud cover was less frequent [42]. This time of year was also similar to when our reference data were collected for the study area in the form of aerial imagery and field measurements. The MOD13 and MOD09 products used in the study were from the Terra satellite, while the MYD13 and MYD09 products were from the Aqua satellite. The 16-day MOD13/MYD13 NDVI products [43] were based on daily MOD09/MYD09 atmospherically corrected reflectance data [44]. The MOD13/MYD13 NDVI products were originally derived using a constrained view angle/NDVI maximum value compositing technique described by [43].  [21,[45][46][47] and the Phenological Parameters Estimation Tool (PPET) [48,49]. The authors of [29] provide a description of CONUS vegetation phenology products output by PPET. For this project, TSPT was used to perform most of the temporal processing of the annual NDVI time series data from MOD13 and MYD13 data, and PPET was used to further refine the fused MOD13/MYD13 NDVI data after TSPT processing (Figure 3).  The same workflow also generally applies to the NDVI change products derived directly from daily MOD09/MYD09 products, except that the MOD09/MYD09 change products used daily NDVI data as inputs to temporal data processing.
TSPT was used initially to read the MOD13 and MYD13 NDVI values and quality assurance (QA) data, including image acquisition dates for each pixel, which were then disaggregated from 16-day composites into daily data. MOD13 and MYD13 QA data were also used to remove pixels affected by clouds and cloud shadows. NDVI noise artifacts were removed using a spike detection and removal algorithm, and data with extremely oblique viewing geometry (e.g., data acquired with view zenith angles exceeding 50 degrees) were removed. Once cleaned, the MODIS Aqua and Terra data was merged via maximum value compositing of NDVI for a given date of daily data. Data voids in the NDVI time series were then filled using an interpolation technique that included temporally filtering via curve fitting with the Savitzky Golay (SG) algorithm [50]. After interpolation, NDVI was reaggregated from daily into 8-day NDVI composites. In processing the MOD13/MYD13 data, the . Workflow used to generate MODIS percent NDVI change maps from MOD13/MYD13 NDVI data. MOD13 is the 16-day composite data from the Terra satellite; MYD13 is from the Aqua satellite. The same workflow also generally applies to the NDVI change products derived directly from daily MOD09/MYD09 products, except that the MOD09/MYD09 change products used daily NDVI data as inputs to temporal data processing.
TSPT was used initially to read the MOD13 and MYD13 NDVI values and quality assurance (QA) data, including image acquisition dates for each pixel, which were then disaggregated from 16-day composites into daily data. MOD13 and MYD13 QA data were also used to remove pixels affected by clouds and cloud shadows. NDVI noise artifacts were removed using a spike detection and removal algorithm, and data with extremely oblique viewing geometry (e.g., data acquired with view zenith angles exceeding 50 degrees) were removed. Once cleaned, the MODIS Aqua and Terra data was merged via maximum value compositing of NDVI for a given date of daily data. Data voids in the NDVI time series were then filled using an interpolation technique that included temporally filtering via curve fitting with the Savitzky Golay (SG) algorithm [50]. After interpolation, NDVI was reaggregated from daily into 8-day NDVI composites. In processing the MOD13/MYD13 data, the TSPT upward and downward spike detection threshold was set to 0.05, the SG frame size was assigned to 23 days and the SG full width half maximum (FWHM) was set to 16 days. In addition, all negative NDVIs were assigned to 0 and the positive NDVIs multiplied by 100.
The resulting NDVI products from TSPT processing were then used in PPET to generate 8-day cumulative integral (CI) NDVI products, yielding 46 total products per year. In doing so, the CI NDVI products were subjected to segmented curve fitting to provide additional smoothing. Individual dates of 8-day NDVI were calculated from the CI NDVI products by subtraction of the previous date of CI NDVI from a current date of CI NDVI [51]. After PPET processing, the 8-day NDVIs were reaggregated into 24-day NDVIs (via maximum value compositing) with each 24-day product date interleaved by 8-days with the preceding and subsequent date [30]. This 24-day temporal composite aggregation is employed for near-real time regional forest disturbance monitoring by the USDA Forest Service ForWarn system in which cloud free NDVIs are needed for regional forest monitoring across the CONUS. For this study, the 24-day interval of 5-28 August was chosen for analyzing NDVI change and MPB forest mortality within the study area.

MOD09/MYD09 Data
We also evaluated products derived from daily MODIS imagery to test if such data could produce a more refined, accurate NDVI change product for estimating percent tree mortality. Daily MOD09 Terra and MYD09 Aqua reflectance data were temporally processed using TSPT software with parameter settings for the noise reduction and void interpolation steps optimized for these data. In addition, daily NDVI was computed from the MOD09/MYD09 reflectance data. TSPT settings included the upward and downward spike detection threshold set to 0.05 maximum slope between dates, the SG frame size set to 23 days, and the SG full width half maximum set to 16 days. Data were only processed when view zenith angles were ≤ 50 degrees. For TSPT to effectively process the daily reflectance data into daily NDVI products, the data collection period included the date of comparison, plus approximately one month of data was used from before and after the targeted date of observation. These extra days of data were needed to interpolate NDVIs for dates with data voids in the time series.

Percent NDVI Change Maps
After preliminary processing of MODIS NDVI data, percent NDVI change products were calculated with the following formula: Percent NDVI change = ((((NDVI current − NDVI historical)/NDVI historical) × 100) + 101), where current is the current year (2008) and historical is either 2000 or 2003. This change metric is comparable to that used in the derivation of forest change product resident to the ForWarn system [30,31]. MODIS NDVI change maps generated for the project included:  Figure 4). Such change products have values that range from −100 to +100%. For this project, we added 101 to rescale the output from 1 to 201. Rarely occurring extreme negative change outliers (those less than −100% change) were set to 1 in the final output and rarely occurring extreme positive change outliers (those greater than +100% change) were set to 201 in the final output. The no-change value was set to 101. This rescaling procedure was performed so that the change products could be output as unsigned 8-bit data files.

Reference Data Preparation
The primary reference data for the project was a 2.4-m resolution forest cover and condition map derived by classification of aerial multispectral imagery acquired on 13 August 2008 ( Figure 5). The multispectral data used in the classification was originally acquired at a spatial resolution of 30 cm and was aggregated to 2.4-m resolution prior to the classification [34]. This map included three forest condition classes (live forest (with green overstory canopies), dead forests (with either dead red and dead gray trees) in addition to herbaceous vegetation (e.g., meadows), bare soil, and shadow/water classes. This map was assessed by [34] using field measurements from 400-m 2 circular

Reference Data Preparation
The primary reference data for the project was a 2.4-m resolution forest cover and condition map derived by classification of aerial multispectral imagery acquired on 13 August 2008 ( Figure 5). The multispectral data used in the classification was originally acquired at a spatial resolution of 30 cm and was aggregated to 2.4-m resolution prior to the classification [34]. This map included three forest condition classes (live forest (with green overstory canopies), dead forests (with either dead red and dead gray trees) in addition to herbaceous vegetation (e.g., meadows), bare soil, and shadow/water classes. This map was assessed by [34] using field measurements from 400-m 2 circular   To compare to MODIS percent NDVI change products, the aerial forest condition map was further processed into a series of fractional abundance maps, one for each mapped class. In each case, a given forest condition class was recoded to 100 and then spatially aggregated into a 240-m fractional abundance product that reported the percent of the pixel occupied by the featured condition (e.g., dead red forest). To derive the aerial percent forest mortality map, fractional abundance images for the dead red and dead gray tree categories were added together into a composite percent mortality map (Figure 5d). Non-mapped pixels outside the effective map area were assigned a null value and were disregarded for subsequent processing. The output was then reprojected to match the MODIS spatial resolution and projection.

Regression Analysis
The project assessed relationships between a "reference" aerial percent forest mortality map and various MODIS-based NDVI change maps. The aerial imagery-based percent forest mortality map was used as the dependent variable "y" for multiple regression analyses that used a given NDVI change map as the independent variable "x". The results of these regression analyses are summarized in Table 2   The effect of extreme spurious outliers on regression results was explored by removing data points with ± three or more standard deviations from a the mean value of a given tested MODIS percent NDVI change product. After outlier removal, the regression analyses were rerun and summarized as Table 3. For the MODIS product with the highest estimated r and r 2 result from linear regression analysis, an r 2 was also computed using a nonlinear exponential equation in order to compare to the linear regression results. For the best linear regression result, percent tree mortality map were computed from the applicable MODIS NDVI change map using linear and nonlinear regression equations. In addition to the scatter plot and r 2 analysis, the MODIS-based forest mortality maps (from the linear and nonlinear models) were compared visually to the aerial-based percent mortality map.

MODIS NDVI Change versus Percent Tree Mortality Map Using All Pixels
The daily reflectance-based NDVI change map for 2008 with the 2003 baseline was compared to the reference aerial image-based percent tree mortality map, yielding an r 2 of 0.544. In contrast, the 16-day MOD13/MYD13 NDVI change product using the 2003 baseline yielded a lower explained variance (r 2 = 0.362). Using 2000 as the baseline instead of 2003 decreased the r 2 for MOD09/MYD09 data (r 2 = 0.425), but the explained variance for MOD13/MYD13-based data compared to reference data increased slightly using the 2000 baseline (r 2 = 0.393). The MOD09/MYD09 data for each of the two baseline years showed higher r 2 values than the best result from MOD13/MYD13 data. The MOD09/MYD09 NDVI change product with the 2003 baseline yielded the lowest standard error (14.577), though the amount was still higher than desired. In contrast, the MOD13/MYD13 change product with the 2003 baseline yielded the highest standard error (17.249).

MODIS NDVI Change versus Percent Tree Mortality Map with Outliers Removal
After outlier removal, processed daily MOD09/MYD09 NDVI data using the 2003 baseline and the reference data yielded an r 2 of 0.570; with the 2000 baseline, r 2 = 0.479 (Table 3). Upon outlier removal, processed MOD13/MYD13 NDVI data with outlier removal produced higher r 2 values with the reference data using the 2003 (r 2 = 0.523) and 2000 baseline (r 2 = 0.406). After outlier removal, the MOD13/MYD13 NDVI change product using the 2003 baseline yielded a higher r 2 value than the MOD09/MYD09 NDVI change product using a 2000 baseline (Table 3). With outlier removal applied, the best temporally processed MOD09/MYD09 data produced better r 2 results (more of the variance explained) with the reference aerial forest mortality map (r 2 = 0.570) compared to that with the best temporally processed MOD13/MYD13 NDVI data with the reference data (r 2 = 0.523). For the NDVI change products with outlier removal, the standard error was somewhat reduced, compared to the same change maps with extreme outlier removal (Tables 2 and 3). With outlier removal applied, the standard error was more similar for the MOD09/MYD09 and MOD13/MYD13 NDVI change maps using the 2003 baseline (14.150 versus 14.888). Although the standard error was reduced for all the regressions with outlier removal applied, there was still a higher than desired standard of error for all tested products with outlier removal. Figure 6 shows the scatter plot of the NDVI change product and the aerial image mortality for the linear regression model with the best r 2 (MOD09/MYD09 data with 2003 baseline). This figure visually depicts the relationship between the MODIS NDVI change map and aerial tree mortality map, including linear and nonlinear (exponential) trend lines overlain onto the scatter plot. Based on the r 2 value, the linear model indicates moderate statistical explanatory power for this NDVI change map overall. The scatter plot shows that a linear relationship occurs between a given MODIS NDVI change map and aerial forest mortality map data. However, a range of NDVI change values can occur for a given percent mortality level from the aerial data. Based on the linear trend line, the percent mortality appears to peak at~80%, though the maximum percent mortality on the aerial map is 85%.
This suggests that the linear model may under predict very high percent mortality levels. The linear model also suggests that the MODIS-based NDVI change of~+3% (on x-axis) is roughly equivalent to 0% mortality on the aerial data (on y-axis). However, the 0% mortality level on the aerial map is more frequently represented by NDVI change values that are negative. In some cases, NDVI change values may be influenced by sub-lethal MPB impacts to attacked but still living trees [30,31]. change map and aerial forest mortality map data. However, a range of NDVI change values can occur for a given percent mortality level from the aerial data. Based on the linear trend line, the percent mortality appears to peak at~80%, though the maximum percent mortality on the aerial map is 85%. This suggests that the linear model may under predict very high percent mortality levels. The linear model also suggests that the MODIS-based NDVI change of~+3% (on x-axis) is roughly equivalent to 0% mortality on the aerial data (on y-axis). However, the 0% mortality level on the aerial map is more frequently represented by NDVI change values that are negative. In some cases, NDVI change values may be influenced by sub-lethal MPB impacts to attacked but still living trees [30,31]. The nonlinear model produced a higher r 2 (0.593) than the linear model (0.570) ( Figure 6). The nonlinear model appears to reduce NDVI change predicted percent mortality values in the mid-range of the percent mortality values represented on the aerial-based forest mortality map. In addition, the nonlinear model may over-estimate the highest percent mortality values, given that the aerial data's maximum percent mortality of 85%. The nonlinear model's trend line does appear to more effectively follow data points in the highest and the lowest percent mortality levels. However, fewer data points are present in the highest percent mortality levels. Figure 7 compares percent tree mortality maps derived from aerial ( Figure 7a) and MODIS data (Figure 7c,d). The MODIS tree mortality maps were computed using the MOD09/MYD09 NDVI change map for 2008 versus 2003 in conjunction with regression equations (linear and exponential) that used the relationship between the rescaled NDVI change map and the aerial mortality map (equations given in Figure 6) to derive predictive MODIS-based mortality maps, applying both linear and exponential The nonlinear model produced a higher r 2 (0.593) than the linear model (0.570) ( Figure 6). The nonlinear model appears to reduce NDVI change predicted percent mortality values in the mid-range of the percent mortality values represented on the aerial-based forest mortality map. In addition, the nonlinear model may over-estimate the highest percent mortality values, given that the aerial data's maximum percent mortality of 85%. The nonlinear model's trend line does appear to more effectively follow data points in the highest and the lowest percent mortality levels. However, fewer data points are present in the highest percent mortality levels. Figure 7 compares percent tree mortality maps derived from aerial ( Figure 7a) and MODIS data (Figure 7c,d). The MODIS tree mortality maps were computed using the MOD09/MYD09 NDVI change map for 2008 versus 2003 in conjunction with regression equations (linear and exponential) that used the relationship between the rescaled NDVI change map and the aerial mortality map (equations given in Figure 6) to derive predictive MODIS-based mortality maps, applying both linear and exponential models. A visual comparison of the MODIS-and aerial-based mortality maps showed many areas with good correspondence in general. In addition, the linear model (Figure 7c) appears to more closely approximate the patterns evident in the aerial forest mortality map (Figure 7a). However, the  (Figure 7d) shows lower mortality levels more comparable to those on the aerial mortality map (Figure 7a) that correspond with areas that have a mixture of live green forests and sometimes non-forest green vegetation on the aerial false color RGB (Figure 7b). Consequently, the exponential model (Figure 7d) may depict the lower mortality level areas more effectively than the linear model (Figure 7c). The scatter plot shows a high amount of variability, though the trend lines appear to dissect the data cloud in an expected manner.

Discussion
Our assessment of the effect of baseline year indicated that the MODIS NDVI change products using the 2003 baseline yielded higher r 2 values than those using 2000. Differences could be related to the reduced data availability of MODIS data in 2000. Only one MODIS instrument was in operation as of 2000 (Terra), and frequent cloud cover occurred in the study area during August 2000, based on QA

Discussion
Our assessment of the effect of baseline year indicated that the MODIS NDVI change products using the 2003 baseline yielded higher r 2 values than those using 2000. Differences could be related to the reduced data availability of MODIS data in 2000. Only one MODIS instrument was in operation as of 2000 (Terra), and frequent cloud cover occurred in the study area during August 2000, based on QA data and visual inspection of reflectance data. The full collection of Aqua data began in 2003, which meant that twice as much daily MODIS data were available for producing the 2003 baseline data set compared to 2000, thereby allowing for more data with greater cloud free coverage and increased frequency of quality NDVIs. Our regression analyses also showed that the temporally processed daily MOD09/MYD09 data produced slightly better r 2 values than the 16-day MOD13/MYD13 products (Tables 2 and 3). The results suggest that the temporal processing of MODIS time series data may be improved by processing daily instead of pre-composited 16-day NDVI data.
The MODIS-based mortality maps shown in Figure 7 were also compared visually to understand differences with aerial mortality data products. For the linear model, the differences were most apparent for the lowest and highest percent tree mortality levels. Some differences between the aerial and MODIS products appeared to occur in areas with small or narrow patches of non-forest and live non-attacked forest that were intermingled with or were adjacent to areas of forest mortality as seen on the 2.4-m aerial imagery. Similarly, edges between larger patches of differing land cover types and conditions also tended to be more spatially mixed on MODIS 250-m data. Differences between the test and reference maps could be due to several reasons: (1) underestimation of mortality for lower NDVI change levels; (2) overestimation of mortality at higher levels of negative NDVI change; (3) mixed pixel effects (forest and non-forest classes); (4) uncertainties in co-registration of the aerial and MODIS images; and resampling of the MODIS and aerial data to a common map projection; (5) the coarser 250-m spatial resolution of the MODIS data relative to the aerial data; and (6) the relative low frequency of low percent mortality areas in the study area. Although infrequent, some residual spurious data values (e.g., from unmitigated atmospheric contamination) in the MODIS change products may occasionally occur. Another possible factor may be classification error in the spatially averaged 240-m aerial percent mortality classification product, though this effect is likely to be a very minor or negligible issue given that the aerial classification of six classes at 2.4-m resolution produced an overall accuracy of 90% compared to ground reference data [34]. The integration of the MODIS and aerial data was another potential factor, given that the two data sources had different spatial resolutions and map projections.
Our study found that the inclusion of infrequent extreme outliers in the regression analyses lowered the r 2 value, especially with respect to the MOD13/MYD13-based forest change products ( Table 2). By pruning the most extreme outliers to exclude values exceeding mean MODIS NDVI change ±3 standard deviations (SDs), the regression results improved, especially for the MOD13/MYD13-based product using the 2003 baseline (Table 3). This refinement provided a means for assessing performance that considers most of the data, only excluding the most extreme yet rare outlier data considered to be "noise". It is possible that the r 2 values may be further improved by using a more stringent threshold for pruning outliers, such as the mean NDVI change plus a 2 SD instead threshold instead.
In terms of r value (as opposed to r 2 ), all the tested MODIS products (without and with outlier removal) exceeded values of −0.60 (Tables 2 and 3), indicating moderate strength between the compared variables [52]. The best MOD09/MYD09 and MOD13/MYD13 results obtained with outlier removal yielded r correlation coefficients greater than −0.7, which [52] considers to be of high strength.
In this study, the best linear regression for using MODIS NDVI data to estimate forest mortality resulted in an r 2 of 0.57, which exceeded the best r 2 of 0.39 reported by [33], which used MODIS 16-day MOD13 NDVI to estimate forest mortality in pine plantations due to Ips bark beetle attacks. Our project's higher r 2 values may be due to additional temporal processing of the MODIS data, the inclusion of MODIS Aqua data along with MODIS Terra data, differences in landscape characteristics (e.g., patch sizes, cover types, and tree mortality levels), plus our study used a much higher number of observations used to compute the correlations. In addition, our study aimed to estimate percent forest mortality as the fraction of forest canopy cover that exhibited forest mortality. In contrast, [33] defined forest mortality in terms of number of dead trees per pixel. It may be more difficult to map the number of dead crowns, given that crown size can vary according to tree species, height, age, and growing conditions (i.e., site).
Environmental factors may also influence the correlation between MODIS and aerial imagery products. The level of cloud cover can significantly influence the number of quality data points available. We used Collection 5 MODIS data for the project, given it was the level of MODIS data that was available at the time the study was conducted. It is possible that using Collection 6 MODIS data may result in improved percent forest mortality maps, given that calibration coefficients, product generation methods, and data quality assessment algorithms have been updated. Correlations may also be improved by including other meaningful explanatory environmental variables in a multiple regression analysis, such as digital elevation data and detailed forest cover type maps from, e.g., the Gap Analysis Project [53] or Landscape Fire and Resource Management Planning Tools Project (LANDFIRE) [54]. With bark-beetle induced forest mortality, the contrast between live and dead forest canopies (regardless of being dead red or gray) is distinct to the visible eye. The comparisons of MODIS-based percent forest mortality maps to reference data may also be aided using IDS data to identify general locations of recorded forest mortality areas, taking into account that the IDS data can have classification and locational errors.
The results presented here pertain to a mature predominately lodgepole pine forest attacked by MPBs in which most of the overstory trees were either in a dead red or gray state as of 2008. MODIS 250-m NDVI change products from both prior to and after MPB attack provide information on the relative severity or intensity of the damage, as a function of forest canopy greenness decline. MODIS NDVI change imagery typically shows relative change in canopy greenness due to extensive, severe forest disturbances, such as MPB-induced forest mortality, especially when the damage is sufficiently extensive and severe within a given pixel. However, the ability of MODIS NDVI forest change maps to be converted into effective percent forest mortality maps may also depend in part on the mixture of land cover types within the MODIS pixel, as well as the amount of reduction in greenness occurring in MPB-killed forests.
We used 250-m MODIS NDVI data, in part because it uses finer spatial resolution spectral band data compared with the 500-m resolution MODIS spectral bands. It is conceivable that other spectral indices could improve MODIS-based estimations of MPB mortality, especially if employing additional temporal processing techniques (e.g., perhaps using a Bidirectional Reflectance Distribution Function (BRDF) as a step). However, some of the MODIS 500-m resolution bands can have data quality issues, such as data dropouts in MODIS Aqua band 6 (1.6-micron short wave infrared [55] and misregistration of short-wave infrared to visible and near infrared bands [56]. Furthermore, MODIS 250-m products have four times finer spatial resolution than the 500-m data, and this finer 250-m resolution may be an important factor for discriminating changes in canopy greenness in smaller patches than observed at the 500-m scale. The best MODIS-based percent mortality map from our project yielded a r 2 of 0.59. However, the standard error for the tested products was high (best case for linear regression was 14.150) and we acknowledge there is room for improvement. We think the performance of the percent forest mortality models may improve with additional research into the following: (1) alternative methods for processing the MODIS data, including methods for noise reduction and co-registration of NDVI and reference data on forest mortality; (2) use of MODIS Collection 6 data which may improve the quality of acquired MODIS data products; (3) use of other explanatory environmental variables in the model in addition to percent change in NDVI, such as digital elevation models and perhaps meteorological data; (4) use of alternative regression methods, such as logistic regression; (5) incorporation of detailed vegetative cover type geospatial data; and (6) use of percent tree canopy geospatial data in refining maps of forest mortality, which in particular may help improve the development of more viable percent forest mortality maps in semi-open forests and woodlands, such as pinyon juniper woodlands. In addition, other relevant research needs to be conducted, including: (1) use of alternatives to aging MODIS systems (e.g., OLCI and VIIRS) for producing percent forest mortality maps at the regional scale; and (2) use of alternative spectral indices to NDVI that may yield improved maps of percent forest mortality.
MODIS percent NDVI change products are currently used in the ForWarn near real time on-line system for viewing, detecting, assessing, and tracking regional forest disturbances across CONUS [30,31]. ForWarn presently uses MODIS Collection 6 near real time NDVI data from the Global Inventory Monitoring and Modeling Studies (GIMMS) Global Agricultural Monitoring (GLAM) system, along with historical baseline NDVIs from temporally processed MOD/MYD13 data, to compute forest disturbance monitoring products that are displayed in the system's geospatial data viewer. The methods we discuss here are amenable to automation within ForWarn in that regression equations, once constructed, could be used to generate MODIS NDVI-based forest mortality maps in near real time, especially if combined with other relevant geospatial data and GIS-based decision rules. The use of MODIS NDVI data for estimating percent forest mortality is a value-added product that can be used to aid the ForWarn system as well as part of an overall effort to synoptically monitor forest conditions with a multi-scale, multi-technology approach involving ground, aerial, and satellite observations. Such a multi-stage sampling approach to forest health monitoring and assessment was discussed by [9].

Conclusions
In this study, we performed regression analysis to compute and assess r and r 2 values for different MODIS 250-m NDVI change products compared to a previously produced, highly accurate percent forest mortality map derived from high spatial resolution aerial multispectral imagery. For the best apparent change map, we applied regression equations from the earlier analysis on the MODIS NDVI change map to compute percent mortality maps. The study resulted in a method using MODIS imagery that could be employed to help scale up forest mortality mapping from a local to regional scale. Our study demonstrated that moderate to high strength correlations (r values) occurred for various MODIS NDVI change products compared to a percent mortality map from aerial reference data. We found that the best r 2 values were obtained with temporally processed daily NDVI computed from MODIS surface reflectance data, though similar results were obtained from temporally processed 16-day NDVI products. The highest r 2 values occurred when using the daily MODIS MOD09/MYD09-based NDVI products. The project also used the best apparent MODIS-based data in conjunction with aerial data to compute mortality maps that compare favorably with a reference mortality map derived from only aerial data. Based on results from this study, our technique should work comparably for other lodgepole pine forests in the region similarly affected by a multiyear MPB outbreak. By 2008, the overstory lodgepole pine forests in the study area were for the most part severely impacted by MPB. While our results are promising, additional studies need to be conducted on other areas with bark beetle damaged forest, including forests less damaged by MPB, to further assess the potential of MODIS data for this application. Other sources of satellite data should also be considered, given that the MODIS sensors are nearing the end of their life span.
The development and demonstration of this MODIS technique could not have been performed without the availability of an accurate high-spatial resolution forest mortality map for a sufficiently sizeable area (i.e., small region). However, now that regression models have been built, the best apparent models (linear and nonlinear) produced from this project could be used to compute MODIS-based mortality maps for other lodgepole pine forested areas. Based on project results, the nonlinear model in conjunction with a MODIS NDVI change map may be able to detect lower levels of mortality than the linear model, which could be useful in initially mapping early stages of MPB outbreaks. For other types of forests with bark beetle outbreaks, it may be necessary to build new regression models for using MODIS data to produce forest mortality maps. Where high spatial resolution aerial imagery-based forest mortality maps do not typically exist as reference data, it may be possible to construct alternative reference forest mortality maps from such sources as National Agriculture Imagery Program (NAIP), Landsat and Sentinel-2 satellite, higher resolution commercial satellite, aerial IDS, and/or in-situ data.
The MODIS NDVI change products used in this study were generated to aid ForWarn and other near real time forest disturbance detection and monitoring systems. This study provided evidence that useful forest mortality mapping products can be generated from temporally processed MODIS 250-m NDVI data from the ForWarn system. Such mortality maps could be integrated into the ForWarn system to help end-users (e.g., resource managers) in assessing regionally extensive forest mortality due to multiyear bark beetle outbreaks. MODIS-based percent forest mortality products can also potentially aid forest salvage and restoration efforts, fire fuel mapping, and forest mortality risk mapping (e.g., during drought years). The project also provided evidence that use of daily MOD09/MYD09 data may improve generation and utility of near real time NDVI change products resident to the US Forest Service's ForWarn system. A similar finding was made by [21] regarding the mapping of gypsy moth defoliation with MODIS data.
While the mortality mapping technique from this project showed potential for MPB impacted lodgepole pine dominated forests, it does not determine or predict causal agents and does not summarize results according to mapped cover type. More effort is needed to develop these additional capabilities. Additional work is also required to further assess the application with respect to other bark beetle damage and other forest disturbances, cover types, and locations. More studies are needed to assess potential of MODIS for estimating tree mortality in other forest and woodland types in the western US, such as high elevation spruce-fir forests and pinyon-juniper woodlands. Additional research is also needed to further improve and refine models of percent forest mortality from MODIS and comparable satellite systems. If this MODIS-based method proves to be effective for other forest types, bark beetle species, mortality levels, and locations, then such a validated procedure could be applied for forest mortality mapping at regional and national scales.

Funding:
The project's geospatial data processing and analysis was conducted with funding from the USDA Forest Service Threat Assessment Centers (EFETAC and WWETAC) to NASA Stennis Space Center, as part of a multi-federal agency collaborative multiyear project to develop and deploy national and regional forest threat early warning system tools. In doing so, personnel from the Computer Sciences Corporation (CSC), Inc. and ASRC Research and Technology Solutions (ARTS) performed project data processing and analysis under contract NNS10AA35C. This manuscript was completed with funding from a USDA Forest Service grant to Leidos (contract HHSN316201200034W, award AG4568D170050).