Surface Elevation Changes Estimation Underneath Mangrove Canopy Using SNERL Filtering Algorithm and DoD Technique on UAV-Derived DSM Data

: Estimating surface elevation changes in mangrove forests requires a technique to ﬁlter the mangrove canopy and quantify the changes underneath. Hence, this study estimates surface elevation changes underneath the mangrove canopy through vegetation ﬁltering and Difference of DEM (DoD) techniques using two epochs of unmanned aerial vehicle (UAV) data carried out during 2016 and 2017. A novel ﬁltering algorithm named Surface estimation from Nearest Elevation and Repetitive Lowering (SNERL) is used to estimate the elevation height underneath the mangrove canopy. Consequently, DoD technique is used to quantify the elevation change rates at the ground surface, which comprise erosion, accretion, and sedimentation. The signiﬁcant ﬁndings showed that region of interest (ROI) 5 experienced the highest volumetric accretion (surface raising) at 0.566 cm 3 . The most increased erosion (surface lowering) was identiﬁed at ROI 8 at − 2.469 cm 3 . In contrast, for vertical change average rates, ROI 6 experienced the highest vertical accretion (surface raising) at 1.281 m. In comparison, the most increased vertical erosion (surface lowering) was spotted at ROI 3 at − 0.568 m. The change detection map and the rates of surface elevation changes at Kilim River enabled authorities to understand the situation thoroughly and indicate the future situation, including its interaction with sea-level rise impacts.


Introduction
In the mangrove forest region, the environment is very dynamic because of dense canopy vegetation and the tidal inundation that causes the surface elevation to be partially submerged, especially during high tide and the highest astronomical tides (HAT). The physical measurement approach is challenging because of the rough terrain conditions full of mangrove roots and muddy soil that complicate the physical work [1,2]. Total stations, terrestrial laser scanners (TLS), electronic distance measurement (EDM), and other surveying equipment that needs a tripod are hardly set up in this terrain. The alternative way to monitor erosion in the mangrove environment is to use the shifting of the mangrove boundary/vegetation line, as discovered by [3][4][5]. However, the study only looks at the surface of the mangrove canopy and ignores the hidden features beneath it. Any physical changes on the ground are hidden by mangrove vegetation, especially from an aerial view, and these obstacles have inspired this study. To quantify surface elevation changes on sites with small-scale areas, a few terrestrial measurements such as physical surveys, erosion pins, and the surface elevation table-marker horizon (SET-MH) method are commonly used.
However, aerial monitoring is still convenient compared to terrestrial in terms of mobility, time-consuming, data abundance, and the ability to cover a large area in a short time. Besides, tidal inundation in mangrove areas that often becomes a problem for researchers (if using a physical survey method) could be avoided using aerial monitoring [6,7].
Most previous surface extraction research has focused on surface elevation and surface change monitoring of the bare earth topography. The studies [8][9][10][11][12][13] use satellite imagery to monitor erosion in the coastal region. Meanwhile, [14][15][16][17][18][19][20] used UAV imagery to assess erosion in bare earth areas. The previous scholar observes the erosion based on the visible changes in topography based on the bare earth's spatial and temporal analysis. Previous researchers had made fewer attempts to evaluate the surface elevation changes caused by erosion, accretion (deposition), and sedimentation rates at riverbank areas covered by canopies or trees. At the canopy-covered riverbank area, the canopy has obstructed the surface, especially from an aerial view. The mangrove canopy covering the riverbank should be removed using the vegetation filtering method to monitor erosion in the canopycovered area.
For decades, researchers have developed many filtering methods to remove aboveground objects such as vegetation, buildings, and other human-made structures. Filtering approaches such as LAStools triangulated irregular networks (TIN) densification, the typical method in Agisoft Photoscan, and others have been used to filter point clouds by utilising an airborne laser scanner or light detection and ranging (LiDAR) [44][45][46]. Similar to UAV imagery data, it uses several algorithms, for instance, a progressive morphological filter (PMF), a simple morphological filter (SMRF), or a cloth simulation filter (CSF), and a structural filter, Caractérisation de Nuages de Points (CANUPO), for vegetation filtering purposes [47]. Unlike LiDAR, a UAV photogrammetric survey could not penetrate partly through the vegetation layer and, therefore, might have failed to generate actual ground points below a vegetated surface. The aerial photograph could not capture the information below dense vegetation using a typical sensor (e.g., a Red, Green, and Blue (RGB) sensor) unless the captured RGB images are analysed using a specific method, like an excessive greenness vegetation index (ExG-VI), for surface extraction in dense forest areas [45].
Previous studies related to vegetation filtering using LiDAR have already been conducted by other researchers, which could lead to various filtering algorithms being invented and improvised for better output accuracy [44][45][46]. However, less study was performed by other researchers to utilise UAS products, especially the DSM model, in estimating surface elevation underneath mangrove-covered areas since it only captured the top surface of mangrove trees instead of the ground surface. Anders et al. (2019) [45] evaluated the comparison of DEM accuracy using common filtering algorithms such as the standard built-in method in Agisoft Photoscan, colour-based filtering using vegetation index (VI), iterative surface lowering (ISL), triangulated irregular networks (TIN), and a combination of iterative surface lowering and the VI method (ISL + VI). Since they revealed the potential of the UAS optical sensor in performing the task just like LiDAR, it inspired this study to use the UAS approach with a novel SNERL filtering algorithm for generating surface elevation underneath mangrove-covered areas. For that reason, capturing the top of mangrove surface information is a sufficient way before proceeding to image processing using the structure from the motion-multi-view stereo (SfM-MVS) method.
The existing problem in the application of UAV in mangrove canopy change monitoring is the difficulty of using the optical photogrammetry sensor (i.e., RGB-Complementary Metal Oxide Semiconductor (CMOS) in estimating elevation changes underneath the mangrove canopy. Optical sensors usually detect the mangrove canopy as the surface height model, not as the maximum height above the mangrove area. The hidden mangrove surface height below the canopy is normally obstructed because it is hard to quantify using a conventional UAV photogrammetry survey, unless a certain filtering algorithm is used to reduce the mangrove canopy height to ground level. Once the hidden mangrove surface height is uncovered, the monitoring of surface elevation changes can be carried out using the GCD technique.
Hence, this study focuses on improving the vegetation filtering algorithm for extracting the surface elevation underneath the mangrove canopy to form a novel algorithm called Surface estimation from Nearest Elevation and Repetitive Lowering (SNERL). This filtering algorithm removes the mangrove canopy from the terrain level. It then generates a DEM whose accuracy is almost comparable to the physical measurement, using total station, levelling, or the Global Navigation Satellite System (GNSS). The vertical datum underneath the mangrove canopy is based on the nearest identified surface level, such as a riverbank or opening surface in the middle of a dense vegetation area. The digital surface model (DSM), which represents elevation at the mangrove canopy height, has been collected from UAV data that comprise two different flying durations, during low and high tide. UAV raw data (in raw aerial images) is processed based on the SfM-MVS method using commercial software such as Agisoft Metashape. Although it can generate orthophoto, contours, and three-dimensional (3D) models, DSM is the most crucial output in this study. The other results, such as orthophoto, serve as supporting data. Since the DSM image from the UAV aerial survey only captured the top of topographic features such as vegetation and buildings, it needs a filtering method to filter and reduce it until it reaches terrain level (ground surface). This study might require additional data, such as physical topographic data from real-time kinematic (RTK) observations, to validate the accuracy of the filtered DEM.

Study Area
The Kilim River in Langkawi Island, Kedah, Malaysia, has been chosen as the study site ( Figure 1). This study area was designated as a UNESCO Kilim Karst Geoforest Park (KKGP) in 2007. The Kilim River has a semi-diurnal tide because it is located at 6 • 21.518 -6 • 26.093 N and 99 • 51.159 -99 • 51.159 E [4]. Although the Kilim River is approximately three kilometres from the coastal area, the tidal effect still exists and affects the Kilim River. This area has been granted Geopark status by the United Nations Educational, Scientific, and Cultural Organization (UNESCO) because of its impressive rock formations (limestone or karst), the richness of mangrove forests, and its unique geological significance. The study area displays natural mountain valley topography with an average slope of 30 • , and some areas are up to 80 • , making it quite challenging for a drone pilot to fly over the entire region. The study area is also surrounded by various geomorphological structures such as muddy soil at the riverbank, dense mangrove forest, limestone, karst, and different kinds of rock formations in hilly areas. This region was chosen as the study area because it involves additional land cover, including high-density vegetation, low vegetation, shrubs, and different geomorphological features. The Kilim River was selected due to the previous study by [3], which discovered significant mangrove dynamic changes in this area.

UAV Data Acquisition
The unmanned aerial vehicle (UAV) employed in this study was a Da-Jiang Innovations (DJI) Phantom 4 Advanced model manufactured by SZ DJI Technology Corporation Limited (Shenzhen, China). The camera model attached has a focal length of 3.61 mm, a resolution of 4000 × 3000, and pixel sizes of 1.56 × 1.56 µm. DJI GO software was used to plan the flight path and flight parameters. It includes the flight altitude, frontal and side overlap, and the size of the coverage area. In 2016, the flying altitude during low tide was 149 m.
Meanwhile, in 2017, the flying altitude during low tide was 228 m. The reason for choosing a flight altitude of 228 m above the ground is the topographic conditions of the study area that is surrounded by hills. During UAV data acquisition, the flight path was set for auto-deploy of the aircraft with the self-propelled piloting system. To perform a smooth take-off and landing, a flight plan should be made with thorough consideration of all aspects by the pilot. It is to avoid any unwanted incident, especially a collision with a hill. In this study area, the maximum elevation of a hill was indicated at 200 m at most. For that reason, 228 m of flying altitude is a safe precaution to minimise the collision risk during the UAV data acquisition process.
According to flight data for both epochs, each flight captured a different number of aerial images. UAV data collection in 2016 captured more images and more camera stations than in 2017 because of its high-flying altitude. For the coverage area, 0.661 km 2 was the total coverage area during low tide in 2016, while 0.688 km 2 was during low tide in 2017.

GNSS Data Observation for GCPs
For this study, eight artificial targets and a rubber mat with the "X" mark (with dimensions of 2 m × 2 m) were set up at distributed locations along the Kilim River (as displayed in Figure 2). One-hour observation of the GNSS survey is carried out using a Topcon Trimble GR-5 model to measure horizontal (latitude and longitude) and vertical data (ellipsoidal height) in GDM2000 at each GCP. The accuracy of the Topcon GR-5 model when using the static technique is 3.0 mm + 0.5 ppm (horizontal) and 5.0 mm + 0.5 ppm (vertical), whereas it is 10 mm + 1 ppm (horizontal) and 15 mm + 1 ppm (vertical) when implementing the Real-Time Kinematic (RTK) technique [48]. In this study, a static method with 1-h observation was chosen over the RTK method for GCP data acquisition. The main reason is the accuracy of the observed GCP that reaches 3.0 mm + 0.5 ppm (horizontal) and 5.0 mm + 0.5 ppm (vertical) using the static method, while the RTK method could reach 10 mm + 1 ppm (horizontal) and 15 mm + 1 ppm (vertical). Not only is it used for GCP data acquisition, but a similar station is also used as the base station for GNSS vertical accuracy assessment in order to verify the quality of the filtered DEM after passing through the SNERL filtering algorithm. The other reason the static method was used is the requirement for great accuracy in ellipsoidal height from GNSS observations in order to determine the orthometric height. To obtain the highest accuracy of orthometric height, the highest accuracy of ellipsoidal height through a static method is necessary while the geoidal height is determined from the Malaysian Geoid model (MyGeoid). The orthometric height is identified by subtracting the ellipsoid height (GNSS survey) from the geoid height (MyGeoid model) to make sure the vertical data (orthometric height) reaches the global mean sea level (MSL). Hence, the final altitude (z) was inserted along with the coordinates (latitude and longitude) during image processing, and it was based on orthometric height (height above mean sea level). After the processing using Agisoft Metashape, the UAV output, especially the orthomosaic and DSM models, resulted in a great GSD value since the inserted value of vertical value (altitude) was referred to as orthometric-based height.
As a local state coordinate system that specifically applies to the Kedah and Perlis regions, the Geocentric Datum of Malaysia 2000 (GDM 2000) was used. To ensure that the vertical data (orthometric height) reaches the global mean sea level, the orthometric height is calculated by subtracting the ellipsoid height (GNSS survey) from the geoid height (MyGeoid model). Hence, the vertical data were inserted as orthometric height (height above mean sea level) during the SfM-MVS process, along with the coordinates (latitude and longitude) of GCPs.

Generating UAV-Derived DSM through SfM-MVS Method
The SfM-MVS photogrammetry technique was used to create three-dimensional point clouds (3D), which allows for the reconstruction of topography from randomly dispersed and orientated images from uncalibrated cameras. SfM-MVS was chosen to reconstruct the three-dimensional (3-D) structure from a series of overlapping images. This is due to its analytical methods being based on a computer-aided solution of mathematical algorithms and digital softcopy of photogrammetry-based data, eliminating the need for opto-mechanical hardware. The images were processed using the commercially available SfM software, Agisoft Metashape v1.6.4, which uses a typical SfM-MVS progression. The data processing was performed using Intel Core i7 3.6 GHz CPUs, NVIDIA Quadro K80 GPUs, and 32 GB of RAM.

Segregation of Vegetation and Non-Vegetation Areas Using Vegetation Indices (VI)
There are numerous vegetation indices for various applications, but the most effective of which are derived using green, red, red-edge, near-infrared, and/or infrared bands, such as the Normalised Difference Vegetation Index (NDVI) and the Enhanced Vegetation Index (EVI) (EVI). In lightweight UAVs, such as the one used in this study (DJI Phantom 4 Advanced), the maximum payload is limited; hence, the types of sensors onboard are restricted by their weight. While multi-spectral cameras are becoming more compact and some sensors can be fitted into small aircraft (e.g., the Parrot Sequoia), many UAVs already implement these new capabilities. These cameras have CMOS sensors that are sensitive to visible light until near-infrared (NIR) (400-900 nm), along with RGB channels and filters that block radiation over 650 nm.
There are vegetation indices solely based on visible light bands, such as the Visible Vegetation Index (VVI) and Excessive Greenness (ExG, Equation (2)), distinguishing vegetation from bare land. VVI is a scale from 0 to 1, with lower values indicating bare ground and higher vegetation values, as shown in Equation (1).
R, G, and B are the red, green, and blue channel values, respectively. RGB 0 is the vector of the reference green colour (in this case, R 0 = 40, G 0 = 60, and B 0 = 10), and W is a weight component (here 1) [49]. As indicated in Equation (2), ExG is a continuous index with higher values representing bright green reflecting surfaces such as vegetation and lower values representing bare ground or obstructions.
The initial stage in extracting the Vegetation Index (VI) is to generate the visible vegetation index (VVI), which is then followed by ExG parameterisation (Figure 3). ExG is a continuous index where larger values represent bright green reflecting surfaces such as vegetation, and lower values represent bare terrain or obstructions. The R, G, B values were then normalised, calculated, and visually compared to determine how much vegetation was detected. It was discovered that VVI was useful for distinguishing minor variations between vegetation types and that ExG provided a greater contrast between vegetation and barren land. ExG was the indicator of selection because the purpose was to classify vegetation points and ground points. During the filtering process, any point in the point cloud with an ExG greater than 0.1 was omitted. The output for this method is the filtered vegetation from DSM that consists only of vegetation with greenness featured and excludes bare or terrain surfaces.
This process utilised VI as the classifier of vegetation and non-vegetation features. Using the Normalised Differential Vegetation Index (NDVI) algorithm, the features with the higher VI were considered the vegetation area and vice versa. VI chooses the features that are higher in the green band (Band 2) as the vegetation area, while the red band (Band 1) and blue band (Band 3), which are low in VI, are non-vegetation ( Figure 2).
Then, Excessive Greenness (ExG) was performed using the Binary Thresholding Function to classify between vegetation and non-vegetation using the binary numbers 0 (low) and 1 (high) as the index. So, any features with ExG > 0.1 in the NDVI transformed raster were considered vegetation and kept while the rest were removed. The ExG raster was converted into an integer separated by the binary number (0, 1) using an integer (spatial analyst) before being converted back into a polygon shapefile, as shown in Figure 2. The last process in this sub-section is to select only features with ExG > 0.1 to be overlayed with UAV DSM data and thus clipped to remove any unwanted Region of Interest (ROI) within the whole DSM data. Hence, the output in this process is the clipped DSM that is overlaid on the UAV DSM data that will be used in further processing. The exact process was repeated for the second epoch of data (Low tide 2017).

Vegetation Filtering Using SNERL Algorithm
The SNERL algorithms originated from a linear prediction by [50], which later developed into iterative surface lowering (ISL), before being combined as Iterative Surface Lowering (ISL) + Vegetation Index (ISL + VI), as explained by [45]. The SNERL algorithms operate as a smoothing algorithm for irregular heights of vegetation while considering the nearest elevation height as the reference. By using the nearest surface elevation height (negative residual) that is not obstructed by vegetation canopy, the repetitive filtering is performed for point cloud data, so that it can be reduced until it achieves the vertical datum level. As a reference vertical datum, any surface height that is not covered by vegetation is considered ground surface, which is also known as "negative residual".
The surface height or topography has a higher probability of having negative residuals than the vegetation point, which has a higher likelihood of having small negative or positive residuals. Hence, H i as vegetation height, represented by a point cloud or DSM profile, is considered a positive residual in this algorithm. Eliminating residuals based on the calculated weight is repetitively carried out until it removes all vegetation points to reach a negative residual level ( Figure 4). The repetitive process of eliminating residuals depends on the weight value, H a,b , which is computed using the proposed SNERL algorithm in Equation (3) as follows (example on the 1st iteration): The difference between this algorithm and other filtering methods is that an exposed bare surface in the middle of the mangrove forest acts as the reference elevation. In the concept of SNERL filtering, the exposed bare surface reveals the elevation height, and there will be negative residuals. Such an approach could guide mangrove surface lowering towards a specific height, and the filtering process would eliminate residuals, iteratively towards a higher accuracy of the hidden ground surface. The higher accuracy of the filtered ground surface beneath the mangrove canopy will be helpful in further study of mangrove geomorphology and, thus, provide an alternative for researchers to use a remote platform for advanced related research instead of only physical techniques.
The SNERL algorithm scripting was performed in R software, as shown in Figure 3. The unfiltered surface height in DSM was clustered into a few categories of height, in which each type represented a particular feature. For example, each feature has an obvious height range. For example, bare land ranged from 0 to 1 m, short grass ranged from 20 to 25 cm, tall grass ranged from 2.1 m tall with roots extending down into the soil to 1.8 m, subshrub ranged from 1 to 2 m, small shrubland ranged from 1 to 2 m, and tall shrubs ranged from 2 to 5 m. The SNERL algorithm will filter the height of DSM until it reaches this 3rd iteration as the final and lowest DSM cluster. With a height ranging from 15 to 20 m, this category was considered the highest DSM cluster.
In the first iteration, the filtering process only affects the surface height of the water between 15 and 20 m. At a height ranging from 10 to 15 m, this category was considered the medium DSM cluster. Then, in the second iteration, the filtering process only affects the surface height between 10 and 15 m, while for another category of height, it will not affect the filtering process. The other type of height was between 5 and 10 m, and the 3rd iteration will affect features included in this group of height. The remaining DSM surface, which has no more extended height above 5 m, was considered to have reached terrain level and no longer has vegetation covering the earth's surface. In this state, the DSM no longer has vegetation or shrubs and can be a Digital Elevation Model (DEM). Hence, the final DEM will be analysed for its accuracy assessment, and if it meets the desired accuracy, it can proceed to the next stage, known as GCD processing.

Quantification of Surface Elevation Changes Using Geomorphological Change Detection (GCD) Method
GCD is a technique to measure geomorphological change processes on the surface that comprises erosion, accretion (deposition), and sedimentation derived from temporal topographic surveys. Because of the qualitative and spatially explicit results, GCD is quickly becoming a more prevalent river monitoring and restoration method. The volumetric change is determined by comparing the surface elevation of the digital elevation model (DEM) obtained from a periodic topographic survey with that of the DEM. Each DEM has an ambiguous surface representation, but GCD could identify changes over several temporal surveys. DEM uncertainty comprises two situations: Complicated and straightforward. The minimum level of detection (MinLoD) is used for simple uncertainty, while spatial variability through a fuzzy inferencing system (FIS) is used for complicated uncertainty. For the Difference of DEM (DoD) analysis, GCD provides tools for quantifying and propagating uncertainties in each DEM [51].
DoD analysis have been used to detect changes in the soil's topography over time and quantify the eroded or deposited volumes of sediment [52]. The digital elevation models obtained from each survey are differenced to provide each DoD result. DoDs are used to calculate the net volumetric change in a reach as it progresses over time. From a geomorphic standpoint, it describes the shift in storage in terms of a sediment budget (due to erosion and deposition). As indicated in Equation (4), DEMs from multiple periods are subtracted from one another to form a morphological change raster, where t1 is the initial time and t2 is the subsequent time of DEM acquisition.
The process starts with selecting the filtered DEM as the input for the DEM survey, and both data are loaded into the GCD software. At least two DEM datasets were required as input, with both datasets being either new or old DEM based on the epoch: New DEM (low tide DEM of 2017) and old DEM (low tide DEM of 2016). The intersection of new and old surface data extents functioned as the area of interest during the change detection analysis of both datasets. The algorithm then selects one of three approaches for evaluating uncertainty: (i) Superficial minimum level of detection (minLoD), (ii) propagated errors, and (iii) probabilistic thresholding. Although this is identical to using the raster calculator in ArcGIS, the GCD software helps apply a threshold to the data. This threshold is helpful since the most commonly used method for managing DEM uncertainties specifies a minimum detection threshold (minLoD) to identify actual surface changes from inherent noise. In this study, minLoD was chosen over propagated errors and probabilistic thresholding because of the simple DEM input using the same threshold, which results in similar errors in the DEM surface. The value for the MinLoD threshold for all ROIs is 0.1 m, according to the value of RMSE for UAV-DSM versus the static GNSS accuracy assessment. Figure 4 displays the flow of the method for generating DEM change detection using GCD stand-alone software.

SfM-MVS Result
All epochs have been analysed for RMSE error to determine the accuracy of each result based on the previous GCP table. For low-tide observations in the year 2016, Table 1 quantified the total for N error as 1.380 cm, E error is 1.003 mm, H error is 0.263 mm, and NE error is 1.723 mm. Meanwhile, for the year 2017, the low-tide observation displayed N error at 1.309 mm, E error at 1.463, H error at 0.232 mm, and NE error at 1.921 mm for all 8 GCPs (Table 1). Last but not least, Table 1 also summarises and displays the side-by-side accuracy assessment of GCP for both epochs. For low-tide observations in the year 2016, the table shows the minimum error at 0.339, the maximum error at 2.612 mm, the ME error at 1.697 mm, and the SDE error at 0.686 mm. Meanwhile, for low-tide observations in the year 2017, 8 GCPs were analysed and showed the minimum error at −0.131 mm, the maximum error at 2.372 mm, the ME error at 1.660 mm, and the SDE error at 0.652 mm.  Table 2 also shows the assessment of camera calibration parameters for UAV data based on the focal length (F), principal point coordinates (C x and C y ), radial distortion polynomial coefficients (K1, K2, and K3), and tangential distortion coefficients (P1, P2, P3, and P4). The

SNERL Filtering Result
The filtering process was comprised of UAV DSM data in 2016 and 2017. For the epoch of 2016, eight regions of interest (ROI) were established in the study area. In ROI 1, the average height before filtering was between 1.5 m and 25 m, and after the filtering process, the average height was reduced to an average of between −3.7 m and 6.3 m. As shown in Figure 6, the average height before filtering was between −5 m and 29 m, and after the filtering process, the average height was reduced to an average of −1.0 m and 6.0 m. For ROI 3, the average height before filtering was between −10 m and 67 m but decreased to −1.0 m and 5.6 m after the filtering process. Meanwhile, in ROI 4, the average height changed from an average of −37 m to 48 m before filtering into −5.0 m and 6.7 m after filtering. Before filtering, the average height in ROI 5 was between −74 m and 50 m, but it was reduced to between −24 m and 15 m after filtering. In ROI 6, the height changes from an average of −10 m to 30 m (before filtering) into −4.1 m and 12.5 m (after filtering). The following ROI (ROI 7) decreases in average height from −10 and 50 m to −2.5 m and 12.5 m after the filtering process. The last ROI (ROI 8) showed a decreasing average height from −10 m and 35 m (before filtering) to −2.5 m and 10.0 m (after filtering).
In 2017, a similar number of ROIs were applied to the study area. In ROI 1, the average height before filtering was between −29 m and 50 m, and after the filtering process, the average height was reduced to an average of between −6.7 m and 12.5 m. In ROI 2, the average height before filtering was between −10 m and 30 m, and after the filtering process, the average height was reduced to an average of 0.2 m and 5.0 m, as shown in Figure 6. For ROI 3, the average height before filtering was between 0 m and 80 m but decreased to 0.0 m and 4.8 m after the filtering process. Meanwhile, in ROI 4, the average height changed from an average of 2 m to 60 m before filtering into 0.3 m and 7.0 m after filtering. In ROI 5, the average height before filtering was between −19 m and 50 m and was reduced to −5.0 m and 14.5 m after the filtering process. In ROI 6, the height was changed from an average of 1 m to 40 m (before filtering) to 0.2 m and 12.5 m (after filtering). The following ROI (ROI 7) decreased in average height from −10 m to 29 m in −2.5 m and 7.7 m after filtering. The last ROI (ROI 8) showed a decreasing average height from −4.0 m to 39 m (before filtering) into −10.0 m and 7.5 m (after filtering). The histogram of the average height comparison between the unfiltered and filtered UAV DSM is displayed in Appendices A and B. After the filtered DEM was generated, surface elevation between epochs was compared using the Global Mapper v22.0 software path profile tools. The path profile line was performed simultaneously on both epochs of filtered DEM. The comparison includes the deviation in height between two epochs of filtered DEM and the dominant height of each ROI. Based on Figure 7, ROI 2 showed an increment of DSM (unfiltered) and DEM (filtered) between 2016 and 2017, and a similar situation also happened at ROI 3, 5, 6, and 8. On the other hand, the contrast situations occurred at ROI 1, 4, and 7, where the decrement dominated the surface. The quality of the DEM using the SNERL filtering algorithm was evaluated and quantified based on the vertical accuracy assessment. The vertical accuracy assessment of DSM post-processed using the SNERL algorithm was compared to GNSS data, which resulted in an RMSE error of 0.089 m, a mean of 0.039 m, and a standard deviation (SD) of 0.322 m. So, it can be concluded that the SNERL filtering algorithm accuracy was less than 1 m.  Figure 8 shows the graphical output of surface elevation changes beneath the mangrove canopy, which includes the graphical view of areal changes, volumetric changes, and vertical averages calculated with GCD software. Furthermore, the GCD map was shown with a threshold value of 20 to provide a quick look at erosion (surface lowering) and accretion (surface raising) in each region of interest (ROI).  In ROI 1 (a), the GCD map showed an almost balanced proportion of erosion and accretion for areal changes, with accretion slightly higher than erosion. In comparison, ROI 1 (b) and 1 (c) also displayed more accretion over erosion in their histograms for volumetric changes and vertical changes. The GCD map in ROI 2 (a) and (b) showed an imbalanced proportion of erosion to accretion for areal and volumetric changes, and a similar thing appeared in the change bars in ROI (c). For ROI 3 (a) and (b), erosion appeared more dominant than accretion, and vertical changes also increased to a significant value. A similar situation occurred in ROI 3, 4, 7, and 8, demonstrating that erosion was more dominant than accretion, but this was not the case in ROI 6, where accretion was far more significant than erosion.

GCD Result
The histogram and change bars in Figure 9 showed ROI 1 (a) experienced almost a balanced proportion of erosion and accretion for areal changes, with accretion slightly higher than erosion, while ROI 1 (b,c) also displayed more accretion than erosion in their histograms for volumetric changes and vertical changes. The histogram and change bars of ROI 2 (a,b) in Figure A1 show an imbalanced proportion of erosion to accretion for areal and volumetric changes, and a similar thing appeared on the change bars of ROI 2 (c). For ROI 3 (a,b) rates, erosion appeared more dominant than accretion, and vertical changes also increased to a significant value, as displayed in Figure A2. A similar situation occurred at ROI 4 and 7 in Figures A3 and A5, while ROI 8 in Figure 6 demonstrated that erosion was more dominant than accretion, but this was not the case in ROI 6 ( Figure A5), where accretion was far greater than erosion. Vertical changes directly affected the imbalance between erosion and accretion, which was proved by the change bars in all ROIs (c) that propagated the values of the average depth of lowering, the average depth of raising, and the average total thickness difference.  GCD analysis has quantified the geomorphological change rates, which comprised volumetric, vertical change rate, and percent imbalanced at each ROI on the Kilim River. Based on Table 2

Geomorphological Changes underneath Mangrove Canopy at Kilim River
As shown by the histogram, change bars, and tabulated information from the GCD result, each ROI has a different surface elevation rate based on the filtered output of the SNERL algorithm. GCD software utilised the Differential of DEM (DoD) method in estimating erosion (surface lowering), accretion (surface raising), and vertical averages for each ROI since the region was not connected to each other and eased the GCD processing since the area became smaller. The indicator for an ROI to have either erosion or accretion depends on the total net volume difference. At the same time, regarding the vertical average, it averages the average net thickness of the difference in the area with detectable changes.
The finding in Table 1 shows that only ROI 1, 5, and 6 were dominated by accretion while the rest were dominated by erosion. The GCD map displayed ROI 1, 5, and 6, showing accretion-based colour. The histogram and change bars also show a similar output in which accretion appears most in specific ROIs. The volumetric change rate for ROI 1, 5, and 6 depicted a positive value, which signifies that the area was experiencing accretion. Specifically, based on numerical results, ROI 1 tabulated 0.101 cm 3 while ROI 5 and 6 recorded 0.566 cm 3 and 1.248 cm 3 of volumetric rates. In the meantime, the rest of the ROIs displayed negative volumetric rates, including ROI 2 (−0.640 cm 3 ), ROI 3 (−0.338 cm 3 ), ROI 4 (−0.003 cm 3 ), ROI 7, (−0.499 cm 3 ), and ROI 8 (−2.469 cm 3 ). The finding summarises that ROI 5 experienced the highest volumetric accretion (surface raising), while the highest erosion (surface lowering) was identified in ROI 8.
Another signature of the accretion-dominated region was shown by the areal changes and vertical average rates of ROI 1, 5, and 6. The findings in Table 1  According to the findings, ROI 6 had the highest vertical accretion (surface raising), while ROI 3 had the highest erosion (surface lowering).

Short and Long-Term Impacts of Geomorphological Changes on Mangrove Surface and the Interaction to Sea Level Rise
The balanced proportion of erosion (surface lowering) and accretion (surface raising) in volumetric, areal, and vertical changes is essential for ecosystem equilibrium, especially in dynamic ecosystems like mangroves. Ocean changes through tides directly influence any changes to the land surface in mangrove areas. Like the ocean current, which is constantly changing due to tides, the balance of tidal flux, which keeps the ocean current in an equilibrium state, somehow experiences unusual changes due to sea-level rise. Mean sea levels are rising due to both thermal expansions of seawater as temperatures rise due to climate change and the melting of polar ice caps and other land ice, which adds additional water to the sea [53]. Sea level rise causes an unbalanced flow of low and high water transitions, affecting coastal regions exposed to ocean water.
Mangrove soils may be unable to increase at the same pace as the local rate of sealevel rise, resulting in the loss of trees in the mangrove area's lower parts and along the seaward edge. Mangroves are likely to invade landward areas that are currently within the tidal frame, assuming adequate substrate and topography exist. Because of the deeper water in mangrove areas, waves may penetrate further into the area, resulting in erosion, particularly along the seaward side. Positive and negative feedback between changes in sea level and surface rates may be influenced. An elevation surplus occurs when surface elevation rises faster than the sea level, whereas an elevation deficit occurs when the sea level rises faster than the mangrove surface.
When the sea level rises or land subsides, the volume of accommodation space increases because the height difference between the substrate and the mean sea level rises. If soil inputs are high enough, this volume can now be filled with soil, allowing the soil surface to rise until the newly formed accommodation space is filled. Soil inputs include organic and inorganic sediments, as well as underground roots. Mangroves would remain in their ideal tidal frame, i.e., between the mean sea level and high tide, as the height of the mangrove soil surface rises. Without such an increase in the soil surface height, the mangrove surface could descend below the mean sea level, causing stress and physical harm to the mangrove trees. The relative height of the mangrove surface remains constant within the tidal range if the change in soil surface height closely matches the change in sea level.
According to Malaysia's sea-level rise phenomenon, especially in the Malacca Straits (west coast of Peninsular Malaysia), the trend shows an increasing rate of 4.95 ± 0.15 mm yr −1 based on a study by [54]. A similar study of sea-level rise trends by [55] discovered the mean increase in the Malaysian coastline region at 4.47 ± 0.71 mm yr −1 based on satellite altimetry and corrected tidal data. Although the values seem tiny, the long-term effects, especially at lower elevations or coastal areas, will be significant in the future. Things appear to become worse if the rate of imbalanced surface elevation changes in the mangrove area escalates rapidly in the future. Frequent surface lowering (erosion) causes the elevation to not adapt to tidal inundation during high water intrusion.
In this study, a few ROIs, such as ROIs 2 and 5, showed a balance of erosion and accretion, as illustrated in Table 1. Both ROIs also numerically displayed a fair number of erosion and accretion rates, which signifies an equilibrium in geomorphological changes within one year. Suppose the situation remains the same in the future. In that case, both ROIs will have a stable and healthy ecosystem, even if the ocean's unusual changes such as the sea-level rise or natural events such as typhoons or storms hit the area. This is because, in regions with very high rates of accretion, mangrove soil surfaces may increase faster than the local rate of sea-level rise, causing terrestrial species to infiltrate landward areas and progradation to occur (i.e., new land is formed seaward of the current mangrove area, which mangroves then colonise). This is most likely to happen along the deltas of big rivers that carry enormous amounts of sediment to the coast. Sea level rise rates might well be matched by an increase in mangrove soil surface elevation, allowing mangroves to remain on the same site while potentially colonising new landward areas with suitable substrate and terrain.

Conclusions
In conclusion, UAV photogrammetry data could be used as the primary data for estimating surface elevation underneath a mangrove canopy. A geospatial approach comprised of SfM-MVS, vegetation index (VI) segregation, and an SNERL filtering algorithm is efficient in generating UAV DEM below the mangrove canopy at the closest level to the terrain. The other approach, based on a geomorphological method known as GCD analysis, can evaluate the rate of surface elevation changes using the DoD technique and can discover interactions between mangrove surface elevation changes and sea-level rise. Hopefully, more research will be conducted to improve and enrich our knowledge of the mangrove ecosystem and prepare for a mitigation plan. The GCD map and the rates of surface elevation changes at the Kilim River enabled authorities such as the Langkawi Development Authority (LADA) and the Department of Drainage and Irrigation (DID) to understand the situation much deeper. This allows the preparation of a mitigation plan to avoid unbalanced surface elevation changes that could bring unstable equilibrium to the mangrove ecosystem in the future.          The content of Appendix B was a continuation of Figure 7. Figure A8 was referred to an average height comparison between filtered and unfiltered for ROI 1, while Figure A9 for ROI 3, Figure A10 for ROI 4, Figure A11 for ROI 5, Figure A12 for ROI 6, Figure A13 for ROI 7, and Figure A14 for ROI 8.