Visualizing the Spatiotemporal Trends of Thermal Characteristics in a Peatland Plantation Forest in Indonesia : Pilot Test Using Unmanned Aerial Systems ( UASs )

The high demand for unmanned aerial systems (UASs) reflects the notable impact that these systems have had on the remote sensing field in recent years. Such systems can be used to discover new findings and develop strategic plans in related scientific fields. In this work, a case study is performed to describe a novel approach that uses a UAS with two different sensors and assesses the possibility of monitoring peatland in a small area of a plantation forest in West Kalimantan, Indonesia. First, a multicopter drone with an onboard camera was used to collect aerial images of the study area. The structure from motion (SfM) method was implemented to generate a mosaic image. A digital surface model (DSM) and digital terrain model (DTM) were used to compute a canopy height model (CHM) and explore the vegetation height. Second, a multicopter drone combined with a thermal infrared camera (Zenmuse-XT) was utilized to collect both spatial and temporal thermal data from the study area. The temperature is an important factor that controls the oxidation of tropical peats by microorganisms, root respiration, the soil water content, and so forth. In turn, these processes can alter the greenhouse gas (GHG) flux in the area. Using principal component analysis (PCA), the thermal data were processed to visualize the thermal characteristics of the study site, and the PCA successfully extracted different feature areas. The trends in the thermal information clearly show the differences among land cover types, and the heating and cooling of the peat varies throughout the study area. This study shows the potential for using UAS thermal remote sensing to interpret the characteristics of thermal trends in peatland environments, and the proposed method can be used to guide strategical approaches for monitoring the peatlands in Indonesia.


Introduction
The peatland issues in Indonesia have drawn attention from national and international stakeholders, particularly after the 2015 megafire that burned 0.9 million ha of peatlands, mostly in the southeastern province of Sumatra and southern Kalimantan and Papua [1].It was estimated that 11.3 Tg CO 2 per day was emitted due to the fires in September-October 2015 (i.e., 84% of Indonesia's annual emissions reported in 2002 [1]).Fires and the resulting haze have had numerous effects in Indonesia and the neighboring countries, and public health has been severely impacted.Large-scale land clearing for plantations and illegal logging activities are both major causes of the peatland fires [2,3] that produce smoke haze pollution and carbon emissions in Indonesia [4].The El Niño cycles also play a significant role in peatland fire incidents.Moreover, the area is characterized by uncontrolled drainage due to poor water management, and the groundwater table has become deeper, which has led to the drying and loss of peat.This process accelerates the decomposition of peat and enhances the emissions of greenhouse gases (GHGs) [5].
The sustainable management of the Indonesian peatlands has been a challenge.To overcome this challenge, the government, non-government, and private sectors have undertaken high-level efforts to reduce peatland deforestation and fires and restore degraded peatlands.Although the government has generally focused on developing policies and regulations, many researchers and entities in the private sector have focused on improving peatland management.Studies have shown that appropriate water management is a key to achieving sustainable peatland management [6]; however, achieving the correct management plans and operation is difficult, lacking in monitoring full spatial heterogeneity of hydrological behavior in peatlands [7,8].The heterogeneity of peat areas makes evaluations difficult to conduct through ground observations, and the spatial characteristics of these areas are poorly understood.The importance of the peatland is characterized by multiple functions, but one important factor is that it controls the chemical processes in ecosystem respiration [9].When the peat dome is subsidized by the drainage of the ground water (mainly by anthropogenic activity), there are enhancements of peat surface temperature, which, in turn, lead to enhancements of peat oxidation (Hooijer et al. [10]).Jauhiainen et al. [11] showed that the emission of GHGs from peat environments increases with the peat temperature, and several conditions were examined in non-shaded areas and different degrees of shaded areas, showing that higher shaded areas, where the shadows contribute to lowering the soil temperature, are suppressing the GHG emission more than non-shaded areas.Hilasvuori et al. [12] investigated temperature sensitivity to soil decomposition and found that a higher decomposition rate is seen in oxic layers than in anoxic layers, and the rate responds to an increase in temperature.The authors showed not only the importance of temperature but also the layer of the peats, where older-layer peats (higher depth peat) are more sensitive at higher temperatures than surface layers.From such works, we can understand that revealing the thermal condition of the peat is important.Peat ecosystems do not coincide with administrative/concession boundaries and must be studied at the landscape level [13,14].However, not much effort has been made so far for thermal exploration at a spatial scale relating to peat oxidation.In some cases, satellite information has been used for thermal exploration [15][16][17], with a focus mainly on the effects of land use and land cover (LULC) on surface temperatures or on forest fires.Nevertheless, spatial analyses of temperature information related to the oxidation of the peat are rare.Notably, temperature is an important factor that can lower groundwater levels and enhance the oxidative breakdown of tropical peats as the temperature increases [18].Additionally, such factors influence the oxidation of tropical peats by microorganisms, root respiration, the soil water content, and so forth [10,11,18].Thus, dynamic fluctuations in the temperature of tropical peats can have a positive feedback on the GHG flux and accelerate the effects on climate change [19].Therefore, it is important to identify temperature trends; however, the 1-km resolution of the Moderate Resolution Imaging Spectroradiometer (MODIS) and the 100-m resolution of the Land Remote Sensing Satellite (Landsat) combine all the important thermal information into one pixel.Therefore, based on these data, the heterogeneous peat environment cannot be accurately delineated.Understanding the spatial information associated with a peat environment is critical, and the spatial and temporal extents, required resources, and safety of ground observations are limited.Thus, fine-resolution methods should be considered for in-depth analyses of the region.To achieve such a goal, emerging methods that rely on unmanned aerial systems (UASs) could be effective.
The utilization of UASs has been on trend within the last few years.Great attention has been drawn by multirotor UASs, collecting multiple aerial photos and processing them with a computer vision method called structure from motion (SfM) [20].SfM is an effective method constructing a three-dimensional (3D) model based on multiple images, which can construct ortho imagery through a photogrammetry process to view or analyze a site for various purposes including forestry [20], topography [21], agriculture [22], and so forth.UASs can effectively improve data collection, especially in locations where accessibility is limited or frequent cloud cover limits the data availability of satellite images [23].Moreover, observation flights can be conducted on any occasion, which lets us improve the temporal frequency and spatial resolution of the data [24].The UAS approach could shed light on how to conduct continuous monitoring over a site.Moreover, the technology is now being developed to use thermal data.Recent trends are suggesting the possibility of combining thermal sensors with UASs to monitor various objects.Currently, the common cases are inspections of solar panels for detecting defects and failures on photovoltaics modules [25], facilitating detection and maintenance of difficult-to-observe sites.Another is where thermal sensors are used to observe natural animals for wildlife monitoring and conservation [26], because the thermal information for animals can clearly be detected.Although a thermal UAS has not yet been utilized as an onboard-type of technology, Hare et al. [27] have evaluated the spatial patterning of groundwater discharge by utilizing a high-resolution forward-looking thermal infrared (TIR) camera.The uses of thermal sensors are progressing; however, only a small number of cases have been reported to show their potential.Nevertheless, even with the emerging trends of UASs and photogrammetry, the mapping of spatial extent with thermal sensors is rarely seen.This is the case especially in the peat environment of Indonesia, which could give important perspectives in visualizing the peat temperature to consider the heterogeneity of the area for the conditions of oxidation.The challenges have been identified in the works by Maes et al. [28], who investigated the optimum processing for thermal imageries through SfM and concluded that improvements can be made in the alignments of the images by considering the camera calibration, air temperature correction, and image position estimation aiding with RGB camera alignment.The instrument of thermal sensor is not considered for direct use with UASs; therefore, improvements can be made by utilizing the optimum sensors for such purposes.By achieving a complete mapping of thermal images and applying them for spatial analysis, new perspectives could be achieved that cannot be delineated from only point observations.
In this first pilot examination, the possibility of using UAS multicopter drones and thermal infrared cameras to map the spatiotemporal characteristics of the peatland plantation area is explored to potentially identify significant trends that can be interpreted through aerial remote sensing.

Study Area
The experimental site is located in West Kalimantan, and the site is managed by an industrial plantation company.The area was planted in January 2017 with Acacia crassicarpa as the main commercial species.Each tree can be effectively interpreted through aerial photos, and these photos provide a clear view of the bare soil; therefore, the area was appropriate for the purposes of this study.The area is 34 hectares, and we focused on a small area at the site.The observed area was approximately 400 m from west to east and 200 m from north to south.The annual precipitation in the area is approximately 3140 mm, and the air temperature can range from 23.2 to 36.6 • C depending on the season.The elevation of the study area is approximately 8 m above sea level.Figure 1 shows the site location.

UAS and Thermal Camera Image Processing
On 30 and 31 January 2018, flights were conducted over the study area to collect aerial photos.Two different types of data were collected from the flights: RGB images (i.e., conventional photos) and thermal infrared images.A DJI Mavic Pro multicopter UAS (DJI, Shenzhen, China) was used to collect the RGB images, and a DJI Inspire 1 Pro (DJI, Shenzhen, China) combined with a Zenmuse-XT thermal camera (FLIR Systems, Wilsonville, OR, USA) was used to collect the thermal imagery.The camera was oriented vertically toward the surface, and the side and forward overlap were set to 85% and 90%, respectively, for both the Mavic and Inspire 1 flights.The thermal images were taken three times on 30 January, including the first at 7:00, the second at 14:00, and the final photo at 19:00.On 31 January, only one image was taken, at 7:00.Data were collected four times at each location.The RGB image was taken on 31 January.
The collected photos were processed with the 3D modeling software Agisoft PhotoScan Pro version 1.4.1 (Agisoft LLC, St. Petersburg, Russia).The software performs automatic photogrammetric processing of images to align individual photos and generate mosaicked orthorectified imagery [21].Specifically, this software was also used to generate the mosaicked thermal imagery.The processing parameters were progressed with the "highest" settings for the alignment process.Tie point limits and key point limits were set to the default value of 4000 and 40,000, respectively, and the adaptive camera model fitting settings were checked.After the alignment process, using the gradual selection module, 1/10 of the total points were selected in the reconstruction uncertainty option, and the selected points were deleted.Unreasonable points that seemed to be misaligned were manually deleted from the visual interpretation.For the dense cloud generation process, the "ultra high" option was selected with the depth filtering set to "aggressive".After the dense cloud progress, the mesh generation was selected and processed with the "arbitrary" option with the face count set to "high" and with hole filling option checked.Finally, the orthomosaic was processed using the mesh as the surface data, with the blending mode as "mosaic".For the thermal imagery, the data were collected in the original digital number (DN) value format; therefore, after the mosaic was processed, the following formula was used to convert the data to absolute temperature values: where Tc is the temperature in degrees Celsius and DN is the original value observed in the image.The data were further analyzed using the converted temperature values.A canopy height model (CHM) was generated by subtracting a digital surface model (DSM) from a digital terrain model (DTM) to calculate the tree heights and determine the growth rates in the area [20,29].

UAS and Thermal Camera Image Processing
On 30 and 31 January 2018, flights were conducted over the study area to collect aerial photos.Two different types of data were collected from the flights: RGB images (i.e., conventional photos) and thermal infrared images.A DJI Mavic Pro multicopter UAS (DJI, Shenzhen, China) was used to collect the RGB images, and a DJI Inspire 1 Pro (DJI, Shenzhen, China) combined with a Zenmuse-XT thermal camera (FLIR Systems, Wilsonville, OR, USA) was used to collect the thermal imagery.The camera was oriented vertically toward the surface, and the side and forward overlap were set to 85% and 90%, respectively, for both the Mavic and Inspire 1 flights.The thermal images were taken three times on 30 January, including the first at 7:00, the second at 14:00, and the final photo at 19:00.On 31 January, only one image was taken, at 7:00.Data were collected four times at each location.The RGB image was taken on 31 January.
The collected photos were processed with the 3D modeling software Agisoft PhotoScan Pro version 1.4.1 (Agisoft LLC, St. Petersburg, Russia).The software performs automatic photogrammetric processing of images to align individual photos and generate mosaicked orthorectified imagery [21].Specifically, this software was also used to generate the mosaicked thermal imagery.The processing parameters were progressed with the "highest" settings for the alignment process.Tie point limits and key point limits were set to the default value of 4000 and 40,000, respectively, and the adaptive camera model fitting settings were checked.After the alignment process, using the gradual selection module, 1/10 of the total points were selected in the reconstruction uncertainty option, and the selected points were deleted.Unreasonable points that seemed to be misaligned were manually deleted from the visual interpretation.For the dense cloud generation process, the "ultra high" option was selected with the depth filtering set to "aggressive".After the dense cloud progress, the mesh generation was selected and processed with the "arbitrary" option with the face count set to "high" and with hole filling option checked.Finally, the orthomosaic was processed using the mesh as the surface data, with the blending mode as "mosaic".For the thermal imagery, the data were collected in the original digital number (DN) value format; therefore, after the mosaic was processed, the following formula was used to convert the data to absolute temperature values: where T c is the temperature in degrees Celsius and DN is the original value observed in the image.The data were further analyzed using the converted temperature values.A canopy height model (CHM) was generated by subtracting a digital surface model (DSM) from a digital terrain model (DTM) to calculate the tree heights and determine the growth rates in the area [20,29].

Image Classification
To delineate the LULC classes, we processed the RGB orthoimagery using a conventional supervised classification method to generate a map of the area [23].The original orthoimagery was decomposed into the R (red), G (green), and B (blue) components, and each was represented as an individual band.The RGB bands were also transformed into hue (H), saturation (S), and lightness (L) information, and each of the RGB and HSL bands, as well as the generated indices, were stacked and used in the classification process.Two additional indices were generated using the following formula: where GRVI is the green red vegetation index [30] and G and R are the green and red bands, respectively.NSLDI is the normalized saturation-lightness difference index, which is similar to the normalized saturation-value difference index (NSVDI), except that the HSV color space "value" is converted to the HSL lightness band, where S and L are the saturation and lightness bands, respectively [31].The GRVI results in better enhancement of vegetated areas, and the NSLDI can be used to effectively detect shadow areas.In addition to those noted above, three more images were generated: ratio images of B/G, R/B, and R/G, where B is the blue band.Training data for each LULC class were constructed based on the orthoimages, and all the variables/indices were stacked and used in the classification process.Multilayer perceptron (MLP) neural network classification was implemented [23].The MLP is a type of artificial neural network classifier widely used for modeling complex behaviors and patterns, as well as image classifications.It uses the back-propagation algorithm to the learn the characteristics of the class boundaries.Further details can be found in Riccioli et al. [32].This method yielded a final LULC product with five classes.Each LULC class was used for further analysis with the thermal information.The validation of the classified map was performed with an accuracy assessment generating with an error matrix.Ground truth samples were collected through the aerial photo, and the LULC classes were confirmed by the photos from the ground observations.In total, 300 sample points were generated using the systematic sampling method per LULC class, and 1200 samples were used for computing the overall accuracy of the developed categorical map.Although there are five classes, the shadow class was not considered for the validation process.

Principal Component Analysis (PCA) of Thermal Images
A PCA method was implemented using the four thermal mosaicked images to compress the information, output new band information, and interpret the thermal trends within the area.PCA in the remote sensing field is usually used to compress the information from multiple bands and convert large amounts of information into several new images called components.This procedure is implemented to reduce the number of similar images and make the results simpler to interpret [33].PCA has also been used to analyze the long-term trends of sea surface temperature (SST) and detect anomalies such as El Niño and La Niña events [34].We used the forward T-mode unstandardized method to determine the synoptic spatial patterns [35].TerrSet GIS software (Clark Labs, Worcester, USA) was used for the analysis.All the thermal images from 7:00, 14:00, 19:00, and 7:00 of the next day were selected and processed.Four components were generated, and we selected feature areas from the images and compared the thermal information at each location and each land cover type.To assess the statistical significance and compare the differences between the selected areas and variables, analysis of variance (ANOVA) and Tukey honestly significant difference (HSD) were conducted in R Studio software (Integrated Development for R. Studio, Inc., Boston, MA, USA).

Mapping of the RGB, Thermal, and Classification Maps
The orthoimage of the RGB data resulted from the 425 images collected during the aerial survey.The resolution of the RGB image was approximately 5 cm and that of the DSM and DTM was 10 cm.For the thermal image, 652, 660, and 668 images were used from the 7:00, 14:00, and 19:00 flights on 30 January, respectively, and 704 images were used from the 7:00 flight on 31 January.Each mosaicked image had a resolution of approximately 12 cm.There were some unconstructed areas in the thermal imagery depending on the observation time, especially at the edges of the scene.This issue was caused by the lack of similarities between the paired images, which made it difficult to locate the tie points for each image.This situation occasionally occurred when paired images had small variations within the image scene that made it difficult to locate matching points.Regardless, we were able to generate a large enough flight area to conduct a spatial analysis of the site.Figure 2 shows the generated RGB images and the thermal images from 7:00, 14:00, and 19:00 on 30 January and 7:00 on 31 January 2018.Each flight was conducted for approximately 40 min including battery changes.

Mapping of the RGB, Thermal, and Classification Maps
The orthoimage of the RGB data resulted from the 425 images collected during the aerial survey.The resolution of the RGB image was approximately 5 cm and that of the DSM and DTM was 10 cm.For the thermal image, 652, 660, and 668 images were used from the 7:00, 14:00, and 19:00 flights on 30 January, respectively, and 704 images were used from the 7:00 flight on 31 January.Each mosaicked image had a resolution of approximately 12 cm.There were some unconstructed areas in the thermal imagery depending on the observation time, especially at the edges of the scene.This issue was caused by the lack of similarities between the paired images, which made it difficult to locate the tie points for each image.This situation occasionally occurred when paired images had small variations within the image scene that made it difficult to locate matching points.Regardless, we were able to generate a large enough flight area to conduct a spatial analysis of the site.Figure 2 shows the generated RGB images and the thermal images from 7:00, 14:00, and 19:00 on 30 January and 7:00 on 31 January 2018.Each flight was conducted for approximately 40 min including battery changes.The Mavic aerial photo captured on 31 January 2018 was classified into the following five land cover types: bare soil, water, trees, shrubs, and shadows.The final map of the area is shown in Figure 3. Table 1 shows the error matrix for the accuracy assessment of the classification, and the overall accuracy was 95.67%.The majority of the errors were associated with the misclassification of trees and shrubs because those classes have similar band information in some areas.Another is the case of bare lands and water, which fall into the shadow class, because some areas with casting The Mavic aerial photo captured on 31 January 2018 was classified into the following five land cover types: bare soil, water, trees, shrubs, and shadows.The final map of the area is shown in Figure 3. Table 1 shows the error matrix for the accuracy assessment of the classification, and the overall accuracy was 95.67%.The majority of the errors were associated with the misclassification of trees and shrubs because those classes have similar band information in some areas.Another is the case of bare lands and water, which fall into the shadow class, because some areas with casting shadows have probably had similarity with the darker shadow areas; however, the result is good enough to be interpreted as representative of the area.The original categorical map was processed via 7 × 7-window majority filtering to smooth the results of the clustered pixels.
shadows have probably had similarity with the darker shadow areas; however, the result is good enough to be interpreted as representative of the area.The original categorical map was processed via 7 × 7-window majority filtering to smooth the results of the clustered pixels.

PCA and Spatial Characteristics of the Thermal Temperature
Four components were generated from the four thermal images observed at 7:00, 14:00, and 19:00 and 7:00 on the following day.Components 1, 2, 3, and 4 explained 67.7%, 16.1%, 10.4%, and 5.6% of the total variance, respectively.Figure 4 shows the images generated from components 1 to 4 of the thermal data analysis, and we focused on a total of five areas that exhibited specific characteristics of the site.The values were normalized, and only the visual interpretation of each component is shown.

PCA and Spatial Characteristics of the Thermal Temperature
Four components were generated from the four thermal images observed at 7:00, 14:00, and 19:00 and 7:00 on the following day.Components 1, 2, 3, and 4 explained 67.7%, 16.1%, 10.4%, and 5.6% of the total variance, respectively.Figure 4 shows the images generated from components 1 to 4 of the thermal data analysis, and we focused on a total of five areas that exhibited specific characteristics of the site.The values were normalized, and only the visual interpretation of each component is shown.shadows have probably had similarity with the darker shadow areas; however, the result is good enough to be interpreted as representative of the area.The original categorical map was processed via 7 × 7-window majority filtering to smooth the results of the clustered pixels.

PCA and Spatial Characteristics of the Thermal Temperature
Four components were generated from the four thermal images observed at 7:00, 14:00, and 19:00 and 7:00 on the following day.Components 1, 2, 3, and 4 explained 67.7%, 16.1%, 10.4%, and 5.6% of the total variance, respectively.Figure 4 shows the images generated from components 1 to 4 of the thermal data analysis, and we focused on a total of five areas that exhibited specific characteristics of the site.The values were normalized, and only the visual interpretation of each component is shown.We focus on the locations of Areas A-E.We extracted the temperature information for each land cover type in those areas.Figure 5 shows the temporal temperature trends, and there are unique characteristics in each extracted area.As an overall trend, Area A exhibits high heating and high cooling for each land cover type.The starting temperature at 7:00 is 28 • C on average, and it increases to an average of 40 • C.After the sun begins to fall and the radiation decreases, the cooling of the land surface accelerates, and the temperature drops to 30.5 • C. Additionally, a decrease of approximately 10 • C can be observed from midday to the evening.Conversely, Area B displays a slightly warmer temperature at 7:00 of approximately 28. ).The temperatures of water areas gradually increased but not as fast as the temperatures of bare soil, and water areas slowly released absorbed heat; thus, the water areas usually had the warmest temperature by night.It is clear that water bodies have high capacities to absorb heat energy, as they take time to warm up and release radiation to the atmosphere.The Tukey HSD analysis indicated that each area was significantly different from all other locations (p < 0.05).
Remote Sens. 2018, 10, x FOR PEER REVIEW 8 of 15 We focus on the locations of Areas A-E.We extracted the temperature information for each land cover type in those areas.Figure 5 shows the temporal temperature trends, and there are unique characteristics in each extracted area.As an overall trend, Area A exhibits high heating and high cooling for each land cover type.The starting temperature at 7:00 is 28 °C on average, and it increases to an average of 40 °C.After the sun begins to fall and the radiation decreases, the cooling of the land surface accelerates, and the temperature drops to 30.5 °C.Additionally, a decrease of approximately 10 °C can be observed from midday to the evening.Conversely, Area B displays a slightly warmer temperature at 7:00 of approximately 28.7 °C.The temperature gradually increases to 35.8 °C at 14:00 and slowly decreases to 31 °C.The temperature increases and decreases in Area B are smaller than those in Area A. Area C exhibits the lowest starting temperature compared with the other areas (22.7 °C), and although the temperature does gradually increase similar to that in Area B, it declines to 28.2 °C in the evening.The trends in Areas D and E are similar to those in the other areas, but the temperature in Area D continues to rapidly decrease until the morning of the next day (25.3 °C), and it actually rises in the morning of the following day in Area E (29.1 °C).Of the individual classes, most of the water areas were detected and had the warmest temperatures in the morning compared with other LULC types (average minimum difference of 0.84 °C).The temperatures of water areas gradually increased but not as fast as the temperatures of bare soil, and water areas slowly released absorbed heat; thus, the water areas usually had the warmest temperature by night.It is clear that water bodies have high capacities to absorb heat energy, as they take time to warm up and release radiation to the atmosphere.The Tukey HSD analysis indicated that each area was significantly different from all other locations (p < 0.05).

Tree Height-to-Area Characteristics
Using the information derived from the CHM, the average tree height for each area was computed.Figure 6 shows the CHM image and the boxplot of the tree heights and shrub heights for Areas A-E for the tree and shrub classes.At the top of the boxplot, the results of the Tukey HSD test are provided as letters to indicate the differences between each area.In the CHM image, the west side of the site exhibited some errors associated with the DTM, and these errors resulted in large height values that were inaccurate.However, these areas were not considered in further analyses; therefore, we did not address these errors.As an overall trend, Area D displayed the highest average tree height and shrub height, and Area B exhibited the lowest heights.Areas C and E had similar tree heights, but there were few high values of the shrub height in Area E. We will form an additional hypothesis in the discussion section.

Tree Height-to-Area Characteristics
Using the information derived from the CHM, the average tree height for each area was computed.Figure 6 shows the CHM image and the boxplot of the tree heights and shrub heights for Areas A-E for the tree and shrub classes.At the top of the boxplot, the results of the Tukey HSD test are provided as letters to indicate the differences between each area.In the CHM image, the west side of the site exhibited some errors associated with the DTM, and these errors resulted in large height values that were inaccurate.However, these areas were not considered in further analyses; therefore, we did not address these errors.As an overall trend, Area D displayed the highest average tree height and shrub height, and Area B exhibited the lowest heights.Areas C and E had similar tree heights, but there were few high values of the shrub height in Area E. We will form an additional hypothesis in the discussion section.

Thermal Mapping
The thermal mapping of the peat environment has successfully worked through the SfM process.As Maes et al. [28] showed the possibilities of thermal mapping, this case has also shown but with better precision due to the sensors used, which were more optimized for the utilization together with the UAS.Maes et al. [28] noted that camera position is important in the alignment process of the SfM, and the Zenmuse-XT could clearly be recognized for such parameters, because it was integrated in the UAS system (pitch-yaw-roll information is available).Instead, some weakness in the alignment process was shown in the morning time, where thermal trends in the peat were not showing much variety across the site, giving difficulties in selecting the tie points for each image.Several alignment processes were repeated to establish the correct alignment; however, selecting manual tie points could be one option to aid in providing a reference to make better alignments if it does not well work.
The mapping of the thermal trends across the site has shown the daily cycle of thermal trends.The thermal characteristics of the individual LULC classes were relatively clear in the different areas of the study site.Overall, bare soil exhibited the highest temperatures, and vegetated areas displayed the lowest temperatures.Bare soil can absorb more heat than can other land cover types [36], and vegetation has lower temperatures due to evapotranspiration from the leaves related to the photosynthesis process [17].The characteristics of water indicate how slowly it heats and cools within a day.The patterns of the temporal data indicated that in the morning, the water temperatures were slightly warmer and that they gradually increased throughout the day.The water temperature did reach that of bare soil; however, it began to decline later in the day but retained more heat than the other cover classes and remained warmer than the other classes in the evening.If the surface peat layer was showing the temperature trend as it was seen, then important factors can be discussed relating to the oxidation process.If the oxidation process is continuously increasing with the temperature trends shown by Jauhiainen et al. [11] and Hilasvuori et al. [12], then with respect to the temperature of the peat in the site, in the morning and night, the majority of the area is more than 25° C, and in the afternoon, it can even peak to 50° C. If so, then the progress of the decomposition is accelerating more than expected, and there is a possibility of underestimating the possible emission rates of the peat areas.Even though the plantation trees can be potential objects for shading the peats, the fast-growing acacia can still take approximate one year to make relevant shading areas.This  Areas A-E).The small letters at the top of the box plot indicate the differences in tree heights (computed by the Tukey HSD method (p < 0.05)).

Thermal Mapping
The thermal mapping of the peat environment has successfully worked through the SfM process.As Maes et al. [28] showed the possibilities of thermal mapping, this case has also shown but with better precision due to the sensors used, which were more optimized for the utilization together with the UAS.Maes et al. [28] noted that camera position is important in the alignment process of the SfM, and the Zenmuse-XT could clearly be recognized for such parameters, because it was integrated in the UAS system (pitch-yaw-roll information is available).Instead, some weakness in the alignment process was shown in the morning time, where thermal trends in the peat were not showing much variety across the site, giving difficulties in selecting the tie points for each image.Several alignment processes were repeated to establish the correct alignment; however, selecting manual tie points could be one option to aid in providing a reference to make better alignments if it does not well work.
The mapping of the thermal trends across the site has shown the daily cycle of thermal trends.The thermal characteristics of the individual LULC classes were relatively clear in the different areas of the study site.Overall, bare soil exhibited the highest temperatures, and vegetated areas displayed the lowest temperatures.Bare soil can absorb more heat than can other land cover types [36], and vegetation has lower temperatures due to evapotranspiration from the leaves related to the photosynthesis process [17].The characteristics of water indicate how slowly it heats and cools within a day.The patterns of the temporal data indicated that in the morning, the water temperatures were slightly warmer and that they gradually increased throughout the day.The water temperature did reach that of bare soil; however, it began to decline later in the day but retained more heat than the other cover classes and remained warmer than the other classes in the evening.If the surface peat layer was showing the temperature trend as it was seen, then important factors can be discussed relating to the oxidation process.If the oxidation process is continuously increasing with the temperature trends shown by Jauhiainen et al. [11] and Hilasvuori et al. [12], then with respect to the temperature of the peat in the site, in the morning and night, the majority of the area is more than 25 • C, and in the afternoon, it can even peak to 50 • C. If so, then the progress of the decomposition is accelerating more than expected, and there is a possibility of underestimating the possible emission rates of the peat areas.Even though the plantation trees can be potential objects for shading the peats, the fast-growing acacia can still take approximate one year to make relevant shading areas.This finding can spatially characterize the peat environment and the actual oxidation progress that could be occurring in the region.

PCA Components Relating to Topographic Variables
The PCA method is usually very effective in transforming large quantities of temporal data, such as thermal information, into compressed samples.Unfortunately, our data were collected in only four time periods, which resulted in imagery that was very similar to the original temperature data.However, with more images and data, the PCA results would be more robust.When focusing on each spatial area, two main patterns could be identified.One area with high temperature increases and decreases (Area A) was observed, and another area with low temperature increases and decreases (Area B) was seen.Some differences in other areas could also be observed, such as an area with low temperatures in the morning (Area C), an area where the temperature further declined in the evening (Area D), and an area where the temperature even began to slightly increase the following morning (Area E).The plantation company has provided microtopographical elevation data for the study site (Figure 7).finding can spatially characterize the peat environment and the actual oxidation progress that could be occurring in the region.

PCA Components Relating to Topographic Variables
The PCA method is usually very effective in transforming large quantities of temporal data, such as thermal information, into compressed samples.Unfortunately, our data were collected in only four time periods, which resulted in imagery that was very similar to the original temperature data.However, with more images and data, the PCA results would be more robust.When focusing on each spatial area, two main patterns could be identified.One area with high temperature increases and decreases (Area A) was observed, and another area with low temperature increases and decreases (Area B) was seen.Some differences in other areas could also be observed, such as an area with low temperatures in the morning (Area C), an area where the temperature further declined in the evening (Area D), and an area where the temperature even began to slightly increase the following morning (Area E).The plantation company has provided microtopographical elevation data for the study site (Figure 7).A brief analysis was performed to interpret the thermal characteristics of the heating and cooling trends in the area in relation to topographical variables.Several factors were computed using the digital elevation model (DEM), and the System for Automated Geoscientific Analysis (SAGA) wetness index, catchment area, modified catchment area, catchment slope, slope length and steepness (LS) factor, wind exposure index, slope, terrain ruggedness index, terrain position index, and diurnal anisotropic heating were all computed using SAGA GIS software ver.6.3.0 [37].Using the heating information from 14:00-7:00 and cooling information from 19:00-14:00 as dependent variables, all the topographic indices were incorporated as independent variables, and machine learning (MLP) regression was implemented to determine which variables exhibited a response to the heating and cooling effects (all LULC types were used as dummy variables).As a result, the heating and cooling models displayed values of R 2 = 0.23 and R 2 = 0.36, respectively.The important variables for the heating model were the wetness index, modified catchment area, and wind exposure, and those for the cooling model were the elevation, wetness index, and wind exposure.The heating and cooling of each LULC type can be described by several factors, such as the wind exposure and topographical openness [38]; however, even with the nonlinear regression approach, not all the model's results could be solely explained by surface terrain information, and there remains some uncertainty regarding the thermal trends at the study site.Just for brief interpretation, Figure 8 shows the regression between the PCA components and the topographic variables.We will not go into detail at this moment; however, the direct terrain factors (e.g., DEM), wetness, and modified catchment areas showed responses to the PCA components.A brief analysis was performed to interpret the thermal characteristics of the heating and cooling trends in the area in relation to topographical variables.Several factors were computed using the digital elevation model (DEM), and the System for Automated Geoscientific Analysis (SAGA) wetness index, catchment area, modified catchment area, catchment slope, slope length and steepness (LS) factor, wind exposure index, slope, terrain ruggedness index, terrain position index, and diurnal anisotropic heating were all computed using SAGA GIS software ver.6.3.0 [37].Using the heating information from 14:00-7:00 and cooling information from 19:00-14:00 as dependent variables, all the topographic indices were incorporated as independent variables, and machine learning (MLP) regression was implemented to determine which variables exhibited a response to the heating and cooling effects (all LULC types were used as dummy variables).As a result, the heating and cooling models displayed values of R 2 = 0.23 and R 2 = 0.36, respectively.The important variables for the heating model were the wetness index, modified catchment area, and wind exposure, and those for the cooling model were the elevation, wetness index, and wind exposure.The heating and cooling of each LULC type can be described by several factors, such as the wind exposure and topographical openness [38]; however, even with the nonlinear regression approach, not all the model's results could be solely explained by surface terrain information, and there remains some uncertainty regarding the thermal trends at the study site.Just for brief interpretation, Figure 8 shows the regression between the PCA components and the topographic variables.We will not go into detail at this moment; however, the direct terrain factors (e.g., DEM), wetness, and modified catchment areas showed responses to the PCA components.
. Regression analysis between the principal components and the topographic variables computed from the System for Automated Geoscientific Analysis (SAGA) GIS software.
Al Kayssi et al. [39] noted that the soil temperature differences between day and night are affected by the soil moisture content.In drier conditions, greater temperature increases can be observed at midday, and at night, large decreases occur but wetter conditions suppress the temperature increases, and gradual cooling occurs.This case is similar to that observed in Area A and Area B. If this trend reflects the differences in the moisture content of the soil, then the surface thermal data may provide information regarding the soil moisture or groundwater layer [40].Additionally, the tree heights in Area B were lower than those in other areas.Observations from other areas where the conditions were often wet indicated poor tree growth.Assuming that this trend is applicable throughout the study area, we can hypothesize that both the surface and underground (water table) conditions are wet in these areas.Therefore, groundwater conditions can potentially be determined by monitoring peat using a thermal camera, which can provide important spatial information that relates to the GHG flux from peat [10,11,18].However, for a more in-depth analysis, additional belowground data should be collected, including the moisture content and the groundwater depth below peatland.This potential analysis represents a second phase in studying the peat environment, and there will be a future study on this issue.

Conclusions
For the first time, we have mapped and visualized both the spatial and temporal trends of the thermal characteristics at an extremely high resolution in the peatland plantation forests of West Kalimantan, Indonesia.These trends were visualized by direct measurements of thermal radiation using a UAS and TIR camera and not via interpolation methods.The high spatial resolution of the land surface clearly displays the trend of each land cover class; moreover, the feature areas within the study site exhibit unique trends and thermal patterns.The direct measurements of thermal information can provide strong information regarding the understanding of how oxidation progress is occurring in the region, and decisions could be made for its strategical management.The PCA method was used to successfully compress and generate new images and identify the associated spatial patterns.Additionally, relationships that might be connected to the wetness of the area were identified.Tree growth also displayed differences among the extracted areas; therefore, the differences in site productivity and various environmental factors, including the groundwater and soil nutrient conditions, will be assessed in the future.With the combination of thermal and multispectral sensors, the vegetation health could be one parameter for monitoring.Several results Al Kayssi et al. [39] noted that the soil temperature differences between day and night are affected by the soil moisture content.In drier conditions, greater temperature increases can be observed at midday, and at night, large decreases occur but wetter conditions suppress the temperature increases, and gradual cooling occurs.This case is similar to that observed in Area A and Area B. If this trend reflects the differences in the moisture content of the soil, then the surface thermal data may provide information regarding the soil moisture or groundwater layer [40].Additionally, the tree heights in Area B were lower than those in other areas.Observations from other areas where the conditions were often wet indicated poor tree growth.Assuming that this trend is applicable throughout the study area, we can hypothesize that both the surface and underground (water table) conditions are wet in these areas.Therefore, groundwater conditions can potentially be determined by monitoring peat using a thermal camera, which can provide important spatial information that relates to the GHG flux from peat [10,11,18].However, for a more in-depth analysis, additional belowground data should be collected, including the moisture content and the groundwater depth below peatland.This potential analysis represents a second phase in studying the peat environment, and there will be a future study on this issue.

Conclusions
For the first time, we have mapped and visualized both the spatial and temporal trends of the thermal characteristics at an extremely high resolution in the peatland plantation forests of West Kalimantan, Indonesia.These trends were visualized by direct measurements of thermal radiation using a UAS and TIR camera and not via interpolation methods.The high spatial resolution of the land surface clearly displays the trend of each land cover class; moreover, the feature areas within the study site exhibit unique trends and thermal patterns.The direct measurements of thermal information can provide strong information regarding the understanding of how oxidation progress is occurring in the region, and decisions could be made for its strategical management.The PCA method was used to successfully compress and generate new images and identify the associated spatial patterns.Additionally, relationships that might be connected to the wetness of the area were identified.Tree growth also displayed differences among the extracted areas; therefore, the differences in site productivity and various environmental factors, including the groundwater and soil nutrient conditions, will be assessed in the future.With the combination of thermal and multispectral sensors, the vegetation health could be one parameter for monitoring.Several results have been shown, and overall, this study yielded interesting results that demonstrate the potential of using UASs for monitoring the thermal conditions in peat environments.
Author Contributions: K.I. took the lead in conceptualizing the study, writing the manuscript, and data collection/analysis.K.W., T.K., N.A.P., S.S., T.K., and O.K. all contributed substantially to data collection, methodology development, validation, analysis, and writing and revising/editing the manuscript.
Funding: This research was funded by the Japan Society of the Promotion of Science (JSPS) KAKENHI under grant number numbers JP16K07978 and 16H05747 and by the Research Institute for Humanity and Nature (RIHN) under grant number 14200117.

Figure 1 .
Figure 1.Location of the study site.

Figure 1 .
Figure 1.Location of the study site.

Figure 3 .
Figure 3. Land cover map of the study area.

Figure 4 .
Figure 4. Component images from the principal component analysis (PCA): (a) Component 1; (b) Component 2; (c) Component 3; and (d) Component 4. Areas A-E were extracted from a visual interpretation of areas that displayed different characteristics compared with other locations.The scale values were normalized.

Figure 3 .
Figure 3. Land cover map of the study area.

Figure 3 .
Figure 3. Land cover map of the study area.

Figure 4 .
Figure 4. Component images from the principal component analysis (PCA): (a) Component 1; (b) Component 2; (c) Component 3; and (d) Component 4. Areas A-E were extracted from a visual interpretation of areas that displayed different characteristics compared with other locations.The scale values were normalized.

Figure 4 .
Figure 4. Component images from the principal component analysis (PCA): (a) Component 1; (b) Component 2; (c) Component 3; and (d) Component 4. Areas A-E were extracted from a visual interpretation of areas that displayed different characteristics compared with other locations.The scale values were normalized.
7 • C. The temperature gradually increases to 35.8 • C at 14:00 and slowly decreases to 31 • C. The temperature increases and decreases in Area B are smaller than those in Area A. Area C exhibits the lowest starting temperature compared with the other areas (22.7 • C), and although the temperature does gradually increase similar to that in Area B, it declines to 28.2 • C in the evening.The trends in Areas D and E are similar to those in the other areas, but the temperature in Area D continues to rapidly decrease until the morning of the next day (25.3 • C), and it actually rises in the morning of the following day in Area E (29.1 • C).Of the individual classes, most of the water areas were detected and had the warmest temperatures in the morning compared with other LULC types (average minimum difference of 0.84 • C

Figure 5 .
Figure 5. Temporal characteristics of temperature trends for each land cover type in Areas A-E.The small letters in the upper-right corners indicate that the location is different from all other locations (computed by the Tukey honesty significance difference (HSD) method (p < 0.05)).

Figure 5 .
Figure 5. Temporal characteristics of temperature trends for each land cover type in Areas A-E.The small letters in the upper-right corners indicate that the location is different from all other locations (computed by the Tukey honesty significance difference (HSD) method (p < 0.05)).

15 Figure 6 .
Figure 6.Canopy height model (CHM) of the study area and the heights of the tree and shrub classes at different locations (Areas A-E).The small letters at the top of the box plot indicate the differences in tree heights (computed by the Tukey HSD method (p < 0.05)).

Figure 6 .
Figure 6.Canopy height model (CHM) of the study area and the heights of the tree and shrub classes at different locations (Areas A-E).The small letters at the top of the box plot indicate the differences in tree heights (computed by the Tukey HSD method (p < 0.05)).

Figure 7 .
Figure 7. Digital elevation model (DEM) representing the microtopography of the study site.Visual interpretation can be used to identify some characteristics of the terrain features that exhibit some similarity with the thermal information.

Figure 7 .
Figure 7. Digital elevation model (DEM) representing the microtopography of the study site.Visual interpretation can be used to identify some characteristics of the terrain features that exhibit some similarity with the thermal information.

Figure 8 .
Figure 8. Regression analysis between the principal components and the topographic variables computed from the System for Automated Geoscientific Analysis (SAGA) GIS software.

Table 1 .
Error matrix of the generated categorical map.

Table 1 .
Error matrix of the generated categorical map.

Table 1 .
Error matrix of the generated categorical map.