Evaluation of FORMOSAT-2 and PlanetScope Imagery for Aboveground Oil Palm Biomass Estimation in a Mature Plantation in the Congo Basin

The present study developed methods using remote sensing for estimation of total dry aboveground biomass (AGB) of oil palm in the Congo Basin. To achieve this, stem diameters at breast height (DBH, 1.3 m) and stem heights were measured in an oil palm plantation located in Gabon (Congo Basin, Central Africa). These measurements were used to determine AGB in situ. The remote sensing approach that was used to estimate AGB was textural ordination (FOTO) based upon Fourier transforms that were applied, respectively, to PlanetScope and FORMOSAT-2 satellite images taken from the area. The FOTO method is based on the combined use of two-dimensional (2D) Fast Fourier Transform (FFT) and Principal Component Analysis (PCA). In the context of the present study, it was used to characterize the variation in canopy structure and to estimate the aboveground biomass of mature oil palms. Two types of equations linking FOTO indices to in situ biomass were developed: multiple linear regressions (MLR); and multivariate adaptive spline regressions (MARS). All best models developed yielded significant results, regardless of whether they were derived from PlanetScope or from FORMOSAT-2 images. Coefficients of determination (R2) varied between 0.80 and 0.92 (p ≤ 0.0005); and relative root mean-square-errors (%RMSE) were less than 10.12% in all cases. The best model was obtained using MARS approach with FOTO indices from FORMOSAT-2 (%RMSE = 6.09%).


Introduction
The amassing information on the total dry aboveground biomass (AGB) of plants is essential for understanding and monitoring their contributions to the global carbon cycle, together with their responses to climate change [1,2]. On small parcels of land, AGB can be estimated without too much difficulty from the complete biomass harvest of field sample plots and, subsequently, extrapolated to the unit area [3]. Another more common approach develops functional allometric relationships between field measurements of aboveground tree dry masses and structural parameters, such as DBH (diameter at breast height), stem height, basal area, wood density or crown area [3][4][5]. Over large areas, the financial constraints that are imposed by high costs of data collection, together with the large quantities of data that must necessarily be collected, have always limited the estimation of AGB.
Several studies have demonstrated the potential for using remote sensing as a complementary alternative for overcoming these limitations [5][6][7]. Predictive relationships can be established between structural variables and remotely sensed data [1,3,6,8,9]. These relationships are generally empirical, often have local or regional scopes, and are frequently difficult to transfer in space and time. In the Despite the diversity of remote sensing data available, the FOTO approach has not been used extensively. In particular, it has not been evaluated in combination with MARS or MLR approaches using various images, such as FORMOSAT-2 and PlanetScope, to describe canopy texture and estimate AGB. Studies on oil palms are rather limited or inexistent. The relationships between FOTO indices and structural parameters of oil palms have not been established yet in the tropical zones, except for young oil palms (≤9 years old) AGB [14].
In this paper, we applied the FOTO method on the near infrared bands of FORMOSAT-2 and PlanetScope to describe the textural structure of the canopy and assess AGB of mature oil palms (≥30 years) in the Congo Basin. Several oil palms AGB models were developed with MARS and MLR using the FOTO indices derived from each image in the study area. The models were compared to select those predicting the best results. This study was the first attempt to produce AGB map of mature oil palms in the Congo Basin using the combination of MARS or MLR and FOTO indices.

Description of the Study Area
The study site was located in Makouké, Moyen-Ogooué Province, Gabon. Gabon is a Central African country covered by moist tropical forest typical of the Congo Basin. Its boundaries extend from about 2 • N to 4 • S (Equatorial Guinea and Cameroon in the north, Republic of the Congo in the south) and from 9 • E to 15 • E (Atlantic Ocean in the west, Republic of the Congo in the east). The site is an oil palm plantation block that was located between 10 • 24 27" E and 10 • 24 57" E, and 0 • 30 6" S and 0 • 30 16" S ( Figure 1). The property is owned by Olam Palm Gabon (OPG), which is a company that has several agro-industrial oil palm plantations. This is the oldest plantation in the country, having been in production for 39 years. It covers 5700 ha, with mature oil palms occupying an area of 1500 ha. A watercourse crosses the site (North-South) in its eastern half ( Figure 2). Individual trees in the study consisted of mature palms, which are mainly in drier uplands, but some of them are found in the poorly drained shallows along the watercourse. Regional climate is warm-humid equatorial, with temperatures ranging from 27 to 38 • C, and annual rainfall ranging from 1800 to 2000 mm.

Satellite Data
Satellite images for this study come from PlanetScope (L3Harris Geospatial, Broomfield, CO, USA) and FORMOSAT-2 (Decommissioned, National Space Organization, Taiwan). PlanetScope images were acquired in April 2017 as an Analytic Ortho Tile Product (spatial resolution: 3.125 m × 3.125 m). The data were radiances at the top of the atmosphere and did not undergo further correction during the study. FORMOSAT-2 images are pansharpening products at 2 m × 2 m resolution (i.e., merging of high-resolution panchromatic and lower resolution multispectral images) from February 2011. Geometric corrections were applied to both types of images. Acquisitions were made in four similar spectral bands: 455-515 nm, 500-590 nm, 590-670 nm and 780-860 nm for PlanetScope; and 450-520 nm, 520-600 nm, 630-690 nm and 760-900 nm for FORMOSAT-2.

Satellite Data
Satellite images for this study come from PlanetScope (L3Harris Geospatial, Broomfield, CO, USA) and FORMOSAT-2 (Decommissioned, National Space Organization, Taiwan). PlanetScope images were acquired in April 2017 as an Analytic Ortho Tile Product (spatial resolution: 3.125 m × 3.125 m). The data were radiances at the top of the atmosphere and did not undergo further correction during the study. FORMOSAT-2 images are pansharpening products at 2 m × 2 m resolution (i.e., merging of high-resolution panchromatic and lower resolution multispectral images) from February 2011. Geometric corrections were applied to both types of images. Acquisitions were made in four similar spectral bands: 455-515 nm, 500-590 nm, 590-670 nm and 780-860 nm for PlanetScope; and 450-520 nm, 520-600 nm, 630-690 nm and 760-900 nm for FORMOSAT-2.

In Situ Measurements
Several field measurements were carried out during the study. Different structural parameters of oil palms that could be used for AGB calculations were measured. Given that the plantation was considered mature, all measurements were made on palm trees that were >30 years old. Instruments that were used in the field included a global positioning system (Garmin GPSMAP 64SC), compass (TOPOCHAIX), DBH tape and a vertex IV associated with a transponder 3. The GPS was used to position the sample plots in the field. The compass helped to materialize the directions of the transects. DBH were measured with a tape. The vertex coupled to the transponder allowed us to measure the stem height at two different positions at a distance of 10 m from the considered oil palm tree. Fieldwork was conducted in October 2017 on 40 sample plots that were established in the plantation, based upon systematic sampling. Plots (0.1 ha) were aligned East-West on four belt transects, each of which was 30.8 m wide and about 900 m long. The total area that was sampled was 25 ha. Natural obstacles, such as watercourses or tall grasses, limited accessibility and constrained sampling that was proposed along all survey lines. Figure 2 shows the sampling scheme for the sample plots. Stem height (HT) and diameter at breast height (DBH, 1.3 m) were recorded for all oil palms within each plot, as was their total number (NP). The density (D) values were determined from ratios between NP and individual plot areas. Basal area (BA, m 2 ha −1 ) was calculated [35] as the crosssectional area of each stem measured at 1.3 m above the ground, i.e., (π × DBH 2 )/4. Mean BA was estimated as the sum of all oil palms that were contained in each sample plot. Mean values for these variables are summarized in Table 1.
Mean in situ AGB (kg) was evaluated in each sample plot using the following ln-transformed allometric equation: This equation was derived from a previous study that had been conducted in the same area [36], which proposed different allometric models that were based upon in situ measurements that were collected during destructive sampling. For the rest of the present study, in situ AGB is expressed in tons ha −1 (Table 1). Table 1. Descriptive statistics for measurements taken from 510 oil palms distributed across 40 sample plots. DBH is diameter at breast height, HT is stem height, NP is the number of stems per plot, D is the density of palm trees per hectare, BA is basal area and AGB is total dry aboveground biomass

In Situ Measurements
Several field measurements were carried out during the study. Different structural parameters of oil palms that could be used for AGB calculations were measured. Given that the plantation was considered mature, all measurements were made on palm trees that were >30 years old. Instruments that were used in the field included a global positioning system (Garmin GPSMAP 64SC), compass (TOPOCHAIX), DBH tape and a vertex IV associated with a transponder 3. The GPS was used to position the sample plots in the field. The compass helped to materialize the directions of the transects. DBH were measured with a tape. The vertex coupled to the transponder allowed us to measure the stem height at two different positions at a distance of 10 m from the considered oil palm tree. Fieldwork was conducted in October 2017 on 40 sample plots that were established in the plantation, based upon systematic sampling. Plots (0.1 ha) were aligned East-West on four belt transects, each of which was 30.8 m wide and about 900 m long. The total area that was sampled was 25 ha. Natural obstacles, such as watercourses or tall grasses, limited accessibility and constrained sampling that was proposed along all survey lines. Figure 2 shows the sampling scheme for the sample plots.
Stem height (H T ) and diameter at breast height (DBH, 1.3 m) were recorded for all oil palms within each plot, as was their total number (NP). The density (D) values were determined from ratios between NP and individual plot areas. Basal area (BA, m 2 ha −1 ) was calculated [35] as the cross-sectional area of each stem measured at 1.3 m above the ground, i.e., (π × DBH 2 )/4. Mean BA was estimated as the sum of all oil palms that were contained in each sample plot. Mean values for these variables are summarized in Table 1. Mean in situ AGB (kg) was evaluated in each sample plot using the following ln-transformed allometric equation: This equation was derived from a previous study that had been conducted in the same area [36], which proposed different allometric models that were based upon in situ measurements that were collected during destructive sampling. For the rest of the present study, in situ AGB is expressed in tons ha −1 (Table 1).

Use of the FOTO Method and Estimation of Aboveground Biomass
Our methodological approach consisted of four parts: (i) the use of Fourier transformations in textural ordination (FOTO) to generate texture indices from the PlanetScope and FORMOSAT-2 images; (ii) development of MLR-type models for estimating AGB; (iii) development of MARS-type models for estimating AGB; and (iv) validation of the different models. The flowchart in Figure 3 summarizes the research approach.

Textural Index Generation Using the FOTO Approach
FOTO was implemented in this study following steps that are described in [37] and [8]. These steps allow gradients of variation in image texture to be characterized by coupling two-dimensional discrete fast-Fourier transformation (FFT-2D) with principal component analysis [8,14,18,20]. On one hand, the approach permits textural characterization of very high-resolution PlanetScope and FORMOSAT-2 images that were acquired over oil palm plantations; on the other hand, we assume textural indices that we generated are related to aboveground biomass of these plants [6,8,11]. In this study, FOTO was performed in MatLab R2016a (MathWorks, Natick, MA, USA). Near-infrared spectral bands of the two images were selected and analyzed for their sensitivity to vegetation [8]. The main steps are explained below.
The first step consisted of masking aspects of the images that were not related to the vegetation cover, i.e., bare ground, clouds and their shadows, among others. Images of the palm trees were then masked. We accomplished this masking using a combination of unsupervised classification with the k-means algorithm on the near infrared spectral band of each of the images in the study and visual examination to isolate unwanted features. Oil palm plantations were distinguished from natural forests and non-forest areas in the classification process. Following this operation, a mask was applied to the infrared band of each original image to retain only areas that were occupied by oil palms, thereby defining the area of application of FOTO in each case.

Textural Index Generation Using the FOTO Approach
FOTO was implemented in this study following steps that are described in [37] and [8]. These steps allow gradients of variation in image texture to be characterized by coupling two-dimensional discrete fast-Fourier transformation (FFT-2D) with principal component analysis [8,14,18,20]. On one hand, the approach permits textural characterization of very high-resolution PlanetScope and FORMOSAT-2 images that were acquired over oil palm plantations; on the other hand, we assume textural indices that we generated are related to aboveground biomass of these plants [8,6,11]. In this study, FOTO was performed in MatLab R2016a (MathWorks, Natick, MA, USA). Near-infrared spectral bands of the two images were selected and analyzed for their sensitivity to vegetation [8]. The main steps are explained below.
The first step consisted of masking aspects of the images that were not related to the vegetation cover, i.e., bare ground, clouds and their shadows, among others. Images of the palm trees were then masked. We accomplished this masking using a combination of unsupervised classification with the k-means algorithm on the near infrared spectral band of each of the images in the study and visual examination to isolate unwanted features. Oil palm plantations were distinguished from natural forests and non-forest areas in the classification process. Following this operation, a mask was applied to the infrared band of each original image to retain only areas that were occupied by oil palms, thereby defining the area of application of FOTO in each case.
The second step consisted of defining window sizes that were needed to calculate the 2D Fourier spectra on the different images. Following most previous work (e.g., [6,14]), we estimate these square window (WS) sizes in this study, such that WS = N*ΔS (N = number of pixels in the X or Y direction; ΔS = pixel size in meters). In a forest context, window sizes of 50 to 150 m are usually selected for applying FFT-2D [37,8,6]. Yet, Proisy et al. [8] demonstrated that window sizes between 75 and 120 m were best suited for characterizing near-infrared data in the IKONOS images that they examined. These sizes improved the results of principal component analysis (i.e., PCA), especially between 75 and 100 m. When choosing window sizes, the best possible compromises should be sought. Smaller sizes may not capture characteristics of the mature canopy, while sizes that are too large tend to  The second step consisted of defining window sizes that were needed to calculate the 2D Fourier spectra on the different images. Following most previous work (e.g., [6,14]), we estimate these square window (WS) sizes in this study, such that WS = N*∆S (N = number of pixels in the X or Y direction; ∆S = pixel size in meters). In a forest context, window sizes of 50 to 150 m are usually selected for applying FFT-2D [6,8,37]. Yet, Proisy et al. [8] demonstrated that window sizes between 75 and 120 m were best suited for characterizing near-infrared data in the IKONOS images that they examined. These sizes improved the results of principal component analysis (i.e., PCA), especially between 75 and 100 m. When choosing window sizes, the best possible compromises should be sought. Smaller sizes may not capture characteristics of the mature canopy, while sizes that are too large tend to represent landscapes more than canopy granularity, thereby reducing the spatial resolution of the FOTO maps [8,37]. Based upon these findings, we adopted three window sizes, i.e., 75 m, 100 m and 120 m. The WS report on ∆S provides the corresponding numbers of pixels N in the X and Y directions for these windows. The respective window sizes that were selected for the PlanetScope image are 24 × 24, 32 × 32 and 38 × 38 pixels. Window sizes (in pixels) for the FORMOSAT-2 image were respectively 38 × 38, 50 × 50 and 60 × 60 pixels.
The third step applied FFT-2D using the selected different sizes of square windows. FFT-2D yielded the radial spectra (or r-spectra), thereby converting the spectral radiance from the space domain to the frequency domain. Hence, a spectral radiance that is expressed in the space domain by I(x,y) of the x-column and y-row image is converted in the frequency domain into a function F(p,q), where p and q are spatial frequencies (or wavenumbers) along both Cartesian directions. These frequencies indicate the number of repetitions of an object within a given radial distance. They are expressed in cycles km −1 , making it easier to compare results that are obtained from different window sizes. For any given window size WS, the Fourier transform is applied over a set of discrete wavelengths, defined as the ratios of WS to wavenumbers [19]. The highest wavenumber corresponds to the Nyquist frequency of WS/2. Wavenumbers p and q are thus constrained to the following intervals: 1 ≤ p ≤ WS/2 and 1 ≤ q ≤ WS/2. Phase information is not very relevant to the texture and, therefore, can be neglected. Thus, periodograms (I pq ) can be calculated for each pair of spatial frequencies (p,q) using Fourier coefficients (a pq and b pq ), such that: The periodogram provides information on the variance of the image, according to the relative contributions of textural categories at various spatial scales. Based on its polar form, it is possible to conveniently characterize and recognize the distribution of surface patterns within the window considered, both in terms of scale and direction. The value of the periodogram in its polar form (I rθ ) reflects the portion of the variance of the image that can be explained by waveforms having a given wavenumber r and direction of motion θ, such that r = p 2 + q 2 and θ = tan −1 (p/q). The r-wavenumber corresponds to the number of times the same waveform, i.e., the same reference pattern, is repeated within the window in the θ direction.
Additional simplifications can be applied to periodograms when the objective is textural analysis. They consist of calculating the average value of periodograms in all possible directions of motion θ to generate an averaged radial spectrum I(r). This spectrum is commonly referred to as the "r-spectrum": where k is the number of periodogram values at spatial frequency r and σ 2 is the variance of the image. The r-spectrum thus makes it possible to summarize the characteristics of textural variation in terms of the state of grain size and fineness, which is very useful for vegetation studies such as oil palm plantations. More details regarding the calculation of r-spectra can be found in Renshaw and Ford [38], Mugglestone and Renshaw [39], Couteron et al. [37], Proisy et al. [8], and Guo and Rees [19]. As previously noted, the current study considers three increasing window sizes for each image and determines which one is the best. The optimal WS (window size) that was used to generate the FOTO textural indices corresponds to the one that exhibits the highest Nyquist frequency among the three. The Nyquist frequency is the maximum frequency that the signal can display without containing errors, i.e., one-half of the sampling frequency. It often takes a wavenumber between 1 and 3. The r-spectra and Nyquist frequencies were obtained for the different windows and summarized in Figure 4. The best window size is 38 × 38 pixels for the PlanetScope image, while 60 × 60 pixels is the best size for the FORMOSAT-2 image. WS is almost the same for the two images, i.e., 119 m for PlanetScope and 120 m for FORMOSAT-2.
The fourth step involves textural ordination. The r-spectra that are produced were incorporated into a two-dimensional table, the rows of which represent the r-spectrum of the given window and the columns are the amplitude I(r) values. Normalization was performed on the columns of the table, after which the normalized spectra were subjected to principal component analysis (PCA). PCA characterizes the variability and measures the dispersion between r-spectra, reducing the dimensionality of the data to a few components (three in the current study). Principal component axis 1 contrasts coarse (negative scores) and fine (positive scores) textures of canopy grain. Coarse textures generally represent large canopy sizes, while fine textures are more consistent with small-to medium-sized tree canopies. Principal component (PC) axis 2 could be interpreted as a gradient of canopy openness, from closed canopy (negative scores) to open canopy (positive scores) conditions, as has been indicated in various studies [6,37,[40][41][42]. By definition, the first three PCs that were extracted explain most of the variation in the data; factor scores for these components were used as textural indices (FOTO indices). Their values were mapped as Red-Green-Blue (RGB) texture images to produce FOTO maps, with a spatial resolution equal to the window size (WS). The fourth step involves textural ordination. The r-spectra that are produced were incorporated into a two-dimensional table, the rows of which represent the r-spectrum of the given window and the columns are the amplitude I(r) values. Normalization was performed on the columns of the table, after which the normalized spectra were subjected to principal component analysis (PCA). PCA characterizes the variability and measures the dispersion between r-spectra, reducing the dimensionality of the data to a few components (three in the current study). Principal component axis 1 contrasts coarse (negative scores) and fine (positive scores) textures of canopy grain. Coarse textures generally represent large canopy sizes, while fine textures are more consistent with small-to medium-sized tree canopies. Principal component (PC) axis 2 could be interpreted as a gradient of canopy openness, from closed canopy (negative scores) to open canopy (positive scores) conditions, as has been indicated in various studies [6,37,[40][41][42]. By definition, the first three PCs that were extracted explain most of the variation in the data; factor scores for these components were used as textural indices (FOTO indices). Their values were mapped as Red-Green-Blue (RGB) texture images to produce FOTO maps, with a spatial resolution equal to the window size (WS).

Estimation of AGB by Multiple Regression
Oil palm structural parameters from the 40 field sample plots were first averaged before proceeding to the AGB model step. To achieve this, field sample plots that were contained in the WS window sizes (120 m × 120 m) of both image types were identified for FOTO analysis. Parameter averages for these plots were then calculated for each window. Depending upon position, the window could contain one or a maximum of three plots. Mean values of oil palm structural parameters that were estimated were examined using scatter plots. Field data and calculations were checked for possible outliers.
Mean values of the principal component (PC) scores for the normalized r-spectra (corresponding to the FOTO indices) were calculated for each window and then combined with the plot averages of the structural parameters. ArcGIS V10.7.1 (ESRI, Redlands, CA) was used for the overlay and spatial processing of the data.
Before developing biomass models using MLR, we examined the relationships between FOTO indices and each palm structural parameter (DBH, HT, BA, D) estimated from field sample plots data in the study. This analysis identified parameters that best explained variations in the textural index, which would be the most influential in the MLR models. Finally, an MLR model was developed to . Nyquist and spatial frequencies as a function of mean r-spectra for window sizes ranging from 24 to 38 pixels for PlanetScope (top row) and 38 to 60 pixels for FORMOSAT-2 (bottom row).

Estimation of AGB by Multiple Regression
Oil palm structural parameters from the 40 field sample plots were first averaged before proceeding to the AGB model step. To achieve this, field sample plots that were contained in the WS window sizes (120 m × 120 m) of both image types were identified for FOTO analysis. Parameter averages for these plots were then calculated for each window. Depending upon position, the window could contain one or a maximum of three plots. Mean values of oil palm structural parameters that were estimated were examined using scatter plots. Field data and calculations were checked for possible outliers.
Mean values of the principal component (PC) scores for the normalized r-spectra (corresponding to the FOTO indices) were calculated for each window and then combined with the plot averages of the structural parameters. ArcGIS V10.7.1 (ESRI, Redlands, CA) was used for the overlay and spatial processing of the data. The values of the structural parameters and the FOTO indices for the study are presented by sample plots in Appendix A.
Before developing biomass models using MLR, we examined the relationships between FOTO indices and each palm structural parameter (DBH, H T , BA, D) estimated from field sample plots data in the study. This analysis identified parameters that best explained variations in the textural index, which would be the most influential in the MLR models. Finally, an MLR model was developed to estimate BA from the PlanetScope and FORMOSAT-2 images, respectively. The general form of the equation is as follows [8]: where a 0 and a C are the MLR coefficients, and FI represents FOTO indices that were obtained from scores for the first three axes of PCA. To develop MLR models for estimating AGB, 75% of the total 40 samples available were used. These data were randomly selected. The remaining 10 samples were held back for the independent validation of the results. Evaluation of various relationships that were developed (FOTO indices vs. structural parameters; AGB vs. FOTO indices) was performed using several statistical metrics. Relationships were deemed satisfactory when coefficients of determination (R 2 ) were greater than 0.5, with p < 0.05. (P-values were calculated at a significance level of 5%). RMSE (Root-Mean-Square Error) and %RMSE (Relative RMSE) were further used to evaluate AGB models with MLR. The lowest values were sought. Ultimately, the choice of the best model in the present study was based more upon RMSE than upon R 2 , following recommendations from previous studies (e.g., [43]).

Estimation of AGB Using MARS
According to various studies that were mentioned in the introduction (Section 1), the MARS approach seems well suited to accounting for non-linearities in the data. It models relationships through summation of a series of so-called "basis functions" [27], which in fact constitute a set of stepwise adaptive linear regressions. We applied this method for the first time in the estimation of oil palm AGB. The generic expression of the MARS model that is based upon the dependent variable (y) and explanatory variable (x) is written as follows [28,32,44]: where ε is the error and f (x) is the unknown regression function, which can be written as: where β 0 is the constant in the model, K m (x) is the m th basis function, β m is the coefficient associated with the function K m (x), and M is the number of basis functions in the model. Each basis function K m (x) consists of a hinge function of the form max(0, x-k) or max (0, k-x), where k is a constant value. This formulation takes into account individual effects of the independent variables plus those resulting from interactions between two or more predictors (x). MARS eliminates less significant and redundant basis functions to create the final model [25]. To do this, basis functions are compared using the Generalized Cross-Validation Criterion (GCV). GCV is the residual Root-Mean-Square Error divided by a penalty determined by the complexity of the model [28]: where

Validation of the Estimates
Validation of the different models that were developed for estimating AGB used 25% of the samples (10 samples), independent of the development or training phase. The different metrics that were considered in performance analysis included R 2 , p-values, RMSE and %RMSE [25,45]. The best model was subsequently used to produce a FOTO map of AGB. The best model was subsequently used to produce a FOTO map of AGB [28,46,47] using ArcGIS software.

Textural Indices Using the FOTO Approach
FOTO was applied to infrared bands of PlanetScope and FORMOSAT-2 images to extract r-spectra for the oil palms. Figure 5a shows the variations in these spectra as a function of the r-wave number. The r-spectra for both images capture the entire gradient of the oil palm canopy cover. Maximum values are reached at r-wavenumber = 1, with low spatial frequencies of 8.4 cycles km −1 and 8.3 cycles km −1 for PlanetScope and FORMOSAT-2, respectively ( Figure 4). These values characterize canopy texture of the oil palm plantation. Beyond these estimates, mean r-spectrum values remain observable at wave numbers 2 and 3 (Figure 5a). The latter appear at spatial frequencies greater than the maxima for both images, as shown in Figure 4. These values further characterize the texture of oil palms for FORMOSAT-2 with a 60-pixel window size, but did not fully characterize the texture for PlanetScope with a 38-pixel window size (Figure 4). Minimum r-spectrum values were obtained at r-wavenumber ≥ 4, with high spatial frequencies (≥33 cycles km −1 ) for both images (Figure 4). These values are more consistent with spatial patterns that do not represent the texture of oil palms.
Remote Sens. 2020, 12, x FOR PEER REVIEW 11 of 27 m), but average stem diameter is smaller (48.8 cm compared to 50.8 cm). This class is localized in lowlands near the watercourse and, thus, potentially subject to growing conditions that are indicative of poorly drained soils.  PCA of PlanetScope and FORMOSAT-2 r-spectra demonstrates that factor scores derived from the first three axes are sufficient to explain most variations in oil palm canopy textures in the images (Figure 5b). For PlanetScope, the three components account, respectively, for 61.3%, 21.9% and 5.4% of the variance (90.8% of total variation). For FORMOSAT-2, they account, respectively, for 63.4%, 19.6% and 9.3% of the variance (90.3% of total variation). Ordinations of the first two axes (PC1 and PC2) are presented in Figure 5c,d for FORMOSAT-2 and PlanetScope, respectively. PC1 and PC2 account for just over 80% of total explained variance. The observed dispersions indicate the lack of correlation (orthogonality) between the two components in each case, consistent with expectation. In the case of FORMOSAT-2, four distinct clusters or texture classes emerge in the ordination (Figure 5c Using factor scores for the first two principal component axes (PC1 and PC2), red-green-blue representations of FOTO canopy texture maps were produced for each of the images (Figure 6). The colors are indicative not only of PC factor scores for both axes, but also textural classes. Yellow denotes the open fine-textured oil palm class. Here, scores are positive on both PC1 and PC2 ( Figure 5). This class occupies the entire central portion of the plantation and extends eastwards. It is dominated by with palms with lower-stature crowns, with an average stem height of around 8.8 m and DBH of 50.8 cm. Oil palms in the coarse-closed class are shown in blue ( Figure 6). This class is characterized by negative scores on both PC axes ( Figure 5). It is found in the western and eastern regions of the plantation. In these areas, palms have high average stem heights (about 9.

Estimation of AGB by Multiple Regression
Correlations between FOTO indices and the various structural parameters of the oil palms were first analyzed. Table 2 summarizes correlations (as R 2 values) with factor scores for the three principal components in both image types. For illustrative purposes, Figure 7 displays the results with PC1. The negative relationships that are obtained between PC1 from FORMOSAT-2 and plant structural parameters are the most significant overall (p < 0.0001). In particular, variation in PC1 not only strongly explains BA (R 2 = 0.84), but also that of stem diameter (R 2 = 0.78). Despite being significant, the weakest correlations are obtained with stem height (R 2 = 0.56) and canopy density (R 2 = 0.54). With PlanetScope, the results follow the same trends, but with weaker correlations in general. Basal area and DBH are the parameters that are best explained by PC1 in this case, with R 2 values of 0.68 and 0.64, respectively (p ≤ 0.0001). Stem height and density are the least explained with PC1 from PlanetScope (R 2 < 0.55). Other than PC1, the other components (PC2, PC3) generally show insignificant relationships with the structural parameters of oil palms, regardless of the satellite data considered ( Table 2).  Structural parameters of oil palms are explanatory variables of the aboveground biomass. A relationship between these parameters and FOTO indices that are from FORMOSAT-2 and PlanetScope, in particular PC1 (Table 2), suggests links between these indices and AGB. Based upon this preliminary analysis, MLR models of AGB were developed by considering the three components (PC1, PC2 and PC3) for each of the two image types. In the case of FORMOSAT-2, the best relationship that was obtained (Model 1) integrates all three indices (Table 3). It is a strong and significant relationship (R 2 = 0.86, p < 0.0001). The associated RMSE was estimated to be 4.03 tons ha −1 . The best MLR model constructed considering PlanetScope FOTO indices (Model 2) uses also the three components (PC1, PC2 and PC3). This model offers comparable metrics to that of FORMOSAT-2 (R 2 = 0.83; p < 0.0001, RMSE = 3.77 tons ha −1 , as indicated in Table 3. Structural parameters of oil palms are explanatory variables of the aboveground biomass. A relationship between these parameters and FOTO indices that are from FORMOSAT-2 and PlanetScope, in particular PC1 (Table 2), suggests links between these indices and AGB. Based upon this preliminary analysis, MLR models of AGB were developed by considering the three components (PC1, PC2 and PC3) for each of the two image types. In the case of FORMOSAT-2, the best relationship that was obtained (Model 1) integrates all three indices (Table 3). It is a strong and significant relationship (R 2 = 0.86, p < 0.0001). The associated RMSE was estimated to be 4.03 tons ha −1 . The best MLR model constructed considering PlanetScope FOTO indices (Model 2) uses also the three components (PC1, PC2 and PC3). This model offers comparable metrics to that of FORMOSAT-2 (R 2 = 0.83; p < 0.0001, RMSE = 3.77 tons ha −1 , as indicated in Table 3. Table 3. Multiple linear regression models developed for estimating aboveground oil palm biomass using indices estimated by textural ordinations on Fourier transforms (FOTO), which were produced from PlanetScope and FORMOSAT-2 images. The constants a, b 1 , b 2 , and b 3 are coefficients of the respective multiple linear regressions (MLR) equations. RMSE (tons ha −1 ) is the Root-Mean-Square-Error and %RMSE is its percentage. p-values are calculated at a significance level of 5%.

Biomass Estimation Using MARS
AGB models were developed with the MARS approach using FOTO indices of the FORMOSAT-2 and PlanetScope images, according to methods that are detailed in Section IV. Several combinations of indices PC1 to PC3 could be used to generate interesting equations. Table 4 summarizes the models that we have developed, and retained for validation purposes and to select the best model for each image type. All three models that were produced from FORMOSAT-2 data have very satisfactory performance criteria (R 2 ≥ 0.89, p < 0.0001; RMSE < 3.6 tons ha −1 ; %RMSE < 5.9%). Among them, Models 3 and 4 incorporate interactions between PC1 and PC3, thereby making additional adjustments compared to Model 5, which rely solely upon basic linear functions. The integration of interaction terms reduces slightly the estimation errors (less than 1%), but increases the degrees-of-freedom that are associated with the equations. Table 4. Aboveground biomass models estimated using multivariate adaptive spline regressions (MARS). Three models were also developed using FOTO indices derived from PlanetScope. Models 6 and 7 incorporate interaction terms, while Model 8 sums piecewise linear functions of the form β i (PC i − k i ), where PC i corresponds to FOTO indices PC1, PC2 or PC3, and β i and k i are constants. As was the case with FORMOSAT-2, models taking interactions into account yield slightly smaller errors, but the differences are non-significant. Overall, the three models developed for PlanetScope appear as strong as those developed with FORMOSAT-2 data, with errors of the same order of magnitude (R 2 ≥ 0.84, p < 0.0001; RMSE < 3.8 tons ha −1 ; %RMSE < 6%).

Validation of Estimates
AGB models that were developed using MLR and MARS approaches, based on the FOTO indices, were validated from independent samples that were not considered during the training or development phase (see Section 4.1). Table 5 summarizes results that were obtained. The validation exercise confirms strong relationships between biomass and FOTO indices for both PlanetScope and FORMOSAT-2, as demonstrated by the summary statistics (R 2 ≥ 0.80; p ≤ 0.0005; RMSE < 5.5 t ha −1 ; %RMSE < 10.2%), regardless of the approach (MLR or MARS) used. In this study, MLR appears to be just as effective as the MARS method, especially when FORMOSAT-2 data are used. Overall, the best validation results are obtained with FORMOSAT-2 data (R 2 ≥ 0.89; p < 0.0001), especially for Models 3, 4 and 5, which yield lowest errors (%RMSE ≤ 6.40%). Model 1 (based on MLR) gives quite similar validation results. Slight elevated errors were obtained with PlanetScope MARS models at the validation stage (on average about 3% above FORMOSAT-2 models in terms of %RMSE). Examination of the various results obtained shows that proposed approaches, when combined with the use of FOTO texture indices, are effective in estimating oil palm aboveground biomass. The use of PlanetScope or FORMOSAT-2 images does not introduce important differences. Of all proposed models that were based upon MARS, Model 3 for FORMOSAT-2 and Model 6 for PlanetScope appear to be the best models. These two models have errors of 6.09% and 8.51%, respectively, on the estimated AGB (Table 5). Model 1 that is based upon the MLR approach with FORMOSAT-2 is just as efficient, with a similar error level of 6.56%.
Following validation, an AGB map was generated with Model 3 using FORMOSAT-2 data for illustrative purposes (Figure 8). Observable spatial variability in the map (Figure 8b) reflects variations in oil palm structure that were observed in the field during sampling. It is also similar to variability in oil palm texture (Figure 6a). Thus, highest AGB (>59.00 t ha −1 ) is found in dry land areas (dark green and medium dark green, Figure 8b). These areas benefit from regular maintenance of the oil palms by OPG. On one hand, they are characterized by an open or closed coarse texture (Classes 2 and 4, Figure 6b), i.e., they are dominated by mature oil palms with large tops. On the other hand, lowest AGB (<59.00 t ha −1 ) is in the poorly drained lowland area surrounding the stream crossing the eastern section of the plantation (Figure 8b, light green and medium light green). Pockets of continuous stagnant water are found here; palms do not benefit from the same silvicultural maintenance. Consequently, these oil palms exhibit lower growth and appear less mature. They are characterized by smaller crowns, as observed during fieldwork. The corresponding textures in the images are rather thin open or closed in these areas (classes 1 and 3, Figure 6b).

Potential of FOTO Indices for Oil Palms
FOTO indices that were obtained using factor scores for the first three axes of the principal component analysis allowed us to capture coarse and fine textures of oil palm crowns. They also made it possible to characterize the degree of canopy openness and closure. Various studies have also shown the ability of FOTO indices to describe the textures of the canopy of forests and oil palms [6,8,14]. The homogeneity and the simple structures of the mature oil palm stand in an agro-industrial plantation context certainly favored the performance of FOTO indices. According to its original applications, the FOTO method is best suited for treating homogeneous stands [18,48]. Overall, the palms that were considered were mature (>30-years-old) and the differences that were observed between them could be mainly attributed to their location, i.e., either in uplands or in low-lying areas, which conferred greater homogeneity to the images. Yet, FORMOSAT-2 data captured texture information better than PlanetScope did. These data resulted from a pansharpening process and showed greater contrast with a fine spatial resolution (2 m × 2 m), thereby favoring a sharper delineation of palm crowns. In addition, radiometric uncertainty of FORMOSAT-2 was less than 5% [49], which helped to reduce noise in the data. To minimize information losses and noise, no radiometric correction was applied to the FORMOSAT-2 image, as recommended in previous work employing the FOTO method [8]. A low noise level favors full coverage of the Nyquist frequency by the spatial frequencies of the r-spectrum over the very first wavenumbers, as was the case in this study with the FORMOSAT-2 image ( Figure 5).
FOTO indices that were derived from PlanetScope also captured variation in the texture of oil palms, but in a less clear-cut manner. By comparing the Nyquist frequencies that were obtained, we observe that spatial frequencies of the PlanetScope r-spectrum contain more noise than those of FORMOSAT-2, especially for the first three wavenumbers. Several factors could explain these differences. The spatial resolution of PlanetScope (3.125 m × 3.125 m), although very good, is lower

Potential of FOTO Indices for Oil Palms
FOTO indices that were obtained using factor scores for the first three axes of the principal component analysis allowed us to capture coarse and fine textures of oil palm crowns. They also made it possible to characterize the degree of canopy openness and closure. Various studies have also shown the ability of FOTO indices to describe the textures of the canopy of forests and oil palms [6,8,14]. The homogeneity and the simple structures of the mature oil palm stand in an agro-industrial plantation context certainly favored the performance of FOTO indices. According to its original applications, the FOTO method is best suited for treating homogeneous stands [18,48]. Overall, the palms that were considered were mature (>30-years-old) and the differences that were observed between them could be mainly attributed to their location, i.e., either in uplands or in low-lying areas, which conferred greater homogeneity to the images. Yet, FORMOSAT-2 data captured texture information better than PlanetScope did. These data resulted from a pansharpening process and showed greater contrast with a fine spatial resolution (2 m × 2 m), thereby favoring a sharper delineation of palm crowns. In addition, radiometric uncertainty of FORMOSAT-2 was less than 5% [49], which helped to reduce noise in the data. To minimize information losses and noise, no radiometric correction was applied to the FORMOSAT-2 image, as recommended in previous work employing the FOTO method [8]. A low noise level favors full coverage of the Nyquist frequency by the spatial frequencies of the r-spectrum over the very first wavenumbers, as was the case in this study with the FORMOSAT-2 image ( Figure 5). FOTO indices that were derived from PlanetScope also captured variation in the texture of oil palms, but in a less clear-cut manner. By comparing the Nyquist frequencies that were obtained, we observe that spatial frequencies of the PlanetScope r-spectrum contain more noise than those of FORMOSAT-2, especially for the first three wavenumbers. Several factors could explain these differences. The spatial resolution of PlanetScope (3.125 m × 3.125 m), although very good, is lower than that of the FORMOSAT-2 product used here. This may possibly affect the discrimination of palm crowns. The PlanetScope image had already undergone radiometric corrections, changing from grey levels to radiances at the top of the atmosphere. Further, radiometric uncertainty of PlanetScope is on the order of 5-6% [50], which is higher than that of FORMOSAT-2. The combination of these different factors can promote loss of information and limit the performance of PlanetScope [50].
Nevertheless, we found that FOTO indices that were derived independently from the two types of images are significantly correlated, except for the third component. Indeed, the correlation between FORMOSAT-2 PC1 and PlanetScope PC1 is r = 0.80, while the correlation between the two PC2 is r = 0.81. PC3, which consisted primarily of noise, is uncorrelated (r = −0.36). These significant correlations indicated that FORMOSAT-2 and PlanetScope captured very similar information in their first two FOTO indices, even though the images were acquired on very different dates (2011 for FORMOSAT-2; 2017 for PlanetScope). This difference in years could initially be viewed as a handicap. However, several studies agree that the growth rate of oil palm in agro-industrial plantations stabilizes after 25 years, and that its height growth gradually slows down towards 30 years [51][52][53]. The stabilization of the growth rate affects several structural parameters (diameter, crown size), which may be important for estimating AGB using FOTO indices [42]. As the study was carried out on a mature plantation, the effect of different dates of image acquisition was not noticed.
Choosing the best window size is a very important step in implementing FOTO. To determine this size in this study, the window size that produced the highest Nyquist frequency was selected for each image. This choice also took into account the decomposition or separation of the principal component scores. This separation was observed in the ordination of factor scores for the two principal components of the r-spectra for each of the images in the study. With this choice, the number of scores that were available in each principal component with each window was also considered sufficient for establishing relationships between the FOTO indices and AGB [8]. These different choices made it possible to facilitate the application of FFT-2D analysis and production of the r-spectra. Some scores in the present study appear very high ( Figure 6). However, as also observed and reported in Couteron et al. [37,40], these high values closely match the fine textured grain, and are not outliers.

Estimation of Aboveground Biomass
The study established different models for estimating AGB based on the MLR approach and the MARS method. Equations that were developed with MLR were very significant overall. The validation process confirmed the robustness of these relationships, especially with FORMOSAT-2 (Model 1, R 2 = 0.85, p = 0.0001, RMSE = 4.18 t ha −1 , %RMSE = 6.56%). With MARS, we developed three models for each type of image. Some of these models used only basic piecewise functions, while others integrated the interactions that occurred between FOTO indices (PC1, PC2 and PC3). Integration of these terms provided additional adjustments that slightly improved model performance. Yet, this improvement frequently was not significant considering the RMSE differences. For the best two models with FORMOSAT-2 data, the difference is <0.2%, while it is less than 0.8% for PlanetScope. With MARS, it is possible to consider several combinations of basic functions. Models 3 to 5 that were developed for FORMOSAT-2 and Models 6 to 8 for PlanetScope were apparently the best compromises in this study. Of these, Model 3 stands out as the most efficient following the validation process, with R 2 = 0.92, p < 0.0001, and %RMSE = 6.09%. For comparison, larger errors, ranging from 28 to 36% were reported in [4] for AGB prediction uncertainties of young oil palms (1 to 5 years), based on models using spectral indices and Normalized Difference Vegetation Index (NDVI). The FOTO PC1 index plays a crucial role in each of the models developed in the present study. This was not very surprising, since PC1 already explains several important structural parameters of the oil palm stand, including basal area (R 2 = 0.68 to 0.84), DBH (R 2 = 0.64 to 0.78), stem height (R 2 = 0.53 to 0.56) and density (R 2 = 0.50 to 0.54) (see Figure 7). These different parameters are critical for biomass estimation. PC2 (and, to a lesser extent, PC3) provided additional information that favors more accurate biomass estimation.
The relatively homogeneous structure of the oil palm plantation is well captured by the FOTO indices, particularly at the FORMOSAT-2 scale, thereby favoring a more accurate estimate of aboveground biomass. In a study in Malaysian Borneo, Singh et al. [14] similarly reported very significant results (R 2 = 0.83) between AGB of young oil palms (<9-year) and FOTO indices that were derived from SPOT 5 imagery acquired at a fine spatial resolution similar to that of FORMOSAT-2 (2.5 m × 2.5 m). The results found in this study are comparable to those reported in several other studies for different vegetation types. Indeed, strong relationships (R 2 = 0.92) have been shown between AGB of mangroves in French Guiana and the FOTO indices (PC1, PC2 and PC3) derived from panchromatic Ikonos image at a spatial resolution of 1 m [8]. In several studies, the FOTO index (PC1) was well correlated to structural parameters of tropical forests [6,14,40,42], as was the case here for oil palms. For example, the study in [40] on tropical forest reported R 2 = 0.57 between FOTO PC1 and stem height, which is similar to the values we found here for oil palms (R 2 = 0.56 with FORMOSAT-2, and 0.53 with PlanetScope). In addition, based on simulations, Barbier et al. [42] reported good relationships between PC1 and forest canopy DBH (R 2 = 0.75), crown size (R 2 = 0.62), but the correlation with density was relatively weak (R 2 = 0.35). Singh et al. [14] found strong relationships between FOTO indices and structural parameters of a mixed forest landscape. The biomass models that were developed in our study with FORMOSAT-2 FOTO indices apparently performed slightly better than those of PlanetScope, as shown during the validation phase ( Table 5). The average difference in terms of %RMSE is in the order of about 3%. Radiometric noise in PlanetScope images (5 to 6%; [50]) was higher than in FORMOSAT-2 (<5%; [49]), as were differences in spatial resolution and the influence of image processing on calculation of the FOTO indices. Aside from the characteristics and the quality inherent to the different types of images, potential measurement errors of oil palm variables in the field may also affect the results found.
Models that were based on the MARS approach are more dependent on the range of data values and their degree of autocorrelation than on pre-processing applied to satellite images [54]. This result underscores the importance of principal component analysis, which reduces correlations in FOTO indices. In all of the results that we obtained, the MARS approach was well suited to estimating AGB in oil palm plantations, both with PlanetScope and FORMOSAT-2 imagery. The performance of MARS is attributable to the method's ability to transform and select important variables, and then to identify non-linear and multidimensional interactions between them [25,33,55]. Therefore, it was not surprising that various previous studies have reported favorable results for the MARS method compared to other statistical methods (logistic regression, MLR, among others) in a range of fields, such as medicine, pedology, hydrology and forestry [23,[31][32][33]. The MARS approach deserves further investigation in remote sensing for AGB estimation, not only of plantations but also of natural forests. Other advanced methods that are based on the concept should also be introduced, including Multivariate Adaptive Regression Spline Differential Evolution (MARS-DE), which has been recently proposed by Al-Sudani et al. [56].
Despite the very significant results that were obtained, the study had some limitations that should be mentioned. The number of sample plots that were considered is relatively small. The use of a sufficient number of sample plots randomly distributed on plantations at various growth stages would improve the robustness of AGB models that are based upon FOTO indices. In this sense, PlanetScope images and other types of fine-resolution optical images (SPOT, PLEADES, SkySat, among others) could be analyzed to better understand the effects of spatial resolution, radiometric noise and pre-processing. The very encouraging results that were obtained here with FOTO indices also permit the use of other texture indices that are derived, for example, from a grey-scale co-occurrence matrix, wavelet decomposition or multiband texture analysis. The addition of auxiliary information, such as topographic, climatic or geomorphological data, could also be considered for improving estimates of AGB [13,57].

Conclusions
This study estimated oil palm AGB with MARS and MLR using PCA scores (components 1, 2 and 3) performed on FORMOSAT-2 and PlanetScope r-spectra. These scores were considered as FOTO textural indices for each image type. The results that were obtained with both approaches are very good overall, especially when FORMOSAT-2 data are considered. Slightly higher errors were obtained with PlanetScope compared to FORMOSAT-2, according to the validation performed. The different models that were developed with MARS, particularly those integrating interaction terms between the FOTO indices (explanatory variables), proved to be more efficient in the study for both type of satellite data. They all integrate the basis functions of the FOTO indices of the three PCA axes and contain their interaction terms. The best model in the study is based on FORMOSAT-2 FOTO indices. It was chosen to produce an AGB map of the study area. Upon further studies, such a model could be used to map AGB of oil palms in the Congo Basin using remote sensing data. Then corresponding sequestrated carbon could be evaluated, with linkage to climate change. The present study focused mainly on a mature oil palm plantation. Subsequent work would have to extend it to larger plantation areas with different stages of palm growth and development, and use a larger number of sample plots. The MARS approach should be further investigated not only with different types of very high-resolution commercial images, but also with free medium-resolution images (e.g., Landsat, Sentinel) for mapping plantation biomass at larger scales.