The accuracy of large-area forest canopy cover estimation using Landsat in boreal region

International Journal of Applied Earth Observation and Geoinformation 53 (2016) 118–127

Contents lists available at ScienceDirect

International Journal of Applied Earth Observation and Geoinformation journal homepage: www.elsevier.com/locate/jag

The accuracy of large-area forest canopy cover estimation using Landsat in boreal region Hadi a,∗ , Lauri Korhonen b , Aarne Hovi a , Petri Rönnholm a , Miina Rautiainen a,c a b c

Aalto University, School of Engineering, Department of Built Environment, PO Box 15800, 00076 Aalto, Finland University of Eastern Finland, School of Forest Sciences, PO Box 111, FI-80101 Joensuu, Finland Aalto University, School of Electrical Engineering, Department of Radio Science and Engineering, PO Box 13000, 00076 Aalto, Finland

a r t i c l e

i n f o

Article history: Received 23 May 2016 Received in revised form 9 August 2016 Accepted 31 August 2016 Keywords: Canopy cover Landsat Tree cover Boreal Forest

a b s t r a c t Large area prediction of continuous field of tree cover i.e., canopy cover (CC) using Earth observation data is of high interest in practical forestry, ecology, and climate change mitigation activities. We report the accuracy of using Landsat images for CC prediction in boreal forests validated with field reference plots (N = 250) covering large variation in latitude, forest structure, species composition, and site type. We tested two statistical models suitable for estimating CC: the beta regression (BetaReg) and random forest (RanFor). Landsat-based predictors utilized include individual bands, spectral vegetation indices (SVI), and Tasseled cap (Tass) features. Additionally, we tested an alternative model based on spectral mixture analysis (SMA). Finally, we carried out a first validation in boreal forests of the recently published Landsat Tree Cover Continuous (TCC) global product. Results showed simple BetaReg with red band reflectance provided the highest prediction accuracy (leave-site-out RMSECV 13.7%; R2 CV 0.59; biasCV 0.5%). Spectral transformations into SVI and Tass did not improve accuracy. Including additional predictors did not significantly improve accuracy either. Nonlinear model RanFor did not outperform BetaReg. The alternative SMA model did not outperform the empirical models. However, empirical models cannot resolve the underestimation of high cover and overestimation of low cover. SMA prediction errors appeared less dependent on forest structure, while there seemed to be a potential for improvement by accounting for endmember variability of different tree species. Finally, using temporally concurrent observations, we showed the reasonably good accuracy of Landsat TCC product in boreal forests (RMSE 13.0%; R2 0.53; bias −2.1%), however with a tendency to underestimate high cover. © 2016 Elsevier B.V. All rights reserved.

1. Introduction Regular large-scale monitoring of forest cover with remote sensing is presently needed more than ever before; the world’s political leaders collectively acknowledge the urgent importance of stopping deforestation and forest degradation in the international climate agreement reached in Paris in 2015. The boreal forests which account for 21% of world’s forest area, while not experiencing serious deforestation threat, are highly vulnerable to global warming which increases disturbance e.g. from fire or insect outbreaks (Bonan, 2008). These disturbances potentially contribute CO2 emissions and further enhance climate warming.

∗ Corresponding author. E-mail address: hadi.hadi@aalto.fi ( Hadi). http://dx.doi.org/10.1016/j.jag.2016.08.009 0303-2434/© 2016 Elsevier B.V. All rights reserved.

Forest cover mapping with remote sensing (RS) has been traditionally done by classification of forest and non–forest areas (or other pre-defined forest classes), which are then upscaled to a larger mapping unit to calculate the proportion of pixels or segments labelled as forest. However, the resulting maps are dependent on the forest definition, i.e. mainly the threshold of tree cover above which an area is classified as forest. The definitional discrepancies have been the main source of inconsistent satellite-based estimates between previous studies (Sexton et al., 2015). This has led to the recommendation of shifting from classification approach toward mapping canopy cover as a continuous field (Hansen et al., 2002; Sexton et al., 2015). In addition to allowing consistent forest area estimates, canopy cover maps would allow for monitoring changes in forest cover at scales smaller than the pixel size. Canopy cover is highly correlated with basal area and thus volume and biomass. Beyond practical forestry context, canopy cover is also an important ecological indicator related to

Hadi et al. / International Journal of Applied Earth Observation and Geoinformation 53 (2016) 118–127

119

Table 1 Summary of study plots. Site

No. of plots

Plot size (m)

Canopy cover (%)

Stand density (trees ha−1 )

Mean tree height (m)

Mean DBHa (cm)

Basal area (m2 ha−1 )

Suonenjoki Koli Tammela Joensuu Rovaniemi Sodankylä Hyytiälä Sotkamo All

98 30 7 7 6 66 24 12 250

25 × 24 30 × 30 rb = 12.5 r = 12.5 r = 12.5 r = 12.5 Variablec r = 12.5 Range Mean ± SD

2.5–96.8 43.7–97.5 24.0–74.0 17.8–87.3 13.5–92.4 6.0–78.2 33.1–96.5 45.2–94.6 2.5–97.5 60.3 ± 21.3

257–15,920 433–2883 552–7630 383–4605 682–12,470 95–17,080 623–8700 965–8815 95–17,080 2575 ± 2446

1.8–29.2 9.8–35.5 2.7–27.6 3.4–21.9 1.3–14.4 0.4–24.4 2.2–33.1 3.4–25.1 0.4–35.5 15.0 ± 7.3

1.4–36.6 10.4–66.5 3.6–35.8 3.8–28.9 0.5–19.1 0.0–34.6 1.1–47.0 2.6–34.6 0.0–66.5 18.5 ± 10.1

0.3–44.8 15.1–43.4 0.8–27.0 2.0–38.0 1.0–26.0 0.0–34.0 0.8–54.1 1.5–32.0 0.0–54.1 19.5 ± 10.7

a b c

Diameter at-breast-height (1.3m). Circular plot with radius r. Mostly 30 × 30m, but also varied from 20 × 30m to 35 × 40m.

spatial heterogeneity of vegetated surface, habitat, rain and light interception, and forest microclimate (Jennings et al., 1999). Canopy cover is defined as the proportion of ground covered by the vertical projection of the tree crowns (FAO, 2004; IPCC, 2003; Jennings et al., 1999), including small gaps inside the crown perimeter (Gschwantner et al., 2009). Unbiased measurement of canopy cover using vertical sighting tubes (Korhonen et al., 2006) is a timeconsuming process. Thus, canopy cover (CC) is seldom included in traditional forest inventory, or is assessed visually either in the field or from aerial photos. Two implications therefore arise: firstly, remote sensing techniques as alternative to laborious field measurements of CC are much needed; and secondly, the accuracy of the remote sensing techniques should be assessed most reliably with unbiased vertical measurements of CC. Regarding the latter, use of aerial photos to generate reference CC data may be inaccurate (McIntosh et al., 2012), while visual estimates of reference CC from high spatial resolution images can vary substantially between human interpreters (Montesano et al., 2009). Currently the most accurate remote sensing technique to estimate CC is the airborne LiDAR given its measurement principle and geometry (e.g., Griffin et al., 2008; Korhonen et al., 2011). However, LiDAR is not suited for continuous, large-scale CC monitoring. The planned NASA satellite LiDAR mission GEDI scheduled for launch in 2018 samples between 50N and 50S latitudes thus leaving the boreal region uncovered. Therefore, passive optical satellites remain an important data source for wall-to-wall CC monitoring. Among the operational optical Earth observation satellites, Landsat is one of the most promising data sources for global CC monitoring as they are freely available, have high spatial resolution and revisit frequency, and provide long term data record (continuity) spanning more than four decades. CC has been globally generated in the past first at 1 km resolution with images from the AVHRR sensor on board NOAA satellite (DeFries et al., 1999), then at 500 m (Hansen et al., 2003) and 250 m (DiMiceli et al., 2011) resolution with MODIS sensor on board Aqua and Terra. More recently, Sexton et al. (2013) further rescaled the MODIS Vegetation Continuous Field (VCF) product to Landsat resolution (30 m) and showed good agreement between the two. The improvement in spatial resolution is a significant progress as previously validation efforts have been lacking due to uneasy sampling for reference ground measurements at MODIS coarse resolution. Importantly, to our knowledge the new Landsat Tree Cover Continuous field (TCC) product has not yet been validated in the boreal region. The aim of this present study was to assess the capability of Landsat image data to estimate CC in boreal forests across wide geographical area with diverse forest structure, stand development, species composition, and site type. We utilized a large ground reference dataset (N = 250) of unbiased vertical measurements of ground reference CC collected in south, central, and north Finland.

Fig. 1. The location of study sites.

The extensive dataset allowed us to evaluate a realistic scenario of predicting CC in a new area by leave-site-out cross validation procedure. We tested empirical regression methods suitable for proportion data (such as CC) and an alternative spectral mixture analysis for CC prediction. We also carried out a first validation of the new Landsat TCC global product in boreal forest using temporally consistent (same year, season) observations. Finally, we analysed the CC prediction errors to identify ways to further improve the prediction accuracy. 2. Materials 2.1. Field measurements Field measurements of CC were carried out at eight forest sites in Finland (Table 1, Fig. 1) in 2005–2009. The 250 plots were located purposely to collect as diverse a dataset as possible from each study site and as a whole. The dataset comprised 135 Scots pine (Pinus sylvestris L.), 95 Norway spruce (Picea abies L. Karst) and 20 birch (Betula spp. L.) dominated plots. The site type varied with 119 xeric,

120

Hadi et al. / International Journal of Applied Earth Observation and Geoinformation 53 (2016) 118–127

Table 2 A summary of Landsat 5 TM images. Site

Field date

Image date

Solar elevation (◦ )

Suonenjoki

June–Aug 05 July 06 May–June 06 May 07 May 07 May 07 July–Sept 07 June–July 08 July 09

09 July 05 28 July 06 19 June 06 02 June 07a 09 Aug 07 02 June 07 04 July 07 02 June 07a 05 Aug 09

48.2 45.1 49.7 49.7 42.1 44.9 45.5 49.7 41.7

Koli Tammela Joensuu Rovaniemi Sodankylä Hyytiälä Sotkamo a

Same scene.

90 mesic, and 41 herb-rich plots. The stand structure (development class) varied with 5 young seedling stands, 45 mature seedling stands, 66 young forest stands, 63 mature forest stands, 66 old growth forest stands, and 5 seed tree stands. Originally not intended for satellite-based estimation, the plot size and shape were variable (Table 1). CC above 1.3 m was measured at all plots in a consistent manner using the vertical sighting (Cajanus) tube by the same person. More details on the measurement protocol can be found in Korhonen et al. (2006). 2.2. Satellite data We used eight scenes (Table 2) from the Landsat 5 TM Surface Reflectance Climate Data Record (USGS, 2016). The scenes were acquired over the sites in the same year as field campaign (except Hyytiälä in which the former year image was used) during June–August to minimize phenological differences. The TM sensor has 6 reflective spectral bands: band 1 (blue, 0.45–0.52 ␮m), band 2 (green, 0.52–0.60 ␮m), band 3 (red, 0.63–0.69 ␮m), band 4 (NIR, 0.76–0.90 ␮m), band 5 (SWIR1, 1.55–1.75 ␮m), and band 7 (SWIR2, 2.08–2.35 ␮m). 2.3. Landsat tree cover (TCC) product The recently published Landsat-derived Tree Cover Continuous field (TCC) global product (Sexton et al., 2013) has been produced by training each Landsat scene with the coincident 250-m MODIS vegetation continuous field (VCF) product (MOD44B, DiMiceli et al., 2011) using ensemble regression tree approach. We downloaded TCC epoch 2005 (made from a mosaic of best quality Landsat pixels from 2003 to 2008) covering our study sites where field measurements were carried out in 2005–2008. To ensure temporal consistency, using the accompanying ‘idx’ layer which informs the input Landsat scene used for each pixel prediction in the TCC product, we selected TCC pixels that correspond (see Section 3.1) to our field plot area and used input scene acquired in the same year as our field measurement, as well as during summer months to minimize possible phenological differences. 3. Methods 3.1. Spectra and TCC extraction To minimize spatial (planimetric) mismatch between field plot area and corresponding area in the Landsat image and TCC, we extracted the spectra (i.e., at-surface bidirectional reflectance factor, BRF) and TCC using exact intersection between field plot and image area, and calculated the BRF and TCC as area-weighted average of the values of the pixels overlapping with the field plot area. Further, to account for residual geometric error in the images and error in field GPS measurements, we applied a systematic shift (x, y) to the field plot location and searched the optimal shift that

Fig. 2. Example correlation surface to identify optimal shift (x, y) applied to spectra extraction area that maximizes correlation between BRFred and CC in Sodankylä.

maximized correlation between CC and BRF values (Fig. 2) in red and SWIR which had highest correlation with CC. The shifts were applied at 2.5 m distance steps up to ±30m. The search distance was extended if the correlation peak was not identifiable. The optimal shift was applied in extracting Landsat TCC from the same study area. 3.2. Analysis of relationship between BRF and CC To assess the dependence of the relationship between CC and BRF on categorical variables site location, image scene, dominant tree species, site fertility type, and stand development class, we calculated the variance in BRF explained by the interaction between CC and each of these factors. We fitted a linear model with response BRF and explanatory variables CC:Site, CC:Scene, CC:Species, CC:Fertile, and CC:Development. The variance accounted by an interaction term which was not accounted by other interaction terms was calculated as the decrease in model R2 from the model with all interaction terms included together, to the model with the interaction term in question removed. 3.3. CC prediction with empirical regression We investigated two statistical models suitable for CC data (proportion, 0–1 interval): the beta regression (BetaReg, Ferrari and Cribari-Neto, 2004) and random forest (RanFor, Breiman, 2001). To be noted that when there are cases of CC = 0 or 100%, the zero-andone inflated BetaReg (Ospina and Ferrari, 2009) can alternatively be used. Landsat-based predictors utilized include BRFs of individual spectral bands, spectral vegetation indices (SVI), and the Tasseled cap (Tass) features. Four 2-band SVIs were assessed: the normalized difference vegetation index (NDVI), reduced simple ratio (RSR), and best (i.e., highest correlation with CC) two-band (i,j) combinations in normalized difference (NDi,j = (BRFi −BRFj )/(BRFi +BRFj )) and simple ratio (SRi,j = BRFi /BRFj ) forms. RSR (Brown et al., 2000) was calculated as follows: RSR =

BRFNIR BRFred

 BRF  SWIR1.max -BRFSWIR1 BRFSWIR1.max -BRFSWIR1.min

(1)

where the minimum (0.031) and maximum (0.254) reflectance values in shortwave infrared (BRFSWIR1.min and BRFSWIR1.max ) were determined from the pooled data (N = 250). We also tested using scene-specific BRFSWIR1.min and BRFSWIR1.max of forest pixels (based on 25m CORINE land cover map by Finnish Environment Institute) which provided lower correlation between RSR and CC.

Hadi et al. / International Journal of Applied Earth Observation and Geoinformation 53 (2016) 118–127

In addition to 2-band SVIs, we also tested 3-band (bands i, j, k) SVIs of various forms as in Verrelst et al. (2015): modified simple ratio mSRi,j,k = (BRFi -BRFk )/(BRFj -BRFk ); modified normalized difference mNDi,j,k = (BRFi -BRFj )/(BRFi +BRFj -2 ∗ BRFk ); three-band spectral indices 3BSIi,j,k = (BRFi -BRFk )/(BRFj +BRFk ); 3BSI Wangi,j,k = (BRFi -BRFj +2 ∗ BRFk )/(BRFi + BRFj -2 ∗ BRFk ); and 3BSI Tiani,j,k = (BRFi -BRFj -BRFk )/(BRFi + BRFj +BRFk ). The Tass procedure transforms the TM data structure into Tass components which correspond to physically meaningful surface characteristics: brightness (TassBright), greenness (TassGreen), and wetness (TassWet). We applied Tass transformation coefficients for reflectance factor data (Crist, 1985) to calculate the three Tass components, as well as angle (TasPowell et al., 2010) and sAngle = arctan(TassGreen/TassBright);  2

distance (TassDist = TassBright +TassGreen2 ; Duane et al., 2010). We first assessed simple univariate (single predictor) prediction with BetaReg. Different response link functions (default logit, probit, cloglog, and loglog) and predictor (X) transformations (log10 (X), sqrt(X), 1/X, (X+X2 +X3 )) were evaluated to account for nonlinearity. We then evaluated multivariate BetaReg by performing exhaustive search for best (lowest RMSECV , see Section 3.5) set of multiple predictors. The alternative RanFor model was tested with all predictors (without X transformations), as well as the exhaustive search approach similar to that applied in multivariate BetaReg.

3.4. CC prediction with spectral mixture analysis The task of estimating CC from remotely sensed data corresponds to estimating sub-pixel fractional (tree) cover, commonly achieved by spectral mixture analysis (SMA). In SMA (Adams et al., 1986), a mixed pixel signal can be modelled as a linear combination of pure spectral signatures of its constituent components (i.e., the endmembers, EM), weighted by their sub-pixel fractional cover. Our SMA implementation assumes the top-of-canopy BRF as a linear mixture of signal from two components: (1) the canopy (which includes tree foliage and living branches following CC definition (Gschwantner et al., 2009)); and (2) the understory: BRFi = fC,i ∗ BRFC,i + fU,i ∗ BRFU,i + εi ; fC,i +fU,i = 1; fC,i ≥ 0,fU,i ≥ 0 (2) where BRFi is the top-of-canopy BRF in TM band i; fC,i is fractional cover of canopy (sunlit and shaded) i.e., the CC to be estimated; fU,i fractional cover of understory (sunlit and shaded); BRFC,i reflectance factor of canopy EM (sunlit and shaded) in band i; BRFU,i reflectance factor of understory (sunlit and shaded) EM in band i; and εi residual (unmodelled spectra) in band i. To obtain physically meaningful cover estimates, the fractions were constrained to sum to unity, and each fraction to be non-negative.



121



a = BRFC,i -BRFU,i and intercept b = BRFU,i , solved with n observations (i.e., number of forest plots in training set). The understory EM spectrum is given by intercept b, while (a+b) gives canopy EM spectrum. 3.5. Accuracy assessment and model comparison We performed a leave-site-out cross validation (CV) i.e., we trained the models using data from all but one site (8-1=7 sites) and predict the independent data in the one site held out for validation (8 CV folds). We then calculated the cross-validated RMSE and bias:



RMSECV =

ni=1 (CCLandsat,i − CCfeild,i )2 n



biasCV =

ni=1 CCLandsat,i − CCfeild,i

(5)



n

(6)

where n is the number of plots, CCLandsat,i is predicted CC, and CCfield,i is reference CC measured in the field for plot i. We also calculated the pseudo R2 as square of the Pearson product-moment correlation coefficient between measured and predicted CC values. RMSECV and biasCV are reported in absolute% cover unit. The significance of bias was tested using t-test comparing predicted and measured CC values. Finally, to determine the significance of difference in accuracy between models, we tested differences in mean and variance of the model residuals (predicted minus measured CC) with ANOVA and Bartlett’s test, respectively. 3.6. Error analysis Finally, we investigated the robustness of the prediction models i.e., if the prediction errors depend on (1) forest structure i.e., mean tree height, DBH, basal area, and stand density; and (2) site/scene, species, and site fertility type. We calculated the proportion of variance in residuals (predicted minus measured CC) that was explained by these factors. Identifying the influential factors may suggest ways to further improve the accuracy. We fitted a multiple linear regression to the residuals with all (1) and (2) above as explanatory variables. Measured CC was also included as explanatory variable, as systematic larger error at low (overestimation) and high (underestimation) CC was observed. The residuals variance accounted by each explanatory variable was calculated using the Lindeman-Merenda-Gold (R2 lmg , Lindeman et al., 1980) metric which handles correlated regressors (as were the forest variables), by averaging the regressor’s sequential R2 contribution over orders of their inclusion and over model sizes. 4. Results and discussion

3.4.1. Endmember estimation Our initial attempt to extract EM spectra from the Landsat image using Minimum Noise Feature (MNF) transformation and Pixel Purity Index indicated lack of pure canopy pixels at Landsat resolution in boreal forest. Thus, we estimated EM spectra using the known fractional cover in training data set of our model crossvalidation procedure (Section 3.5). That is, we estimated BRFC,i and BRFU,i in Eq. (2) by substituting fC,i = field-measured CC and fU,i = (1 − CC): BRFi = CC ∗ BRFC,i + (1-CC) ∗ BRFU,i + εi





Generally similar peak in the correlation surface was observed for BRFred and BRFSWIR1 (not shown). BRFred was selected as it provided the shifts (Table 3) which gave slightly higher overall band-by-band correlation (rmax higher by 0.01–0.05) compared to BRFSWIR1 . 4.2. Relationship between CC and Landsat spectral predictors

(3)

Given two components, rearranging Eq. (3) we obtain: BRFi = BRFC,i -BRFU,i ∗ CC + BRFU,i + εi

4.1. Spectra (BRF) extraction

(4)

which effectively reduces to simple linear regression between CC and BRFi with slope

Relatively strong negative linear correlation was observed between CC and BRF in the visible and SWIR region (highest in red, r = −0.78), while CC appeared uncorrelated with BRFNIR (Fig. 3). BRFred and BRFSWIR1 response to CC saturated at around CC = 70%, and there was a slight increase in BRFSWIR1 above CC 80%. The absence of relationship between boreal forest CC and BRFNIR (as was

122

Hadi et al. / International Journal of Applied Earth Observation and Geoinformation 53 (2016) 118–127

Fig. 3. Relationship between field-measured CC and Landsat at-surface BRF.

Table 3 Optimal shift (x, y) applied to BRF extraction footprint (field plot area) that maximizes correlation (rmax ) between CC and BRFred . Scene (Site)

No. of plots (N)

x, y (m)

rmax

Suonenjoki 2005 Suonenjoki 2006 Koli Tammela Joensuu Rovaniemi Sodankylä Hyytiälä Sotkamo

73 25 30 7 7 6 66 24 12

−32.5, +7.5 −12.5, −12.5 +22.5, −45.0 −25.0, 0.0 −10.0, −7.5 −2.5, −27.5 −10.0, −20.0 +32.5, +10.0 +5.0, +7.5

0.70 0.86 0.68 0.90 0.88 0.88 0.76 0.77 0.40

Table 4 Variance in BRF explained by CC interaction with site location, image scene, dominant tree species, site fertility type, and stand development class. Explanatory variable

CC:Site CC:Scene CC:Species CC:Fertile CC:Development a

Explained variancea (%) in BRFred

BRFNIR

BRFSWIR

0.3 0.0 1.2 0.1 5.5

0.1 0.0 4.8 1.7 20.9

0.2 0.0 2.4 0.1 8.8

See Methods Section 3.2.

also observed for leaf area index e.g., Heiskanen et al. (2011)) could indicate that in this biome the understory (sunlit and viewed by the sensor at low CC) and canopy are equally bright in NIR (especially during growing season), and that stand BRFNIR decreases as the stand matures possibly due to increase in inter-crown shadowing. The relationship between BRFs and CC did not appear strongly dependent on dominant tree species (Fig. 3): the CC-species interaction accounted for <5% variance in BRFs (Table 4). At similar CC, broadleaved birch plots were only slightly brighter than coniferous plots, possibly due to within-plot admixture with other tree species. Spruce was slightly darker than pine. Stand development class was related to BRFs response to CC, especially in NIR where it accounted for 21% variance in BRF. The noticeable 10% higher BRFNIR

observed for some plots (Fig. 3d) were mature seedling stands. This could be due to two reasons: first, the typically short trees (mean tree height <9 m) in these stands allowing a larger proportion of the understory to receive direct sun illumination; and second, the generally abundant green cover over the understory in these stands; both increasing stand NIR reflectance. The BRFs-CC relation did not strongly depend on site, scene, or site fertility (contributing <2% variance in BRFs). Considering spectral predictors derived from multiple individual bands (Fig. 4), NDVI showed large scatter in its relationship with CC at low cover (<60%). For both SR and ND forms, visible BRFred and BRFgreen combination were best correlated with CC. RSR seemed to be the most promising index. All the other 3-band SVIs did not provide stronger correlation with CC, with the best one was 3BSI = (BRFred − BRFgreen )/(BRFblue + BRFgreen ) achieving r = −0.71. Among the Tass components (Fig. 4d–f) and all derived spectral predictors, TassWet had the strongest correlation with CC. TassWet contrasts the sum of reflectance in SWIR (which are sensitive to water content and hence stem volume), with the sum of visible and NIR bands. TassBright decreased as expected with increasing CC from low to moderate CC (typical during early forest succession) as the canopy increasingly obscures the understory, which is brighter than canopy. During this early succession, TassGreen did not show a declining pattern as in the theoretical simulation of successional reflectance trajectory of conifer forest with grass understory (Song et al., 2002). This was a direct consequence of the observed absence of correlation between BRFNIR and CC as TassGreen contrasts NIR reflectance with the sum of other bands. This unique spectral signature of boreal forest makes TassAngle and TassDist (r = 0.67 and −0.33 respectively) unsuitable to monitor CC changes based on the hypothetical opposite direction of TassGreen and TassBright spectral trajectory in this biome (i.e., TassGreen increases while TassBright decreases with increasing CC during forest succession/recovery; the opposite for tree cover removal). Angle and distance between TassBright and TassWet (not shown) had lower correlation with CC (r = 0.67 and −0.61, respectively) compared to TassWet alone. The slight increase in both TassBright and TassGreen at CC > 60% could be because as trees

Hadi et al. / International Journal of Applied Earth Observation and Geoinformation 53 (2016) 118–127

123

Fig. 4. Relationship between field-measured CC and derived Landsat spectral predictors.

continue to grow, at some point they begin to block their own shadows from sensor viewing (Song et al., 2002). 4.3. Prediction accuracy 4.3.1. Empirical regression Among the univariate BetaReg models (Table 5), the best single predictor was individual band BRFred (Fig. 5a) with RMSECV = 13.7%, R2 CV = 0.59, and bias = 0.3% (not significant). The transformations of more than one band into indices and Tasseled cap components did not improve the prediction accuracy. Different BetaReg link functions and individual predictor transformations (e.g., log10 (X)) did not increase accuracy either. The observed adequacy of individual band was quite surprising, but has also been observed in previous studies estimating CC in similar forest type (Smith et al., 2009; Kuusinen et al., 2015). However, previous related studies rarely reported the performance of individual bands against e.g., widely used indices, and therefore it remains a question whether this was in fact a common finding. Further, including more than one predictor in BetaReg (Table 6) did not increase prediction accuracy (Fig. 5b). RanFor did not outperform simpler BetaReg (Table 6, Fig. 5c). RanFor with prior selection of predictor subset by exhaustive search was slightly more accurate than when all predictors were included at once (RanFor with 4 best predictors had 0.7% lower RMSECV than RanFor with all predictors). The residuals (normally distributed, Shapiro-Wilk p > 0.05) were not significantly different in mean (ANOVA p = 0.98) and variance (Bartlett’s p = 0.99) among the empirical models. The lack of improvement by nonlinear model RanFor may suggest that in optical satellite-based CC estimation, the responsepredictor relationship is not highly nonlinear and therefore accounting for the sources of random (unsystematic) error (noise) is more important than adjusting the regression model (curve fitting). We also tested generalized additive model with spline smoothing with no improvement either (not shown). As geometric error was minimized in this study, the unsystematic error sources could be the residual atmospheric noise, the influence of adjacent

(neighboring) pixels on signal from a single pixel, and other forest characteristics that cause spectral variability. We analyzed the latter in Section 4.5. A common limitation of the empirical models was the apparent systematic overestimation at low CC and underestimation at high CC (>70%). This was due to the inherent stand reflectance response to CC variation in boreal forests because empirical models do not explicitly account for the understory signal contribution at low CC, while optical signal simply saturates in dense canopies. The overestimation at low CC is especially a practical issue as internationally CC threshold of 10% is used to classify forest and non-forest areas (FAO, 2004). Despite the large geographical area, multiple image scenes, and very diverse forest plots in our study, compared to previous studies on optically-based CC estimation in boreal forest (Korhonen et al., 2013; Zald et al., 2016) or other conifer forest sites (Pu et al., 2008; Smith et al., 2009; Stojanova et al., 2010), the highest prediction accuracy achieved in this study (RMSECV 13.7%) is within the published range (RMSE 12–18%). The achieved accuracy is also reasonably good considering the error in visual estimation of tree cover by human experts using high resolution imagery was shown to reach 18.7% RMSE (Montesano et al., 2009). 4.3.2. Spectral mixture analysis The estimated EM spectrum (Fig. 6a) shows noticeably higher understory EM BRF in visible and SWIR bands than canopy EM, along with similar values in NIR. The EM BRF agrees well with satellite-measured BRF of near-pure canopy and near-pure understory plots in visible and SWIR, while the EM BRF in NIR was lower than satellite measurements. Little difference in BRFNIR between near-pure canopy and near-pure understory was also observed in satellite measurements. Compared to field measurements of understory spectra (Majasalmi and Rautiainen, 2016), in general the magnitude of the estimated understory BRF was lower than field-measured BRF. Field understory measurements show noticeably higher BRFNIR than BRFSWIR unlike the estimated understory EM here. Across the leave-site-out EM estimation (Fig. 6b), the

124

Hadi et al. / International Journal of Applied Earth Observation and Geoinformation 53 (2016) 118–127

Table 5 Prediction accuracy of univariate BetaReg. Spectral predictors

RMSECV (%)

R2 CV

BiasCV a (%)

BRFblue BRFgreen BRFred BRFNIR BRFSWIR1 BRFSWIR2 NDVI RSR Best SR=BRFred /BRFgreen Best ND=(BRFred -BRFgreen )/(BRFred +BRFgreen ) Best 3BSI=(BRFred -BRFgreen )/(BRFblue +BRFgreen ) TassBright TassGreen TassWet TassAngle TassDist

15.7 14.7 13.7 22.6 14.6 14.9 17.2 15.0 15.4 15.5 15.5 18.6 22.6 15.0 17.5 21.1

0.45 0.52 0.59 0.20 0.53 0.51 0.38 0.51 0.47 0.47 0.47 0.24 0.00 0.51 0.36 0.04

−0.1 0.5 0.3 0.9 1.1 0.5 0.1 1.2 0.7 0.7 0.6 0.2 2.5 0.8 0.4 −0.3

a

All bias (prediction-measured) are not significant (t-test p-value > 0.1).

Fig. 5. Predicted CC against measured CC. Dotted line is linear regression between predicted and measured CC.

canopy EM showed no spectral variability, while understory EM was slightly more variable (also observed in Kuusinen et al., 2015), especially in NIR. Prediction accuracy was lower for SMA (RMSECV = 17.8%) than empirical models (Table 5, Table 7, Fig. 5d). Unmixing without the uncertain EM BRFNIR , or EM BRFblue typically affected by atmospheric noise, did not improve the accuracy. We also tested (not

shown) including shade component which endmember spectra was approximated as deep clear lake or darkest Landsat forest pixel (based on CORINE map) in addition to the estimated canopy and understory components, which did not improve accuracy either. Previously estimating mixed conifer forest CC from ASTER data, Smith et al. (2009) also obtained lower accuracy for

Hadi et al. / International Journal of Applied Earth Observation and Geoinformation 53 (2016) 118–127

125

Table 6 Prediction accuracy of multivariate BetaReg and RanFor with best (exhaustive search) multiple predictors. Model-n best predictor seta

RMSECV (%)

R2 CV

BiasCV b (%)

BetaReg-2 (BRFgreen , best ND) BetaReg-3 (BRFgreen , BRFred , best 3BSI) BetaReg-4 (BRFblue , BRFgreen , best SR, best 3BSI) RanFor-2 (BRFred , best SR) RanFor-3 (BRFred, RSR, best SR) RanFor-4 (BRFred , best SR, best ND, TassDist) RanFor-5 (BRFgreen , RSR, best SR, best ND, TassBright) RanFor (all predictors)

13.6 13.7 13.6 14.2 13.8 13.6 13.6 14.3

0.59 0.59 0.60 0.57 0.59 0.59 0.59 0.55

0.5 0.3 0.9 −0.4 0.8 0.3 0.8 0.6

a b

Based on exhaustive search. All bias not significant (t-test p-value > 0.1).

Fig. 6. (a) The estimated endmember (EM) spectral signature in comparison to BRF of Landsat pixels (averaged to minimize noise) with near-full CC (near-pure canopy) and near-zero CC (near-pure understory). Shown is example for Suonenjoki site i.e., EM spectra were estimated from other sites for leave-site-out cross validation. (b) The endmember spectral variability across the leave-site-out CV (8 fold).

Table 7 Prediction accuracy of spectral mixture model. Input

RMSECV (%)

R2 CV

BiasCV a (%)

Spectra residual

All bands Exclude BRFNIR Exclude BRFblue

17.8 18.1 17.8

0.53 0.55 0.53

1.2 0.3 1.3

0.041 0.010 0.041

a

Bias not significant (t-test p-value > 0.1).

SMA (RMSE = 22.7%) than empirical regressions based on individual bands and indices. The foremost source of mixture modelling error is the uncertainty in EM characterization (Song, 2005). Previous mixture modelling studies have made use of reference EM spectra measured in the field or laboratory. In this study, we focused on image-based EM as for large area application field spectral measurements may not always be available, and reference EM spectra may not match image spectra due to imperfect atmospheric corrections, different illumination conditions, and phenology difference (Song, 2005). It is also possible at the Landsat TM scale and due to significant spectral contribution of the understory layer beneath the canopy to boreal forest reflectance (25–45%, Rautiainen and Lukeˇs, 2015), parts of the errors were due to not accounting for the multiple scattering between canopy and understory i.e., the nonlinear spectral mixture. However, Kuusinen et al. (2015) showed the diffuse understory spectral contribution to forest reflectance is minor (<2%), while the quantitative interpretation of the nonlinear fraction (interaction between EM) in nonlinear SMA is still not well understood (Somers et al., 2011). Another limitation of SMA is it employs a fixed EM spectral signature and therefore does not account for spatial and temporal variability of EM spectra (Somers et al., 2011). This is especially an issue in boreal forest given the substantial spectral variability of the

understory. Therefore, multiple EM SMA approach (Roberts et al., 1998) which allows EM set to vary on per pixel basis may improve the accuracy. We did not estimate EM specific to tree species and/or site fertility (understory) type as the leave-site-out accuracy assessment did not allow reliable EM estimation by species and/or site fertility types from all of the training set realizations. The extrapolation approach required field measurements data with wide CC range not possible if the data were stratified by species and/or site fertility types, especially for the limited number of birch plots. Further, regardless of species and site fertility, spectral signature of tree and understory component also likely varies between stands in different development stages (Song et al., 2002). In the future, the endmember characterization can benefit from more reference CC data with the increasingly available airborne laser scanning data, along with the ongoing efforts in building global library of vegetation spectro-temporal signature from EO-1 Hyperion archive. The increasingly available airborne laser scanning data could serve as an accurate reference CC as alternative to field measurements providing more observations to improve the EM estimation approach in this study. Specifically, with more observations of different tree species with wide CC range, speciesspecific EM can be reliably estimated. Alternatively, very high resolution image would make possible direct EM extraction from the image as a single tree crown (possibly also the sunlit and shaded parts) and the understory can be resolved. 4.4. Landsat TCC validation with field-measured CC A total of 30 plots all from Suonenjoki field measurements in 2005 and 2006 had a concurrent (same year, summer months i.e., July–August) CC estimate from Landsat TCC 2005 epoch (path 187, row 016). Clear underestimation of plots with CC above 70% can be

126

Hadi et al. / International Journal of Applied Earth Observation and Geoinformation 53 (2016) 118–127

Fig. 7. (a) Comparison between Landsat TCC 2005 epoch and field-measured CC (N = 30) from same year and summer months. (b) Leave-site-out prediction for the same validation plots using best model (BetaReg-BRFred ) and local reference CC in this study.

observed (Fig. 7a). The overall error (RMSE = 13.0%) is lower than the validation of the product in a few temperate and tropical sites (RMSE = 17.4%; Sexton et al., 2013), however the linear agreement here is considerably poorer (R2 = 0.8 in Sexton et al., 2013). Compared to leave-site-out prediction using best model BetaReg-BRFred trained with local nationwide accurate ground reference in this study (Fig. 7b), the Landsat TCC product was considerably accurate in terms of overall error (RMSE) when temporally concurrent comparison was made as here, despite the underestimation of high cover. The systematic error (underestimation of high cover, overestimation of low cover) of the Landsat TCC product is a known issue which is inherited from the parent data i.e., MODIS VCF. To correct for the systematic errors, Sexton et al. (2013) proposed a linear post-calibration using the slope and intercept of linear regression between Landsat TCC and reference CC. This was recently applied in Taiga-Tundra ecotone (Montesano et al., 2016). The slope of the linear regression in our validation (0.43, Fig. 7a) is close to the value (0.46–0.49) for pre-calibrated Landsat TCC in Montesano et al. (2016), while the intercept is notably higher here (31.2% > 16–17%). Further validation with a larger number of temporally concurrent observations is however needed to examine whether the degree of systematic error is rather universal, thus allowing reliable postcalibration over large geographical area and possibly also globally. Further, attention must also be given to resolve the still relatively large unsystematic (random) error (as indicated by the low R2 ). Another likely reason for the disagreement was the definition, thus structural measurement difference between CC and the TCC: CC in forestry application includes within-crown small opening as part of the canopy, while TCC was rescaled from MODIS VCF which is defined as component of tree cover that obstructs skylight (Hansen et al., 2003) and thus excludes within-crown gaps from the measured tree cover (i.e., ‘effective’ CC, Rautiainen et al., 2005). However, managed boreal forests dominated by conifer trees with narrow crowns are typically open so that between-crown gaps contribute significantly more than within-crown small gaps to total stand/plot gap fraction. 4.5. Error analysis Forest structure (CC, height, DBH, basal area, stand density) in total explained a large proportion (43%) of the BetaReg model residuals and relatively smaller proportion (22%) of SMA residuals (Table 8). BetaReg residuals were particularly dependent on

Table 8 Variance in model residuals explained by forest structural variables, site location/image scene, dominant tree species, and site fertility class. g is number of groups/categories. Explanatory variables

CC (%) Mean tree height (m) DBH (cm) Basal area (m2 ha−1 ) Stand density (trees ha−1 ) Site/Scene (g = 8) Species (g = 3) Site fertility (g = 3) Sum

%Variance of residuals BetaReg-BRFred

SMA-6 band

24.5 4.4 1.8 9.7 2.5 16.0 6.8 1.6 67.3

5.2 7.6 3.0 1.9 4.3 19.8 13.5 0.6 55.9

the measured CC, followed by basal area which was the structural variable most highly correlated with CC (r = 0.77). This was the systematic error as noted above. This unresolved limitation of empirical models gives motivation for further improvements in alternative methods that account for understory spectral contribution in low CC stands and potential increase in inter-crown shadowing in high CC stands more explicitly e.g., SMA or physically based models. The latter requires intensive input variables often not available, and thus the more parsimonious SMA can be a better alternative for large area application especially because its errors were found to be more independent from forest structure than in empirical models. Site location/image scene accounted for considerable variance in residuals (>16%) in both models (Table 8). One possible reason was that the most northern site Sodankylä is very different from the sites located in central Finland, because the harsh climate conditions and poor soils result in considerably lower mean CC. However, it cannot be inferred exactly whether site or scene was responsible for this residual variance contribution as in our study the scene corresponded to site. The only exceptions to this were 1) the Suonenjoki site where we used two scenes (corresponding to field measurements in 2005, 2006), while 2) Hyytiälä and Tammela were covered in one scene (Table 2). Tree species accounted for relatively high variance in residuals of SMA model (14%). This suggests SMA may be improved with species-specific endmembers. SMA error was also quite dependent on tree height (explaining 7.6% residual variance), which could indicate the effects of different degree of shading by trees given the low sun elevation (<50◦ ) at the study area latitudes. Finally, the prediction errors were independent of

Hadi et al. / International Journal of Applied Earth Observation and Geoinformation 53 (2016) 118–127

site fertility classes possibly because significant spectral variability exists within each site fertility class.

5. Conclusion We demonstrated the feasibility of large area (i.e., across sites and image scenes) canopy cover (CC) prediction in boreal forests from Landsat images, in which simple beta regression with single predictor red reflectance was found optimal (RMSECV 13.7%). However, the differences with the next best spectral predictors (BRFgreen , BRFSWIR1 , RSR, or TassWet) were marginal (1.3% maximum RMSECV difference). Results showed the tendency of Landsat Tree Cover (TCC) global product and in general empirically-based prediction to underestimate high CC (>70%) and overestimate low CC (<30%) in boreal forests. TCC product accuracy was however reasonably good (RMSE 13.0%) compared to our leave-site-out prediction trained with local reference CC data (RMSE 14.0%). Owing to the unresolved saturation issue of empirical models, there is a strong motivation for future efforts to be directed to improvements of CC prediction with spectral mixture analysis, by accounting for endmembers spectral variability of different tree species and possibly different degrees of shading.

Acknowledgements The study was supported by Aalto University (grant 91160122) and the Academy of Finland (grants 13286390, 288142). We thank the two anonymous reviewers for their comments which greatly helped us to improve the analysis and writing in our manuscript.

References Adams, J.B., Smith, M.O., Johnson, P.E., 1986. Spectral mixture modeling: a new analysis of rock and soil types at the Viking Lander 1 Site. J. Geophys. Res. 91, 8098–8112. Bonan, G.B., 2008. Forests and climate change: forcings, feedbacks, and the climate benefits of forests. Science 320, 1444–1449. Breiman, L., 2001. Random forests. Mach. Learn. 45, 5–32. Brown, L., Chen, J.M., Leblanc, S.G., Cihlar, J., 2000. A shortwave infrared modification to the simple ratio for LAI retrieval in boreal forests: an image and model analysis. Remote Sens. Environ. 71, 16–25. Crist, E.P., 1985. A TM tasseled cap equivalent transformation for reflectance factor data. Remote Sens. Environ. 17, 301–306. DeFries, R.S., Townshend, J.R.G., Hansen, M.C., 1999. Continuous fields of vegetation characteristics at the global scale at 1-km resolution. J. Geophys. Res. 104, 16911–16923. DiMiceli, C.M., Carroll, M.L., Sohlberg, R.A., Huang, C., Hansen, M.C., Townshend, J.R.G., 2011. Annual Global Automated MODIS Vegetation Continuous Fields (MOD44B) at 250 m Spatial Resolution for Data Years Beginning Day 65, 2000–2010 Collection 5 Percent Tree Cover. University of Maryland, College Park, MD, USA http://www.landcover.org/data/vcf/. Duane, M.V., Cohen, W.B., Campbell, J.L., Hudiburg, T., Turner, D.P., Weyermann, D.L., 2010. Implications of alternative field-sampling designs on Landsat-based mapping of stand age and carbon stocks in Oregon forests. For. Sci. 56, 405–416. FAO, 2004. Global Forest Resources Assessment Update 2005 (FRA 2005): Terms and Definitions (No 83/E). FRA Working Paper 1, Rome. Ferrari, S., Cribari-Neto, F., 2004. Beta regression for modelling rates and proportions. J. Appl. Stat. 31, 799–815. Griffin, A.M.R., Popescu, S.C., Zhao, K., 2008. Using LiDAR and normalized difference vegetation index to remotely determine LAI and percent canopy cover. In: Proceedings of SilviLaser 2008: 8th International Conference on LiDAR Applications in Forest Assessment and Inventory, Edinburgh, UK, pp. 446–455. Gschwantner, T., Schadauer, K., Vidal, C., Lanz, A., Tomppo, E., di Cosmo, L., Robert, N., Duursma, D.E., Lawrence, M., 2009. Common tree definitions for national forest inventories in Europe. Silva Fenn. 43, 303–321. Hansen, M.C., DeFries, R.S., Townshend, J.R.G., Sohlberg, R., Dimiceli, C., Carroll, M., 2002. Towards an operational MODIS continuous field of percent tree cover algorithm: examples using AVHRR and MODIS data. Remote Sens. Environ. 83, 303–319. Hansen, M.C., DeFries, R.S., Townshend, J.R.G., Carroll, M., Dimiceli, C., Sohlberg, R.A., 2003. Global percent tree cover at a spatial resolution of 500 meters: first results of the MODIS vegetation continuous fields algorithm. Earth Interact. 7, 1–15.

127

Heiskanen, J., Rautiainen, M., Korhonen, L., Mõttus, M., Stenberg, P., 2011. Retrieval of boreal forest LAI using a forest reflectance model and empirical regressions. Int. J. Appl. Earth Obs. Geoinf. 13, 595–606. IPCC, 2003. Good Practice Guidance for Land Use, Land-Use Change and Forestry. Institute for Global Environmental Strategies (IGES), Hayama, Japan. Jennings, S.B., Brown, N.D., Sheil, D., 1999. Assessing forest canopies and understorey illumination: canopy closure, canopy cover and other measures. Forestry 72, 59–74. Korhonen, L., Korhonen, K., Rautiainen, M., Stenberg, P., 2006. Estimation of forest canopy cover: a comparison of field measurement techniques. Silva Fenn. 40, 577–588. Korhonen, L., Korpela, I., Heiskanen, J., Maltamo, M., 2011. Airborne discrete-return LIDAR data in the estimation of vertical canopy cover, angular canopy closure and leaf area index. Remote Sens. Environ. 115, 1065–1080. Korhonen, L., Heiskanen, J., Korpela, I., 2013. Modelling lidar-derived boreal forest canopy cover with SPOT 4 HRVIR data. Int. J. Remote Sens. 34, 8172–8181. Kuusinen, N., Stenberg, P., Tomppo, E., Bernier, P., Berninger, F., Kuusinen, N., Stenberg, P., Berninger, F., Tomppo, E., Bernier, P., 2015. Variation in understory and canopy reflectance during stand development in Finnish coniferous forests. Can. J. For. Res 45, 1077–1085. Lindeman, R.H., Merenda, P.F., Gold, R.Z., 1980. Introduction to Bivariate and Multivariate Analysis. Scott Foresman, Glenview, IL. Majasalmi, T., Rautiainen, M., 2016. The potential of Sentinel-2 data for estimating biophysical variables in a boreal forest: a simulation study. Remote Sens. Lett. 7, 427–436. McIntosh, A.C.S., Gray, A.N., Garman, S.L., 2012. Estimating canopy cover from standard forest inventory measurements in Western Oregon. For. Sci. 58, 154–167. Montesano, P.M., Nelson, R., Sun, G., Margolis, H., Kerber, A., Ranson, K.J., 2009. MODIS tree cover validation for the circumpolar taiga-tundra transition zone. Remote Sens. Environ. 113, 2130–2141. Montesano, P.M., Neigh, C.S.R., Sexton, J., Feng, M., Channan, S., Ranson, K.J., Townshend, J.R., 2016. Calibration and validation of landsat tree cover in the taiga-Tundra ecotone. Remote Sens. 8, 551. Ospina, R., Ferrari, S.L.P., 2009. Inflated beta distributions. Stat. Pap. 51, 111–126. Powell, S.L., Cohen, W.B., Healey, S.P., Kennedy, R.E., Moisen, G.G., Pierce, K.B., Ohmann, J.L., 2010. Quantification of live aboveground forest biomass dynamics with Landsat time-series and field inventory data: a comparison of empirical modeling approaches. Remote Sens. Environ. 114, 1053–1068. Pu, R., Gong, P., Yu, Q., 2008. Comparative analysis of EO-1 ALI and Hyperion, and Landsat ETM+ data for mapping forest crown closure and leaf area index. Sensors 8, 3744–3766. Rautiainen, M., Lukeˇs, P., 2015. Spectral contribution of understory to forest reflectance in a boreal site: an analysis of EO-1 Hyperion data. Remote Sens. Environ. 171, 98–104. Rautiainen, M., Stenberg, P., Nilson, T., 2005. Estimating canopy cover in scots pine stands. Silva Fenn. 39, 137–142. Roberts, D.A., Gardner, M., Church, R., Ustin, S., Scheer, G., Green, R.O., 1998. Mapping chaparral in the Santa Monica Mountains using multiple endmember spectral mixture models. Remote Sens. Environ. 65, 267–279. Sexton, J.O., Song, X.-P., Feng, M., Noojipady, P., Anand, A., Huang, C., Kim, D.-H., Collins, K.M., Channan, S., DiMiceli, C., Townshend, J.R., 2013. Global, 30-m resolution continuous fields of tree cover: landsat-based rescaling of MODIS vegetation continuous fields with lidar-based estimates of error. Int. J. Digit. Earth 6, 427–448. Sexton, J.O., Noojipady, P., Song, X.-P., Feng, M., Song, D.-X., Kim, D.-H., Anand, A., Huang, C., Channan, S., Pimm, S.L., Townshend, J.R., 2015. Conservation policy and the measurement of forests. Nat. Clim. Change, 1–6. Smith, A.M.S., Falkowski, M.J., Hudak, A.T., Evans, J.S., Robinson, A.P., Steele, C.M., 2009. A cross-comparison of field, spectral, and lidar estimates of forest canopy cover. Can. J. Remote Sens. 35, 447–459. Somers, B., Asner, G.P., Tits, L., Coppin, P., 2011. Endmember variability in spectral mixture analysis: a review. Remote Sens. Environ. 115, 1603–1616. Song, C., Woodcock, C.E., Li, X., 2002. The spectral/temporal manifestation of forest succession in optical imagery: The potential of multitemporal imagery 82, 285–302. Song, C., 2005. Spectral mixture analysis for subpixel vegetation fractions in the urban environment: how to incorporate endmember variability? Remote Sens. Environ. 95, 248–263. Stojanova, D., Panov, P., Gjorgjioski, V., Kobler, A., Dˇzeroski, S., 2010. Estimating vegetation height and canopy cover from remotely sensed data with machine learning. Ecol. Inf. 5, 256–266. USGS, 2016. Product Guide: Landsat 4–7 Climate Data Record (CDR) Surface Reflectance Version 6.3. ˜ J., Clevers, J.G.P.W., Verrelst, J., Rivera, J.P., Veroustraete, F., Munoz-Marí, Camps-Valls, G., Moreno, J., 2015. Experimental Sentinel-2 LAI estimation using parametric, non-parametric and physical retrieval methods −A comparison. ISPRS J. Photogramm. Remote Sens. 108, 260–272. Zald, H.S.J., Wulder, M.A., White, J.C., Hilker, T., Hermosilla, T., Hobart, G.W., Coops, N.C., 2016. Integrating Landsat pixel composites and change metrics with lidar plots to predictively map forest structure and aboveground biomass in Saskatchewan, Canada. Remote Sens. Environ. 176, 188–201.