DOES TOPOGRAPHIC NORMALIZATION OF LANDSAT IMAGES IMPROVE FRACTIONAL TREE COVER MAPPING IN TROPICAL MOUNTAINS?

Fractional tree cover (Fcover) is an important biophysical variable for measuring forest degradation and characterizing land cover. Recently, atmospherically corrected Landsat data have become available, providing opportunities for high-resolution mapping of forest attributes at global-scale. However, topographic correction is a pre-processing step that remains to be addressed. While several methods have been introduced for topographic correction, it is uncertain whether Fcover models based on vegetation indices are sensitive to topographic effects. Our objective was to assess the effect of topographic correction on the accuracy of Fcover modelling. The study area was located in the Eastern Arc Mountains of Kenya. We used C-correction as a digital elevation model (DEM) based correction method. We examined if predictive models based on normalized difference vegetation index (NDVI), reduced simple ratio (RSR) and tasseled cap indices (Brightness, Greenness and Wetness) are improved if using topographically corrected data. Furthermore, we evaluated how the results depend on the DEM by correcting images using available global DEM (ASTER GDEM, SRTM) and a regional DEM. Reference Fcover was obtained from wall-to-wall airborne LiDAR data. Landsat images corresponding to minimum and maximum sun elevation were analyzed. We observed that topographic correction could only improve models based on Brightness and had very small effect on the other models. Cosine of the solar incidence angle (cos i) derived from SRTM DEM showed stronger relationship with spectral bands than other DEMs. In conclusion, our results suggest that, in tropical mountains, predictive models based on common vegetation indices are not sensitive to topographic effects.


INTRODUCTION
Landsat satellite images became available for free in 2008, which together with systematic data acquisition plan has fortified Landsat's role as a primary source of information for global land change research (e.g., Wulder et al., 2012). Landsat images also have good data consistency among Landsat missions and historical archive since 1972. The pre-processing methods of Landsat data have reached the level of maturity and preprocessed data sets have become available (e.g. Landsat Climate Data Record (CDR)) similar to moderate resolution data e.g. various MODIS data products (Masek et al., 2006). However, topographic correction remains as one of the pre-processing steps that have not been addressed globally.
Topographical effects in satellite images are caused by illumination differences between the sunlit and shaded slopes. The magnitude of the effect depends on the time of the year because of variations in the sun elevation. Hence, in the topographic normalization, the dependency of reflectance factors on topographic position is removed. Topographic correction has been shown to improve land cover mapping using object-based classification (Moreira & Valeriano, 2014) and land cover classification accuracy (Pellikka, 1996;Vanonckelen et al., 2013).
Topographic correction methods can be grouped into three categories: those based on band ratios, those based on a Hyperspherical Direction Cosine Transformation (HSDC), and those requiring digital elevation model (DEM). DEM-based methods, can be summarized as three broad types: (1) empirical methods, (2) Lambertian methods and (3) non-Lambertian methods (Gao & Zhang, 2009). The most often used methods are * Corresponding author.
Continuous fields of vegetation attributes can be estimated using a multitude of predictors derived from Landsat images. However, the topographic effect on different types of predictors can be different. For example, canopy cover, fractional tree cover (Fcover) and leaf area index (LAI) are typically estimated based on vegetation indices such as Normalized Difference Vegetation Index (NDVI) and Reduced Simple Ratio (RSR) (Majd et al., 2013;Korhonen et al., 2013;Wu, 2011).
NDVI is a simple index used to accentuate vegetation from imagery containing reflectance in the red and NIR portions of the spectrum. RSR is an empirical SWIR modification to the simple ratio (SR) vegetation index (Brown et al., 2000). Both indices attempt to depress background reflectance and improve the accuracy in extracting vegetation information from remotely sensed data.
The correction of satellite data for illumination differences due to Several authors have compared topographic correction methods for Landsat data (e.g., Hantson & Chuvieco, 2011;Moreira & Valeriano, 2014;Vanonckelen et al., 2013). However, these assessments are typically focusing on the removal of the topographic effect, instead of assessing how the accuracy of output products is affected. Although some studies have used topographic correction as a step to improve the land cover classification accuracy, it has been only rarely tested if topographic correction improves the accuracy of continuous variables (Törmä & Härmä, 2003), such as Fcover. Furthermore, classification is typically based on individual bands, but not vegetation indices, which are commonly used in Fcover, LAI and biomass mapping.
In this study, our objective was to assess the effect of topographic correction on the accuracy of Fcover predictions based on Landsat images with and without topographic correction. In particular, we aimed to examine if predictive models based on NDVI, RSR and TC indices are affected differently, and if results depend on the source of DEM.

Study area
The Taita Hills (3°25´S, 38°20´E) are located in the northernmost part of the Eastern Arc Mountains in southeastern Kenya ( Figure  1). The area is characterized by distinct topographical variation and the mountainous hills raise up to 2200 m a.s.l. from the Tsavo Plains at 600−900 m a.s.l. Taita Hills are considered as one of the world's most important biodiversity hotspots. However, due to the favorable climate and edaphic conditions, the indigenous cloud forests of the Taita Hills have suffered substantial deforestation and degradation due to agriculture, grazing and logging of forest for firewood and charcoal manufacturing since the early 1960´s (Pellikka et al., 2009).

Landsat images
Landsat surface reflectance (CDR) is a product of the Landsat Ecosystem Disturbance Adaptive Processing System (LEDAPS) at the National Aeronautics and Space Administration (NASA) Goddard Space Flight Center. The images have been calibrated to radiance, converted to top-of-atmosphere reflectance and then atmospherically corrected using the MODIS/6S methodology (Masek et al., 2006).
The Taita Hills lies on the south east quarter of the Landsat image (WRS path 167 and row 62). Landsat scenes with cloud cover less than 30% in the lower right quarter in 2013 were downloaded and two images less affected by cloud, shadow and missing data due to Scan line corrector (SLC) off and corresponding to the minimum and maximum sun elevation were used in this study. Both images were acquired at 7:32 UTC time in the morning and image processing level was L1T.

Digital elevation models
DEM used for topographic correction were obtained from various sources. Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) Global Digital Elevation Model (GDEM) has been generated from the stereoscopic ASTER satellite images. ASTER GDEM is available at the resolution of 30 m. The ASTER GDEM covers land surfaces between 83°N and 83°S. It is a product of Ministry of Economy, Trade, and Industry (METI) and the National Aeronautics and Space Administration (NASA) (https://www.jspacesystems.or.jp/ersdac/GDEM/E/4.html).
Shuttle Radar Topography Mission (SRTM) 1 arc-second global elevation data has been available since 23 September 2014 (https://lta.cr.usgs.gov/SRTM1Arc). SRTM provides DEM for 80% of the earth's surface (all land areas between 60° N and 56° S). SRTM is a joint project between the National Geospatial-Intelligence Agency (NGA) and NASA. The DEM has been filled to remove small artifacts.
All Landsat images, ASTER and SRTM DEM were downloaded from the United States Geological Survey (USGS) Earth Explorer platform (http://earthexplorer.usgs.gov/).
A regional DEM (TOPO DEM) for the Taita Hills was created using the contour lines of the topographic maps at 1:50 000 scale produced by the Survey of Kenya (Pellikka et al., 2004). DEM has a pixel size of 20 m × 20 m, and it was resampled to 30 m × 30 m pixel size similar to the Landsat images.

Airborne LiDAR data and fractional tree cover
We used Fcover (1-canopy gap fraction) derived from airborne LiDAR as a reference data. Previous studies have shown that LiDAR provides canopy cover and gap fraction estimates with an accuracy comparable to the field measurements (Korhonen et al., 2011;Heiskanen et al., in press). Furthermore, wall-to-wall LiDAR data provided a large sample size and the use of random sampling scheme in the logistically difficult terrain.  The discrete return LiDAR data was acquired 4-5 February 2013 for 10 km × 10 km area (Table 2). Data were pre-processed by the data vendor (Topscan Gmbh) and delivered as georeferenced point cloud in UTM/WGS84 coordinate system with ellipsoidal heights. Vendor also filtered ground returns by using Terrascan software (Terrasolid Oy). Furthermore, we filtered buildings and powerlines. Then, the ground returns were used to generate DEM at 1 m cell size. We also removed the overlap between the adjacent flight lines based on minimum scan angle using lasoverage tool in LAStools software (rapidlasso GmbH).
FUSION software (McGaughey, 2014) was used for computing a proxy of Fcover. We extracted LiDAR returns for 90 m × 90 m sample plots corresponding to 3 × 3 pixel windows of Landsat ETM+ images. Heiskanen et al. (in press) found that all echo cover index (ACI) (e.g., Morsdorf et al., 2006) gave an unbiased estimate of canopy gap fraction using the same LiDAR data. Hence, we used ACI as a proxy of Fcover. ACI was computed as: where, All refers to the returns of all return types (i.e. single, first, intermediate and last returns) and canopy refers to the returns from the forest canopy. Laser returns with height less than 1.5 m from the ground were considered as ground in order to exclude understory vegetation from Fcover.

Topographic correction of Landsat-data
C-correction is a semi-empirical topographic correction method, which consists of a modified cosine correction with the empirical parameter cλ (Reese & Olsson, 2011). cλ (Eq. 4) is calculated for every band (λ) separately based on the linear relationship (Eq. 3) between the spectral data and the cosine of the solar incidence angle (cos i) (Eq. 2). Each band were topographically corrected using Eq. 5.
where, i = solar incidence angle with respect to surface normal sz = solar zenith angle sl = slope az = solar azimuth angle as = aspect ρλ,t = topographically influenced (t) reflectance of band λ bλ = intercept of linear regression mλ = slope of linear regression cλ = c-factor calculated for every band separately ρλ,n = topographically normalized (n) reflectance of band λ cos i, slope and aspect were calculated from each DEM. Resolution of each DEM was made similar to Landsat image. Areas with a slope ≤ 2% were not included in the parameter estimation. Land cover stratification based on NDVI has been successfully used to consider the land cover dependency of cλ   Figure 3). 500 samples from each sub class, i.e. 2500 samples for each class and 7500 samples for each image were selected for linear regression between cos i and spectral bands. cλ was calculated for each NDVI class and each band.

Predictor variables for fractional tree cover
NDVI, RSR and TC indices (Brightness, Greenness and Wetness) were used as predictor variables. NDVI combines information only from NIR and Red bands, but RSR includes additional SWIR bands to reduce the effects of background reflectance (Brown et al., 2000). Since Landsat ETM+ images were in surface reflectance, tasseled cap coefficients for Landsat TM reflectance factor data were used.

Regression analysis
We used simple linear regression to model relationships between Fcover and NDVI, RSR, and TC indices. The results were evaluated based on coefficient of determination (R 2 ) and root mean square error (RMSE).

RESULTS
The various predictor variables explained variation in Fcover differently (Figure 4). Out of all predictors, the highest R 2 and the lowest RMSE were obtained using RSR. NDVI and Greenness showed the weakest relationships with Fcover. NDVI, RSR, Greenness and Wetness showed positive correlation with Fcover as they are affected by the amount of vegetation, while Brightness had negative correlation as it increases with the higher amount of open soil and lower vegetation cover.
The strength of the relationship also varied between the images acquired under different illumination conditions (Table 3). NDVI, RSR, Brightness, Greenness and Wetness explained more variation in Fcover when extracted from the image with higher sun elevation angle (29.9.2013). For example, topographically uncorrected Brightness indices explained only 37.5% of variation in Fcover when extracted from the lower sun elevation angle image in comparison to 45.2% when extracted from the higher sun elevation angle image. NDVI was also affected by the time of image acquisition. This was consistent for all the predictor variables.
The effect of topographic correction on Fcover regression models was not consistent for all predictors (Figure 4 and Table 3). The models based on Brightness and Wetness showed some improvement in terms of R 2 and RMSE but other models were not significantly affected, particularly when based on the larger solar elevation angle image. In the case of NDVI, there was no improvement at all due to topographic correction in either image. The source of DEM also affected the results in different ways ( Figure 4, Table 3), although differences were small like improvements in the model fit due to topographic correction. However, SRTM DEM performed better for Brightness and Greenness indices which were most affected by the topographic correction. Furthermore, when estimating c-factors by regression, the relationships between cos i and spectral bands were the strongest in terms of R 2 (Table 4) when derived from SRTM DEM.

DISCUSSION
According to our result RSR had the strongest linear relationship with Fcover obtained from LiDAR data. NDVI showed weak relationship with Fcover. These results are similar to those obtained by Korhonen et al. (2013) from a boreal forest site. However, our results differ from Majd et al. (2013) and Wu (2011) who found strong relationship between Fcover and NDVI in pecan orchards and savannas ecosystem. Weak relationship between NDVI and Fcover can be explained by the heterogeneous landscape in the study area and saturation of NDVI as the study area contains also patches of mountain rain forests.
However, it is difficult to say if this is due to weaker topographic effects or a result of variations in vegetation phenology. Image with lower sun elevation match better in terms of phenology to LiDAR acquisition.
Furthermore, our results demonstrate that topographic correction has very small or no effect on Fcover models. In the case of NDVI, no improvement was observed. Therefore, topographic correction might not be necessary for achieving the best Fcover prediction accuracy in tropical mountains that have relatively high solar elevation angles. In particular, models based on ratio based vegetation indices seem to be not affected by topographic correction. The fact that NDVI is less sensitive to topographic conditions has also been demonstrated by Matsushita et al. (2007). However, with multiplicative indices (e.g. Brightness), topographic correction is likely to have considerable effect.
Even though artifacts could be clearly observed in the SRTM DEM, it performed well in comparison to other sources of DEM, in agreement with similar previous studies (Vanonckelen et al., 2013;Balthazar et al., 2012). Most noticeably, SRTM DEM provided better R 2 in the C-correction than regional DEM based on the best topographic maps of the study area.
More comprehensive assessments on the effect of topographic correction in modelling vegetation attributes is necessary. Uncertainties related to the use of LiDAR data as a source of reference data needs to be further analyzed. Optimally, field data would be used but LiDAR data allowed us to sample Fcover from all topographical conditions with large sample size. Furthermore, we only considered one topographic correction method (Ccorrection) and results could be different when other topographic correction methods are used. However, C-correction has usually performed well in comparison to other DEM-based topographic correction methods. Furthermore, similar comparisons would be needed across larger latitudinal range because topographic effects increase with increasing latitude and decreasing solar elevations.

CONCLUSIONS
We assessed the effect of topographic correction on the accuracy of Fcover predictions based on Landsat images. Several sources of DEM were evaluated when performing the topographic correction. Among the tested indices, RSR provided the best regression model for Fcover. There were no major difference between the linear models based on topographically corrected or non-corrected NDVI and RSR, which indicates that NDVI and RSR are relatively robust against topographical effects. Therefore, in tropical mountain environments, studies based on NDVI or RSR do not necessarily need to consider the effect of topography. However, topographic correction should be considered for lower sun elevation angle images if models are developed using TC Brightness or Wetness. Finally, SRTM was shown to be the best DEM source for topographic corrections in the tested conditions.