Inverting Aboveground Biomass?Canopy Texture Relationships in a Landscape of Forest Mosaic in the Western Ghats of India Using Very High Resolution Cartosat Imagery

Aboveground Biomass?Canopy Texture Relationships a Landscape Very Resolution Cartosat Abstract: Large scale assessment of aboveground biomass (AGB) in tropical forests is often limited by the saturation of remote sensing signals at high AGB values. Fourier Transform Textural Ordination (FOTO) performs well in quantifying canopy texture from very high-resolution (VHR) imagery, from which stand structure parameters can be retrieved with no saturation effect for AGB values up to 650 Mg · ha − 1 . The method is robust when tested on wet evergreen forests but is more demanding when applied across different forest types characterized by varying structures and allometries. The present study focuses on a gradient of forest types ranging from dry deciduous to wet evergreen forests in the Western Ghats (WG) of India, where we applied FOTO to Cartosat-1a images with 2.5 m resolution. Based on 21 1-ha ground control forest plots, we calibrated independent texture–AGB models for the dry and wet zone forests in the area, as delineated from the distribution of NDVI values computed from LISS-4 multispectral images. This stratiﬁcation largely improved the relationship between texture-derived and ﬁeld-derived AGB estimates, which exhibited a R 2 of 0.82 for a mean rRMSE of ca. 17%. By inverting the texture–AGB models, we ﬁnally mapped AGB predictions at 1.6-ha resolution over a heterogeneous landscape of ca. 1500 km 2 in the WG, with a mean relative per-pixel propagated error <20% for wet zone forests, i.e., below the recommended IPCC criteria for Monitoring, Reporting and Veriﬁcation (MRV) methods. The method proved to perform well in predicting high-resolution AGB values over heterogeneous tropical landscape encompassing diversiﬁed forest types, and thus presents a promising option for affordable regional monitoring systems of greenhouse gas (GhG) emissions related to forest degradation.


Introduction
Earth's forests, which cover 30% of the total land area [1], are dynamic systems which are constantly in a state of change and drive/respond to the changes taking place in our environment.
Forests store and sequester carbon and play an important role in climate change mitigation.Deforestation and forest degradation are the second largest source of anthropogenic carbon emissions, after fossil fuel emissions, which reached a contribution of nearly 20% of the total greenhouse gas (GhG) emissions for the period between 1970 and 2004 [2].Although these estimates have been scaled downwards to 10%-15% of the total GhG emissions in recent studies [3][4][5], it appears that most of the emissions occur in the tropical forests [6].Earths tropical forests cover only about 13% of the total ice free land area, but account for nearly 50% of the total plant biomass on the planet [1,7].Recognizing the importance of reducing GhG emissions from tropical forests, the United Nations Framework Convention on Climate Change (UNFCCC) developed an international program on Reducing Emissions from Deforestation and forest Degradation (REDD), which proposes mechanisms to provide financial and technological support to developing countries for stabilizing/increasing their forest carbon stocks and keeping carbon emissions below a baseline [8].The implementation of this program depends upon our ability to accurately assess current forest carbon stocks, particularly aboveground biomass (AGB) [9] at operational spatial and temporal scales [10].Mapping AGB variation due to conversion of forests to open lands and vice versa (deforestation/reforestation), which creates drastic biomass loss or gain, are nowadays achievable with reasonable accuracy [11].However, assessing changes related to various intensities of forest degradation, i.e., when forest structure changes but does not entail a change in land use, remains a challenge [12].
Though one of the largest contributors, carbon emissions due to deforestation and forest degradation in tropical forests are still poorly quantified because of the use of different data, assumptions and methodologies for estimating carbon stocks [5].Given the extent of the tropical forests, inaccessibility and limited ground measurements, remote sensing has long been seen as a way to increase effectiveness in assessing tropical forest AGB at regional level.Combining forest structure parameters from stand level estimates with extensive forest cover maps derived from remote sensing have been shown to have the potential for providing reasonable biomass estimates at local, regional and global level [13,14].However, from stand level AGB estimates to wall-to-wall extrapolation using different satellite data, the assessment of carbon stocks or AGB for tropical forests is still spoilt by uncertainty [15] thus weakening the results at regional scales.One of the major reasons for uncertainty in AGB estimates from most medium to high spatial resolution space-borne optical and radar platforms is the saturation of remote sensing signals at fairly low levels of biomass (150-200 Mg•ha −1 ), far below the upper range of values found in tropical forests (>500 Mg•ha −1 ) [16].Synthetic aperture radar (SAR), is a promising remote sensing sensor to map forest structure because of the ability of longer wavelength (L-band and P-band) radar backscatter to penetrate forest canopies [13].However, AGB values saturate between 100-150 Mg•ha −1 [13,17] for L-band radar signatures and for up to 300 Mg•ha −1 for P-band data [13].Recently, small-footprint aerial LiDAR data have proven to be a robust tool for forest structure modeling and biomass estimations [18], immune to signal saturation over the complete range of observed biomass [19], thus opening promising perspectives for the monitoring of forest degradation.However, aerial LiDAR covers very small extents and is very expensive to regularly operate as part of a regional monitoring system.In addition, its use also requires flight authorizations, which may delay independent monitoring.Up to now, it has been mostly used at local scales to provide estimates of canopy heights and indirectly AGB, which can be used as a basis for calibration of moderate resolution satellite data from which regional level AGB prediction models are built [18,20].In the absence of airborne data, filling the gap between field data and medium/high resolution imagery is difficult.This is probably one of the main reasons why uncontrolled discrepancies, between ground-inventory-based and large-scale remotely sensed biomass extrapolations, still exist for tropical forests [21,22].Devising and validating mapping methods of forest stand structure parameters and AGB over vast areas of tropical forests that would not need up-scaling to medium resolution is a challenging task.Recent results of canopy texture analysis of very high resolution (VHR) panchromatic images provide encouraging prospects for describing forest stand structures from airborne [23,24] and spaceborne observations [25][26][27][28].Forest science has long used visual interpretation of aerial photographs to characterize and map forest stands, and foresters have long known that the aspect of forest canopy, which directly reflects the size distribution of tree crowns and inter-crown gaps, conveys information on stand structure.This old paradigm has received support from both empirical and theoretical studies dealing with allometries between tree parts, notably crowns and trunks [29].
The increasing availability of VHR digital canopy images pleads for the development of methods for extracting information from canopy aspect.One such method is Fourier Transform Textural Ordination (FOTO), which characterizes the spatial patterns in forest canopies through texture indices that can be linked with stand structure parameters such as AGB.The FOTO method decomposes the spatial variation of pixel values within the canopy image window into frequency domain, i.e., with respect to spatially periodic functions, to compute Fourier spectra [24].One of the advantages of FOTO is that it is unsupervised and requires minimal input from the user, an advantage because of the nature of many forest landscapes in the tropics for which field information is scarce.FOTO gives output often interpretable in terms of size distribution of crowns or gaps that make the canopy aspect.The FOTO method has been in development for the past decade and has been tested on different tropical forest sites across the world [25][26][27]30,31], and those case studies provided robust estimates of forest AGB, even at high biomass levels, based on empirical relationships between canopy texture and stand structure parameters.The method has been tested using a range of different VHR satellite images which are available nowadays such as IKONOS, Quickbird-2, Geoeye-1, Spot-6/-7, and even freely available Google Earth (GE) imagery.FOTO canopy texture indices and their calibration with field plot data proved to be useful for inverting information from VHR images to characteristics of vegetation; spatial patterns in semi-arid vegetation [24,32]; stand structure parameters in the tropical rainforests of French Guiana [30] and oil palm-tropical forest landscape in Malaysia [28]; AGB in wet evergreen forest in the southern Western Ghats (WG), India [26,33]; and mangroves in French Guiana [25].Regional level AGB estimation was achieved over large territories using two different satellite images in the Congo basin forests [27].However, there are still challenges when dealing with multiple images (from same or different sensors) which are needed to cover vast forest tracts.Indeed, acquisition parameters reflecting the geometry of scene, sun and sensor may be different [16], and thereby result in unavoidable inconsistencies in texture measures.
Calibrating remotely sensed AGB estimates also relies on the accuracy of field estimates of biomass.Plot level biomass estimates are rarely obtained directly [34].Estimates are instead derived from already established allometric equations that convert easily measurable tree attributes such as diameter at breast height (DBH), tree height and wood density into AGB [35].Mapping regional level biomass from remotely sensed data primarily consists of two steps: (1) field AGB estimates obtained from sample plots; and (2) correlating the field estimates of AGB with the remote sensing derived forest canopy indices to obtain biomass estimates at the scale of the output pixels derived from texture analysis.Strong errors may be associated with these estimates [35,36] and plot level AGB uncertainties have been investigated [36][37][38].However, the way these errors propagate up to the AGB inversion of the remote sensing signal has rarely been assessed (but see [39]), especially in the context of the AGB inversion from texture-derived indices from VHR optical imagery.
In this paper, we used the FOTO method to derive canopy texture indices from Cartosat-1a VHR images, covering thousands of hectares of dry deciduous (DD), moist deciduous (MD) and small patches of wet evergreen (WE) tropical forests in the Western Ghats (WG) of India.The goals of the study were to explore the potential of the FOTO method to devise a texture-based biomass model at regional level across different forest types and to propagate the field-based errors up to the final AGB estimates for the output pixels.Though the FOTO method has been developed with satellite images of primarily metric level resolution (pixel size ≤ 1 m), a secondary objective was to test the method on Cartosat-1a imagery having 2.5 m spatial resolution.Lastly, we intend to assess the overall error on AGB prediction integrating both uncertainties on field plot measurements and errors deriving from FOTO predictions of plot AGB values.

Study Area
The study area (Figure 1) is located in the Yellapur Forest Division, Uttara Kannada District, Western Ghats of Karnataka, India.The area comprises nearly 1500 km 2 , with Yellapur (14 • 57 44.108 N, 74 • 42 22.12 E) roughly at the center, with a mean elevation of 457 m.The study area is a transitional zone between young Deccan trap rock formations and older Archean shield crystalline rocks of the Indian Peninsula [40].The southwest (SW) monsoon occurring between the months of June to October is responsible for most of the rainfall received.Mean annual precipitation is 2500 mm and monthly mean temperatures range between 25 • to 33 • Celsius.The Western Ghats forms an orographic barrier, which exacerbates the SW monsoon and is responsible for the steep bioclimatic gradients which are found in the region [41].The study area comprises of highly mixed, heterogeneous forest landscape, and has a continuum from dry deciduous forests in the lower rainfall zones to the northeast, transiting to moist deciduous and high rainfall zones of secondary wet evergreen forests to the south and south west with small patches of primary wet evergreen forests [42].The tropical wet/semi evergreen, moist deciduous and dry deciduous forest types in the study area have been described by Champion and Seth [43] while Pascal [42] mapped them as part of a general map of forest types at a scale of 1:250,000.The forests in the study area were observed to be frequently disturbed and in different stages of degradation due to anthropogenic stresses [44].

Study Area
The study area (Figure 1) is located in the Yellapur Forest Division, Uttara Kannada District, Western Ghats of Karnataka, India.The area comprises nearly 1500 km 2 , with Yellapur (14°57′44.108′′N,74°42′22.12′′E)roughly at the center, with a mean elevation of 457 m.The study area is a transitional zone between young Deccan trap rock formations and older Archean shield crystalline rocks of the Indian Peninsula [40].The southwest (SW) monsoon occurring between the months of June to October is responsible for most of the rainfall received.Mean annual precipitation is 2500 mm and monthly mean temperatures range between 25° to 33° Celsius.The Western Ghats forms an orographic barrier, which exacerbates the SW monsoon and is responsible for the steep bioclimatic gradients which are found in the region [41].The study area comprises of highly mixed, heterogeneous forest landscape, and has a continuum from dry deciduous forests in the lower rainfall zones to the northeast, transiting to moist deciduous and high rainfall zones of secondary wet evergreen forests to the south and south west with small patches of primary wet evergreen forests [42].The tropical wet/semi evergreen, moist deciduous and dry deciduous forest types in the study area have been described by Champion and Seth [43] while Pascal [42] mapped them as part of a general map of forest types at a scale of 1:250,000.The forests in the study area were observed to be frequently disturbed and in different stages of degradation due to anthropogenic stresses [44].Forest cover in the study area (grey, forest; and white, non-forest), Yellapur division, Uttara Kannada District, Karnataka, Western Ghats, India.The binary map is prepared from a mosaic of Cartosat-1a base imagery.Crosses represent location of 21 1-ha field plots.

Datasets
Cartosat-1 is a stereoscopic earth observation satellite launched aboard ISRO's Polar Satellite Launch Vehicle (PSLV) on 5 May 2005 [45].It orbits the Earth in a circular, sun-synchronous (near polar) orbit at an altitude of 624 km.Cartosat-1 satellite carries two panchromatic sensors, a PAN-F (PAN Forward-pointing camera, with a forward tilt of 26 • ahead of nadir) and PAN-A (PAN Aft-pointing camera with a aft tilt of −5 • behind nadir).In this study, we use the images from the near nadir PAN-A camera, hereafter designated as Cartosat-1a images.For canopy texture analysis, a mosaic of two Cartosat-1a panchromatic images with a spatial resolution of 2.5 m was used to define the limits of the study area.The two Cartosat-1a scenes were acquired and orthorectified using UTM coordinate system.The two images were acquired on the same day under similar acquisition conditions (Sun Azimuth, 136.93 • and 136.69 • ; Sun Elevation, 51.70 • and 51.83 • ).For computing the normalized difference vegetation index (NDVI) for the study site, a mosaic of two orthorectified and atmospherically corrected LISS-4 multispectral images (spatial resolution = 5.8 m) was used and resized to Cartosat-1a limits.Both Cartosat-1a and LISS-4 images were acquired by National Remote Sensing Center (NRSC), Hyderabad, India.Characteristics of both sensors are specified in Table 1.Vegetation maps (scale, 1:250,000) of the Western Ghats of India established in 1982 [42] and updated in 1997 [46], were used as reference and are freely available as shapefiles [47].Shapefiles of rivers, water-bodies, location of major settlements and road network were extracted from OpenStreetMap ® online database [48].Water-body and river shapefiles were used to create a water mask.After masking out water-bodies and rivers, results from preliminary runs of FOTO (see below) were used to create non-forest area mask (which included villages, agricultural land, and open areas within forests) on the basis of their textural signature.ASTER Digital Elevation Model (30 m resolution) was used to compute the slope for the study area and areas with slope >25% were discarded as to avoid the illumination artifacts on tree canopy crowns [49,50].Outputs from all three steps were combined to create non-forest mask and delineate the workable forest area (Figure 1).The same non-forest mask was used to subset the LISS-4 image.NDVI image was computed from LISS-4 multispectral imagery, and classified into two broad forest classes viz.dry zone forests (dry deciduous forest) and wet zone forests (moist deciduous and wet evergreen forests), based upon the histogram of observed NDVI values.

Field Plot Establishment, Tree Inventory
Field sampling design, for identifying locations of ground control plots, integrated preliminary information on canopy texture gradient (fine to coarse grain) observed from preliminary FOTO runs of Cartosat-1a images, forest type information from IFP vegetation maps, topography, road and river networks and inputs from local forest officials.We established 21 permanent 1-ha plots, during two field missions in 2014 and 2015, covering the complete gradient of canopy texture and different forest types present in the study area.Of the 21 1-ha plots, 6 were established in dry deciduous (DD) forest, 8 in moist deciduous (MD) forest and 7 in wet evergreen (WE) forest.The 21 permanent plots were established following the RAINFOR [51] protocols and recommendations from Condit et al. [52].Each plot was a 100 m × 100 m perfect square marked using Trimble Differential Global Positioning System (DGPS) and Leica Total Station to ensure north-south orientation of the plots and to minimize geo-location errors.In each plot, stem DBH was measured at 1.3 m height from the ground or above the buttress if any, using calibrated measuring tapes.A total of 8913 individual stems ≥ 10 cm DBH were recorded in all the 21 1-ha plots.Total tree height (H) was measured for randomly selected trees in all DBH classes using a NIKON Forestry Pro laser rangefinder.Botanical nomenclature was recorded for all the trees, identified up to species or genus level, with reference nomenclature to the Herbarium of the French Institute of Pondicherry (HIFP).One hundred thirty-two tree species were recorded in the field belonging to 43 botanical families.As we did not get permission for collecting wood core samples, wood density (ρ) values for the recorded tree species were compiled from the Global Wood Density (GWD) database [53].Where the wood density of the species was not available, genus level mean was assigned and where even genus level mean was not available, mean wood density of the plot was assigned to that particular species/genus.

Plot-Level Biomass Estimation and Error Propagation
Plot AGB estimation and error propagation were implemented by using the R BIOMASS package [54].The AGB of each tree is estimated from its DBH, total height and wood density using Chave's global pan-tropical allometric model (i.e., model 4 in [35]), that is: where AGB tree is the tree level AGB estimate in Mg, ρ is the specific wood density in g/cm 3 , DBH is the diameter at breast height in centimeters, and H is the total tree height in meters.As H was measured for a subset of total trees in each plot, a Michaelis-Menten model [55] was used to develop plot-specific H-DBH allometric relationships so to predict the missing H values.For a particular plot with no height information (YP07), allometric model coefficients from the plots within same forest type and similar species composition were used to estimate tree heights.Plot level AGB estimates (in Mg•ha −1 ) were obtained by summing individual tree AGB estimates.
For uncertainty analysis we assumed no errors were made on DBH measurements.For uncertainty in wood density values, a mean standard deviation (sd) was associated to each tree at the appropriate taxonomic level i.e., the one at which the mean wood density was calculated for the tree, using the approach in BIOMASS package.A plot-level mean sd was assigned to individuals for which wood density was not assigned at the species or genus level.For tree height error assessment the residual standard error (RSE) of the plot-specific H-DBH models was used.The uncertainties associated with plot-based AGB estimates were obtained as described in [54], i.e., through a Markov Chain Monte Carlo (MCMC) scheme, where the errors related to diameter, wood density and H-DBH model were propagated up to the plot-level with 1000 iterations.This step yielded probability distribution of the tree AGB values for each plot separately from which the mean AGB and 95% confidence interval were calculated along with other stand structure parameters, such as plot basal area (G) and quadratic mean diameter (QMD), reported in Table 2.

Fourier Transform Textural Ordination (FOTO)
FOTO [25] relies on 2-D spectral analysis by Fourier transform of panchromatic very high-resolution satellite images.The applications of 2-D Fourier transform to explain spatial patterns in digital images were described in detail in [56,57] and first applied to monitor vegetation patterns in aerial photographs [24,58].The method was further tested on tropical forest canopies [26,27,30,31].Using multivariate ordination, i.e., Principal Component Analysis (PCA) of Fourier spectra, it characterizes the canopy images with respect to canopy grain and allows quantification of patterns as present in tropical forest canopy.The first step is to divide the satellite image into unit windows of reasonable size which would allow for sufficient repetitions (3)(4)(5) of the largest patterns expected in the study area, here the crowns of large canopy trees (maximum of 25-30 m in crown diameter).For Cartosat-1a image the optimal window size was estimated to be 125 m × 125 m (50 × 50 pixels) after running the FOTO method with different window sizes and comparing the variance as explained by the eigenvalues of the PCA.The windows of relevant size were individually subjected to 2D Fast Fourier Transform (2D-FFT) so as to decompose the spatial variation of the spectral radiance within the image window with respect to frequency domain (for details refer to [24,30]).
Any digital image is an n × n array of pixels, depicting a range of digital numbers (DN) (e.g., 0-255), which express the radiance of each pixel.The 2D-FFT, for each unit window, computes the Fourier coefficients (a 2 pq , b 2 pq ) on mean corrected DN values.Each pair of coefficients corresponds to sine and cosine waveforms which are defined by harmonic wavenumbers, p and q along the two Cartesian XY directions in the geographical space.Wavenumbers p and q, quantify the spatial frequency by expressing the number of times each wave repeats itself in the window along a particular direction.The squared amplitude of Fourier coefficients is the 2D periodogram or power spectrum [56].It is a square array which represents the decomposition of the variance in image spectral radiance according to all possible integer pairs (p, q) of wavenumbers (with 1 ≤ p ≤ n/2 and 1 ≤ q ≤ n/2), where n is the window size in pixels.The periodogram can also be expressed in polar form by averaging the values with respect to bins of harmonic wavenumbers values (r) and orientations (θ).Within the scope of this study, we compute the radial spectrum or r-spectrum for which periodogram values were averaged by bins of harmonic wavenumbers (r = p 2 + q 2 ) over all possible orientations (θ = tan −1 p/q).Each of the window r-spectra provides the proportion of window variance explained by successive wavenumbers.R-spectra thus represent the spatial structure per sampling window in terms of spatial frequencies, which can be related to the canopy textural properties in terms of coarse to fine grain.
A Fourier r-spectrum was computed for each 125 m × 125 m window, all r-spectra being further assembled into a single matrix with individual windows as rows and the harmonic spatial frequencies of the Fourier r-spectrum as columns.However, high redundancy was regularly found at higher frequencies, i.e., shorter wavelengths, with more than half of the harmonic Fourier components with wavelengths shorter than 10 m.The issue was resolved by modifying the algorithm as presented in [59], where the size of the sampling window is doubled in each direction (zero-padding) [59,60] and r-spectra are linearly interpolated in wavelength space.This step allows for a more regular sampling between low and high spatial wavelengths and improves the characterization of large tree crowns and gaps [59].The number of columns retained after zero-padding was 13 Fourier components.These spatial wavelengths are used to detect the dominant pattern sizes as observed in the unit window, which are in our case expected to relate to the apparent crown sizes of the canopy trees or to canopy gaps.Although the two Cartosat-1a images had very similar scene geometries, r-spectra computed from the two images were inter-calibrated as introduced in [27,31], thanks to a overlapping area between the two images.Standardized major axis (Model II) linear regressions were performed for each Fourier harmonic component of the image r-spectra tables corresponding to the overlap region [29].The r-spectra of one of the images were then transformed using slopes and intercepts of the aforementioned regression model and concatenated with the r-spectra of the other image.The combined r-spectra table was then subjected to standardized PCA and window scores on the most prominent ordination axes were used as texture indices.

Texture-AGB Model Selection, Validtion and Inversion
Texture-AGB models were calibrated using multiple linear regressions between FOTO-derived texture indices and field-derived mean plot AGB values (after error propagation).Model performance was evaluated based on R 2 , mean root mean square error (RMSE) and mean residual standard error (RSE) computed through a Leave One Out Cross Validation (LOOCV) scheme [61].Models were adjusted considering: (i) all 21 plots together; and (ii) with a stratification scheme based on forest ecological classes (dry zone and wet zone forests) based on NDVI segmentation from LISS-4 image.The best model for inversion was selected by comparing FOTO-derived AGB values (i.e., predictions of the above multiple linear models) to plot-derived mean AGB values (i.e., estimations of the pan-tropical allometry).
To propagate errors associated with field-based AGB estimates up to the AGB map, we built independent texture-AGB regressions using the 1000 simulated field-based AGB estimates provided by the MCMC error propagation method.For each model, a LOOCV scheme was applied to assess the associated residual standard error (RSE).A random error was then assigned independently to the AGB prediction for each pixel using a normal distribution with a mean of 0 and a standard deviation equal to the model-specific RSE.This resulted in 1000 AGB estimations per output pixel from which we calculated the per-pixel coefficient of variation (CV), computed as the ratio of standard deviation to the mean of AGB estimations associated with the pixel.The final step is to project the above estimated values to come up with regional level maps of mean predicted AGB, and CV.

Landscape Scale Mapping of Forest Canopy Texture
After delineation of non-forest areas, the Cartosat-1a panchromatic image was subjected to FOTO texture analysis.Each image was divided into 50 pixel windows (1 pixel = 2.5 m) and an r-spectrum is computed for each of the unit windows.The r-spectra tables from the two images were then inter-calibrated and the complete resulting table was subjected to standardized PCA which yielded three prominent PCA axes explaining about 57% of the total variability (Figure 2). Figure 2a shows the scatter-plot of all 125 m × 125 m Cartosat-1a image windows on the first two PCA axes along with the sample plot windows illustrating the canopy texture from fine to coarse texture.Canopy windows corresponding to ground control points (black dots) having fine-grained canopy texture (e.g., YP09, YP10, and YP20) are displayed on the extreme negative side of PCA1 axis, while those having coarse-grained canopies (e.g., YP16 andYP17) are on the positive extreme of PCA1.The windows with intermediate canopy texture were more near the center of the PCA plane.Windows with high scores on the positive side of PCA2 displayed a mixed, coarse to intermediate texture gradient, signifying the simultaneous presence of smaller and larger spatial patterns.The correlation circle (Figure 2c), shows that the first principal axis (PCA1) had higher loadings for shorter spatial wavelengths (≤6, i.e., pattern sizes ≤15 m, fine textures) and larger spatial wavelengths (≥10, i.e., pattern sizes ≥25 m, coarse textures).The r-spectra of selected canopy image windows with fine, intermediate and coarse grain canopy textures are shown in Figure 3, illustrating that the r-spectra for finer textures are skewed towards low spatial wavelengths and the r-spectra for coarse texture canopies are skewed towards larger spatial wavelengths (i.e., with lower repetition of patterns per unit window).Using the geographic coordinates associated with each unit window, PCA scores on the first three axes were mapped as RGB composite image (Figure 4b), which shows the highly mixed and frequently heterogeneous nature of the forests present in the study area.The blue color in the PCA image (low scores on PCA1) corresponds to fine canopy texture, whereas purple to pink color gradient correspond to intermediate to coarse canopy texture.Note that the mask of non-forest areas has seemingly increased in the PCA image as compared to the original Cartosat-1a image because the 2D-FFT computation in the FOTO method with zero-padding excludes all the windows at the mask edges (Figure 4a,b).An NDVI image was computed from LISS-4 multispectral imagery using histogram thresholding (Figure 4c).It showed bimodality with, the NDVI values from −0.23 to 0.35 for dry zone forests and 0.35 to 0.75 for the wet zone forests (Figure 4c).The values for wet zone forests were within the range specified in [62] (i.e., 0.2 to 0.8 with a mean of 0.5) for the moist deciduous forests in the Yellapur area.The threshold of 0.35 was thus used to create two distinct forest masks, one for dry   An NDVI image was computed from LISS-4 multispectral imagery using histogram thresholding (Figure 4c).It showed bimodality with, the NDVI values from −0.23 to 0.35 for dry zone forests and 0.35 to 0.75 for the wet zone forests (Figure 4c).The values for wet zone forests were within the range specified in [62] (i.e., 0.2 to 0.8 with a mean of 0.5) for the moist deciduous forests in the Yellapur area.The threshold of 0.35 was thus used to create two distinct forest masks, one for dry An NDVI image was computed from LISS-4 multispectral imagery using histogram thresholding (Figure 4c).It showed bimodality with, the NDVI values from −0.23 to 0.35 for dry zone forests and 0.35 to 0.75 for the wet zone forests (Figure 4c).The values for wet zone forests were within the range specified in [62] (i.e., 0.2 to 0.8 with a mean of 0.5) for the moist deciduous forests in the Yellapur area.The threshold of 0.35 was thus used to create two distinct forest masks, one for dry zone forests and one for wet zone forests.From dry zone NDVI values we can infer that most of the dry deciduous forests are located on the northeastern portion of the study area, which conforms with the vegetation maps established by the IFP based on intensive field survey [42,46] and with our own observations during the two field missions.The NDVI masks are then used for creating two separate PCA images corresponding to dry zone and wet zone forests which are further used for AGB inversion.

Canopy Texture-AGB Relationships
We fitted, for the 21 1-ha field plots in Table 2, multiple linear regression models between field-derived mean plot AGB values and FOTO-PCA scores on the first three axes.PCA images corresponding to dry zone and wet zone forests were used to fit two separate linear regressions models of mean AGB on PCA scores.Table 3 gives a summary of model parameters and errors related to mean coefficients after LOOCV.Texture-AGB model applied without any distinction based on texture classes or forest zones gave relatively poor, with mean relative root mean square error (rRMSE) of 26.7% and R 2 = 0.53 (Figure 5a).For wet zone ground control plots a single texture-AGB model was validated with rRMSE = 17.1% and R 2 = 0.76.The model is based on a gradient of canopy texture with high biomass plots corresponding to coarse grain canopies (large tree crowns) and low biomass plots corresponding to fine grain canopies (of smaller tree crowns).The predicted mean AGB values for the wet zone plots ranged between 141 Mg•ha −1 and 486 Mg•ha −1 without any signs of saturation.A separate texture-AGB model was fitted for dry zone ground control plots which predicted low biomass values ranging between 159 Mg•ha −1 to 228 Mg•ha −1 , with mean rRMSE of 0.9% and R 2 = 0.98.The dry zone plots reflected a gradient of canopy texture from fine to coarse grain, but with coarse grain canopies corresponding to plots with canopy gaps.With a mean absolute error < 11% between the field-derived AGB and FOTO-derived AGB estimates (13.6% for wet zone forests and 1.4% for dry zone forests), the FOTO canopy texture indices have proven to be robust in predicting AGB at plot level.The maximum absolute error for a plot was 29% (110 Mg) for the moist deciduous plot YP07.Comparison between predicted AGB values from two separately fitted models for dry and wet zone forests and the estimated AGB values gave a mean rRMSE of 16.7% and R 2 = 0.82 (Table 3, Figure 5b).Figure 5 shows the correlation between predicted AGB values obtained and field-derived AGB for both single texture/forest class model and the model implemented after grouping into ecological forest zones.We can see significant improvement in the correlation once we clustered the field plots based on ecological zones.Although we are concentrating only on AGB in this study, we also tested the potential of texture indices to predict other forest stand structure parameters.Stand basal area estimates were closely correlated to canopy texture indices (R 2 = 0.58), whereas quadratic mean diameter had weak correlation with texture indices (R 2 = 0.29), for the fifteen wet zone plots.Comparison between predicted AGB values from two separately fitted models for dry and wet zone forests and the estimated AGB values gave a mean rRMSE of 16.7% and R 2 = 0.82 (Table 3, Figure 5b).Figure 5 shows the correlation between predicted AGB values obtained and field-derived AGB for both single texture/forest class model and the model implemented after grouping into ecological forest zones.We can see significant improvement in the correlation once we clustered the field plots based on ecological zones.Although we are concentrating only on AGB in this study, we also tested the potential of texture indices to predict other forest stand structure parameters.Stand basal area estimates were closely correlated to canopy texture indices (R 2 = 0.58), whereas quadratic mean diameter had weak correlation with texture indices (R 2 = 0.29), for the fifteen wet zone plots.

Mapping Forest Aboveground Biomass
Texture indices, derived from Cartosat-1a and calibrated with plot-level AGB estimates using the forest zone-wise models in Table 3 were used to predict and map AGB values throughout the study area along with the associated relative error of prediction per pixel or CV (Figure 6).

Mapping Forest Aboveground Biomass
Texture indices, derived from Cartosat-1a and calibrated with plot-level AGB estimates using the forest zone-wise models in Table 3 were used to predict and map AGB values throughout the study area along with the associated relative error of prediction per pixel or CV (Figure 6).The AGB map shows a gradient from high to low biomass from southern part of the study area towards the north, which follows the rainfall gradient present in the study area.Wet zone forests are dominated by patches of wet evergreen forests with high biomass with predicted AGB exceeding 500 Mg•ha −1 and low relative error per-pixel (<16%), depicted by blue color in the CV image (Figure 6b).Secondary wet evergreen forests and moist deciduous forests constituted intermediate AGB range, whereas areas with degraded formation of wet evergreen forests and old forest plantations in the moist deciduous forests exhibited low biomass values in the prediction map (green pixels).Dry zone forests in the north-east are dominated by lower to intermediate predicted biomass values ranging between 75 Mg•ha −1 and 310 Mg•ha −1 and higher per-pixel relative error (>30%).

Predictive Power of FOTO in Estimating Stand Structure Parameters
In combination with the field information, forest canopy texture may provide us with information regarding status and type of forests present in the study area [27].In this study, we demonstrated that canopy texture analysis of VHR images does provide robust estimates of AGB, even in the presence of diverse tropical forest types in a heterogeneous landscape, as frequently observed in the WG of India.Most previous studies [25,26] reached positive results in the context of close canopy evergreen tropical forests and mangroves, and even succeeded [27] in dealing with a mosaic of forest types in Democratic Republic of Congo (Central Africa) when different inversion The AGB map shows a gradient from high to low biomass from southern part of the study area towards the north, which follows the rainfall gradient present in the study area.Wet zone forests are dominated by patches of wet evergreen forests with high biomass with predicted AGB exceeding 500 Mg•ha −1 and low relative error per-pixel (<16%), depicted by blue color in the CV image (Figure 6b).Secondary wet evergreen forests and moist deciduous forests constituted intermediate AGB range, whereas areas with degraded formation of wet evergreen forests and old forest plantations in the moist deciduous forests exhibited low biomass values in the prediction map (green pixels).Dry zone forests in the north-east are dominated by lower to intermediate predicted biomass values ranging between 75 Mg•ha −1 and 310 Mg•ha −1 and higher per-pixel relative error (>30%).

Predictive Power of FOTO in Estimating Stand Structure Parameters
In combination with the field information, forest canopy texture may provide us with information regarding status and type of forests present in the study area [27].In this study, we demonstrated that canopy texture analysis of VHR images does provide robust estimates of AGB, even in the presence of diverse tropical forest types in a heterogeneous landscape, as frequently observed in the WG of India.Most previous studies [25,26] reached positive results in the context of close canopy evergreen tropical forests and mangroves, and even succeeded [27] in dealing with a mosaic of forest types in Democratic Republic of Congo (Central Africa) when different inversion models were used for fine, intermediate and coarse canopy texture classes.However, our study is the first case for which both evergreen and moist deciduous stands (i.e., wet zone forest plots) were considered in a single inversion model.Texture-AGB relationships in the wet zone (R 2 = 0.76) are consistent with the results obtained in other tropical forests of the world [25][26][27]30,31], with AGB values ranging from 27 Mg•ha −1 in some degraded Congo basin forests in Africa [27] to more than 650 Mg•ha −1 in a protected forest in the Western Ghats of India [26].In addition, we were able to characterize a texture-AGB model separately for the dry zone plots (R 2 = 0.98).This demonstrated the robustness of the texture based indices in consistently predicting forest biomass estimated from the pan-tropical allometric model when applied either to wet zone or to dry zone forests.Also, the FOTO texture indices appeared successful in other studies in predicting the stand structure parameters such as apparent crown diameter [63], and total stand characteristics such as tree density, quadratic mean diameter, mean canopy height and basal area [30,33].For our study, texture indices were found to be closely correlated with basal area (R 2 = 0.63) but had a weak correlation with quadratic mean diameter (R 2 = 0.36) inside the wet zone forests. Te correlation of FOTO indices with total stem density was poor and this may be because of the heterogeneous forest structure i.e., a mosaic of forest succession stages from highly degraded and disturbed patches of secondary evergreen and moist deciduous forests to small patches of undisturbed evergreen forests, encountered in the study area.Another reason is the expected limitation of canopy images regarding the understory trees, as they are not visible in the canopy.For this reason, it is not expected that FOTO provides acceptable predictions of stem density systematically.However, this problem was not serious for AGB predictions because canopy trees are known to account for a large share of stand AGB [64,65].
For the present study, we benefited from two VHR images, nearly synchronous and acquired in very similar scene geometric conditions.However, a major issue for applying FOTO to vast areas through the use of multiple images (from same or different sensors), is that the images may display heterogeneity in scene geometry due to varying acquisition parameters which results in inconsistency of texture measures [16].Acquisition geometries can affect the textural properties in the image and can blur texture-structure relationships between different images [16,59].However, Barbier and Couteron [59] showed from sensibility analysis based on extensive set of simulated images that the problem is only serious to deal with scenes imaged in backwards mode.Different approaches for correcting texture indices have been in development lately, such as partitioned standardization of Fourier r-spectra [16,28] and standardized major axis (Model II) linear regression of r-spectra frequency bins, which is possible if overlap between images exist [27,59].The latter solution, i.e., regression of r-spectra based on frequency bins, is the most efficient and is particularly promising for diachronic studies that would be a part of degradation monitoring (e.g., to track logging impacts [59]).In the present study, we applied this principle for r-spectra corrections, thanks to the overlap between the two images for which good correlations were observed between the two r-spectra tables.However, images acquired with extreme backward viewing conditions close to hotspot mode are not liable to be corrected and should not be used [59].

Effect of Different Forest Types Present in the Study Area
In stark contrast to the homogeneous, continuous and closed canopy forests encountered in previous studies [25,26], forests in Yellapur study area are highly fragmented as seen from the Cartosat-1a masked image (Figure 4a,c).This is representative of most locations in the Western Ghats of India, and to an increasing share of tropical forest landscapes across the world.Thus, it becomes very difficult to find large homogeneous patches of evergreen forests.Fine texture classes mostly corresponded to the old growth plantations and degraded forest patches in wet zone forests, characterized by small crown size, few gaps and low AGB.Coarse texture classes corresponded to areas with large crown size and high biomass in wet zone forests.In dry zone forests with fairly open and low biomass stands, the coarse texture is determined by large canopy gaps.It was because of this inverse relationship (coarse texture due to canopy gaps-low biomass) exhibited by few dry zone plots (YP04 and YP22) that the regression model based on single texture/forest class gave poor correlation (R 2 = 0.53).Thus, texture information is not unequivocal in the study area and there is a gradient of textures, from fine to coarse, in both dry zone and wet zone forests of contrasting meaning.Because of this, texture-AGB models behave differently, with respect to stand structure parameters, in dry zone forests with fairly open canopies and plots with canopy gaps are classified as high biomass areas and hence high prediction error is observed when a single model is fitted across all the forest types (Figure 6b).Combination of reflectance information (NDVI) from LISS-4 image improved the overall texture-structure relationships by allowing discrimination of forest types into different zones and enabling us to fit a separate model for each of the two forest zones.It can be safely stated that we need to develop two different approaches for closed canopy and open canopy forests.In the case of forest stands with canopy gaps, it may be useful to increase the size of sampling plots above 1-ha for better quantification of canopy gaps.In addition, a lacunarity based statistic, found to be strongly associated with gap volume [23], may also be explored as a complement to FOTO in order to better quantify canopy gaps and to distinguish between gaps vs. crowns induced coarse grained textures.
Field measurement errors were propagated up to the plot level using an MCMC scheme in a Bayesian framework with uninformative priors [54] and were depicted by the 95% CI of the AGB estimates at the plot-level.The texture-AGB predictions at the plot level did not show any apparent saturation with a global rRMSE of 16.7%.The errors were then mapped up to the regional level at the scale of output pixel size and we found high relative errors in the dry zone forests (due to limited number of field plots) and lower relative errors in the wet zone forests.However, a very high correlation was observed for texture-AGB model for the six plots in the dry zone forests (R 2 = 0.98).Increasing the size of sampling plots in those fairly open, low biomass forest stands or increasing the number of sample plots, to reach accurate AGB estimations, would be affordable, in terms of cost and time, for operational services such as forest departments.

Potential of Cartosat-1a Imagery
Cartosat-1a images have provided encouraging results in characterizing forest stand structures through a systematic analysis of canopy features by Fourier textural ordination.The texture indices extracted from Cartosat-1a image did not show any saturation for high levels of biomass present in the study area, even more than 500 Mg•ha −1 .The aforementioned range of errors with Cartosat-1a images was comparable to those exhibited by higher resolution Quickbird-2/Geoeye-1 images [27] and IKONOS/GE images [33].We may also note that the same range of errors was also observed from airborne LiDAR based biomass predictions in [66].Cartosat-1a provides a comprehensive and large database of VHR satellite images at affordable cost, which can be exploited for large scale and multi-temporal analysis of variations in structure, composition and biomass of vast tracts of tropical forests.In addition to Cartosat-1a, other spaceborne systems such as SPOT-6/-7 [67] have a spatial resolution of 1.5 m and large swaths (60 km, which limits the need for inter-image corrections) and are of higher affordability than VHR images of higher resolutions (i.e., pixel size between 0.5-1 m).For the present study, though the discrepancies in sun and sensor geometries and image acquisition parameters were limited, we applied the recommended r-spectra corrections, as image overlap was available.However, in the case of regional level studies, successive Cartosat-1a images acquired in along-track mode, have overlapping areas from which models for correction of Fourier indices can be devised as in [27].With the growing demand for low-cost VHR datasets, Cartosat-1a offers an acceptable alternative albeit at a spatial resolution above 2 m.To the best of our knowledge, Cartosat-1a images have not previously been used for characterizing forest structure and predicting AGB of a large heterogeneous forest landscape.

Conclusions
In this study, we utilized a consistently evolving and robust technique (i.e., FOTO) to characterize forest canopy texture in terms of repetition of patterns in very high spatial resolution images which are readily available and affordable.FOTO texture analysis has proven to be efficient in different tropical forests across the world.Inclusion of reflectance information based on classical indices like NDVI, allows us to gain complimentary information to texture as to efficiently characterize and distinguish different forest types present in the study area.Texture indices were calibrated with the ground based AGB estimates, and the overall error resulting from the complete upscaling process leading from individual tree allometries to regional level FOTO predicted estimates of AGB has been computed and mapped.This is the first example of such an error assessment for pixel-level AGB predictions from texture-derived indices.
The FOTO method is a powerful tool in predicting stand structure parameters from canopy texture analysis of VHR satellite images.It remained relevant even when applied to slightly lower resolution images (2.5 m instead of 1 m spatial resolution) and heterogeneous forest landscapes.This study is a good promotion for the use of very high-resolution Cartosat-1a images, and can be used for large scale and multi-temporal analyses of forest above ground biomass.Other sources of images of similar resolution and large swath (e.g., SPOT-6/-7) can similarly contribute to improved characterization and mapping of tropical forests stand parameters and related AGB.

Figure 1 .
Figure 1.Forest cover in the study area (grey, forest; and white, non-forest), Yellapur division, Uttara Kannada District, Karnataka, Western Ghats, India.The binary map is prepared from a mosaic of Cartosat-1a base imagery.Crosses represent location of 21 1-ha field plots.

Figure 1 .
Figure 1.Forest cover in the study area (grey, forest; and white, non-forest), Yellapur division, Uttara Kannada District, Karnataka, Western Ghats, India.The binary map is prepared from a mosaic of Cartosat-1a base imagery.Crosses represent location of 21 1-ha field plots.

Figure 2 .
Figure 2. PCA on Fourier r-spectra of 50 × 50 pixels (125 m × 125 m) forest canopy windows extracted from Cartosat-1a panchromatic image collected over forests in Yellapur division of the Western Ghats of India.(a) Window scores of first two PCA axes (PCA1 and PCA2) in grey with the black dots representing canopy windows over the ground control plots (as supplementary individuals).Selected canopy image windows are displayed for illustration.The arrow at the bottom describes the general direction illustrating canopy texture from fine to coarse grain along PCA1; (b) eigenvalues showing the percentage of variance in the r-spectra matrix explained by each PCA axis; and (c) correlation circle between apparent spatial wavelengths (i.e., number at the end of arrows is the number of 2.5 m pixels) and PCA axes.Each arrow (radial direction) gives the PCA loadings of the particular spatial wavelength with the two PCA axes.

Figure 2 .
Figure 2. PCA on Fourier r-spectra of 50 × 50 pixels (125 m × 125 m) forest canopy windows extracted from Cartosat-1a panchromatic image collected over forests in Yellapur division of the Western Ghats of India.(a) Window scores of first two PCA axes (PCA1 and PCA2) in grey with the black dots representing canopy windows over the ground control plots (as supplementary individuals).Selected canopy image windows are displayed for illustration.The arrow at the bottom describes the general direction illustrating canopy texture from fine to coarse grain along PCA1; (b) eigenvalues showing the percentage of variance in the r-spectra matrix explained by each PCA axis; and (c) correlation circle between apparent spatial wavelengths (i.e., number at the end of arrows is the number of 2.5 m pixels) and PCA axes.Each arrow (radial direction) gives the PCA loadings of the particular spatial wavelength with the two PCA axes.

Figure 3 .
Figure 3. Fourier r-spectra of the three selected image windows corresponding to Fine, Intermediate and Coarse grain canopy texture.

Figure 4 .
Figure 4. (a) Input Cartosat-1a panchromatic image (grey-levels), non-forest area masked out (black); (b) RGB image of the window scores of the first three principal PCA axes on Fourier r-spectra (see Figure2).The blue color in the PCA image (low scores on, PCA1) corresponds to fine canopy texture, whereas purple to pink color gradient correspond to intermediate to coarse canopy texture; and (c) NDVI image of the study area with dry zone forests (dry deciduous) in yellow and wet zone forests (moist deciduous and wet evergreen) in green.

Figure 3 .
Figure 3. Fourier r-spectra of the three selected image windows corresponding to Fine, Intermediate and Coarse grain canopy texture.

Figure 3 .
Figure 3. Fourier r-spectra of the three selected image windows corresponding to Fine, Intermediate and Coarse grain canopy texture.

Figure 4 .
Figure 4. (a) Input Cartosat-1a panchromatic image (grey-levels), non-forest area masked out (black); (b) RGB image of the window scores of the first three principal PCA axes on Fourier r-spectra (see Figure2).The blue color in the PCA image (low scores on, PCA1) corresponds to fine canopy texture, whereas purple to pink color gradient correspond to intermediate to coarse canopy texture; and (c) NDVI image of the study area with dry zone forests (dry deciduous) in yellow and wet zone forests (moist deciduous and wet evergreen) in green.

Figure 4 .
Figure 4. (a) Input Cartosat-1a panchromatic image (grey-levels), non-forest area masked out (black); (b) RGB image of the window scores of the first three principal PCA axes on Fourier r-spectra (see Figure2).The blue color in the PCA image (low scores on, PCA1) corresponds to fine canopy texture, whereas purple to pink color gradient correspond to intermediate to coarse canopy texture; and (c) NDVI image of the study area with dry zone forests (dry deciduous) in yellow and wet zone forests (moist deciduous and wet evergreen) in green.

Figure 6 .
Figure 6.Model inversion and error estimation: (a) aboveground biomass (AGB) predicted using the model coefficients for dry and wet zone forest plots in Table 3; and (b) relative error (CV) in percentage.Pixel Size of the output maps: 125 m × 125 m.

Figure 6 .
Figure 6.Model inversion and error estimation: (a) aboveground biomass (AGB) predicted using the model coefficients for dry and wet zone forest plots in Table 3; and (b) relative error (CV) in percentage.Pixel Size of the output maps: 125 m × 125 m.

Table 2 .
Stand structure parameters of 21 1-ha forest plots inventoried during two field missions in 2014-2015 and 2015-2016.QMD, Quadratic mean diameter of the plot; Mean WD, Mean Wood Density; D max , Maximum diameter at breast height (DBH); H max , Maximum tree height; AGB, mean Aboveground biomass estimated from field measurements; CI, Confidence Interval of AGB estimates computed using Markov Chain Monte Carlo (MCMC) error propagation scheme (see text).

Table 3 .
Parameters of AGB estimation and errors computed from Leave One out Cross Validation (LOOCV).a, b, and c are coefficients of PCA1, PCA2, and PCA3, respectively; RMSE, root mean square error; rRMSE, relative root mean square error; RSE, residual standard error; rRSE, relative residual standard error.