Estimating and Examining the Sensitivity of Different Vegetation Indices to Fractions of Vegetation Cover at Different Scaling Grids for Early Stage Acacia Plantation Forests Using a Fixed-Wing UAS

Understanding the information on land conditions and especially green vegetation cover is important for monitoring ecosystem dynamics. The fraction of vegetation cover (FVC) is a key variable that can be used to observe vegetation cover trends. Conventionally, satellite data are utilized to compute these variables, although computations in regions such as the tropics can limit the amount of available observation information due to frequent cloud coverage. Unmanned aerial systems (UASs) have become increasingly prominent in recent research and can remotely sense using the same methods as satellites but at a lower altitude. UASs are not limited by clouds and have a much higher resolution. This study utilizes a UAS to determine the emerging trends for FVC estimates at an industrial plantation site in Indonesia, which utilizes fast-growing Acacia trees that can rapidly change the land conditions. First, the UAS was utilized to collect high-resolution RGB imagery and multispectral images for the study area. The data were used to develop general land use/land cover (LULC) information for the site. Multispectral data were converted to various vegetation indices, and within the determined resolution grid (5, 10, 30 and 60 m), the fraction of each LULC type was analyzed for its correlation between the different vegetation indices (Vis). Finally, a simple empirical model was developed to estimate the FVC from the UAS data. The results show the correlation between the FVC (acacias) and different Vis ranging from R2 = 0.66–0.74, 0.76–0.8, 0.84–0.89 and 0.93–0.94 for 5, 10, 30 and 60 m grid resolutions, respectively. This study indicates that UAS-based FVC estimations can be used for observing fast-growing acacia trees at a fine scale resolution, which may assist current restoration programs in Indonesia.


Introduction
Quantitative assessments of the green vegetative covers of terrestrial environments are essential for understanding ecosystem dynamics. The functions of green environments (e.g., vegetation, forests) provide important benefits to ecosystems, such as controlling air quality through photosynthesis, generating an energy supply from woody biomass, preventing soil erosion, improving water quality and balancing the heat fluxes of the earth [1][2][3][4][5]. The worldwide terrestrial environment is currently showing rapid changes in particular regions from anthropogenic activities, causing land degradations that engulf the natural environment [6]. Some areas show transitions away from green areas, which 2 of 18 results in substantial impacts on local to global ecosystems, sociocultural and economic impacts [7,8]. To mitigate these impacts, known cooperative agencies and organizations are planning actions for the recovery of green areas [9,10]. Consistent monitoring of the rapid environmental changes in green vegetative coverages is important for conservation and maintaining the sustainability of the natural environment.
Indonesia's landmass includes approximately 24 million hectares (Mha) of peatlands, which represents approximately 83% of the peatlands found in Southeast Asia. Peatlands in Indonesia are distributed mainly among the four large islands of Sumatera (9.2 Mha), Kalimantan (4.8 Mha), Sulawesi (0.06 Mha) and Papua (6.6 Mha) [11]. One of the common uses of peatlands in Indonesia is for Industrial Forest Plantations (IFP). The area of IFP concessions in Indonesia, which are located on peatlands, is 2 Mha [12], or 54.79% of the total IFP area in 2006, which reached 3.65 Mha [13], and Acacia crassicarpa is a fast-growing species that has been developed as a staple plant for most IFPs on peatlands [14].
Indonesia peatland forests provide important local and global benefits. However, their drainage and conversion into agricultural lands without well-planned management has caused considerable and irreversible environmental, social and economic damage. The catastrophic 2015 fires in Indonesia [15] drew national and international attention. That event reinforced the commitments of the Indonesian government to both reduce peatland deforestation and fires and to rehabilitate and restore degraded peatlands via reforestation. Strategic and operational approaches for monitoring the peat ecosystem together with the conditions of the green vegetation are crucial.
Various researchers and institutions have performed related studies of quantitative analyzations of both local and global vegetation coverage. The products of MODIS Vegetation Continuous Fields or the fCover (fraction of vegetation cover, hereafter denoted as FVC) [16] were used, while other researchers utilized the MODIS reflectance data provided by Land Processes Distributed Active Archive Center (LP DAAC) for developing improved FVC data [17]. The remote sensing techniques for FVC development utilize the multispectral information observed from space and validates its product with ground truth information (e.g., field surveys). The estimation methods can vary depending on the model type used, including simple empirical models [18], linear spectral models [19], decision tree method [20], machine learning techniques [21], and so on. Although the input information is rather simple, remote sensing can use various spectral information or the computed vegetation indices (Vis) for its estimation, although correctly delineating the FVC for various regions of the world is still a challenge.
Many studies have indicated the capability of space-borne remotely sensed data for mapping and/or monitoring of regional to global vegetation cover. However, depending on the products or specific locations used, there can be constraints and challenges in processing or accurately estimating a fractional cover. One of the conventional issues observed is the cloud cover, which blocks otherwise available information for analyzing a terrestrial environment [22]. The missing information can be aided with a gap-filling technique [17] for including continuous land information. However, this technique can result in large differences in the spectral information of the area, and then the possibility of an incorrect FVC estimation rises. In tropical regions including Indonesia, obtaining sufficient land information from areas with lesser cloud cover can be challenging; even when cloud removal and gap filling are performed, correct land information for a certain period of time can be lacking. If a large area of peatlands in Indonesia is observed in IFPs with fast-growing trees, then even a small temporal gap of land information may show erroneous assumptions of the vegetation coverage ( Figure 1). To accurately and effectively detect green areas, considering different platforms or sensors is an important step for effectively monitoring the tropical areas. This is especially true in Indonesia where rapid land transitions are occurring. Figure 1. Example of temporal differences for fast-growing species. The left image shows the fastgrowing Acacia trees in its early stages in August, 2018, while the right image shows its rapid growth in October, 2018. Even with small temporal differences, the situation of the land area would change dramatically.
In recent years, many studies using unmanned aerial systems (UASs) were carried out. The UAS platform provides alternatives to space-borne platforms since optical data can be observed in a clearhigh spatial/temporal resolution for the region of interest [23]. This technique has been used in research on ecology [24], precision agriculture/forestry [25,26] and even analyses for estimating vegetation cover [27][28][29]. UAS was successful [27] in clearly estimating vegetation fractions and flower fractions in crop fields with the changing VIs, and work by Chen et al. [28] showed that utilizing UAS-captured imagery may clearly detect grassy vegetation covers due to its highresolution data. Riihimäkia et al. [29] showed that the UAV-derived information can be aided by satellite-observed information in FVC estimations. As Indonesia is exposed to high and frequent cloud coverage nationwide, obtaining clear satellite information is often difficult. Even if this information is collected, radiometric corrections for both atmospheric and topographic data are mandatory, which is a difficult task [30]. Riihimäkia et al. [29] recently showed an approach for estimating the FVC at arctic vegetation using UAS and satellite data through multiple spatial scaling's and different indices. They have expressed that there is a strong correlation between the UAS-based FVC for validation data that can be used to bridge with the satellite data and noted that the sensitivity of VIs was better when using Red-edge or Short Wave Infrared (SWIR) information. The prior study of Riihimäkia et al. [29] shows a relationship analysis between the VIs and FVC that is based on only a single class that classified the area into vegetation/non-vegetation. Depending on regions where heterogeneous land use/land cover (LULC) types are seen, there may be more classes requiring further analysis and how those classes affect the VI response. Minimal research has been conducted in rapid changing environments such as Indonesia for estimating the fractional cover of green vegetation by utilizing UASs, especially in rapidly growing industrial forest plantations (IFP). Higher spatial/temporal resolution imagery may have a high potential to analyze where the changes in green vegetative cover are occurring.
The objective of this study is to develop a method for retrieving the FVC by utilizing UAS and multispectral sensors for the fast-growing Acacia plantation forests in Indonesia. Several VIs are computed using the raw band information to compare the sensitivity of the VIs to FVC, moreover, the result is also compared at different spatial resolutions and with other LULC types. The developed product is compared with the existing method for computing FVC by using satellite imagery, and it examines how the UAS observed product can compensate for conventional space-borne products. This work mainly focus on if UAV-based FVC can be obtained in the forested area, while it is out of the scope at the moment for generalizing the result which could be utilized for global estimations. This study may present advances in UAS research in developing FVC estimations and the possibility of utilizing the platform for collecting ground truth information to bridge airborne sensing with space-borne sensing. Example of temporal differences for fast-growing species. The left image shows the fast-growing Acacia trees in its early stages in August, 2018, while the right image shows its rapid growth in October, 2018. Even with small temporal differences, the situation of the land area would change dramatically.
In recent years, many studies using unmanned aerial systems (UASs) were carried out. The UAS platform provides alternatives to space-borne platforms since optical data can be observed in a clear-high spatial/temporal resolution for the region of interest [23]. This technique has been used in research on ecology [24], precision agriculture/forestry [25,26] and even analyses for estimating vegetation cover [27][28][29]. UAS was successful [27] in clearly estimating vegetation fractions and flower fractions in crop fields with the changing VIs, and work by Chen et al. [28] showed that utilizing UAS-captured imagery may clearly detect grassy vegetation covers due to its high-resolution data. Riihimäkia et al. [29] showed that the UAV-derived information can be aided by satellite-observed information in FVC estimations. As Indonesia is exposed to high and frequent cloud coverage nationwide, obtaining clear satellite information is often difficult. Even if this information is collected, radiometric corrections for both atmospheric and topographic data are mandatory, which is a difficult task [30]. Riihimäkia et al. [29] recently showed an approach for estimating the FVC at arctic vegetation using UAS and satellite data through multiple spatial scaling's and different indices. They have expressed that there is a strong correlation between the UAS-based FVC for validation data that can be used to bridge with the satellite data and noted that the sensitivity of VIs was better when using Red-edge or Short Wave Infrared (SWIR) information. The prior study of Riihimäkia et al. [29] shows a relationship analysis between the VIs and FVC that is based on only a single class that classified the area into vegetation/non-vegetation. Depending on regions where heterogeneous land use/land cover (LULC) types are seen, there may be more classes requiring further analysis and how those classes affect the VI response. Minimal research has been conducted in rapid changing environments such as Indonesia for estimating the fractional cover of green vegetation by utilizing UASs, especially in rapidly growing industrial forest plantations (IFP). Higher spatial/temporal resolution imagery may have a high potential to analyze where the changes in green vegetative cover are occurring.
The objective of this study is to develop a method for retrieving the FVC by utilizing UAS and multispectral sensors for the fast-growing Acacia plantation forests in Indonesia. Several VIs are computed using the raw band information to compare the sensitivity of the VIs to FVC, moreover, the result is also compared at different spatial resolutions and with other LULC types. The developed product is compared with the existing method for computing FVC by using satellite imagery, and it examines how the UAS observed product can compensate for conventional space-borne products. This work mainly focus on if UAV-based FVC can be obtained in the forested area, while it is out of the scope at the moment for generalizing the result which could be utilized for global estimations. This study may present advances in UAS research in developing FVC estimations and the possibility of utilizing the platform for collecting ground truth information to bridge airborne sensing with space-borne sensing.

Study Area
The study site is located in West Kalimantan and is the same area used in Iizuka et al. [31] (Figure 2). The large plantation area is managed by an industrial plantation company. The area was planted in January 2017 with Acacia crassicarpa as the main commercial species, which is one of the fast-growing species that can grow from saplings to up to a few meters in one year. Usually, the plantation site has a cycle of planting to logging in four-year intervals, which is a dramatic change rate. One section of the compartment site is considered for the test. Brief details of the area are explained in Iizuka et al. [31].

Study Area
The study site is located in West Kalimantan and is the same area used in Iizuka et al. [31] ( Figure  2). The large plantation area is managed by an industrial plantation company. The area was planted in January 2017 with Acacia crassicarpa as the main commercial species, which is one of the fastgrowing species that can grow from saplings to up to a few meters in one year. Usually, the plantation site has a cycle of planting to logging in four-year intervals, which is a dramatic change rate. One section of the compartment site is considered for the test. Brief details of the area are explained in Iizuka et al. [31].

Materials and Methods
The overall workflow is shown in Figure 3. The data from the UAS were collected for RGB imagery and multispectral data. These data were used to produce a LULC map of the study area and to compute the fractions of the LULC classes. The multispectral data were used to further compute multiple VIs. Different scaled polygons were generated and within each grid, the fractional coverage and VI values were analyzed for their relationships. The data were further divided into training and validation sets, and an estimation of FVC was implemented for different spatial resolutions.

Materials and Methods
The overall workflow is shown in Figure 3. The data from the UAS were collected for RGB imagery and multispectral data. These data were used to produce a LULC map of the study area and to compute the fractions of the LULC classes. The multispectral data were used to further compute multiple VIs. Different scaled polygons were generated and within each grid, the fractional coverage and VI values were analyzed for their relationships. The data were further divided into training and validation sets, and an estimation of FVC was implemented for different spatial resolutions.

Study Area
The study site is located in West Kalimantan and is the same area used in Iizuka et al. [31] ( Figure  2). The large plantation area is managed by an industrial plantation company. The area was planted in January 2017 with Acacia crassicarpa as the main commercial species, which is one of the fastgrowing species that can grow from saplings to up to a few meters in one year. Usually, the plantation site has a cycle of planting to logging in four-year intervals, which is a dramatic change rate. One section of the compartment site is considered for the test. Brief details of the area are explained in Iizuka et al. [31].

Materials and Methods
The overall workflow is shown in Figure 3. The data from the UAS were collected for RGB imagery and multispectral data. These data were used to produce a LULC map of the study area and to compute the fractions of the LULC classes. The multispectral data were used to further compute multiple VIs. Different scaled polygons were generated and within each grid, the fractional coverage and VI values were analyzed for their relationships. The data were further divided into training and validation sets, and an estimation of FVC was implemented for different spatial resolutions.

GCP Collection Using Low Cost GNSS
For improving the 3D modeling and georeferencing, ground control points (GCPs) were placed and coordinates were collected before the flight campaign. The Reach (Emlid, Hong Kong, China) global navigation satellite system (GNSS) was utilized for collecting the XYZ geographical coordinates (Figure 4). A total of two stations were used. One was set on the tip of a long envy pipe and was fixed at 4 m high off the ground so the surrounding obstacles would not block the GNSS signals; this station was used as a base station. The other receiver was placed on a tripod and utilized as a rover station during the collection of coordinates at each GCP target. A total of six GCPs were observed and the GNSS signals at each GCP were recorded for 5 min. The signals were further processed with a post-processing kinematic (PPK) method to enhance the XYZ coordinate precision. Depending on the satellite signals and processing, the coordinate data can enhance the precision using the PPK method up to a centimeter-level error; when using non-processed GNSS data, the precision can result in a 10 m error (when only using the GPS [32]). The Reach GNSS was set to observe GNSS signals from GPS, QZSS, Galileo and Beidou at a logging frequency of 5 Hz.

GCP Collection Using Low Cost GNSS
For improving the 3D modeling and georeferencing, ground control points (GCPs) were placed and coordinates were collected before the flight campaign. The Reach (Emlid, Hong Kong, Hong Kong) global navigation satellite system (GNSS) was utilized for collecting the XYZ geographical coordinates (Figure 4). A total of two stations were used. One was set on the tip of a long envy pipe and was fixed at 4 m high off the ground so the surrounding obstacles would not block the GNSS signals; this station was used as a base station. The other receiver was placed on a tripod and utilized as a rover station during the collection of coordinates at each GCP target. A total of six GCPs were observed and the GNSS signals at each GCP were recorded for 5 min. The signals were further processed with a post-processing kinematic (PPK) method to enhance the XYZ coordinate precision. Depending on the satellite signals and processing, the coordinate data can enhance the precision using the PPK method up to a centimeter-level error; when using non-processed GNSS data, the precision can result in a 10 m error (when only using the GPS [32]). The Reach GNSS was set to observe GNSS signals from GPS, QZSS, Galileo and Beidou at a logging frequency of 5 Hz.

PPK Processing for Precise GCP Coordinates
For the postprocessing of the GNSS signals, the open source RTKLIB software [33] was used. RTKLIB is a program package for standard and precise positioning with GNSS and can perform a PPK analysis by using Receiver Independent Exchange Format (RINEX) files. First, using the RTKLIB, the log files for both the base and rover stations were opened and examined for the GNSS signal quality. The signals were checked for each satellite, and the only satellites that seemed to receive continuous L1 frequency signals during the recording were included for further processing, and the other satellites were omitted ( Figure 5). Each observation of the GCP targets was postprocessed with the "static" option in the positioning mode, "combined" in the filter type, "Fixed and Held" for the integer ambiguity resolution and 15 degrees for the elevation mask, and 35 dB for SNR (signal-tonoise ratio) filtering were used. After the coordinate of each GCP was computed, the coordinate logs were averaged from the fixed solution (Q1) and the ratio factor of ambiguity validation ≧500. However, if there were no logs for Q1 then a float solution (Q2) was used.

PPK Processing for Precise GCP Coordinates
For the postprocessing of the GNSS signals, the open source RTKLIB software [33] was used. RTKLIB is a program package for standard and precise positioning with GNSS and can perform a PPK analysis by using Receiver Independent Exchange Format (RINEX) files. First, using the RTKLIB, the log files for both the base and rover stations were opened and examined for the GNSS signal quality. The signals were checked for each satellite, and the only satellites that seemed to receive continuous L1 frequency signals during the recording were included for further processing, and the other satellites were omitted ( Figure 5). Each observation of the GCP targets was postprocessed with the "static" option in the positioning mode, "combined" in the filter type, "Fixed and Held" for the integer ambiguity resolution and 15 degrees for the elevation mask, and 35 dB for SNR (signal-to-noise ratio) filtering were used. After the coordinate of each GCP was computed, the coordinate logs were averaged from the fixed solution (Q1) and the ratio factor of ambiguity validation 500. However, if there were no logs for Q1 then a float solution (Q2) was used.

Fixed-Wing RGB Acquisition and Processing
On 17 October, 2018, a flight campaign was conducted at one of the compartment areas in the plantation site. A Firefly6 Pro (BirdsEyeView Aerobotics, New Hampshire, US) fixed-wing VTOL (vertical takeoff and landing) UAS was utilized to collect the aerial photography of the area. For the photo shoot, the SONY α6000 camera was embedded in the gimbal attached beneath the fixed-wing UAS and collected throughout the flight path. The flight path of the UAS was set to a south-north direction with the camera time-lapse set to two seconds, corresponding to a forward (side) overlap of 75% (70%), and flying at 140 m at altitude with an approximate 3 cm ground sampling distance (GSD). The flight course did not consider the wind direction. The collected aerial images were further processed with the software Photoscan Version 1.4.5 (Agisoft, St. Petersburg, Russia) for 3D modeling and mosaicking of the whole study area. The Photoscan parameter was considered with a "High" alignment with default tie and key points, and a dense point cloud of "High" with "Mild" depth filtering. A digital surface model (DSM) was generated using the dense-point cloud and an orthoimage was computed using the DSM. The first routine was processed without adding the GCP information. After the first orthoimage was generated, the GCPs were additionally added and the model was reconstructed by adjusting the point clouds with the GCP information. Furthermore, DSM and the orthoimage was reconstructed based on the GCP-provided point clouds. This processing method makes placing manual GCPs throughout the software easier, especially when a large quantity of images is used.

Multispectral Data Collection and Processing
On 17 October, 2018, the flight campaign was conducted at the compartment area in the plantation site, which is the same as the corresponding area for the RGB acquisition. Right after the RGB collection flight, another flight was conducted to collect multispectral images of the area. The SlantRange 3p (SlantRange, California, US) sensor was embedded in the gimbal of the VTOL UAS and images were collected throughout the flight path. The path was the same as the RGB flight. However, the forward (side) overlap was set to 40% (40%), and the flying altitude as 150 m, corresponding to an approximate 6 cm GSD. The SlantRange 3p can collect four bands, namely green, red, red-edge and near infrared. The ambient illumination sensor is placed on top of the VTOL UAS during the shooting ( Figure 6) and the information is later used for calibrating each image for solar illuminations. The metadata of each image was processed with the SlantView software, which calibrates each image for solar illumination and includes coordinate information. The preprocessed

Fixed-Wing RGB Acquisition and Processing
On 17 October, 2018, a flight campaign was conducted at one of the compartment areas in the plantation site. A Firefly6 Pro (BirdsEyeView Aerobotics, New Hampshire, US) fixed-wing VTOL (vertical takeoff and landing) UAS was utilized to collect the aerial photography of the area. For the photo shoot, the SONY α6000 camera was embedded in the gimbal attached beneath the fixed-wing UAS and collected throughout the flight path. The flight path of the UAS was set to a south-north direction with the camera time-lapse set to two seconds, corresponding to a forward (side) overlap of 75% (70%), and flying at 140 m at altitude with an approximate 3 cm ground sampling distance (GSD). The flight course did not consider the wind direction. The collected aerial images were further processed with the software Photoscan Version 1.4.5 (Agisoft, St. Petersburg, Russia) for 3D modeling and mosaicking of the whole study area. The Photoscan parameter was considered with a "High" alignment with default tie and key points, and a dense point cloud of "High" with "Mild" depth filtering. A digital surface model (DSM) was generated using the dense-point cloud and an orthoimage was computed using the DSM. The first routine was processed without adding the GCP information. After the first orthoimage was generated, the GCPs were additionally added and the model was reconstructed by adjusting the point clouds with the GCP information. Furthermore, DSM and the orthoimage was reconstructed based on the GCP-provided point clouds. This processing method makes placing manual GCPs throughout the software easier, especially when a large quantity of images is used.

Multispectral Data Collection and Processing
On 17 October 2018, the flight campaign was conducted at the compartment area in the plantation site, which is the same as the corresponding area for the RGB acquisition. Right after the RGB collection flight, another flight was conducted to collect multispectral images of the area. The SlantRange 3p (SlantRange, San Diego, CA, USA) sensor was embedded in the gimbal of the VTOL UAS and images were collected throughout the flight path. The path was the same as the RGB flight. However, the forward (side) overlap was set to 40% (40%), and the flying altitude as 150 m, corresponding to an approximate 6 cm GSD. The SlantRange 3p can collect four bands, namely green, red, red-edge and near infrared. The ambient illumination sensor is placed on top of the VTOL UAS during the shooting ( Figure 6) and the information is later used for calibrating each image for solar illuminations. The metadata of each image was processed with the SlantView software, which calibrates each image for solar illumination and includes coordinate information. The preprocessed data were exported and Remote Sens. 2019, 11, 1816 7 of 18 further processed via Photoscan to generate the overall view of the study area. The parameter slightly changes from RGB processing, where the alignment is "Very High", the dense point cloud is "Ultra high," while the other process stays the same. The multispectral data again performed its first routine without GCP and the second process was added with GCP calibration.
Remote Sens. 2019, 11, x FOR PEER REVIEW 7 of 17 data were exported and further processed via Photoscan to generate the overall view of the study area. The parameter slightly changes from RGB processing, where the alignment is "Very High", the dense point cloud is "Ultra high," while the other process stays the same. The multispectral data again performed its first routine without GCP and the second process was added with GCP calibration.

Image Classification
The study area is a small section of the compartment area of the Acacia plantation. However, there is a variety of LULC types in the area, as not only acacias but also water bodies, bare soils and grassy shrubs can be observed. Using the RGB image and multispectral imagery, a conventional supervised image classification was performed to develop a categorical map that shows the distribution of the general LULC types of the site. The Multilayer-Perceptron (MLP) classification [31] was implemented to classify each LULC class, which is a classifier that can handle even nonlinear trends between variables. A total of four classes were classified: Acacia trees, bare soil, water bodies and grass/shrubs. For the validation of the classified map, an accuracy assessment was additionally performed. Using the developed RGB orthoimage and ground information for the study area, 600 evenly distributed sampled points per class were manually collected throughout the scene, and the error matrix was computed. The LULC map was further used to compute the fraction of LULC classes within each gridded location. Other authors indicate that UAS information is highly correlated with traditional ground-based hemispherical photography, which can be used as reference for groundtruth information [34]. Therefore, the UAS-based classification result will be further used likewise as ground-truthing.

Various Vegetation Indices
Three different VIs were computed from the multispectral bands collected using SlantRange. Namely, the Normalized Difference Vegetation Index (NDVI), Green NDVI (GNDVI) and Red-edge NDVI (ReNDVI) were used, where each index can be computed using the following formula: ReNDVI IR Re IR Re , where IR is the infrared band, R is the red band, G is the green band and Re is the red-edge band. The indices all have similar functions. However, they use slightly different characteristics of the electromagnetic spectra. The IR and R are used often and characterize the vegetation function where

Image Classification
The study area is a small section of the compartment area of the Acacia plantation. However, there is a variety of LULC types in the area, as not only acacias but also water bodies, bare soils and grassy shrubs can be observed. Using the RGB image and multispectral imagery, a conventional supervised image classification was performed to develop a categorical map that shows the distribution of the general LULC types of the site. The Multilayer-Perceptron (MLP) classification [31] was implemented to classify each LULC class, which is a classifier that can handle even nonlinear trends between variables. A total of four classes were classified: Acacia trees, bare soil, water bodies and grass/shrubs. For the validation of the classified map, an accuracy assessment was additionally performed. Using the developed RGB orthoimage and ground information for the study area, 600 evenly distributed sampled points per class were manually collected throughout the scene, and the error matrix was computed. The LULC map was further used to compute the fraction of LULC classes within each gridded location. Other authors indicate that UAS information is highly correlated with traditional ground-based hemispherical photography, which can be used as reference for ground-truth information [34]. Therefore, the UAS-based classification result will be further used likewise as ground-truthing.

Various Vegetation Indices
Three different VIs were computed from the multispectral bands collected using SlantRange. Namely, the Normalized Difference Vegetation Index (NDVI), Green NDVI (GNDVI) and Red-edge NDVI (ReNDVI) were used, where each index can be computed using the following formula: where IR is the infrared band, R is the red band, G is the green band and Re is the red-edge band. The indices all have similar functions. However, they use slightly different characteristics of the electromagnetic spectra. The IR and R are used often and characterize the vegetation function where the leaves' chlorophyll absorbs the red spectra and reflects infrared spectra [35]. Using this information one can sense the vegetation activity (greenness) of the area. GNDVI utilizes the green band instead of the red band. The GNDVI is considered an alternative index to NDVI, as it has a wider dynamic range and a higher sensitivity to chlorophyll [36]. ReNDVI also uses the slightly longer wavelength from the red channel instead of R or G, which corresponds to vegetation status information, especially water stress [37]. These three indices were examined with the fractions of the LULC classes to understand the sensitivity of the VIs.

Various Grid Scaling for Relationship Analysis
Four differently scaled square gridded polygons were continuously constructed along the study area with 5, 10, 30 and 60 m grid sizes. The 5 m resolution grid corresponds to satellites such as Rapideye, while the 10 m resolution grid corresponds to Sentinel-2 and the 30 m resolution corresponds to the satellites from the LANDSAT series. The 60 m grid is an additional size used to examine the viewability when the resolution is upscaled to a coarser grid size, which is almost the maximum size that can be determined in the current study area for a relevant analysis. Within each grid, the percentage of LULC classes and the average values of each VI were extracted and examined for the relationship analysis.

Estimating FVC Using UAS Data and Comparison with Satellite-Based FVC
The sample dataset that was used for relationship analysis is divided into sets of training and validation data. For the 5, 10, 30 and 60 m grid data, the samples were randomly divided into 50% for training and the remaining 50% for validating the estimated UAS-derived FVC. The model evaluation is considered by computing the R 2 and the root mean square error (RMSE) and the relative RMSE for each grid scale model [38].
Furthermore, Sentinel-2 imagery was downloaded via the webpage and processed through the Sentinel Application Platform (SNAP) software ver. 6.0.5 to compute the satellite-based FVC. The imagery was observed on 14 September, 2018, just one month before the UAS observation. The biophysical processor module within the SNAP software is used for computing biophysical parameters. While the UAS-based FVC uses a simple empirical model, the SNAP FVC is considered via the neural network method [39]. The satellite-based FVC is later compared with the result from the UAS-based FVC.

PPK Processing for Coordinate Precision
Using the base station and rover GNSS log collected through a small GNSS unit, postprocessing of the GNSS signals was performed using the RTKLIB. A total of six points were processed. Most of the targets are showing fixed positioning of the GNSS equipment (<1 cm precision in XYZ direction), while some GCPs give slightly larger errors when using the float solution. However, even though not all observations give a clear fixed solution, the PPK method can nonetheless give much higher precision than just a single GNSS observation, even with an extremely lightweight, small unit GNSS. The geolocation (XYZ coordinates) computed from the PPK process was used for each GCP target location.

Orthoimagery from UAS Flights for RGB and Multispectral Data
From the flight using the VTOL UAS, totals of 567 and 1204 images were obtained for RGB and multispectral data, respectively. Figure 7 shows the same location of the mosaicked images for the RGB and multispectral data. Due to the irrelevant overlap from the multispectral sensor, there are some gaps in the southern part of the compartment site. However, the procedure succeeded overall in developing an orthorectified image. Gap areas are omitted from further analysis. some gaps in the southern part of the compartment site. However, the procedure succeeded overall in developing an orthorectified image. Gap areas are omitted from further analysis.

LULC Map of the Study Area and Its Errors
Utilizing the RGB imagery and multispectral data, image classification was performed to develop a categorical map. A total of four classes were generated that consisted of the general landscape of the plantation site. Figure 8 shows the result of the mapping. For the validation of the map, an accuracy assessment was performed, and the error matrix is computed in Table 1. The result indicates the overall accuracy of the categorical map as 90.07%. The user accuracies (Producers accuracy) for Acacia forest, bare land, water bodies and grass/shrubs were 83% (96%), 96% (83%), 91% (95%) and 94% (89%), respectively.

LULC Map of the Study Area and Its Errors
Utilizing the RGB imagery and multispectral data, image classification was performed to develop a categorical map. A total of four classes were generated that consisted of the general landscape of the plantation site. Figure 8 shows the result of the mapping. For the validation of the map, an accuracy assessment was performed, and the error matrix is computed in Table 1. The result indicates the overall accuracy of the categorical map as 90.07%. The user accuracies (Producers accuracy) for Acacia forest, bare land, water bodies and grass/shrubs were 83% (96%), 96% (83%), 91% (95%) and 94% (89%), respectively.

LULC Map of the Study Area and Its Errors
Utilizing the RGB imagery and multispectral data, image classification was performed to develop a categorical map. A total of four classes were generated that consisted of the general landscape of the plantation site. Figure 8 shows the result of the mapping. For the validation of the map, an accuracy assessment was performed, and the error matrix is computed in Table 1. The result indicates the overall accuracy of the categorical map as 90.07%. The user accuracies (Producers accuracy) for Acacia forest, bare land, water bodies and grass/shrubs were 83% (96%), 96% (83%), 91% (95%) and 94% (89%), respectively.    Figure 9 shows the relationship analysis between the VIs and fractions of LULC covers for different scaling resolutions. Each VI is examined for three different LULC types: Acacia, grass/shrubs and non-vegetation (i.e., bare soil + water bodies). Each degree of correlation can be seen within the figure.
In general, for the Acacia class, a clear exponential relationship occurs between the FVC and VIs. At the 5 and 10 m grid scales, ReNDVI shows the highest correlation, with GNDVI and NDVI showing lower R2 values. With the resolution becoming coarser, the sensitivity of the VIs to FVC shows a lesser error among all VIs, and for the 60 m grid, all VIs show a strong correlation (R 2 > 0.9). The grass/shrub class, which is classified as another vegetation type, does not show a clearer trend than the Acacia trees. The trend line for the grass/shrub type is shown with a second order polynomial for reference. The polynomial appears to not show many relevant relationships, but it does seem to show a combined trend of the area that has different land conditions. One example is a segment that shows a positive correlation between FVC and VIs, which is shown at areas such as on the west side of the study site, along with the track of the bare soils and water body that crosses from north to south direction. This area shows a simple positive correlation between vegetation and the VIs. Another segment is found within the compartment of the plantation area that shows a negative correlation between the FVC and VIs. In such areas, larger fractions of acacias are found when the grass/shrub class reduces, resulting in higher VIs. Other than the vegetative classes, non-vegetative areas were also examined for the sensitivity of VIs. Interestingly, at all grid scales, a clear negative correlation is observed between the fraction of non-vegetation cover (denoted hereafter as FnVC) and VIs. Similar to the Acacia class, the errors obviously reduce with coarser grids. However, instead of ReNDVI being more sensitive for acacias, the NDVI shows a higher correlation throughout different VIs and at different grid scales for FnVC. Table 1. Accuracy assessment and the error matrix of the generated LULC map. Error C is the error of commission, and Error O is the error of omission. The overall accuracy of the categorical map results is 90.7%.

Acacia
Bare Soil Water Body  Figure 9 shows the relationship analysis between the VIs and fractions of LULC covers for different scaling resolutions. Each VI is examined for three different LULC types: Acacia, grass/shrubs and non-vegetation (i.e., bare soil + water bodies). Each degree of correlation can be seen within the figure. In general, for the Acacia class, a clear exponential relationship occurs between the FVC and VIs. At the 5 and 10 m grid scales, ReNDVI shows the highest correlation, with GNDVI and NDVI showing lower R2 values. With the resolution becoming coarser, the sensitivity of the VIs to FVC shows a lesser error among all VIs, and for the 60 m grid, all VIs show a strong correlation (R2 > 0.9). The grass/shrub class, which is classified as another vegetation type, does not show a clearer trend than the Acacia trees. The trend line for the grass/shrub type is shown with a second order polynomial for reference. The polynomial appears to not show many relevant relationships, but it does seem to show a combined trend of the area that has different land conditions. One example is a segment that shows a positive correlation between FVC and VIs, which is shown at areas such as on the west side of the study site, along with the track of the bare soils and water body that crosses from north to south direction. This area shows a simple positive correlation between vegetation and the VIs. Another segment is found within the compartment of the plantation area that shows a negative correlation between the FVC and VIs. In such areas, larger fractions of acacias are found when the grass/shrub class reduces, resulting in higher VIs. Other than the vegetative classes, non-vegetative areas were also examined for the sensitivity of VIs. Interestingly, at all grid scales, a clear negative correlation is observed between the fraction of non-vegetation cover (denoted hereafter as FnVC) and VIs. Similar to the Acacia class, the errors obviously reduce with coarser grids. However, instead of ReNDVI being more sensitive for acacias, the NDVI shows a higher correlation throughout different VIs and at different grid scales for FnVC.   respectively. The visual interpretation shows that the higher resolution image gives more detail of the area. However, the 30 m resolution data can characterize the overall trend of the site, as the bare soils on the west side can be seen clearly, and the difference between the north compartment with the lower FVC and the southern compartment with the higher FVC is acknowledged. The 60 m resolution can also characterize the area to some degree, however within the spatial extent of the study area, the information is much more aggregated compared to the other, finer resolutions.

FVC Estimation and Satellite-Derived FVC vs. UAS-Derived FVC
The FVC computed from Sentinel-2 is shown in Figure 10e. The evaluation of the Sentinel-2 FVC is compared with the UAS FVC (Figure 10b). From the visual interpretation, the Sentinel-2 FVC shows a more smoothened result than the sharper UAS FVC. Since the Sentinel-2 FVC uses the neural network method for all different bands, the resampled information of the bands appears in the image (RGB is 10 m, however IR and SWIR are 20-30 m and the aerosol/water vapor bands are 60 m). The model validation shows overestimations for lower FVCs and underestimations for higher FVCs for the Sentinel-2 data. Two potential explanations can be made for the FVC difference. One is the error from the neural network approach, and the other is the temporal difference between the satellite and UAS. This subject will be discussed in the next section. respectively. The visual interpretation shows that the higher resolution image gives more detail of the area. However, the 30 m resolution data can characterize the overall trend of the site, as the bare soils on the west side can be seen clearly, and the difference between the north compartment with the lower FVC and the southern compartment with the higher FVC is acknowledged. The 60 m resolution can also characterize the area to some degree, however within the spatial extent of the study area, the information is much more aggregated compared to the other, finer resolutions.
The FVC computed from Sentinel-2 is shown in Figure 10e. The evaluation of the Sentinel-2 FVC is compared with the UAS FVC (Figure 10b). From the visual interpretation, the Sentinel-2 FVC shows a more smoothened result than the sharper UAS FVC. Since the Sentinel-2 FVC uses the neural network method for all different bands, the resampled information of the bands appears in the image (RGB is 10 m, however IR and SWIR are 20-30 m and the aerosol/water vapor bands are 60 m). The model validation shows overestimations for lower FVCs and underestimations for higher FVCs for the Sentinel-2 data. Two potential explanations can be made for the FVC difference. One is the error from the neural network approach, and the other is the temporal difference between the satellite and UAS. This subject will be discussed in the next section. respectively. The visual interpretation shows that the higher resolution image gives more detail of the area. However, the 30 m resolution data can characterize the overall trend of the site, as the bare soils on the west side can be seen clearly, and the difference between the north compartment with the lower FVC and the southern compartment with the higher FVC is acknowledged. The 60 m resolution can also characterize the area to some degree, however within the spatial extent of the study area, the information is much more aggregated compared to the other, finer resolutions.
The FVC computed from Sentinel-2 is shown in Figure 10e. The evaluation of the Sentinel-2 FVC is compared with the UAS FVC (Figure 10b). From the visual interpretation, the Sentinel-2 FVC shows a more smoothened result than the sharper UAS FVC. Since the Sentinel-2 FVC uses the neural network method for all different bands, the resampled information of the bands appears in the image (RGB is 10 m, however IR and SWIR are 20-30 m and the aerosol/water vapor bands are 60 m). The model validation shows overestimations for lower FVCs and underestimations for higher FVCs for the Sentinel-2 data. Two potential explanations can be made for the FVC difference. One is the error from the neural network approach, and the other is the temporal difference between the satellite and UAS. This subject will be discussed in the next section.

Discussion
We will first discuss the trends of FVC to the sensitivity of VIs and the resolution. The resolution (grid scale) shows a straightforward trend, where a coarser resolution shows a stronger correlation. This also coincides with the findings of Riihimaki et al. (2019), which explains that when the finer data is aggregated to a larger unit, the variation in the data decreases, resulting in a higher correlation [29,40]. Starting from the 5 m up to the 60 m resolution, the R 2 improves for the explanatory power and the estimated FVC shows a much lower RMSE for coarser resolutions. This indicates that utilizing UAS data can improve the efficiency and quality of collected ground-truth information for validating coarser imagery FVCs. This approach may outperform expectations by using high-resolution Google Earth scenes [30] or direct observations in the field [28]; this dynamic remains in issues such as the temporal difference of scenes, coarse resolution images that make delineating fractional covers difficult and sampling scale/registration differences between field plots and imagery [29,30]. Since these issues can be controlled for UAS data, new possibilities in bridging to larger-scale earth observations may occur [41].
The sensitivity of VIs to FVC and FnVC was also clear. For the Acacia trees, an exponential increase of FVC occurred with increasing VIs, while for the FnVC, the decreasing trend of FVC with increasing VIs was seen throughout different grid scales. The grass/shrub class showed almost no trends when an overall trend is shown for the whole study area. However, their characteristics can be seen at two different segments, where the pure grass/shrub cover is dominant or if it is mixed within the Acacia plantation area. Many authors have indicated the issues caused by the mixing reflectance from vegetation or the background response to the spectral variability that is observed within the grid [42][43][44]. The background soil reflectance brings errors in observing correct vegetation signatures [42,44], and the mixture of vegetation types can lead to errors due to the spectral difference, especially for woody vegetation and green grass [43]. In our study, Acacia and grass/shrubs were mixed within the compartment area. When the FVC as separated clearly by species, it gave a strong correlation for woody vegetation between the VIs taken from UAS. As expressed by Xiao and Moody [43], the mixture of the vegetation types will result in uncertainties in the estimation. When both the Acacia and grass/shrub classes are aggregated, the R 2 of the exponential model was reduced from a minimum of 0.008 to a maximum of 0.086 depending on the grid scales (for ReNDVI). On the other hand, the correlation between the aggregated FVC (Acacia + Grass/Shrub) and VIs was the strongest for NDVI compared to other the VIs for all grid scales. Estimating single species fractions via satellite data also seems to show stronger correlations when using ReNDVI. However, when considering all vegetation types, NDVI is preferred, although the overall correlation would decrease if all vegetation types were considered. However, depending on different biomes this dynamic may also show trend changes [17]. Therefore, the relationship found here may be suitable mainly for IFP sites within Indonesia for the Acacia species. Notably, estimating the FnVC seems to also be possible, as it has a superior explanatory power compared to directly estimating the vegetation cover. The prospects of estimating the vegetation cover by inversely estimating FnVC look probable, although a large-scale analysis would need to be made for further confirmation.
The comparison between the satellite FVC and UAS FVC gave interesting results. The satellite FVC was based on the neural network approach, which generally tends to show a high accuracy [21,45]. However, from our observation, the lower boundary of the FVC was overestimated and the higher FVC was underestimated. The neural network approach needs to be treated with care during the training phase, as it can produce poor results if poor information is used during the training [46]. Therefore, the neural network approach might also show strong results in other areas, but this was not the case here. Another potential issue is the temporal difference between the satellite data and UAS. The observation period was valid for similar months. However, due to clouds, the data were unusable, which forced us to use only September's data. Although it is only a 1-month difference, the fast-growing acacias may show differences in the FVC result, hence, the validation between the UAS-FVC and satellite FVC could have had shown such a trend. This is the conventional issue of satellite data, where clouds and haze limit data collection on preferable observation dates. We would like to emphasize that this result does not conclude that the FVC estimation via regression is better than the neural network, but the conventional platforms could limit the desired information of the environment at various scenes. The UAS data can overcome this issue and may collect data at any occasion, time and resolution [31,47]. Thus, the UAS FVC may become a potential method for future FVC developments. The potential obstacles posed to future studies would be aviation law and engineering issues [47]. For example, the available flight time per UASs would limit the spatial extent, and the aviation law that limits the flight altitude would also limit the spatial extent. Even collecting land information at a 5 m resolution can provide sufficient green cover information, but expanding the observation area with a limited flight time would be constrained by the flying altitude, as it is mainly regulated by the law. Further investigation will be performed for more precise comparison between the satellite and UAS data in the future by coinciding both data sets.

Conclusions
This work demonstrated utilizing UAS to observe the RGB and multispectral data for analyzing the LULC of the plantation area, and further examined the relationship between the fraction of each LULC class with comparison to different VIs at multiple spatial resolutions. UAS observation was successful for RGB and multispectral data. An LULC map of the area was developed using the orthorectified RGB and multispectral data, and the fraction of each LULC type was extracted within different resolution grids. A comparison of the fractional vegetation cover (FVC) to NDVI, GNDVI and ReNDVI was performed. The result indicates the possibilities of computing the vegetation cover of terrestrial environments through UAS-derived data at IFP in Indonesia. The differences between the FVC and the grid scale and VIs were clear. ReNDVI showed a stronger correlation at a finer resolution for Acacia classes and NDVI was superior for estimating non-vegetation classes. Methods to delineate or monitor the land environment through UAS tended to show promising results and the possibility exists of expanding their potential applications to various fields and approaches. The emerging trend of UAS remote sensing is increasing globally. While airborne and spaceborne sensing are still developing, we hope that the new technology of UAS can integrate with the conventional sensing technology for new findings and developments in the near future. We believe our methods can also be useful for related restoration programs that are currently in progress in Indonesia. We will continue with the approach next to generalize the UAS result and modify the relationships using, for example, radiative transfer models.
Author Contributions: K.I. took the lead in designing, processing the research and writing the manuscript; T.K. arranged for the field surveys, research grants and revising the paper; S.S. and A.Y.S. contributed in field surveying, data collection, data process and revising/editing the manuscript; O.K. contributed in field surveys, research grants and revising of the paper; and all authors have substantially contributed to the work.