This would be better explained in the next section when the difference of digital elevation is calculated.

## 4.2 Digital Terrain Models

Regarding the comparisons between the DJI Terra DTM, the Agisoft Metashape DTM, and the EAGLE FVG, the analysis showed that it was favourable to use the DJI Terra DTM to compare with the EAGLE FVG DTM. Regarding the comparisons between the DJI Terra DTM, the Agisoft Metashape DTM, and the EAGLE FVG, the analysis showed that it was favourable to use the DJI Terra DTM to compare with the EAGLE FVG DTM. This was due to having little difference between the Agisoft Metashape and DJI Terra DTM, but the Agisoft Metashape DTM interpolated areas that were not acquired in the point cloud. Also, the DJI Terra DTM had an advantage of containing the data from all the points instead of a subsampled point cloud.

Consequently, a Difference of Digital Elevation (DoD) was performed between the DJI Terra DTM of the year 2024 and the EAGLE FVG DTM of the year 2017.

<figure><img src="image_70.png" /><figcaption>Figure 2: Difference of Digital Elevation (DoD) Between DJI Terra 2024 and EAGLE FVG 2017 DTMs.</figcaption></figure>

Clearly, there were differences between the two digital terrain models particularly highlighted in the pinnacles of the terrain and tops of the rocks in that inclined area. When highlighting the $μ + 3σ$ and $μ - 3σ$ areas we see that areas contoured in blue show an elevation decrease whereas areas with a red contour

signify an elevation increase. The areas that remained in teal blue and showed no contour experienced minimal change. Rockfall can be tied to the high difference between 2024 and 2017. Based on figure, the areas that lost elevation ended up depositing the elevation lower in the basin which is typical of rockfall situations. There were some areas shown in yellow that were trees not cleaned properly with the DJI Terra software and was considered as ground points while clearly, they were not.

## 4.3 Land Cover Classification Using a Deep Learning Model

Upon completion of the model training, the model began to recognize patterns and identify the classes found in the control images based on their classified masks. After this step, the Rio Cucco orthophotos were introduced. Both datasets which refer to the Rio Cucco 2024 orthophoto and the Rio Cucco 2017 orthophoto were tested to check the effectiveness of the model at identifying the classes in an unintroduced image. The semantic segmentation with this deep learning model in its current form produced mixed results in terms of its classification as it detects vegetation, and bare land well but struggles with accurately detecting roads and buildings in the new image. Based on Figure 3, it is evident that the model works very well at predicting the classes based on the images that it has gotten used to (its original dataset). In Figure 4 however, it can be observed that the model recognizes the shapes and is able to detect well the land (represented in light purple) while mistaking vegetation as water (orange instead of purple). The model produced mixed results with buildings as it was able to recognize their location but did not have proper edge detection for them. The testing of this model yielded mixed results in terms of its segmentation abilities.

<figure><img src="image_71.png" /><figcaption>Figure 3: Deep Learning Model Training Results.</figcaption></figure>
<figure><img src="image_72.png" /><figcaption>Figure 4: Deep Learning Model Results with Rio Cucco Orthophoto.</figcaption></figure>

<figure><img src="image_73.png" /><figcaption>Figure 5: Rockfall in Sloped Areas in Rio Cucco (2024).</figcaption></figure>

# 5 DISCUSSIONS

With respect to the Rio Cucco basin, comparing past data with present data is quintessential to understand the differences both spatially and temporally. Based on the findings of the data when comparing the data of the survey in 2024 and the FVG 2017 data taken from EAGLE FVG, the differences that existed were evident especially when conducting the DoD. The DoD pointed to an increase of elevation in certain areas and a decrease in elevation in others. The increase of elevation was found further down the slope and the decrease of elevation was found higher in the slopes. Several explanations could be tied to why this variation of elevation occurred. The variation could have been tied to possible rockfall events that happened and lead to sediment deposition by removing debris from higher up the slope to be deposited further down the slope. However, the variations could have also been due to errors in the model itself. To clarify which case most likely caused this variation between the two DTM, an examination of the drone images for the June 2024 survey was done. After careful examination, rock debris was noticed downstream in the areas consistent with the blue and red colour contours. Fig.5 shows an image of rocks deposited near a slopped area in the basin. Based on the location of the rock deposition, it can be inferred that the elevation differences in the DTM were most likely due to a rockfall event.

Based on the history of the area and the nature of its geology, it is expected that rockfall events will continue to occur in the area especially because of the calcareous and dolomitic set of rocks that are subject to thermoclastism and cryoclastism. With respect to the deep learning model, the testing showed mixed results in terms of its power at classifying unintroduced images. Other steps to increase the effectiveness of the model would include data augmentation techniques and introducing the digital

terrain models or the digital surface models in the training phase so that the model can segment the area with respect to certain elevation patterns differentiating between water, road, buildings, vegetation and other classes.

# 6 CONCLUSIONS

In summary, the Rio Cucco basin is subject to several changes over the years and needs to be monitored in terms of natural hazard studies. The geological parameters as well as the morphology makes the Rio Cucco basin particularly liable to the effects of climate change. The nature of its rock formations and slopes which are prone to weathering due to abrupt changes in temperature highlights the importance of consistent monitoring. With the advancement of geomatic techniques, artificial intelligence (AI) models, and internet of things (IoT) sensors, perhaps a predictive approach could be pursued to model where and when natural hazards in the area would most likely occur. Consequently, preventative measures could be taken to reduce the impact of these events that are likely to become more common with the emergence of climate change.

# ACKNOWLEDGEMENTS

This study was carried out within the «GeoMORPHoMetry through Scales for a resilient landscape» project – funded by European Union – Next Generation EU within the PRIN 2022 program (D.D. 104 - 02/02/2022 Ministero dell'Università e della Ricerca). This manuscript reflects only the authors' views and opinions and the Ministry cannot be considered responsible for them.
A special thanks is to be extended to our team members at CNR-IRPI Padova who aided us with the information regarding the Rio Cucco basin as well as helping us with the June 2024 survey.

## REFERENCES

Agisoft LLC. (n.d.). *Agisoft Metashape*.
https://www.agisoft.com

Belcore, E., Piras, M., & Pezzoli, A. (2022). Land Cover Classification from Very High-Resolution UAS Data for Flood Risk Mapping. *Sensors*, 22(15), 5622. https://doi.org/10.3390/s22155622

Belcore, E., Piras, M., & Wozniak, E. (2020). Specific Alpine Environment Land Cover Classification Methodology: Google Earth Engine Processing for Sentinel-2 Data. *The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences*, XLIII-B3-2020, 663–670. https://doi.org/10.5194/isprs-archives-XLIII-B3-2020-663-2020

Cavalli, M., Trevisani, S., Comiti, F., & Marchi, L. (2013). Geomorphometric assessment of spatial sediment connectivity in small Alpine catchments. *Geomorphology*, 188, 31–41. https://doi.org/10.1016/j.geomorph.2012.05.007

DJI. (n.d.). *DJI Terra - Make the World Your Digital Asset* - DJI. https://enterprise.dji.com/it/dji-terra

Girardeau-Montaut, D. (n.d.). *CloudCompare - Open Source project*. https://www.danielgm.net/cc/

Huggel, C., Clague, J. J., & Korup, O. (2012). Is climate change responsible for changing landslide activity in high mountains? *Earth Surface Processes and Landforms*, 37(1), 77–91. https://doi.org/10.1002/esp.2223

Indian Space Research Organization. (2023, May 31). *Land cover Classification: Bhuvan Satellite Data*. Kaggle. https://www.kaggle.com/datasets/khushiipatni/satellite-image-and-mask

International Society for Photogrammetry and Remote Sensing. (2024, March 20). *Urban Segmentation - ISPRS*. Kaggle. https://www.kaggle.com/datasets/aletbm/urban-segmentation-isprs

Kaggle. (n.d.). *Kaggle: your machine learning and data science community*. https://www.kaggle.com/

Roia Foundation. (2020, May 29). *Semantic segmentation of aerial imagery*. Kaggle. https://www.kaggle.com/datasets/humansintheloop/semantic-segmentation-of-aerial-imagery

Sangelantoni, L., Gioia, E., & Marincioni, F. (2018). Impact of climate change on landslides frequency: The Esino river basin case study (Central Italy). *Natural Hazards*, 93(2), 849–884. https://doi.org/10.1007/s11069-018-3328-6

Schneiderbauer, S., Fontanella Pisa, P., Delves, J. L., Pedoth, L., Rufat, S., Erschbamer, M., Thaler, T., Carnelli, F., & Granados-Chahin, S. (2021). Risk perception of climate change and natural hazards in

global mountain regions: A critical review. *Science of The Total Environment*, 784, 146957. https://doi.org/10.1016/j.scitotenv.2021.146957

Semnani, S. J., Han, Y., Bonfils, C. J., & White, J. A. (2025). Assessing the impact of climate change on rainfall-triggered landslides: A case study in California. *Landslides*. https://doi.org/10.1007/s10346-024-02428-0

Zhang, C., Sargent, I., Pan, X., Li, H., Gardiner, A., Hare, J., & Atkinson, P. M. (2019). Joint Deep Learning for land cover and land use classification. *Remote Sensing of Environment*, 221, 173–187. https://doi.org/10.1016/j.rse.2018.11.014
SHORT PAPERS
<figure><img src="image_74.png" /></figure>
# Differential Interferometric Synthetic Aperture Radar-Based Landslide Monitoring: A Case Study of Wayanad, India

Manvi Kanwar<fnref n="1" /> ID a and Samsung Lim<fnref n="2" /> ID b

**Keywords:** Rainfall-Induced Landslides, D-Insar, Sentinel-1, Remote Sensing.

**Abstract:** On July 30th, 2024, the Wayanad region of Kerala, India, experienced a devastating landslide triggered by intense monsoon rainfall. This event highlighted the region's need for effective landslide monitoring and early warning systems. Differential Interferometric Synthetic Aperture Radar (DInSAR) is a powerful and efficient remote sensing technique for detecting ground deformation and monitoring landslides. This study utilizes Sentinel-1 SAR images to monitor the Wayanad landslide by processing and analyzing SAR data using the DInSAR technique to identify ground displacement patterns. The results demonstrate the effectiveness of DInSAR in capturing pre- and post-event deformations, offering valuable insights into the landslide dynamics. This study underscores the potential of SAR-based real-time landslide monitoring and risk mitigation in landslide-prone regions.

## 1 INTRODUCTION

Wayanad, a district in Kerala's Western Ghats, India, is highly susceptible to landslides, a hazard exacerbated by both natural and anthropogenic factors. The region's topography, characterized by elevated slopes and a soil cover consisting mainly of boulders, colluvium, and laterite, combined with the region's propensity for heavy monsoon rainfall, makes it particularly vulnerable to landslide events (Kuriakose et al., 2009). On July 30th, 2024, a massive rainfall event which brought over 14 cm of rain in a single day — 493% above the average — triggered multiple landslides in Chooralmala, Mundakkai, and Vellarimala villages (Mishra, 2024). These landslides are not isolated; Wayanad has a history of similar events, such as the over 3,000 landslides reported during the 2018 floods and the multiple landslides in Puthumala in 2019 (National Remote Sensing Centre, 2023).

While natural factors such as topography and rainfall primarily contribute to landslides, anthropogenic activities have increasingly significantly exacerbated these events. Studies following the 2018 and 2019 floods identified extensive land-use changes, particularly the

establishment of tea and rubber plantations and the construction of houses and other infrastructure in landslide-prone areas (Premlet B, 2019). These activities obstruct natural drainage systems and saturate the soil, making it more susceptible to landslides during heavy rains. Moreover, the velocity of rivers such as Iruvazhinji, which runs through steep terrain before descending rapidly in elevation, can intensify the impact of floods and landslides.

Kerala's Western Ghats, including Wayanad, have long been recognised as vulnerable. In 2011, the Madhav Gadgil-led Western Ghats Ecology Expert Panel recommended that a sizeable portion of the Western Ghats, including areas in Wayanad, be designated as Ecologically Sensitive Zones (ESZ) (Gadgil Madhav, 2011). This designation would have restricted industrial activities, mining, and construction in these areas. However, these recommendations were not implemented, exposing the region to the ongoing risks of landslides and other environmental disasters.

In areas affected by excessive rainfall, using optical images to determine landslide extent can be severely hindered by cloud cover and adverse weather conditions. However, microwave-based Synthetic

<footnote n="1">¹School of Civil and Environmental Engineering, University of New South Wales, Australia</footnote>

<footnote n="2">²Associate Professor, School of Civil and Environmental Engineering, University of New South Wales, 2052, Australia
s.lim@unsw.edu.au</footnote>

<footnote>a https://orcid.org/0009-0006-6680-6068</footnote>

<footnote>b https://orcid.org/0000-0001-9838-8960</footnote>
Aperture Radar (SAR) images offer a significant
advantage in such scenarios. Unlike optical imagery,
SAR can penetrate clouds and provide consistent data
irrespective of weather conditions, making it an
invaluable tool for monitoring large and inaccessible
areas quickly and reliably. This advantage is evident
in recent studies (Debevec Jordanova et al., 2024; Lau
et al., 2024), demonstrating DInSAR's utility in
mapping and monitoring landslides. Similarly, this
research applies DInSAR to assess the Wayanad
landslide, a rainfall-induced event, by leveraging the
phase differences between SAR images captured at
different times. While other studies have used
advanced techniques like SBAS or integrated in-situ
monitoring, this study focuses on rapid assessment of
ground deformation using Sentinel-1 SAR data,
which is freely accessible. This approach offers
valuable insights into hazard assessment and disaster
management, allowing response teams to prioritise
the most vulnerable areas and overcome the
limitations of traditional optical methods.

## 2 STUDY AREA

Wayanad District, located in the northeastern part of Kerala, India, is known for its diverse geography and significant natural beauty. The district is characterized by its hilly terrain, with elevations ranging from 700 to 2,100 meters above sea level and is part of the Western Ghats Mountain range. Wayanad's landscape includes dense forests, lush green tea and coffee plantations, and numerous rivers and waterfalls. The region experiences a tropical monsoon climate, with an average annual rainfall of 2,500 to 3,000 millimetres (India Meteorological Department, 2024), contributing to its rich biodiversity and dense vegetation. It is bordered by the districts of Kannur to the west, Kozhikode to the south, and Malappuram to the southwest. To the east, Wayanad shares its boundaries with the state of Karnataka. The district's geographical and climatic conditions make it highly susceptible to landslides, particularly during the monsoon season, which highlights the importance of effective monitoring and disaster management in the region. The present study focuses on one sub-swath of a Sentinel-1 SAR image, which encapsulates approximately 20,825 km² of area, as shown in Figure 1.

<figure><img src="image_75.png" /><figcaption>Figure 1: Study Area.</figcaption></figure>
<table><caption>Table 1: Data Characteristics of the Pre-Event SAR Images Using Stack Overview in Snap.</caption><thead><tr><th>Acquisition date</th><th>Track</th><th>Orbit</th><th>Temporal Baseline (days)</th><th>Perpendicular Baseline (m)</th><th>Modelled Coherence</th></tr></thead><tbody><tr><td>8 July 2024 (reference)</td><td>165</td><td>54662</td><td>0</td><td>0</td><td>1</td></tr><tr><td>20 July 2024 (secondary)</td><td>165</td><td>54837</td><td>12</td><td>-123</td><td>0.89</td></tr></tbody></table>

<table><caption>Table 2: Data Characteristics of the Post-Event SAR Images Using Stack Overview in Snap.</caption><thead><tr><th>Acquisition date</th><th>Track</th><th>Orbit</th><th>Temporal Baseline (days)</th><th>Perpendicular Baseline (m)</th><th>Modelled Coherence</th></tr></thead><tbody><tr><td>1 August 2024 (reference)</td><td>165</td><td>55012</td><td>0</td><td>0</td><td>1</td></tr><tr><td>13 August 2024 (secondary)</td><td>165</td><td>55187</td><td>12</td><td>-142</td><td>0.87</td></tr></tbody></table>

# 3 METHODOLOGY

## 3.1 Data Acquisition

This research leverages SAR datasets obtained from the Sentinel-1A and 1B missions, launched in 2014 and 2016. These satellites are equipped with C-band SAR technology, which offers open-access data capable of capturing Earth's surface imagery under all weather conditions, day and night. Sentinel-1 provides data well suited for medium-resolution applications, with a revisit period of 12 days, allowing for extensive area coverage. The satellites operate in the Interferometric Wide (IW) swath mode, utilizing Terrain Observation by Progressive Scans (TOPS) to achieve a ground coverage of approximately 250 km. The geometric resolution of the SAR products is 5 m by 20 m in range and azimuth, respectively.

For the application of the D-InSAR method, two images were acquired before the landslide event happened, and two were after the Wayanad landslide to demonstrate the changes in the landform. Data processing was performed using SeNtinel Applications Platform (SNAP) version 10.0.0, a free software provided by the European Space Agency (ESA). The SAR data was sourced from the Alaska Satellite Facility (ASF) and provided by the European Space Agency (2023) as part of the Sentinel-1 mission. Each Sentinel-1 Single Look Complex (SLC) product consists of three Interferometric Wide Swaths (IW1, IW2, and IW3) and nine bursts. The IW3 sub-swath was selected for this study as it encompasses the area of interest. The acquired images were in both VH and VV polarization modes. Tables 1 and 2 present the features of the pre-event and post-event datasets, respectively.

The perpendicular baseline was kept within 200 meters to maintain data quality. This restriction helps ensure accurate interferometric measurements and reliable results in SAR data processing. The perpendicular baseline is a critical parameter in SAR interferometry, affecting the interferograms' spatial resolution and coherence (Li & Bethel, 2008). Keeping it within this limit helps minimize errors and maintain the overall quality of the interferometric data (Zhang et al., 2005).

## 3.2 Preprocessing

Preprocessing of the Sentinel-1 SAR data involves several steps to ensure the accuracy of the DInSAR analysis:

**Orbit File Application:** Accurate orbit data are applied to the SAR images to correct for satellite position errors. Precise orbit files provided by ESA were utilized to ensure accurate geolocation and alignment of the Sentinel-1 imagery.

**Radiometric Calibration:** This step involves converting the SAR image data from digital numbers (DN) to radar backscatter intensity, which is essential for quantitative analysis.

**Multilooking:** The SAR images are multilooked to reduce speckle noise inherent in SAR data by averaging adjacent pixels. This study used a multi-look factor of 4 in range and 1 in azimuth, reducing noise while maintaining spatial resolution suitable for analyzing landslide deformation.
**Co-registration:** Pre and post-event SAR images are precisely aligned to ensure accurate phase difference calculations. This step minimizes phase errors due to misalignment.

**Interferogram Generation:** An interferogram is generated by computing the phase difference between the co-registered SAR images. This interferogram represents the ground displacement in the radar line of sight (LOS).

## 3.3 DInSAR Processing

The core of the DInSAR technique involves the following steps:

**Phase Unwrapping:** The interferometric phase is wrapped between -π and π, necessitating a phase unwrapping process to retrieve the actual displacement values. The Statistical-cost Network-flow Algorithm for Phase Unwrapping (SNAPHU) algorithm is commonly used. This study employed SNAPHU with a coherence mask to restrict unwrapping to areas with a coherence threshold of 0.3, ensuring accurate displacement retrieval in high-coherence regions.

**Topographic Phase Removal:** The influence of topography on the interferometric phase is removed using a Digital Elevation Model (DEM), such as the Shuttle Radar Topography Mission (SRTM) DEM, to isolate the phase related to ground deformation. A 30 m resolution SRTM DEM was used in this study to ensure the precise removal of the topographic component in the Wayanad region, characterized by hilly terrain.

**Temporal and Spatial Filtering:** Temporal and spatial filtering techniques reduce noise and enhance the signal related to ground deformation. With an adaptive window size of 32 x 32 pixels, the Goldstein filter preserves deformation-related signals while reducing noise.

**Displacement Map Generation:** The unwrapped phase is converted into displacement maps representing the ground movement in the radar LOS direction. These maps are then geocoded to produce spatially accurate deformation maps. This study projected the displacement maps onto the WGS84 coordinate system for spatial analysis and visual interpretation of ground deformation pattern.

# 4 RESULTS AND DISCUSSIONS

The results obtained after the dataset's preprocessing demonstrate the landform change before and after the event. The LOS displacement and coherence maps were generated. This analysis helps us showcase the land movement that happened towards or away from the radar line of sight. P et al. (2025) used machine learning models to determine the landslide susceptibility in the Wayanad district and found that the western part was most susceptible to landslides. Thus, the coherence and displacement values of these areas are discussed.

The observation from the Maps in Figures 2 to 5 from the estimated coherence suggests several important conclusions:

## 4.1 Coherence Analysis

**Before the Landslide:** The pre-landslide surface stability of the Wayanad region is evident from the coherence map in Figure 2, where 57.8% of the pixels exhibit coherence values above 0.6. This indicates predominantly stable ground. 25.3% of the region displays high coherence values (0.8–1.0), typically associated with undisturbed surfaces like rocky terrain or dense vegetation. Such areas reflect minimal displacement or deformation.

Conversely, regions with coherence values below 0.4, constituting 15.7% of the pixels, primarily occur in riverbeds and forested areas, as seen from the histogram in Figure 3. These zones experience temporal decorrelation due to vegetation growth, soil moisture variations, or minor ground movements, which are unstable precursors. Some areas' relatively lower coherence values (0–0.53) suggest potential surface instability or subtle deformations, even before the landslide event.

**After the Landslide:** Following the landslide triggered by intense monsoon rainfall, significant displacement and surface disruption are observed, particularly in areas near **Mundakkai Hills, Punchrimattom, Mundakkai**, and **Chooralmala** (Figure 4). These locations were selected for monitoring due to their proximity to steep slopes, historical landslide occurrences, and the presence of vulnerable settlements and infrastructure.

The coherence values in these regions decreased notably, reflecting ground instability and displacement. At **Mundakkai Hills**, coherence dropped from an average of 0.75 (pre-event) to 0.42 (post-event). Similarly, **Punchrimattom** experienced
a reduction from 0.81 to 0.46, and **Mundakkai** decreased from 0.84 to 0.48. The most severe displacement occurred near **Chooralmala**, where coherence dropped from 0.77 to 0.32, indicating significant material flow and surface disruption.

Figure 5 demonstrates only 34.2% of the pixels-maintained coherence above 0.6 after the landslide, compared to 57.8% before the event. Regions with coherence values in the range of 0.2–0.4 increased significantly, from 9.8% (pre-event) to 28.6% (post-event), indicating widespread ground displacement. Additionally, areas with coherence below 0.2, representing 12.3% of the pixels, correspond to zones of extreme surface disturbance caused by material displacement during the landslide.

The spatial analysis of displacement highlights the utility of SAR-based techniques in mapping and quantifying ground movement. The coherence maps provide a reliable measure of surface stability, with changes in coherence values reflecting the severity and extent of the landslide. These findings align with observations by (Ferretti et al., 2001), who emphasized the effectiveness of SAR coherence analysis in monitoring ground deformation following natural disasters.

From all above, it can be observed that the coherence analysis underscores the substantial alterations in surface properties caused by the landslide, with significant changes in coherence values reflecting both the extent of the ground displacement and the areas of severe surface disturbance. These results are critical for understanding the landslide's overall impact, as they provide detailed insights into the spatial distribution of ground instability and the transformation of surface conditions. The observed increase in coherence in certain regions post-event suggests the exposure of more stable and homogeneous surfaces, likely due to the compaction and redistribution of material during the landslide. This analysis highlights the event's severity and establishes a reliable foundation for developing future monitoring strategies, enabling better risk assessment, mitigation efforts, and informed planning for disaster-prone regions.

## 4.2 Displacement and Surface Stability

Figures 6 and 7 illustrate the LOS displacement maps of the Wayanad region before and after the landslide event, respectively. These maps provide critical insights into the extent and movement of the landslide as detected by the radar system. The LOS displacement, measured in the radar line-of-sight

direction, highlights surface movement both towards and away from the radar sensor.

**Before the Landslide:** The displacement map before the event, shown in Figure 6, indicates stable surface conditions with relatively uniform displacement values across the region. Displacement values predominantly range between -0.0014 m and 0.0048 m, suggesting minor ongoing ground movement potentially related to environmental or geological activity. Areas such as Mundakkai Hills, Punchrimattom, and Chooralmala exhibit minimal displacement, underscoring the absence of significant pre-event deformation.

**After the Landslide:** The post-event displacement map, shown in Figure 7, demonstrates substantial changes in surface movement triggered by the landslide. Negative displacement values (presented in blue) represent movement away from the radar sensor and are prominently observed in the areas most impacted by the landslide, such as Chooralmala. Maximum negative displacement values reach approximately -0.0118 m, reflecting the material's downward and outward movement from the hillside during the landslide's initiation. In contrast, high positive displacement values (depicted in yellow) correspond to regions where the surface moved towards the radar sensor, likely indicating debris accumulation. This phenomenon is evident in the towns of Mundakkai and Chooralmala, where positive displacement values of up to 0.0118 m are observed, suggesting debris flow and subsequent compaction. The displacement maps effectively capture the dynamic nature of the landslide. The pre-event SAR images reveal small-scale displacements indicative of ongoing ground activity, while the post-event deformation map highlights pronounced changes resulting from the landslide. The stark contrast between the negative and positive displacement values represents the landslide's impact, illustrating the initiation of material flow from the hillside and the accumulation of debris in surrounding areas, including Mundakkai and Chooralmala. These results underscore the utility of InSAR-derived LOS displacement data in detecting and quantifying surface changes. The analysis demonstrates the radar system's capability to track the landslide's initiation and the subsequent debris flow, providing valuable insights into the event's impact and progression.
<figure><img src="image_76.png" /><figcaption>Figure 2: Coherence monitored before the landslide event.</figcaption></figure>

<figure><img src="image_77.png" /><figcaption>Figure 3: Pixel-based analysis of coherence before the landslide event.</figcaption></figure>

<figure><img src="image_78.png" /><figcaption>Figure 4: Coherence monitored after the landslide event.</figcaption></figure>
<figure><img src="image_79.png" /><figcaption>Figure 5: Pixel-based analysis of coherence after the landslide event.</figcaption></figure>

<figure><img src="image_80.png" /><figcaption>Figure 6: Displacement towards radar LOS before the landslide event.</figcaption></figure>

<figure><img src="image_81.png" /><figcaption>Figure 7: Displacement towards radar LOS after the landslide event.</figcaption></figure>
# 5 CONCLUSION

This study demonstrated the successful application of DInSAR using Sentinel-1 SAR data to monitor the Wayanad landslide of July 30, 2024. By analyzing four SAR images captured before and after the event, coherence and LOS displacement maps were generated, providing a detailed assessment of the landslide's extent and dynamics. The pre-event coherence map revealed stable surface conditions, with 57.8% of pixels maintaining coherence above 0.6. In contrast, the post-event coherence map indicated significant ground deformation, with coherence values above 0.6 dropping to 34.2% and regions with coherence below 0.2 increasing to 12.3%. These changes highlight the substantial disturbance caused by the landslide.

The LOS displacement map further quantified surface movement, with negative values (up to -0.011 m) reflecting downward displacement in the affected areas and positive values indicating debris accumulation. This displacement pattern revealed the spatial dynamics of the landslide, from initial material displacement to subsequent debris deposition in the Mundakkai, Punchrimattom, and Chooralmala regions.

Triggered by slope instability due to intensified rainfall, the Wayanad landslide demonstrates the utility of DInSAR for rapid and accurate assessment of rainfall-induced landslides. The technique effectively detected and quantified ground deformation, providing critical insights into the event's mechanisms and extent. These results emphasize the importance of integrating SAR-based remote sensing into landslide monitoring systems, enabling timely assessments for emergency response and long-term hazard mitigation. Future studies could incorporate advanced InSAR techniques to monitor the annual rate of slope movement and improve early warning systems in regions prone to recurring landslides. This research highlights the potential of DInSAR to enhance our understanding of landslide dynamics and support disaster management efforts in vulnerable areas.

European Space Agency (ESA). SNAP-Sentinel Application Platform, v.10.0.0. https://step.esa.int/main/toolboxes/snap

Ferretti, A., Prati, C., & Rocca, F. (2001). Permanent scatterers in SAR interferometry. *IEEE Transactions on Geoscience and Remote Sensing, 39*(1), 8–20. https://doi.org/10.1109/36.898661

Gadgil Madhav. (2011). *Report of the Western Ghats Ecology Expert Panel Part I.*

India Meteorological Department. (2024). *Cumulative Rainfall For Kerala-Sw Monsoon 2024.*

Kuriakose, S. L., Sankar, G., & Muraleedharan, C. (2009). History of landslide susceptibility and a chorology of landslide-prone areas in the Western Ghats of Kerala, India. *Environmental Geology*, 57(7), 1553–1568. https://doi.org/10.1007/s00254-008-1431-9

Lau, R., Seguí, C., Waterman, T., Chaney, N., & Veveakis, M. (2024). InSAR-informed in situ monitoring for deep-seated landslides: insights from El Forn (Andorra). *Natural Hazards and Earth System Sciences*, 24(10), 3651–3661. https://doi.org/10.5194/nhess-24-3651-2024

Li, Z., & Bethel, J. (2008). *Image Coregistration In Sar Interferometry.*

Mishra, A. (2024, August 1). *Explained with maps: How topography, human density led to devastation in Kerala's Wayanad.* The Print. https://theprint.in/the print-essential/explained-with-maps-how-topography-human-density-led-to-devastation-in-keralas-wayanad/ 2202631/

National Remote Sensing Centre, I. S. R. O. (2023). *Landslide Atlas of India (Mapping, Monitoring and R&D studies using Remote Sensing data)*. https://www.nrsc.gov.in

P, L., C, M., Mathew, A., & Shekar, P. R. (2025). Machine learning and deep learning-based landslide susceptibility mapping using geospatial techniques in Wayanad, Kerala state, India. *HydroResearch, 8*, 113–126. https://doi.org/10.1016/j.hydres.2024.10.001

Premlet B. (2019). *Landslides In Kerala 2018 Chair, Educational activities IEEE Kerala Section* (Vol. 23).

Zhang, H., Zeng, Q., Liu, Y., Li, X., & Gao, L. (2005). *The Optimum Selection of Common Master Image for Series of Differential SAR Processing to Estimate Long and Slow Ground Deformation. Master Image for Series of Differential SAR Processing to Estimate Long and Slow Ground Deformation.* http://dx.doi.org/10.1109/IGARSS.2005.1526688

# REFERENCES

Debevec Jordanova, G., Popović, Z., Yastika, P. E., Shimizu, N., Oštir, K., & Verbovšek, T. (2024). SBAS DInSAR and in situ monitoring of the Šumljak landslide (SW Slovenia) dynamics driven by rainfall and piezometric-level fluctuation. *Landslides*. https://doi.org/10.1007/s10346-024-02408-4
# Preservation and Protection of Cultural Heritage in High Tourism Areas Using GIS Technology: A Case Study of the Medieval City of Rhodes

Foteini-Pelagia Leventi<fnref n="1" />, Lemonia Ragia<fnref n="1" /> and Dorina Moullou<fnref n="1" />,<fnref n="2" />

**Keywords:** Medieval City of Rhodes, Cultural Heritage, Overtourism, Climate Change, GIS, Protection and Preservation.

**Abstract:** The Medieval City of Rhodes is one of the most famous destinations in Greece, as it attracts more and more tourists, especially during the summer season. Given that this particular site is included in the UNESCO list, its protection and preservation are of utmost importance. Through a literature review, an effort is made to analyse the phenomenon of overtourism in the area. Both quantitative and qualitative data available on the topic were used from various sources in order to accurately frame the issues and implications of the phenomenon. Finally, several proposals are provided using a mix of strategies and mitigation measures relevant to the issue, emphasizing the importance of using QGIS.

## 1 INTRODUCTION

The Medieval City of Rhodes, a UNESCO World Heritage Site, is a remarkable example of medieval architecture including structures for residential, military, commercial, public and religious use (UNESCO, 1988). Specifically, it is particularly known for the interesting dialogue between these different kinds of buildings, as it offers a glimpse into both antiquity and the complex history of the Crusades and subsequent periods of influence, with ancient remains, primarily temples, and its impressive Medieval and Ottoman structures. However, the city faces increasing challenges that threaten both its physical integrity and its cultural heritage (UNESCO, 1972) (UNESCO, 2016).

Due to its location and history, the Medieval City of Rhodes attracts an increasing number of visitors, which, despite the economic benefits, may lead to overtourism (Avdikos, 2011). As the city is filled with narrow streets and ancient and medieval buildings that were not designed to have large numbers of tourists, their deterioration is likely. Additionally, overcrowding impacts visitors' experience and leads to physical degradation of the landscape (ICOMOS, 2001-2002), (García-Hernández, Dela Calle-Vaquero, & Yubero, 2017).

The Medieval City of Rhodes, apart from overtourism, is facing crucial environmental threats due to climate change. Higher temperatures are causing cracks and erosion in the stone structures, while rising sea levels and storms threaten the city's fortifications. Furthermore, earthquakes in the region heighten the risk of damage. These threats require immediate action to protect the cultural heritage of Rhodes (European Union, 2019).

In order to address these challenges, the use of modern technology is necessary. Geographic Information Systems (GIS), particularly QGIS (Quantum Geographic Information System), provide powerful tools for analyzing spatial data and managing complex urban environments. Some of its advantages pertain to the capability of mapping, visualizing, monitoring and analyzing various aspects of the city's physical and environmental conditions, providing researchers with the chance to propose mitigation strategies and protection measures (National Geographic, n.d.).

The aim of this report is to examine the various challenges facing the Medieval City of Rhodes with a focus on the dual threats of overtourism and environmental risks driven by climate change. By understanding the important vulnerabilities of this historic region, the report seeks to explore how Geographic Information Systems (GIS) can be used

<footnote n="1">¹School of Applied Arts, Hellenic Open University, Patras, Greece</footnote>

<footnote n="2">²Ephorate of Antiquities of East Attica, Hellenic Ministry of Culture, Athens, Greece</footnote>

<footnote>leventi.felia@gmail.com, ragia.lemonia@ac.eap.gr, tmoullou@culture.gr</footnote>
to monitor, analyze and develop strategies to mitigate
these risks. Ultimately, the purpose is to propose
sustainable solutions to protect both the physical
integrity and the cultural heritage of the Medieval
City of Rhodes for future generations.

## 2 CASE STUDY

Rhodes is the largest island in the Dodecanese and is known for its long and important history dating back to ancient times. Its strategic location in the eastern Aegean Sea made it a crossroads between East and West, as well as North and South. Archaeological evidence indicates that the island of Rhodes has been inhabited since the New Stone Age, and by the 11th century BC, it was settled by the Dorians who established the city-states of Ialysos, Kamiros, and Lindos. Another important section of island history is the fact that Rhodes continued to play an important role in the ancient world, particularly through its colonies, such as Gela in Sicily (UNESCO, 1988).

The Medieval City of Rhodes, located at the northern tip of the island, became significant during the Crusades, when the Knights Hospitaller took control of the island in 1309. In particular, the island was initially developed into a fortified stronghold that used to provide medical services to Crusader knights, creating its medieval character. Specifically, the Knights were responsible for the city's architectural development, constructing fortifications, gates and other characteristic buildings that still dominate the urban landscape. Despite numerous challenges, including sieges, natural disasters and Ottoman occupation, Rhodes maintained its historical significance and integrity (UNESCO, 1988) (Luttrell, 2003).

The historical and architectural value of the Old Town of Rhodes is the reason why UNESCO included the site in its World Heritage List. For instance, its fortifications, influenced by Gothic and Ottoman architecture, are remarkable examples of medieval military construction and offer deep insight into these type of structures in the eastern Mediterranean. Therefore, it is important to preserve the original building materials of buildings and the urban organization in order to maintain the authenticity of the city. However, modern pressures, such as overtourism and environmental impacts, pose threats to the city's integrity, necessitating mitigation strategies that could aid in the protection and preservation of its historical and cultural significance (García-Hernández, De la Calle-Vaquero, & Yubero, 2017).

<figure><img src="image_82.png" /><figcaption>Figure 1: The Medieval City of Rhodes (UNESCO, 1988).</figcaption></figure>

### 2.1 Threats

The island of Rhodes, like many Mediterranean regions, is characterized by a great number of threats that are connected to its geographical location and the effects of climate change that can be observed, especially in recent years. In general, the results of climate change are apparent both worldwide and specifically in Rhodes impacting the natural and cultural history of the island (European Commission, n.d.) (Gruber, 2011) (Sesana, Gagnon, Ciantelli, Cassar, & Hughes, 2020).

#### 2.1.1 Temperature

Extreme weather events, including heatwaves, are a significant threat. The decade from 2011 to 2020 was the warmest on record for the island with the data showing that Rhodes is getting warmer, with 2023 being the warmest year at 20.4°C. The increase in temperature deteriorates heritage materials exposed to sunlight and promotes invasive species like the lionfish Pterois miles, affecting local marine life and tourism (Lindsey & Dahlman, 2024) (Meteoblue, n.d.).

#### 2.1.2 Wildfires

Higher temperatures contribute to the risk of wildfires. The EU Environment Program predicts a 14% increase in wildfires by 2030 and 50% by 2100. In summer 2023, Rhodes faced severe wildfires, with around 175,000 acres burned, including protected areas and agricultural land (Jones, Burton, Kelley, & Doerr, 2023).

#### 2.1.3 Coastal Erosion

Coastal erosion is exacerbated by rising sea levels and tectonic activity. The island is in a tectonically active
region, impacting its coastline dynamics. Since 1880, global sea levels have increased by 21-24 cm, affecting Rhodes' coastal zones and historical sites, including the Medieval City (Vandarakis, et al., 2021).

### 2.1.4 Precipitation / Rainfalls

Rhodes experiences scarce and heavy rainfalls. In 2023, it had only 332 mm of rain, one of the driest years in recent history, something that can lead to wildfires. On the other hand, heavy rains can cause flooding and landslides, as seen in 2013 and 2023 (Meteoblue, n.d.).

### 2.1.5 Earthquake

Rhodes has a history of significant earthquakes. In fact, during the medieval period, some significant earthquakes occurred in Rhodes affecting island's structures. For instance, in 1957, an earthquake caused a plethora of damages in structures in the Old Town (Papadopoulos, 2014, p. 181).

## 2.2 The Negative Impact of Overtourism

Overtourism is one of the most important threats that the Old City of Rhodes faces. In fact, increased tourism is strictly connected with overcrowding. This phenomenon is mainly visible during summers that are the peak tourist seasons. In particular, the influx of visitors leads to congestion in the city's narrow streets that originally designed for limited pedestrians. However, overcrowding not only diminishes the visitors experience, but also accelerates the wear and tear on historic buildings and roads (Alamineh, Hussein, Mulu, & Taddesse, 2023) (ICOMOS, 2001-2002).

The constant foot traffic can lead to the deterioration of old stone facades and vibrations in cobblestone pavements, threatening the structural integrity of these historical structures. Moreover, overtourism can result in environmental degradation, which is another critical issue. Specifically, overcrowding leads to pollution, including air, water, and noise pollution, as, for example, the demand for natural resources, such as water, rises significantly, often exceeding the area's capacity (Esteban-Cantillo, Menendez, & Quesada, 2024).

Another important negative impact of overtourism pertains to the commercialization of the Old City, which impacts its authenticity, leading to the loss of its original purpose and meaning. Commercialization

is also connected with the constant reuse of medieval buildings of the town in order to create more hotels, restaurants, shops or anything else that pertains to tourists' needs (Giannakopoulos, et al., 2022).

Last but not least, overtourism negatively affects the local residents. Overcrowding, noise pollution, and rising rental prices make it increasingly difficult for locals to maintain their quality of life, especially during the tourist season. This imbalance highlights how overtourism affects not only the integrity of the Medieval City of Rhodes but also the well-being and the everyday-life of its inhabitants (Shahzalal, 2016).

## 2.3 Mitigation Strategies

The Medieval City of Rhodes faces significant challenges from increasing tourism, environmental degradation, and natural disasters. To address these issues in an effective way, QGIS (Geographic Information System) provides a powerful tool for creating data-driven strategies to mitigate risks and protect the city's cultural and natural heritage (National Geographic, n.d.). Here are some key mitigation strategies utilizing QGIS:

1. **Real-Time Visitor Flow Monitoring**

Through QGIS, authorities can identify the most overcrowded areas of the city by tracking the movements of tourists in real time. This provides the opportunity of control visitors flow and propose redirections to tourists in order to congested roads be avoided reducing the pressure of specific areas of the town.

2. **Predicting Impact on Vulnerable Areas**

The identifications of vulnerable areas to wear and tear because of overcrowding is crucial and the QGIS software can help. Specifically, data analysis through this tool may aid in predictions regarding these areas and taking appropriate measures, like installing protective barriers or limiting visitor access during peak times.

3. **Scenario Planning for Visitor Management**

This software has the ability to model various visitor flow scenarios, something that helps city planners manage the impact of tourists, including entry time regulations or proposing alternative routes.

4. **Visualizing Environmental Impacts**

Through QGIS environmental data are also available. For instance, the identifications of alterations in coasts can be visualized, enabling authorities to understand the
development of the phenomenon and take relevant mitigation measures in order to protect the integrity of the historical site.

## 5. Creating Pedestrian-Only Zones

As the use of QGIS provides city planners with the knowledge of areas of the town with heavy tourist activity they have the ability to propose pedestrian-only zones so as to protect the cobblestone roads and improve both the visitors' experience and the everyday life of locals. Thus, in this map, an effort is made to present the proposed Vehicle Ban Section. This specific layer is depicted with a dotted black line, indicating sections where vehicle access might be restricted in the future to preserve the historical integrity of the area. This map is a valuable tool for managing tourist flows and urban planning within the Medieval City. By distinguishing historical buildings, open spaces, and proposed vehicle-free areas, it helps planners balance preservation with modern demands of overtourism and overcrowding.

In conclusion, QGIS plays a critical role in creating and implementing mitigation strategies for the Medieval City of Rhodes. By integrating data on visitor flow, environmental impacts, and heritage vulnerability, QGIS enables authorities to make informed decisions that protect the city's historical integrity while ensuring sustainable tourism practices.

## 2.4 Using QGIS for Visualization

The preservation and management of historical sites like the Medieval City of Rhodes require not only traditional methods of conservation but also modern technological tools. Geographic Information System (GIS) software, particularly QGIS, plays a crucial role in this process. QGIS is increasingly becoming an invaluable resource for urban planners, archaeologists, and conservationists alike, providing a way to monitor, document, and manage complex environments like Rhodes (Petrescu, 2007).

<figure><img src="image_83.png" /><figcaption>Figure 2: The Medieval City of Rhodes (created using QGIS).</figcaption></figure>
The current state of the Medieval City of Rhodes is presented in detail through the use of QGIS. The use of QGIS is particularly important as it contributes both to a better understanding of the area and to its protection, whether from the significant impacts of natural disasters or from the effects of increasing tourism in the region. What it especially helps with is the improved management of the archaeological site. This is mainly achieved due to the fact that QGIS allows for the creation of accurate and detailed maps, including historical buildings, walls, roads or anything else needed. Essentially, these specific maps contribute to the continuous monitoring of the condition of the research points within the archaeological site, while making it easier to understand any changes occurring in different zones, helping to preserve the integrity of the site (Petrescu, 2007).

Of course, it is worth noting that a QGIS map is particularly helpful in monitoring the medieval city concerning natural disasters and the impacts of climate change. In fact, this island is especially exposed to several natural disasters, such as earthquakes, extreme temperatures or floods. Thus, the application of the software can assist in

visualizing specific areas that are more vulnerable and ultimately create a city plan, helping researchers prioritize the issues that arise. In addition, the creation of maps with QGIS for the Old Town of Rhodes is important for the overall management of the city. For example, through the different layers, historical buildings can be easily distinguished from newer constructions, or even the most frequented by tourists, streets can be highlighted, aiming for a more effective management of overtourism on the island. In conclusion, the use of QGIS for mapping and visualization of the Medieval City of Rhodes is highly important for the optimal preservation and protection of the site, as it enables successful urban planning, monitoring, and overall management of the archaeological site.

As one of the most important steps aimed at addressing various threats and issues related to the Medieval City of Rhodes is the accurate understanding of the area and its significant structures, the mapping process through QGIS is particularly useful. Thus, in this case, the mapping of the archaeological site was mainly done to distinguish certain buildings, structures, and roads. Specifically,

<figure><img src="image_84.png" /><figcaption>Figure 3: The Medieval City of Rhodes with proposed section of pedestrian only zone (created using QGIS).</figcaption></figure>
the buildings, structures, and roads are represented on the map using different colours and shapes, forming distinct layers:

*   **Buildings (Historical & Religious):** Represented in blue. These are significant structures within the medieval city, likely indicating historical or religious buildings. They are marked with numbers for reference and are scattered throughout the city. Their distribution offers insight into the concentration of important buildings within the medieval city.

*   **Squares and Outdoor Sites:** Depicted in pink with cross-hatching, these areas indicate open spaces such as public squares or other significant outdoor sites. Their placement on the map highlights the spatial relationship between the built environment and open spaces, critical for understanding the urban area structure of the city.

*   **Gates:** Shown in green, the map identifies several gates along the city walls. These gates would have been points of entry and exit from the medieval city. In fact, the gates are part of the city's historic fortifications and are vital for controlling access to the city. Their positioning is crucial for analyzing how people entered and exited the city, both historically and in modern times, particularly in relation to tourism flow. Walls: The city walls are represented by a yellow-orange cross-hatched line surrounding the city. These walls form a protective barrier around the city and are a key feature of its medieval fortifications, emphasizing its historical strategic importance.

As mentioned above, the map of the Old Town of Rhodes created through QGIS facilitates the precise identification of each building or open space, which are numbered with details such as their names and uses. It is worth noting that the visualization of these current structures contributes to highlighting the city's architecture, while also providing information on how the urban design can impact daily life and tourist activity within or outside the city.

### 3 CONCLUSIONS

In conclusion, the important phenomenon of overtourism has become a significant global concern that impacts more and more destinations worldwide. This is primarily because, along with the rise of the

phenomena of climate change affecting our planet, it may lead to horrible effects, both regarding the preservation of global cultural and natural heritage and regarding people's everyday-lives (Nugroho & Hardilla, 2020). This basically includes the daily lives of local communities that are disrupted by overtourism in affected areas, as well as the experiences of tourists, who eventually do not enjoy their trips due to this serious phenomenon. Local authorities and countries, especially those that have already affected the consequences of overtourism, should take specific mitigation measures in order to protect cultural heritage and themselves. In order to achieve this, the use of modern technology and especially the use of GIS are essential.

In Greece, a country renowned for its rich cultural attractions, the Medieval City of Rhodes stands out as a prime example of how overtourism along with the effects of climate change can impact a historically significant site. The city's narrow streets and medieval architecture are particularly vulnerable to the pressures of heavy tourist traffic. Without appropriate measures, the integrity of this archaeological site is at risk.

In conclusion, the future of the Medieval City of Rhodes depends on a multi-faceted approach that combines traditional conservation methods with modern technological innovations. Geographic Information Systems, specifically QGIS, provide a crucial platform for analysing, managing, and protecting this historical site. The lessons learned from this case study can serve as a model for other historical cities worldwide that face similar challenges in balancing preservation with tourism and environmental pressures.

### REFERENCES

Alamineh, G., Hussein, J., Mulu, Y., & Taddesse, B. (2023, June 17). *The negative cultural impact of tourism and its implication on sustainable development in Amhara Regional State*. Retrieved from tandfonline.com: https://www.tandfonline.com/doi/full/10.1080/23311983.2023.2224597

Avdikos, V. (2011). *Local hegemonic blocs: The case of tourism in Rhodes*. Retrieved from academia.edu: https://www.academia.edu/561874/Local_hegemonic_blocs_The_case_of_tourism_in_Rhodes

Esteban-Cantillo, O., Menendez, B., & Quesada, B. (2024, April 15). *Climate change and air pollution impacts on cultural heritage building materials in Europe and Mexico*. Retrieved from sciencedirect.com: https://www.sciencedirect.com/science/articlle/pii/S0048969724010842
European Commission. (n.d.). *Causes of climate change*. Retrieved from climate.ec.europa.eu: https://climate.ec.europa.eu/climate-change/causes-climate-change_en

European Union. (2019). *European framework for action on cultural heritage*. Retrieved from op.europa.eu: https://op.europa.eu/en/publication-detail/-/publication/5a9c3144-80f1-11e9-9f05-01aa75ed71a1

García-Hernández, M., De la Calle-Vaquero, M., & Yubero, C. (2017, August 4). *Cultural Heritage and Urban Tourism: Historic City Centres under Pressure*. Retrieved from mdpi.com: https://www.mdpi.com/2071-1050/9/8/1346

Giannakopoulos, D., Karekou, Z., Menegaki, E., Tsilimantou, E., Ioannidis, C., Maistrou, E., Moropoulou, A. (2022, August 1). *Reuse of Historic Buildings in the Medieval City of Rhodes to Comply with the Needs of Sustainable Urban Development*. Retrieved from mdpi.com: https://www.mdpi.com/2073-445X/11/8/1214

Gruber, S. (2011). *The Impact of Climate Change on Cultural Heritage Sites: Environmental Law and Adaptation*. In The Impact of Climate Change on Cultural Heritage Sites: (pp. 209-219). Lexxion Verlagsgesellschaft

ICOMOS. (2001-2002). *Heritage at Risk from tourism*. Retrieved from icomos.org: https://www.icomos.org/public/risk/2001/to urism.htm

Jones, M., Burton, C., Kelley, D., & Doerr, S. H. (2023, July 28). *Greece wildfires: how climate change is involved, and what we can do about it*. Retrieved from theconversation.com: https://theconversation.com/greece-wildfires-how-climate-change-is-involved-and-what-we-can-do-about-it-210404

Lindsey, R., & Dahlman, L. (2024, January 18). *Climate Change: Global Temperature*. Retrieved from climate.gov: https://www.climate.gov/news-features/understanding-climate/climate-change-global-temperature

Luttrell, A. (2003). *The Town of Rhodes 1306-1356*. Rhodes: City of Rhodes Office for the Medieval Town (Ministry of Culture- Archaeological Receipts Fund- City of Rhodes Contract).

Meteoblue. (n.d.). *Climate Change Rhodes*. Retrieved from meteoblue.com: https://www.meteoblue.com/en/climate-change/rhodes_greece_400666

National Geographic. (n.d.). *GIS (Geographic Information System)*. Retrieved from education.nationalgeographic.org: https://education.nationalgeographic.org/resource/geog-raphic-information-system-gis/

Nugroho, A. C., & Hardilla, D. (2020, January). *The Importance of Cultural Heritage Conservation in Society: A Review and Prospect for Future Cities, with Bandar Lampung as Cased Study*. Retrieved from researchgate.net: https://www.researchgate.net/publication/338351963_T he_Importance_of_Cultural_Heritage_Conservation_in

*Society_A_Review_and_Prospect_for_Future_Cities_ with_Bandar_Lampung_as_Cased_Study*

Papadopoulos, G. (2014). *Rhodes: Earthquakes and Tsunamis from Antiquity to the Present Day*. Athens: Aselotos Publications.

Petrescu, F. (2007, October). *The Use of GIS Technology in Cultural Heritage*. Retrieved from researchgate.net: https://www.researchgate.net/publication/228636422_T he_use_of_GIS_technology_in_cultural_heritage

Sesana, E., Gagnon, A. S., Ciantelli, C., Cassar, J., & Hughes, J. J. (2020, March 31). *Climate change impacts on cultural heritage: A literature review*. Retrieved from wires.onlinelibrary.wiley.com: https://wires.onlinelibrary.wiley.com/doi/epdf/10.1002/wcc.710

Shahzalal. (2016). *Positive and Negative Impacts of Tourism on Culture: A Critical Review of Examples from the Contemporary Literature*. Retrieved from https://core.ac.uk/download/pdf/234696865.pdf

UNESCO. (1972, November 16). *Convention Concerning the Protection of the World Cultural and Natural Heritage*. Retrieved 2024, from whc.unesco.org: https://whc.unesco.org/en/conventiontext/

UNESCO. (1988). *Medieval City of Rhodes*. Retrieved from https://whc.unesco.org/en/list/493/

UNESCO. (2016). *The Protection of Cultural Heritage*. Retrieved 2024, from en.unesco.org: https://en.unesco.org/creativity/policy-monitoring-platform/protection-cultural-heritage

Vandarakis, D., Panagiotopoulos, I. P., Loukaidi, V., Hatiris, G.-A., Drakopoulou, P., Kikaki, A., Kapsimalis, V. (2021, August 7). *Assessment of the Coastal Vulnerability to the Ongoing Sea Level Rise for the Exquisite Rhodes Island (SE Aegean Sea, Greece)*. Retrieved from mdpi.com: https://www.mdpi.com/2073-4441/13/16/2169
# ADAclassifier: Trying to Ascertain Why the Ground Is Moving

José A. Navarro<fnref n="1" /> Ⓓª, Anna Barra<fnref n="1" /> Ⓓᵇ, María Cuevas-González<fnref n="1" /> Ⓓᶜ, Pablo Ezquerro<fnref n="2" /> Ⓓᵈ, Silvia Bianchini<fnref n="3" /> Ⓓᵉ, Marta Béjar-Pizarro<fnref n="2" /> Ⓓᵉ, Riccardo Palamà<fnref n="1" /> Ⓓᵉ, Oriol Monserrat<fnref n="1" /> Ⓓᵉ

**Keywords:** Ground Motion, Ground Motion Classification, Active Deformation Areas.

**Abstract:** The large availability of ground deformation measurements generated using Multi-Temporal Synthetic Aperture Radar (MT-InSAR), further increased by the contribution made by the European Ground Motion Service (EGMS), made displacement maps an increasingly common tool. However, their analysis is a complex task due to the large volume of Measurement Points (MPs) provided. ADAfinder, a tool within the ADAtools suite, allows for a reduction in the volume of information to be analyzed by identifying Active Deformation Areas (ADAs), i.e., areas of the terrain that actually move in a coherent and perceptible manner, including an estimate of the reliability of said identification. From here, the natural next step is identifying the reason why these areas are moving. This work presents ADAclassifier, another tool included in the ADAtools suite, still under development, aimed at evaluating up to five possible causes of ground movement, such as subsidence, landslide, uplift, sinkhole, and constructive settlement.

## 1 INTRODUCTION

In recent years, the widespread availability of displacement maps over large areas, generated using Multi-Temporal Synthetic Aperture Radar (MT-InSAR) satellite interferometry techniques, has significantly grown. The launch of the Copernicus Sentinel-1 satellites in 2014, which provide global and consistent data acquisitions under an open-access and free data distribution policy, represented a major shift in both their use and application.

A key result of this advancement has been the creation of regional, national, and continental ground motion services, offering detailed displacement maps that deliver in-depth insights into both human activities and natural events. Since 2022, the European

Ground Motion Service (EGMS) (EGMS, 2017) has been freely offering billions of displacement Measurement Points (MP), updated annually, with coverage spanning almost all of Europe. These maps are known for their high accuracy and precision.

Although these maps hold great potential for territorial management and risk assessment, they remain underutilized due to the difficulties related to their interpretation; thus, automated tools that can streamline and speed up the processes of data extraction, analysis, and interpretation are needed. The creation of derived and simplified maps is essential to increase the adoption of said products by both expert and non-expert InSAR users. In this context, we introduce the ADAtools (Active Deformation Area tools (Barra et al., 2017; Tomás et al., 2019; Navarro et al., 2020).

The ADAtools include several tools; among them, it is important to highlight ADAfinder and ADAclassifier. The first one, ADAfinder, facilitates the automatic extraction and selection of the most significant Active Deformation Areas (ADAs) (Barra et al., 2017), serving as a critical first step in transforming a large number of individual MPs into a more manageable set of polygons for further analysis or appli-

<footnote n="1">¹Centre Tecnològic de Telecomunicacions de Catalunya-CERCA, Castelldefels, Spain</footnote>

<footnote n="2">²Geohazards InSAR Laboratory and Modelling Group (InSARlab), Geohazards and Climate Change Department, Geological and Mining Institute of Spain (IGME-CSIC), Madrid, Spain</footnote>

<footnote n="3">³Università Degli Studi di Firenze (UNIFI), Dipartimento di Scienze della Terra, Firenze, Italy
{jose.navarro, anna.barra, maria.cuevas}@cttc.cat, p.ezquerro@igme.es, silvia.bianchini@unifi.it, m.bejar@igme.es,
{riccardo.palamà, oriol.monserrat}@cttc.cat</footnote>

<footnote>ªhttps://orcid.org/0000-0001-7877-1516</footnote>

<footnote>ᵇhttps://orcid.org/0000-0001-6254-7931</footnote>

<footnote>ᶜhttps://orcid.org/0000-0002-4988-5669</footnote>

<footnote>ᵈhttps://orcid.org/0000-0001-8667-5030</footnote>

<footnote>ᵉhttps://orcid.org/0000-0003-2724-5641</footnote>

<footnote>fhttps://orcid.org/0000-0001-7449-4048</footnote>

<footnote>ghttps://orcid.org/0000-0001-6121-9485</footnote>

<footnote>hhttps://orcid.org/0000-0003-2505-6855</footnote>
cation.

Figure 1 shows the process performed by ADAfinder. At the top, all existing MPs in some area can be observed. It can be seen that the volume of available information makes the analysis a difficult task. At the bottom, the ADAs detected by ADAfinder are displayed. It is important to emphasize that the amount of information the expert needs to analyze is significantly lower. Additionally, together with the polygons describing the areas of interest, an assessment of the reliability associated with the polygons is included, which establishes a scale of four values, ranging from 1, “very reliable” to 4 “very unreliable.” This quality index let the expert focus on those findings that have a greater likelihood of being accurate.

<figure><img src="image_85.png" /><figcaption>Figure 1: ADAfinder: from a plethora or measurement points to a few polygons depicting the Active Deformation Areas (ADAs). The quality index stands for the reliability of the results. ADAs in class 4 (“very unreliable”) are not shown. Source: (Barra et al., 2017).</figcaption></figure>

ADAfinder is a well-seasoned tool, used in numerous projects since 2017, including U-Geohaz (Solari et al., 2020) and MOMPA (Gasc-Barbier et al., 2021) among others. One of the most recent projects (Navarro et al., 2022; Navarro et al., 2023) aimed to publish ADAs for almost all of Europe on the web, processing MPs published by the EGMS (see Figure 2). (Navarro, 2024b) includes a very detailed description of the ADAfinder tool.

Once moving areas or ADAs have been identified, ADAclassifier—this work’s primary focus—enables an initial assessment of the processes likely responsible for the displacement. By integrating auxiliary data—such as inventories of several kinds and a Digital Terrain Model (DTM)—each ADA is given a preliminary characterization with regard to five different ground movement processes: landslides, subsidence, sinkholes, construction settlements, or uplift. Like ADAfinder, ADAclassifier also provides an

<figure><img src="image_86.png" /><figcaption>Figure 2: The European ADA web map (https://groundmotionadas.com). The ADAs were identified by the ADAfinder tool.</figcaption></figure>

assessment of the reliability of the results obtained, again with the goal of helping the expert to analyze said results. This assessment is once again summarized into four levels of reliability, from best to worse.

Figure 3 shows the classification results produced by the original version of ADAclassifier for the example area shown in Figure 1, considering only three processes: landslides, subsidence, and constructive settlement. The reliability of the results has been highlighted using different colors: red, yellow and green (from positive to negative certainty).

<figure><img src="image_87.png" /><figcaption>Figure 3: ADAclassifier (original version): for each ADA, its potential correspondence to one of the different types of analyzed ground movements is assessed. Source: (Navarro et al., 2020).</figcaption></figure>

The first implementation of the ADAtools dates to 2017. Since then, they have been used in numer-
ous projects and have undergone a series of improve-
ments leading to the current version. ADAclassifier
was included in the ADAtools since their very first
version. However, unlike ADAfinder, its initial re-
lease had many limitations, which have recently been
addressed. This work discusses how.

## 2 ADAclassifier

As just mentioned, the original version of ADAclas-
sifier had several flaws, making it highly advisable to
replace it with a new implementation that addressed
these issues. This section describes these flaws and
explains how the new version of the tool addressed
them.

## 2.1 The First Version of the Tool

The core problem lay in the simplicity of the algo-
rithms used to determine whether a ground movement
corresponded to one phenomenon or another. Those
algorithms were an *initial approach* to solving the
problem. In fact, each of the tests aimed at determin-
ing whether a ground movement was due to a specific
type of process (i.e., landslide or settlement) relied on
a *unique* decision tree (one per process). For exam-
ple, Figure 4 depicts the tree used by the landslide
test. Additionally, the number of variables these al-
gorithms considered, such as the slope value or hor-
izontal ground displacement, was relatively limited.
These algorithms were still far from modeling reality
reasonably, as they can be considered a first and early
attempt at approaching the problem.

Furthermore, these algorithms largely depended
on the availability of inventories related to the type of
process being evaluated, such as geological or con-
struction settlement inventories. Their absence in-
evitably led to results that never confirmed that the
ground movement in question was caused by the
tested process. In other words, if a landslide inven-
tory was not available, the ground movement would
never be classified as a landslide; at best, it would be
labeled as a "potential landslide". See again Figure 4.

The simplicity of the algorithms and the lim-
ited availability of inventories explain the low con-
fidence in ADAclassifier's results. These algorithms
only consider a small portion of the reality of ground
movements, overlooking many factors. Additionally,
the lack of inventories forces the algorithms into a
more conservative mode, hindering their ability to
confidently affirm the occurrence of a specific phe-
nomenon. This undermines trust in the results.

<figure><img src="image_88.png" /><figcaption>Figure 4: The unique algorithm (decision tree) used by the old version of ADAclassifier to decide whether a ground movement was produced by a landslide. The threshold values are provided by the user. Excerpted from (Tomás et al., 2019).</figcaption></figure>

Another significant issue plaguing the first ver-
sion of ADAclassifier was performance. As previ-
ously mentioned, inventories played a crucial role in
the implementation of the old version of the tool. In
fact, they were used (when available) by the algo-
rithms implementing every ground motion detection
process. From ADAclassifier's perspective, an inven-
tory is a file containing polygons that define the areas
of the terrain where a specific phenomenon occurs.
Some inventories may also include additional data re-
lated to each polygon; the geological inventory is a
clear example: a value states the type of lithology re-
lated to the polygon. The typical "inventory check"
consisted of computing the percentage of intersection
of the polygon defining the ADA with those included
in said inventory. Due to the very poor implementa-
tion of this process, ADAclassifier's performance was
far from satisfactory, making its use on projects with
territorial coverage beyond the local level impossible.

The team responsible for both the development
(and on many occasions also for the exploitation) of
ADAclassifier were aware of the need for a redesign
of the tool that could be relied upon. These are the
reasons that led to the decision to create the new ver-
sion of ADAclassifier that is described in the next sec-
tion.

## 2.2 The New ADAclassifier

This section analyzes the changes undergone by
ADAclassifier from different perspectives: its inputs,
algorithms, outputs, performance and reliability.
### 2.2.1 The Inputs

From a user standpoint, the new version of ADAclassifier is virtually identical to the original. In fact, almost all of the input files are the same in both cases. The following is a brief list of these files, indicating when applicable which ones are needed only in the new version.

1. ADAs (polygons) and the points used to construct them.

2. Digital Terrain Model (DTM) or, alternatively, slope and aspect maps (the last two being exclusive to the new version).

3. Inventories: landslides, sinkholes, subsidence, constructive settlements, geological.

4. Horizontal displacement and, only for the new version, vertical displacement.

From the previous list, only the files from points 1 and 2 are mandatory, and, in fact, in point 2, it is possible to choose between the DTM or slope and aspect data. If the DTM is provided, then ADAclassifier computes both the slope and aspect; if not, the user must provide these two files as input instead of the DTM. The rest of the files are optional, as in the original version.

Actually, there are two additional input files in the latest version of ADAclassifier (scores and score-to-classes defaults). However, it is very rare for the end user to have to manipulate them, as they include default values that should not normally be changed unless you are an expert in the subject. The content of these files will be discussed later, since they are crucial to making ADAclassifier a very flexible tool.

### 2.2.2 The Algorithms

The most significant changes have been made to ADAclassifier's algorithms. Within the framework of the RASTOOL project (Montserrat et al., 2024; RASTOOL project, 2024), a major effort has been made to characterize the reasons behind different types of ground movements. Specifically, for landslides, constructive settlements, sinkholes, subsidence and, new for this version of ADAclassifier, also uplifts.

Several changes have affected the algorithms. Firstly, their complexity has increased. Unlike the initial version, where a single decision tree based on a few parameters was used to check a specific ground movement process, now multiple trees have been incorporated. More specifically, these are 4 for landslides, 6 for constructive settlements, 5 for sinkholes, 5 for subsidence, and 4 for uplifts. This sums up to

a total of 24 checks, compared to the original 5. Figure 5 shows two of the four decision trees implementing the landslide detection process. Note the THLAnn (thresholds) and SCLAnn (scores) labels in said Figure. These are values input by the user controlling the behavior of the algorithms and will be discussed below.

This increase in the number of decision trees used by the tool makes the detection processes for the different kinds of ground movements more exhaustive. Additionally, the set of physical magnitudes being checked has been expanded, including horizontal and vertical displacement velocities, line of sight velocity, aspect, slope, the percentage of intersection of the ADA with different inventories, and the $r^2$ coefficient of fit of the deformation time series and an inverse exponential function. Some relationships between these magnitudes are also checked. This is expected to result in greater reliability of the obtained results.

Together with the notable improvement resulting from the increased number of decision trees and physical magnitudes involved, another significant key improvement factor is the way in which all these elements (trees, magnitudes, relationships) contribute to the final outcome. In the first version of ADAclassifier, all decision trees reached a definitive conclusion about whether a given ground movement could be labeled one way or another. This was materialized with three possible outcomes for each kind of process: "Not X", "Potential X" or "Is X", where "X" stands for the different process types being checked (such as sinkhole or subsidence). The user could only control the values of the thresholds used in the decision trees (see Figure 4).

Now, the procedure (for each type of process) consists of adding or subtracting points to a cumulative total or score. Each leaf of the decision trees contributes a specific number of points, either positive or negative, to the process evaluation.

These points are represented by the labels named SCLAnn in the leaves of the decision trees in Figure 5, where SC stands for "Scoring", LA indicates that this score refers to the landslide process, and nn is a numeral to distinguish one score from another. These labels change depending on the detected process; for example, in the case of sinkholes, the pattern is SCsinn, where SI denotes the said sinkhole process. Note that the thresholds for some magnitudes follow a similar pattern; in this case these start with the letters TH instead of SC.

Each of the leaves of the decision trees involved in the detection contributes a certain number of points to the final score, which will subsequently be used to decide the status of the ground movement with respect
<figure><img src="image_89.png" /><figcaption>Figure 5: The four decision trees used by the new version of ADAclassifier to determine if a ground movement might be a landslide. Thresholds and scores are denoted by labels like THLAnn and SCLAnn.</figcaption></figure>

to the process being identified. For instance, in the case of subsidence, there are five decision trees; consequently, the final score will consist of the sum of the points obtained in each of them.

The sum of both improvements (better decision trees plus scoring) makes the detection of the type of ground movement in which an ADA is involved more reliable than before, since more situations are now evaluated. In fact, it can be stated that the new version of ADAclassifier is more resilient, as the score assigned to a given process is obtained through a combination of algorithms (the decision trees) and not from a single one as before. The unavailability of some data does not prevent an attempt to reach a conclusion.

It is worth noting that the values of thresholds and the number of points rewarded are customizable through the tool's configuration files. This offers experts a high degree of flexibility, allowing them to fine-tune those parameters without requiring code changes. For example, to accurately model real-world phenomena, threshold values such as slopes or horizontal velocities may need adjustment. With regard to scoring, changing the number of points assigned to the leaves makes possible giving more or less importance to the different cases depicted by the decision trees. Obviously, having greater flexibility implies greater complexity: only experts are capable of changing said values without breaking the application.

### 2.2.3 From Scores to Classes

Despite the increased flexibility and reliability provided by the scoring system, it is necessary to offer a final, simple categorization of the results. This implies that these must be expressed through a reduced number of classes. ADAclassifier performs this process, translating the points obtained for each process into a value that can be one of the following: "It is not X", "It may be X", "It should be X", and finally, "It is X, where X, once again, represents each of the processes being checked."

The key here is how the points are translated into classes. This is done through a configuration file that the expert can modify to fine-tune the system. In this file, for each process, the ranges of points corresponding to each class are specified. In this way, the tool becomes even more flexible, allowing the expert to make adjustments without needing to modify the tool itself.

### 2.2.4 The Output

The unique output file is an extension of the input file with ADAs. There, the original attributes are preserved, such as the polygon defining the perimeter of the ADA or the Quality Index (QI) stating its reliability. The interested reader may check the complete
list of attributes in (Navarro, 2024b) or in (Navarro, 2024a).

The extra fields added by ADAclassifier for each ADA in the data set, which are related to the results of the classification process, are:

* The mean values of slope and aspect for the ADA,

* as well as an extra flag stating whether the value of the aspect above is good for some tests,

* the $r^2$ coefficient, measuring how good is the fit of the ADA's deformation time series and an inverse exponential function,

* the percentage of intersection between the ADA's polygon and those in every available inventory or geological maps,

* the score obtained and the corresponding class for each of the five processes tested and, finally,

* a summary of the results, stating which is (or are) the predominant process(es) making the ADA move.

As can be observed in the list above, ADAclassifier not only generates an answer regarding which process(es) possibly originate(s) the ADAs, but also exposes all the information used to reach that conclusion. This not only serves to justify the answer provided but also allows the expert to analyze the results in case they do not correspond to reality. Having this information will allow the adjustment of algorithms, thresholds, and scores, if necessary.

### 2.2.5 Performance Assessment. A Nationwide Project

As stated in section 2.1 the low performance of the first version of ADAclassifier was a serious problem, caused, mainly, but among other issues, by the poor implementation of the inventory checks. This new release has solved all the problems the authors were aware of, which has produced a noticeable boost in performance that makes it possible to tackle not only regional, but even nationwide projects.

In fact, to put the tool to the test, a total of 661 ADA datasets have already been processed, covering the entire surface of Spain, which is 506,000 km² and includes 21,526 ADAs. The ADAs were already available since they were computed using ADAfinder for the project described in (Navarro et al., 2022; Navarro et al., 2023). To this data, we added a DTM, as well as geological, landslide, and constructive settlement inventories covering the entire territory. The sinkhole inventory was only available for a very small region measuring 17,274 km², so its impact on performance is negligible. The slope and aspect maps for all

Spain were calculated only once using the DTM. This way, ADAclassifier did not need to regenerate them when processing each data set.

The characteristics of the hardware follow:

* Processor (64-bit): 2 × Intel® Xeon® Silver 4309Y CPU @ 2.80 GHz, 2801 MHz, 8 Core(s), 16 Logical Processes.

* RAM: 384 Gb.

* Disk: 18 Tb, Magnetic (not SSD).

* Operating system: Microsoft Windows Server 2022 Standard.

* Year: 2020.

Processing the complete datasets (all Spain) took just 8 hours, that is, less than 0.06 seconds per km², 1.34 per ADA or 43.6 per data set. Note that the server described above was *not* fully dedicated to ADAclassifier for it is shared by about a dozen users launching CPU-consuming tasks. Therefore, results could be still better if no concurrent processes competed for the server's resources.

Looking at these results, the authors believe that the tool is well suited to process nationwide projects. Said results also seem to indicate that it could be used at a continental level.

### 2.2.6 Reliability Assessment

The reliability of the results has been continuously verified since the first release of the new version of ADAclassifier using some of the already available datasets. The results of these initial analyses have served to fine-tune the threshold and scoring files.

Thanks to the availability of EGMS data and thus the processing of the entire Spanish datasets, we now have a much larger volume of information to perform this verification. Obviously, this is a large task due to the overabundance of data (661 datasets), but it is precisely such abundance that will allow us to determine with a reasonable level of certainty whether the conclusions of ADAclassifier are correct or if, in some cases, the system needs further tuning.

In fact, and this is because the determination of the causes of ground movements is not a settled issue, some inconsistencies have already been detected in a very small number of cases: apparently, a few ADAs seems to be moving, according to ADAclassifier, because of both landslide (down) and uplift (up) processes, which is completely contradictory.

At the time of writing this work, the reason for this problem was being analyzed. These cases will be used to understand the algorithm behavior, thus to adjust some thresholds and scores or to modify
some of the decision trees that control the affected processes—thus modifying the application.

Except for these rare exceptions, most ADAclassifier's results are reliable. However, verification by an expert will always be recommended to ensure that situations like the one just described do not go unnoticed.

## 3 CONCLUSION

Despite being available since 2018 as an integral part of the ADAtools, the original version of ADAclassifier has not been regularly used to attempt to determine the causes of ground movement due to this tool's performance issues but, primarily, because of the low reliability of the results offered, which were obtained by means of overly simplistic algorithms.

Recently, and within the framework of the RAS-TOOL project, a great effort has been made to identify the reasons why the five analyzed ground movement processes occur. Ultimately, this effort was aimed at improving the understanding of these processes, which has resulted in a set of detection algorithms closer to reality, materialized in the aforementioned decision trees.

ADAclassifier is not yet a mature tool, but today it is much closer to being able to be exploited regularly in expert environments not necessarily involved in research. It is true that notable defects have still been detected, such as the already mentioned landslide/uplift dichotomy, but this is one more step to be solved like the many others that have been overcome to reach the current state of ADAclassifier. The availability of a large volume of EGMS data making possible the processing of ADAs at a national level opened the doors to an extensive validation task, thanks to which the authors are reasonably satisfied with the reliability of the results produced by tool.

Perhaps it is adventurous to advance a date on which ADAclassifier can be considered a production tool (in expert environments), but the authors consider that, despite the aforementioned defects, that moment is not far off. Probably this will happen by the end of 2024 or the beginning of 2025.

ADAclassifier is part of the ADAtools, which can be downloaded for free at https://adatools.cttc.es.

## ACKNOWLEDGMENTS

This work is part of the Spanish grant SARAI, PID2020-116540RB-C21, funded by MCIN/AEI/10.13039/501100011033.

## REFERENCES

Barra, A., Solari, L., Béjar-Pizarro, M., Monserrat, O., Bianchini, S., Herrera, G., Crosetto, M., Sarro, R., González-Alonso, E., Mateos, R. M., Ligüerzana, S., López, C., and Moretti, S. (2017). A methodology to detect and update active deformation areas based on Sentinel-1 SAR images. *Remote Sens.*, 9(10).

EGMS (2017). European Ground Motion Service (EGMS) - A proposed Copernicus service element. White paper, European Ground Motion Service.

Gasc-Barbier, M., Barra, A., Buxó, P., Trapero, L., Crosetto, M., Colell, X., Fabregat, I., Echeverria, A., and Marturia, J. (2021). Monitoring deformations related to geological risks with InSaR data - the MOMPA project. *IOP Conference Series: Earth and Environmental Science*, 833(1):012142.

Montserrat, O., Barra, A., Béjar-Pizarro, M., Rivera, J., Galve, J. P., Guardiola, C., Cuevas-González, M., Mateos, R. M., Ezquerro, P., Azañón, J. M., Shahbazi, S., Navarro, J., Crosetto, M., and Luzi, G. (2024). ADAtools: herramientas gratuitas y fáciles de usar para extraer y analizar semiautomáticamente mapas de desplazamiento basados en interferometría multitemporal. aplicaciones al European Ground Motion Service (EGMS). In *Proceedings of XX Congreso de la Asociación Española de Teledetección*.

Navarro, J., García, D., and Crosetto, M. (2023). New challenges in the implementation and exploitation of a low-cost web map of the Active Deformation Areas across Europe. In *Proceedings of the 9th International Conference on Geographical Information Systems Theory, Applications and Management - Volume 1: GISTAM*, pages 192-199.

Navarro, J. A. (2024a). ADAclassifier's user guide. Included in the ADAtools software package. Available at: https://adatools.cttc.es.

Navarro, J. A. (2024b). ADAfinder's user guide. Included in the ADAtools software package. Available at: https://adatools.cttc.es.

Navarro, J. A., García, D., and Crosetto, M. (2022). An in-house, self-developed low-cost WebGIS relying on open-source and self-made tools to show the Active Deformation Areas across Europe. In *Proceedings of the 8th International Conference on Geographical Information Systems Theory, Applications and Management*.

Navarro, J. A., Tomás, R., Barra, A., Pagán, J. I., Reyes-Carmona, C., Solari, L., Vinielles, J. L., Falco, S., and Crosetto, M. (2020). ADAtools: Automatic detection and classification of Active Deformation Areas from PSI displacement maps. *ISPRS Int. J. Geo-Inf.*, 9:584.

RASTOOL project (2024). RASTOOL. https://rastool.cttc.es. Accessed on October 25, 2024.

Solari, L., Bianchini, S., Franceschini, R., Barra, A., Monserrat, O., Thuegaz, P., Bertolo, D., Crosetto, M., and Catani, F. (2020). Satellite interferometric data for landslide intensity evaluation in mountainous regions. *Int. J. Appl. Earth Obs. Geoinf.*, 87.
Tomás, R., Pagán, J. I., Navarro, J. A., Cano, M., Pastor, J. L., Riquelme, A., Cuevas-González, M., Crosetto, M., Barra, A., Monserrat, O., López-Sánchez, J. M., Ramón, A., Iborra, S., del Soldato, M., Solari, L., Bianchini, S., Raspini, F., Novali, F., Ferreti, A., Constantini, M., Trillo, F., Herrera, G., and Casagli, N. (2019). Semi-automatic identification and pre-screening of geological-geotechnical deformational processes using Persistent Scatterer Interferometry datasets. *Remote Sens.*, 11(14).
# Unsupervised Image Classification Algorithms Applied to Fire-Prone Area Detection

I. Rahimi<fnref n="1" />,<fnref n="2" /> IDª, L. Duarte<fnref n="1" />,<fnref n="2" /> IDᵇ and A. C. Teodoro<fnref n="1" />,<fnref n="2" /> IDᶜ

**Keywords:** Remote Sensing, Wildfire, Classification, Zagros Mountains.

**Abstract:** Remote sensing data has become critical in identifying fire-prone areas, providing essential insights through satellite imagery and various geospatial inputs. These data sources allow for real-time monitoring, mapping fire susceptibility, and assessing factors such as vegetation, fuel moisture, land use, and environmental conditions. Numerous supervised and unsupervised models combined with remote sensing data have shown great potential in predicting fire-prone regions, offering accurate and timely information for early warning systems and resource allocation. This study focuses on applying two unsupervised methods—PCA, and K-means—using inputs like Sentinel-2 imagery, elevation, and the Zagros Grass Index (ZGI) to identify fire-prone areas in the Kurdo-Zagrosian forests, an area increasingly vulnerable to wildfires. Among the two methods evaluated, PCA demonstrated superior performance in predicting fire-susceptible areas, accurately classifying 80% of the burned regions from 2021 to 2023 as moderate to high-risk zones.

## 1 INTRODUCTION

Forest fires' increasing frequency and intensity worldwide is an escalating concern, driven by natural and human-induced factors such as extreme weather conditions, shifting land use patterns, and rapid urban expansion (Zema, 2020; Bowman, 2017). These factors, especially under the growing influence of climate change, exacerbate the risk of wildfires. This significant loss of forest cover highlights the urgent need for effective wildfire monitoring and prevention strategies.

In recent years, integrating geospatial and remote sensing (RS) data/technologies has provided invaluable insights into wildfire risk factors (Teodoro, 2013). Through RS data, researchers can monitor and analyze variables such as land cover, temperature, and vegetation phenology (Duarte, 2018). This spatial data, combined with Geographic Information Systems (GIS), allows for the continuous monitoring of large areas and provides timely information about the likelihood of fire incidents (Mishra, 2024). RS data enables researchers to track environmental changes in real time and analyze

critical variables such as fuel moisture content, temperature trends, and human activity patterns, all of which contribute to the increased risk of fires.

A key application of this technological advancement is the development of Forest Fire Susceptibility Maps (FSMs). These maps are crucial for identifying areas at high risk of wildfires, enabling authorities to allocate resources efficiently and implement mitigation strategies in advance. By locating and assessing fire-prone regions, FSMs play a central role in reducing the vulnerability of ecosystems and communities while supporting informed ecological and urban planning decisions (Ghorbanzadeh, 2019). The effectiveness of FSMs has been significantly enhanced by integrating GIS, RS, and image classification algorithms, allowing researchers to predict fire susceptibility better and improve early warning systems (Rihan, 2023).

The growing popularity of image classification algorithms in fire susceptibility mapping is due to their ability to process large-scale, high-dimensional datasets and model complex non-linear relationships more effectively than traditional statistical methods (Kantarcioglu, 2023). These models perform well at

<footnote n="1">¹Department of Geosciences, Environment and Spatial Planning Faculty of Sciences of the University of Porto, Portugal</footnote>

<footnote n="2">²Institute of Earth Sciences, FCUP Pole, Porto, Portugal</footnote>

<footnote>ªID https://orcid.org/0009-0002-7411-8637</footnote>

<footnote>ᵇID https://orcid.org/0000-0002-7537-6606</footnote>

<footnote>ᶜID https://orcid.org/0000-0002-8043-6431</footnote>
integrating diverse data types, from climatic conditions and vegetation indices to human activity metrics, which allows for more precise fire risk assessments. In recent years, the development of FSM has been advanced through both single-model applications and ensemble methods, frequently showing improved performance over more conventional techniques (Piao, 2022; Saha, 2023). These advancements underscore FSM's potential to surpass traditional frameworks, particularly by leveraging methods that incorporate broader environmental and spatial analysis techniques.

Unsupervised classification methods, such as Principal Component Analysis (PCA) and K-means clustering, have proven valuable for extracting meaningful patterns from satellite imagery. PCA reduces large datasets into key components, capturing the most significant variations (Sharma, 2021), which is crucial for identifying fire susceptibility factors like vegetation, moisture, and land use patterns. K-Means, on the other hand, clusters satellite image data based on spectral similarity, enabling the identification of regions with similar environmental characteristics (Theodoridis, 2009). These methods are especially effective in handling large-scale, heterogeneous datasets from RS data (Celik, 2009).

The objective of this study is to apply two unsupervised image classification techniques—PCA and K-means—using inputs like Sentinel-2 imagery, elevation, and the Zagros Grass Index (ZGI) (Rahimi, 2024) to identify fire-prone areas in the Kurdistan Zagrosian forests. This region, vulnerable to increasing wildfire incidents, presents a complex environmental landscape where these methods can highlight patterns indicating higher fire susceptibility. By utilizing unsupervised algorithms, the study aims to provide a data-driven approach to wildfire susceptibility mapping, improving prevention and mitigation strategies in fire-prone regions.

# 2 METHODS

## 2.1 Study Area

The method is designed for areas with Semi-Arid (SA) and Semi-Mediterranean (SM) climates, exemplified by the forests of Marivan and Sarvabad in the Kurdistan Province, western Iran, bordering Iraq. These regions have experienced frequent fire incidents over the past decade. Geographically, the study area spans between longitudes 45°57'50" E and 46°46'41" E and latitudes 35°1'1" N and 35°49'51" N. Located in the northern Zagros mountains, Marivan

sits at an average elevation of 1,287 meters above sea level and features varied topography, from mountains to valleys, shaped by its cross-border proximity with Iraq (Rahimi, 2023).

<figure><img src="image_90.png" /><figcaption>Figure 1: Study area: Marivan and Sarvabad, Kurdistan province, Iran.</figcaption></figure>

The climate in this region is SM, characterized by cold winters and hot summers. Annual rainfall averages 991 mm, with considerable variability (±235 mm). Humidity averages around 54%, which is insufficient to sustain green pastures during summer, leaving behind dry grass that elevates fire risks (Rahimzadeh, 2008).

## 2.2 Data Collection and Preparation

Table 1 summarizes the data used in this study. Sentinel-2 satellite imagery and related products provided high-resolution multispectral data essential for environmental analysis and fire risk assessment.

<table><caption>Table 1: Data sources used.</caption><thead><tr><th>Data Type, Projection System, Spatial Resolution (m)</th><th>Time Period</th><th>Data Source</th></tr></thead><tbody><tr><td>Sentinel 2 (13 bands), UTM<sup>1</sup>, 10</td><td>13/05/2020 and 15/09/2020</td><td rowspan="2">ESA<sup>3</sup></td></tr><tr><td>Sentinel Burned Area Product, UTM, 10</td><td>2021 - 2023</td></tr><tr><td>ZGI, UTM, 10</td><td>2020</td><td>(Rahimi, 2024)</td></tr><tr><td>DEM/SRTM<sup>2</sup>, UTM, 10</td><td>-</td><td>USGS<sup>4</sup></td></tr></tbody></table><ol><li>Universal Transverse Mercator</li><li>Digital Elevation Model/Shuttle Radar Topography Mission</li><li>European Space Agency</li><li>United States Geological Survey</li></ol>

## 2.3 Methodology

The collected data layers were first stacked using Python 3.12, ensuring spatial alignment and consistency across all layers. To integrate Sentinel-2 and ZGI layers, the datasets were clipped to the boundaries of the study area using a defined polygon.
The clipped layers were then aligned to ensure
consistent spatial resolution and extent. Finally, the
processed layers were stacked into a multi-band raster,
enabling the combination of spectral information and
vegetation index data for subsequent fire susceptibility
analysis. This approach ensured a precise and
comprehensive dataset for modeling and assessment.
These layers were integrated into a matrix consisting
of all raster layers, which included Sentinel-2 bands
(13 bands), elevation, and the ZGI. The resulting
matrix had dimensions of 15×8020×6557, meaning
that 15 different values represented each pixel in the
study area. Figure 2 illustrates the methodology and
process flow of this approach.

<figure><img src="image_91.png" /><figcaption>Figure 2: Methodological framework.</figcaption></figure>

This study employed PCA and K-means to detect
fire-prone areas effectively. To evaluate the
effectiveness of these models, the Leave-One-Out
(LOO) method was implemented, allowing for a
thorough assessment of how different data layers
contributed to fire susceptibility mapping. The
validation of the results with post-2020 data (from

2021–2023) demonstrated the reliability of these
methods for accurately predicting fire-prone areas
over time, showcasing their potential for long-term
wildfire risk assessment and management.

K-means and PCA were chosen as they are
common unsupervised methods suitable for studies
with limited field data. Unlike supervised learning
approaches, which require extensive datasets and
computational resources, these methods provide a
practical and efficient framework for fire
susceptibility classification. The components were
retained to classify fire susceptibility into high,
average, and low categories. This decision reflects the
study area's vegetation and land cover characteristics
alongside the limited field data available. This
approach balances capturing sufficient data variance
with the need for interpretability, ensuring robust and
practical classification.

To validate the results, the post-2020 burned area
was used. They derived from sentinel 2 by applying
the Normalize Burned Ratio (NBR) on data from
2021 to 2023 (Giddey, 2021). The NBR index proved
to be a valuable tool for identifying burned areas by
leveraging spectral differences between vegetation
and charred surfaces. This methodology ensured
reliable validation of the predictive models and
reinforced the importance of integrating remote
sensing indices into fire susceptibility research.

3 RESULTS

Figure 3 shows the results of applying PCA and K-means. The results revealed varying levels of effectiveness (High, Average, and Low) in predicting fire-prone areas in the study region using 2020 data. The black polygons are the burned areas calculated from Sentinel 2 data from 2021 to 2023, using the NBR index.

These patterns emphasize the spatial variability of
fire-prone areas, which is influenced by factors such
as vegetation type, fuel density, and topographic
conditions, as suggested by earlier studies (Rahimi et
al., 2024). The classification levels in this study were
established based on the number of components
selected for each of the unsupervised methods, PCA
and K-means, with three components defined for
optimal model performance in both cases. The trade-
off between preserving critical information and
minimizing computational complexity informed this
selection process for the number of components. By
reducing the dimensionality of the data, both PCA
and K-means enhanced interpretability while
maintaining the integrity of key predictive features
(Giddey, 2021).
<figure><img src="image_92.png" /><figcaption>Figure 3: Maps generated by a) PCA and b) K-means according to data in 2020. The black polygons demonstrate the burned areas from 2021 to 2023.</figcaption></figure>

The visual representation in Figure 3 highlights distinct patterns of fire susceptibility zones, with clear differentiation between regions of varying risk levels.

Among the methods, PCA produced the most accurate and reliable results, as shown in Figure 3a,

compared to K-means (Figure 3b). The red areas show high-fire susceptible areas, while the yellow and green colors show the average and low-fire susceptible areas. The stark contrast between the red zones and the remaining areas underscores the ability of PCA to highlight regions requiring urgent attention. This distinction aligns with practical applications in fire management, where prioritizing high-risk zones can significantly enhance resource allocation and mitigation efforts.

The validation process compared high-susceptibility areas identified by the models with actual burned areas in subsequent years. PCA exhibited a higher overlap between its predicted high fire-susceptible zones and the areas affected by post-2020 wildfires, with approximately 41%. In contrast, this value was 31.3% for the K-means (Figure 4). This significant difference in overlap suggests that PCA's ability to identify subtle patterns in the data gives it a clear edge over K-means. The higher accuracy of PCA aligns with its proven effectiveness in environmental modeling and its capacity to extract critical components from high-dimensional datasets. Regarding High and Average fire-prone areas together, both methods offered almost the same results. This indicates that PCA outperformed the K-means models in accurately detecting fire-prone regions, highlighting its robustness in mapping fire susceptibility. However, the similar performance of the two methods in identifying combined High and Average susceptibility areas suggests that K-means may still hold utility in broader-scale applications where granular accuracy is less critical. This underscores the need for selecting methods based on specific use-case requirements.

<figure><img src="image_93.png" /><figcaption>Figure 4: Detection of fire-prone areas for the post-2020 years by 2020's data using PCA and K-means.</figcaption></figure>

Moreover, the analysis indicated that certain Sentinel-2 spectral bands could be excluded without substantially reducing the model's accuracy. This
finding highlights the efficiency of PCA in dimensionality reduction, where less critical spectral bands are discarded while retaining the most informative features. Such optimization reduces computational demands and facilitates faster data processing, which is crucial for large-scale studies (Kantarcioglu, 2023). This suggests that PCA effectively reduced dimensionality, preserving only the most informative features while maintaining high predictive accuracy. The ability to maintain predictive accuracy while reducing data complexity makes PCA an attractive option for fire susceptibility mapping, particularly in resource-constrained settings where computational efficiency is essential. The fact that the PCA model consistently provided reliable predictions across multiple years demonstrates its potential for long-term application in fire susceptibility mapping.

This consistency underscores PCA's robustness and adaptability, especially in evolving environmental conditions and fire patterns. Its reliability across diverse temporal datasets reaffirms its suitability for integration into long-term fire management frameworks.

## 4 CONCLUSIONS

In conclusion, this study highlights the effectiveness of using RS data combined with unsupervised image classification techniques to identify fire-prone areas in the Kurdo-Zagrosian forests. Among the two methods tested, PCA performed better in accurately predicting fire-susceptible zones, with 80% of the burned areas from 2021 to 2023 correctly classified as moderate to high-risk. Its ability to reduce dimensionality while preserving critical information, especially when excluding less relevant Sentinel-2 bands, further solidifies PCA as a robust tool for long-term fire risk mapping. In contrast, K-means showed moderate success, identifying around 50% of the burned areas. These findings emphasize the potential of PCA for improving fire management strategies, while K-means may require further refinement to achieve similar predictive accuracy.

In the future, more unsupervised classification methods will be compared to further evaluate their effectiveness and explore potential improvements in fire risk prediction.

## ACKNOWLEDGMENTS

The work is funded by national funds through FCT – Fundação para a Ciência e Tecnologia, I.P., in the framework of the UIDB/04683/2020 (https://doi.org/10.54499/UIDB/04683/2020) and UIDP/04683/2020 633 (https://doi.org/10.54499/UIDP/04683/2020) – Instituto de Ciências da Terra programs.

## REFERENCES

Bowman, M.J.S.D.; Moreira-Munoz, A.; Kolden, C.A.; Chavez, R.O.; Munoz, A.A.; Salinas, F.; Gonzalez-Reyes, A.; Rocco, R.; de la Barrera, F.; Williamson, G.J.; et al. Human- Environmental Drivers and Impacts of the Globally Extreme 2017 Chilean Fires. *Ambio* 2018, doi:10.1007/s13280-018-1084-1

Celik, N. T. (2009). Unsupervised Change Detection in Satellite Images Using Principal Component Analysis and $sk$-Means Clustering. *IEEE Geoscience and Remote Sensing Letters*, *6*(4), 772–776. https://doi.org/10.1109/lgrs.2009.2025059

Duarte, L., Teodoro, A., Monteiro, A.T., Cunha, C., & Gonçalves, H., "PhenoMetrics: An open source software application to assess vegetation phenology metrics," *Computers and Electronics in Agriculture* 148, 82–94 (2018); https://doi.org/10.1016/j.compag.2018.03.007.

ESA European Space Agency (ESA), "Copernicus Open Access Hub, Sentinel-2 Level-2A products," https://scihub.copernicus.eu.

Ghorbanzadeh, O.; Valizadeh, K.K.; Blaschke, T.; Aryal, J.; Naboureh, A.; Einali, J.; Bian, J. Spatial prediction of wildfire susceptibility using field survey GPS data and machine learning approaches. *Fire* 2019, *2*(3), 43. https://doi.org/10.3390/fire2030043.

Giddey, B. L., Baard, J. A., & Kraaij, T. (2021). Verification of the differenced Normalised Burn Ratio (dNBR) as an index of fire severity in Afrotemperate Forest. South African Journal of Botany, 146, 348–353. https://doi.org/10.1016/j.sajb.2021.11.005

Kantarcioglu, O.; Schindler, K.; Kocaman, S. Forest fire susceptibility assessment with machine learning methods in north-east Türkiye. *Int. Arch. Photogramm. Remote Sens. Spatial Inf.* Sci. 2023, XLVIII-M-1-2023, 161–167. https://doi.org/10.5194/isprs-archives-xlviii-m-1-2023-161-2023.

Mishra, M.; Guria, R.; Baraj, B.; Nanda, A.P.; Santos, C.A.G.; Da Silva, R.M.; Laksono, F.A.T. Spatial analysis and machine learning prediction of forest fire susceptibility: A comprehensive approach for effective management and mitigation. *Sci. Total Environ.* 2024, 926, 171713. https://doi.org/10.1016/j.scitotenv.2024.171713.

Piao, Y.; Lee, D.; Park, S.; Kim, H.G.; Jin, Y. Forest fire susceptibility assessment using Google Earth Engine in Gangwon-do, Republic of Korea. Geomatics Nat.
*Hazards Risk* 2022, 13, 432–450. https://doi.org/
10.1080/19475705.2022.2030808.

Rahimi, I., Duarte, L., & Teodoro, A. C. (2024). Zagros Grass Index—A new vegetation index to enhance fire fuel mapping: A case study in the Zagros Mountains. *Sustainability*, 16(10), 3900. https://doi.org/10.3390/su16103900

Rahimzadeh Bajgiran, P., Darvishsefat, A., Khalili, A., & Makhdoum, M., (2008). "Using AVHRR-based vegetation indices for drought monitoring in the Northwest of Iran," *Journal of Arid Environments* 72, 1086–1096

Rihan, M.; Bindajam, A.A.; Talukdar, S.; Shahfahad, Naikoo, M.W.; Mallick, J.; Rahman, A. Forest fire susceptibility mapping with sensitivity and uncertainty analysis using machine learning and deep learning algorithms. *Advances in Space Research* 2023, 72(2), 426–443. https://doi.org/10.1016/j.asr.2023.03.026.

Saha, S.; Bera, B.; Shit, P.K.; Bhattacharjee, S.; Sengupta, N. Prediction of forest fire susceptibility applying machine and deep learning algorithms for conservation priorities of forest resources. Remote Sensing Applications: Society and Environment 2023, 29, 100917. https://doi.org/10.1016/j.rsase.2022.100917.

Sharma, D., Sharma, T., & Singhai, J. (2021). Extraction of Water Bodies from Visible Color Satellite Images Using PCA Feature Map. *2021 IEEE International India Geoscience and Remote Sensing Symposium (InGARSS)*, 1–4. https://doi.org/10.1109/ingarss515 64.2021.9791857

Teodoro, A.C., Duarte, L., (2013).Forest fire risk maps: a GIS open source application a case study in Northwest of Portugal, *International Journal of Geographical Information Science* 27(4), DOI:10.1080/13658816.20 12.721554.

Theodoridis, S., & Koutroumbas, K. (2009). Clustering Algorithms III: Schemes based on Function Optimization. *In Elsevier eBooks* (pp. 701–763). https://doi.org/10.1016/b978-1-59749-272-0.50016-5

USGS EROS Archive - Digital Elevation - Shuttle Radar Topography Mission (SRTM), https://www.usgs.gov/centers/eros/science/usgs-eros-archive-digital-elevation-shuttle-radar-topography-mission-srtm-1.

Zema, D.A.; Nunes, J.P.; Lucas-Borja, M.E, (2020). Improvement of Seasonal Runoff and Soil Loss Predictions by the MMF (Morgan-Morgan-Finney) Model after Wildfire and Soil Treatment in Mediterranean *Forest Ecosystems*. *Catena*, 188, 104415, doi: 10.1016/j.catena.2019.104415
# Urban Growth in Metropolitan Regions Using Dynamic Modeling by Cellular Automata: A Comparative Analysis Between Brazil and Portugal

Elizabeth Maria Feitosa da Rocha de Souza<fnref n="1" /> IDª, Antonio Alberto Teixeira Gomes<fnref n="2" /> IDᵇ
and Vandre Soares Viegas<fnref n="1" /> IDᶜ

**Keywords:** Cellular Automata, Spatial Dynamic Models, Urban Areas, Metropolitan Regions.

**Abstract:** This study synthesizes the outcomes of land use changes obtained through the implementation of dynamic modeling by cellular automata across two metropolitan regions in Portugal and Brazil. The purpose is to analyze the primary findings acquired, considering the particularities of each nation and evaluate the potentialities of the used data. The study examined the metropolitan regions of MRRJ (Rio de Janeiro, Brazil) and AMP (Porto, Portugal). Modifications were implemented in the DinamicaEgo software to the fundamental data representing static and dynamic variables for each context. The findings revealed a substantial increase of urban areas in the MRRJ, and the modeling demonstrated its applicability across the two contexts, considering the requisite modifications for the data accessible in each country.

## 1 INTRODUCTION

Cellular models in urban modelling became implicit in early computer models of the 1960s, e.g., the work of Chapin and Weiss for North Carolina (Chapin and Weiss, 1968) Waldo Tobler's model for Detroit in 1970 (Tobler, 1970), and the model developed by Couclelis for Los Angeles in 1989 (Couclelis, 1989).

Cellular automata (CA) have emerged as the prevailing architectural design employed in spatial simulation models. Cellular automata, as defined by Stephen Wolfram, are mathematical representations of physical systems in which time and space are discrete (Wolfram, 1983). These representations take the form of a regular uniform grid, where each location contains a discrete variable. In accordance with a predetermined set of local rules and the values of the variables in their immediate vicinity at the previous time increment, the variables in each cell are updated simultaneously. CA models have the capability to emulate various concerns associated

with land use changes. They can be used to create self-modifiable models for high-resolution urban land use dynamics, prototypes to simulate land conversion via integrated Geographic Information Systems (GIS) (Clarke, et al 1997; White et al, 1997; Wu, 1998) and have numerous applications for simulating scenarios involving urban dynamics (Barredo et al, 2003; White, and Engelen, 1993). The use of dynamic models can support the definition of environmental public policies that involve the assessment of land use and land cover as a fundamental component to understand changes that may result from combined biophysical and socioeconomic factors, with both short-term and long-term impacts (Meyer and Turner, 1992).

In this context, the research aimed to identify potential scenarios of land use and land cover change until the year 2050 in the Metropolitan Area of Porto (AMP) and the Metropolitan Region of Rio de Janeiro (MRRJ), specifically within the boundaries of urbanized territories, while observing changes between classes and specificities.

<footnote n="1">¹Federal University of Rio de Janeiro, Department of Geography, Graduate Program in Geography Geosciense Institute, Rio de Janeiro, Brazil</footnote>

<footnote n="2">²University of Porto, Faculty of Arts and Humanities, Porto, Portugal
{elizabethmfr, albgomes, v.vegas18}@gmail.com</footnote>

<footnote>ªID https://orcid.org/0000-0002-2672-8169</footnote>

<footnote>ᵇID https://orcid.org/0000-0002-1764-0659</footnote>

<footnote>ᶜID https://orcid.org/0000-0002-2880-8314</footnote>
## 2 STUDY AREA

The study area encompasses two Metropolitan Regions: Metropolitan Area of Porto (AMP) and the Metropolitan Region of Rio de Janeiro (MRRJ - Figure 1). Although referred to by different acronyms—AMP for the Metropolitan Area of Porto in Portugal and MRRJ for the Metropolitan Region of Rio de Janeiro in Brazil—both regions fulfil analogous roles within their respective national contexts. In 2021, the Porto Metropolitan Area (AMP), composed of 17 municipalities, had a population of approximately 1,7 million inhabitants distributed across an area of 2,040 km². Meanwhile, the Rio de Janeiro Metropolitan Region (RMRJ), consisting of 22 municipalities, was home to around 12 million inhabitants in an area of 6,700 km², with significantly higher population density due to the concentration of people in municipalities like Rio de Janeiro, São Gonçalo, and Duque de Caxias. These regions represent important urban hubs in their respective countries, facing distinct challenges in terms of planning and urban development.

Each area functions as a logistical hub for regional planning, economic development, and the coordination of urban services. Both regions represent densely populated areas that support integrated infrastructure alongside public policy frameworks. This functional alignment underscores the significance of metropolitan regions as foundational elements in urban planning and governance in each country.

<figure><img src="image_94.png" /><figcaption>Figure 1: Study area - Areas considered in the study for comparison purposes</figcaption></figure>

Situated in southeastern Brazil, the Metropolitan Region of Rio de Janeiro (MRRJ) is an extremely urbanized, multifaceted, and diverse region. This region, which includes the capital city, Rio de Janeiro, and significant neighbouring municipalities including São Gonçalo, Duque de Caxias, and Niterói, is renowned for its geographical, historical, and economic importance. Based on a survey carried out

in 2022 by the Brazilian Institute of Geography and Statistics (IBGE, 2022), the Metropolitan Region of Rio de Janeiro is the second most populous metropolitan area in Brazil, comprising 74% of the total population of the state of Rio de Janeiro, with an estimated 12 million inhabitants.

The Metropolitan Area of Porto is the second largest in Portugal, encompassing 17 municipalities around the city of Porto. According to data from the Land Use and Occupation Map - COS2022 (Instituto Nacional de Estatística, 2022).

The AMP has emerged as a strategic area for advancing mobility and sustainability projects, essential for the region's future.

## 3 MATERIALS AND METHODS

The CA model used in this work is implemented in the Dinamica EGO (Environment for Geoprocessing Objects) platform 7.4, a software developed by researchers from Federal University of Minas Gerais (UFMG). The methodology used is presented according to the scheme shown in Figure 2.

<figure><img src="image_95.png" /><figcaption>Figure 2: Workflow (*WoE – Weights of Evidence).</figcaption></figure>

Simulation models utilize two types of explanatory variables, divided between Static and Dynamic. The former remains unchanged at each temporal step of the model, while the latter changes with each iteration based on the distance from the central axis being analyzed.

The first step in the modeling process is parameterization using the Weight of Evidence method, following the Bayesian approach. This
involves assigning weights to different input variables based on their relevance to the specific problem.

Historical transition matrices are calculated. The transition matrix describes a system that changes in discrete increments of time. In the Dinamica EGO platform, matrices are calculated using a Markovian model, which combines the Markov Chains technique with Cellular Automata. Dinamica EGO also employs Markov chains to determine the amount of change, as well as cellular automata to reproduce patterns of these changes from probability maps, which are calculated using the Weight of Evidence statistical method (Soares Filho *et al*, 2002). Cellular automata (CA)-based models have been widely used due to their ability to simulate dynamic spatial patterns. The choice of the CA-Markov model in this study is justified by its integration of Markov transition matrices with dynamic spatial modeling, enabling the capture of land use changes with high temporal and spatial granularity.

For the desired period—a range of years—two types of matrices are generated: a global matrix, representing the transition rates for the entire training period, and a multistep matrix, which reflects annual changes. The global matrix aggregates all transitions across the specified period, while the multistep matrix allows for more granular modeling by representing changes on an annual basis.

It is important to clarify which transitions are being modeled, as this directly influences the simulation outcomes. In DINAMICA, transitions are managed through two Cellular Automata (CA) algorithms: the Patcher, which simulates the aggregation of land patches, and the Expander, which models the expansion of existing areas. These mechanisms ensure that spatial patterns align with the observed dynamics, providing a robust framework for projecting future scenarios.

In Dinamica EGO, a set of sub-regional functors is used to process the data separately in each sub-region. In this model, transitions are simulated annually for each of the subdivisions, dividing a map into parts for separate data processing and then combining the results. This allows the modeler to define operations that should be applied only to specific sub-regions or different parameters and coefficients. The result is a model that respects the regional context. In both cases, the transitions of interest focus on changes to urban, vegetation, and cultivated areas.

The dynamic variables considered in this study are exclusively related to land use and land cover in the years specified for each region: AMP and MRRJ.

These variables reflect temporal changes that directly influence the projected scenarios and are updated in each iteration of the model. The static variables, on the other hand, remain constant throughout the modeling process and represent potential factors influencing the observed and projected land use changes. These variables include hydrography, transportation systems, and topography, which play a fundamental role in defining spatial patterns of urban expansion and enable the simulation of future scenarios. The dynamic variables include distances from existing urban areas, which are updated at each iteration of the model, while static variables, such as elevation and the hydrographic network, remain unchanged. Dinamica EGO uses Markov matrices to determine transition rates and spatial patterns based on weights of evidence calculated using the Bayesian method.

In summary, the calibration was conducted using the Expander and Patcher functions, which simulate, respectively, the expansion of existing patches and the formation of new urban patches. To address the specificities of the study area, the model was regionalized, dividing the territory into subareas with distinct characteristics and adjusting parameters for each region. This process allowed for greater accuracy in annual simulations, respecting regional contexts and producing results tailored to the geographic complexity of each studied area. Validation was based on spatial similarity between simulated and observed maps from 2021 (RMRJ) and 2019 (AMP), using fuzzy methods and moving windows. The simulation was regionalized to account for local characteristics of the territory, such as differences between rugged terrain areas and lowlands. The model annually simulates transitions, adjusting parameters to reflect local dynamics and generate accurate predictions of urban expansion.

## 3.1 Metropolitan Area Porto – Portugal - Model

The Metropolitan Area of Porto (AMP) was modeled utilizing the subsequent variables: land use and occupation, roads and railways, hydrography, elevation and slope.

The elevation and slope values were derived using the USGS-provided NASADEM Digital Elevation Model (30m spatial resolution).

At a 1:10,000 scale, the road network is an integral component of the National Motorway Network, while the railway network is an extension of the National Railway Network.
The hydrographic network is a 1:200,000 scale representation of the principal channels and rivers extracted from the Hydrography Network of Continental Portugal. Only the main rivers and channels were considered. While pertinent to the simulation process, this class is regarded as limiting in the modeling process for the expansion of the urbanized class.

In the study, the Land Use and Occupation Map (COS) for 2015 and 2018 was considered for the preliminary analysis, focusing on the base period 2015-2018. The primary categories analyzed were vegetated cover, agricultural regions, and urbanized areas. To project future scenarios, the model incorporated simulations for specific years: 2030, 2040, and 2050. A cross-tabulation analysis was conducted using the converted matrix dataset, maintaining the same spatial resolution to assess changes over time and validate the model's outputs.

## 3.2 Metropolitan Region in Rio De Janeiro/Brasil - Model

In MRRJ, the initial and final maps come from the time series provided by Mapbiomas (collection 7.1), including the years 2016 and 2021, that is, a period of 05 years, con-sidered sufficient to identify areas of urban expansion.

The MapBiomas 7.1 collection includes annual land use and land cover data for the period 1985 to 2021. The secondary data taken from the platform have some spatial inconsistencies in the thematic classifications, often due to noise and mainly in classes and transversal themes, e.g., agriculture and pasture. Currently, there are 30 classes available, from macro-classes 1 - Forest; 2 - Non-Forest Natural Formation; 3 - Farming; 4. Non vegetated area; 5 - Water. The static variables were obtained from official Brazilian agencies (hydrography, road system, and topography).

The first stage of adjusting the Mapbiomas data consists of correcting incorrect transitions that normally occur in the form of isolat-ed pixels and have no correspondence in the real world, for example, a pixel classified as urban infrastructure should not transition to the water class.

The established rule was that no pixel originally classified as urban infrastructure (class 24) - Areas with significant density of buildings and roads - should transition to other classes. After the correction, the other subclasses other than urban infrastructure were converted to the “Others” class, maintaining binary use maps.

Finally, the last phase involves data calibration, validation, and scenario observation. The maximum value of the Fuzzy similarity index is considered, using exponential decay for window sizes of 3x3, 5x5, 7x7, 9x9, and 11x11 as a validation method. The index varies from 0 to 1, where 1 indicates perfect spatial agreement.

Kernel Density Estimation was applied in the study to enhance the identification of change-prone areas.

# 4 RESULTS

Regarding the parameterization of the transition algorithms, Expander and Patcher, the model achieved better validation results with a higher percentage assigned to the Expander algorithm: 0.92% (RMRJ) and 0.87% (AMRJ).

Considering the threshold mentioned in (Soares Filho et al, 2002), where values close to 0.6 in the similarity index with exponential decay indicate adequate spatial congruence between simulated and real maps, it is reasonable to conclude that, in terms of suitable similarity, the period from 2015 to 2018 (AMP) and from 2016 to 2021 (MRRJ) showed satisfactory adjustments, especially in the 3x3 window.

## 4.1 Metropolitan Region in Rio De Janeiro/Brasil

To Metropolitan Region of Rio de Janeiro Based on the model, the following areas are identified as particularly promising for urban expansion in 2050: Situated west of the Capital, on the boundary between the Santa Cruz and Guaratiba communities. Additional noteworthy localities within the Capital comprise Campo Grande, Vila Militar, Campo dos Afonsos, Vargem Grande, and Recreio dos Bandeirantes. These are areas with available spaces for urban occupation and increased construction. Furthermore, the Baixada Fluminense encompasses municipalities Duque de Caxias, Magé, Nova Iguacu, and Seropedica, all of which exhibit substantial prospects for development.

In the East Fluminense, Maricá, Niterói, and São Gonçalo are the municipalities that are undergoing the most significant urban development. In the Mountain subregion, Petrópolis municipality observes noteworthy urban expansion in the vicinity of the central area, which encompasses the Vila Militar and Valparaíso neighborhoods. In the southern region, notable residential areas include
Independência and Quitandinha. Notable are the regions of Mosela to the west and Bonsucesso and Itaipava to the north. The aforementioned regions symbolize the anticipated epicenters of urban growth as predicted by the model for the year 2050 (Figure 3).

<figure><img src="image_96.png" /><figcaption>Figure 3: Urbanization by subregion through 2050, represented graphically by absolute and relative values.</figcaption></figure>

Kernel density estimation (Figure 4) was employed to depict urban expansion in 2050. This interpolator enables the determination of the magnitude of a particular occurrence throughout an entire region (Freire, 2015). Druck *et al*, (2004) describe how the Kernel density estimator creates a surface whose value is proportional to the sample intensity per unit area and fits a two-dimensional function to the considered events in order to estimate the point intensity of the process across the entire study region. Kernel Density estimates the intensity of events (in this case, urbanization) per unit of area, generating a continuous surface that highlights regions with higher concentrations of change. This technique complemented the analysis based on transition matrices, enabling the mapping of regions with significant urban change. The use of this methodology was essential due to the high variability in land use in the MRRJ, providing a more detailed

<figure><img src="image_97.png" /><figcaption>Figure 4: Urban Expansion Scenario for 2050 in MRRJ.</figcaption></figure>

view of urbanization hotspots and enhancing future scenarios for urban planning and environmental management.

Municipalities outside the capital (Rio de Janeiro) highlight a current trend of urban expansion and demand for housing and new developments.

## 4.2 Metropolitan Area Porto - Portugal

The Figure 5 highlights the results obtained after modeling for the years 2030, 2040, and 2050 with a detailed analysis of the evolution.

<figure><img src="image_98.png" /><figcaption>Figure 5: Future Landuse Scenarios for AMP in 2030, 2040, and 2050.</figcaption></figure>

A 2.66 percent increase in urbanized territories is observed across the entire study area from 2030 to 2040, and a 2.55 percent increase from 2040 to 2050. The projected increases for agricultural regions are 1.13% and 1.08%, correspondingly. On the contrary, a decline in forested regions is noted, fluctuating by approximately 1.8% during both time periods.

The Figure 6 delineates the regions where anticipated alterations occur during the examined time frames. There is a total of 2274 polygons. From 2023 to 2040, there is a projected total of 1097 regions undergoing changes, namely involving the conversion of brush areas into urbanized zones and the transformation of brush areas into agricultural land.
<figure><img src="image_99.png" /><figcaption>Figure 6: Regions with alterations in landuse during 2023 – 2040 and 2040 and 2050.</figcaption></figure>

There is a significant concentration of expansion near the cities of Vila Nova de Gaia, Matosinhos, and Gondomar, which are already densely urbanized areas close to the metropolitan center of Porto. This growth suggests a continuous trend of urbanization in peri-urban areas and an expansion into zones with existing infrastructure.

In the period from 2040 to 2050, growth appears more dispersed, with new urbanization points emerging in regions such as Paredes, Santo Tirso, and Oliveira de Azeméis, indicating a possible gradual decentralization of the observed expansion. This pattern suggests a potential saturation of central areas and a search for new development zones on the periphery of the AMP, leading to a more balanced growth between municipalities in the northern and southern parts of the Metropolitan Area. These data highlight the need for integrated urban planning to meet the growing demands for infrastructure and services in these emerging regions.

Figure 7 highlights the projection for the Serra do Monte area (Northern Region of the AMP), with three representations: the 2050 projection showing urbanized area growth; a 2023 GE image as a

reference; and a combination of both, indicating urban area densification and anticipated changes.

<figure><img src="image_100.png" /><figcaption>Figure 7: Changes in 2050 for Serra do Monte, Northern Region of the AMP.</figcaption></figure>

Figure 8 demonstrates the predominant vegetation scenario near the A28 highway, which may be impacted by increased urban construction by 2050. The imagem also highlights an on-site view of the area. The red arrow indicates the direction of the photo's viewpoint.

<figure><img src="image_101.png" /><figcaption>Figure 8: Changes in 2050 for Póvoa de Varzim.</figcaption></figure>

Vila Nova de Gaia, located south of Porto, is expected to expand residential growth in areas with exposed soil and available land, as highlighted in Figure 9.

The results presented for the Metropolitan Regions of Rio de Janeiro (RMRJ) and the Metropolitan Area of Porto offer distinct perspectives regarding urban expansion projections. In the RMRJ, there is a notable focus on specific expansion areas, particularly in sub-regions with significant growth, such as the Baixada Fluminense and the East Fluminense. In the case of the Metropolitan Area of Porto, direct changes in land use stand out, detailing transformations in agriculture and urbanized territories. The conversion of areas from scrubland and agriculture into urbanized territories is emphasized. In both contexts, the results converge by
<figure><img src="image_102.png" /><figcaption>Figure 9: Changes in 2050 for Vila Nova de Gaia.</figcaption></figure>

highlighting challenges related to urban growth, such as inadequate infrastructure and environmental impacts, underscoring the need for sustainable public policies and careful urban planning in both regions.

# 5 CONCLUSIONS

This study presents a synthesis of the outcomes achieved by employing cellular automata to perform dynamic modeling in multiple metropolitan regions of Portugal and Brazil. The research centers on the primary discoveries, taking into account the specific characteristics of each nation and the data utilized, while analyzing the metropolitan areas of Rio de Janeiro (RMRJ) and Porto (AMP). The DinamicaEgo software was modified with respect to the fundamental data that represented static and dynamic variables in each respective context. The findings unveiled a significant rise in RMRJ, thereby illustrating the model's versatility in different settings, contingent upon the data requirements of each country.

The findings pertaining to urban expansion projections for the Metropolitan Region of Rio de Janeiro (RMRJ) and the Metropolitan Area of Porto (AMP) present unique and discernible viewpoints. RMRJ places considerable emphasis on particular expansion zones, particularly in rapidly developing subregions like Baixada Fluminense and East Fluminense. Particular attention is paid to explicit alterations in land use within AMP, which pertain to

agricultural and urbanized regions. Both contexts emphasize the transformation of agricultural and scrubland regions into urbanized ones. The findings converge in that they emphasize the environmental impacts and insufficient infrastructure associated with urban expansion, thereby emphasizing the necessity for sustainable public policies and meticulous urban planning in both areas.

In future stages of the research, new variables are intended to be included in the modeling process. The boundaries of protected areas in the AMP may improve the delineation of urban areas for future decades.

# ACKNOWLEDGEMENTS

We would like to express our sincere gratitude to CNPq (National Council for Scientific and Technological Development) and FAPERJ (Foundation for Research Support of the State of Rio de Janeiro) for their financial support, which was essential to the successful development of this research. Their contributions have played a crucial role in enabling us to carry out this work.

# REFERENCES

- Barredo, J.; Kasanko, M.; McCormick, N.; Lavalle, C. (2003) *Modelling dynamic spatial processes: simulation of urban future scenarios through cellular automata*. Landscape and Urban Planning, 64, 145-160.

- Chapin, J.R.; Weiss, S.F. (1968). *A probabilistic model for residential growth*. Transportation Research, v.2 p 375-390.

- Clarke, K.; Hoppen, S.; Gaydos, L. (1997) *A self-modifying cellular automaton model of historical urbanization in the San Francisco Bay area*. Environment and Planning B, 24, 247-261.

- Couclelis, H. (1989) *Macrostructure and Microbehavior in a Metropolitan Area*. Environment and Planning B, v.16, p.141-154.

- Druck, S., Carvalho, M. S., Câmara, G., & Monteiro, A. V. M. (2004). *Análise Espacial de Dados Geográficos*. Brasília: Embrapa.

- Freire, F. H. M. A. (2015). *Introdução à Estatística Espacial*. Recuperado de http://www.observatoriodas metropoles.ufrj.br/download/aulasanalise-espacial.pdf/

- Instituto Brasileiro de Geografia e Estatística (IBGE). (2022). *Censo results 2022*. Acesso https://censo2022.ibge.gov.br/panorama/

- Instituto Nacional de Estatística (INE). (2022). *Land Use and Occupation Map (COS2022)*. Lisboa: INE.
Meyer, W. B., & Turner, B. L. (1992). *Human population growth and global land-use/cover change*. Annual review of ecology and systematics, 23(1), 39-61.

Soares-Filho, B. S., Cerqueira, G. C., & Pennachin, C. L. (2002). *DINAMICA—a stochastic cellular automata model designed to simulate the landscape dynamics in an Amazonian colonization frontier*. *Ecological Modelling*, 154(3), 217-235.

Tobler, W. R. (1970). *A computer movie simulating urban growth in the Detroit region*. Economic Geography, 46(2), 234–240.

White, R.; Engelen, G.; Uljee, I. (1997) *The use of constrained cellular automata for high-resolution modeling of urban land use dynamics*. Environment and Planning B, 24, 323–343.

White, R.; Engelen, G. (1993). *Cellular automata and fractal urban form: a cellular modelling approach to the evolution of land use patterns*. Environment and Planning A, 25, 1175–1199.

Wolfram, S. (1983) *Statistical mechanics of cellular automata*. Review of modern physics, v. 55, p.601-643.

Wu, F. (1998) *SimLand: A prototype to simulate land conversion through the integrated GIS and CA with AHP-derived transition rules*. International Journal of Geographical Information Systems, 12(1), 63–82.
# Navigating Boundary Discrepancies in SAD69-Based Delimitation:
A Case Study and Practical Guidelines

L. R. Costa<fnref n="1" /> IDª, H. S. Delabary<fnref n="1" /> IDᵇ and R. Z. Araujo<fnref n="2" /> IDᶜ

**Keywords:** Geoprocessing, Environmental Forensics, SAD69, QGIS, ProGriD.

**Abstract:** The delimitation of the 84,130-hectare Serra do Tabuleiro State Park, a fully legally protected Conservation Unit in Santa Catarina, Brazil, relies on shapefiles provided by the managing institution and based on the SAD69 Coordinate Reference System. However, user-defined parameters when handling these shapefiles may result in up to three slightly different polygon representations, each affecting the perception of boundaries shared with adjacent territories. This study investigates these polygon discrepancies and assesses which representation most accurately reflects the intended delimitation. Although a definitive solution is not reached, the authors provide valuable recommendations for public authorities and GIS users to standardize interpretations and improve boundary accuracy.

## 1 INTRODUCTION

The last change to the boundaries of a Conservation Unit, a legally protected area henceforth referred to only as a **PARK** took place in 2009, by means of a State Law, which re-evaluated and redefined the already existing Park's boundaries (Figure 1). The law's documentation includes the approximate plane coordinates of the points that define the Park's polygons and their associated areas (SANTA CATARINA, 2009). The law also establishes that the shapefile published by the public institution responsible for managing the Park would be the reference for the purposes of interpreting and delimiting the Park's boundaries (SANTA CATARINA, 2010).

*'Paragraph 3 For the purposes of interpreting and demarcating the defined boundaries... the delimitation in shapefile format prepared by .... is established as the reference.'*

The shapefile made available by the management institution has served its purpose and has become the basis for society's interpretation of the Park's

boundaries. The establishment of administrative or criminal procedures in the event of unauthorized intervention in the protected area depends on the correct and precise interpretation of the boundaries, an interpretation that is a routine activity in forensic examinations involving alleged crimes against the environment.

The Environmental Military Police, municipal environmental bodies with powers of protection and inspection, the Scientific Police, technical assistants, lawyers and society in general all rely on the Park's demarcation polygons to argue their case in proceedings and fines involving non-authorized interventions inside that Protection Unit (Park).

Unfortunately, two versions of the shapefile were published, each one using a different Coordinate Reference System (SAD69 and SIRGAS2000). This paper aims to show that the officially available shapefiles of the Park's delimitation can lead to three polygons that are slightly offset from each other (in the order of meters), which has caused doubts and confusion in cases where the areas under discussion are near the Park's borders.

The possible alternatives that may have produced versions with slightly different borders will be

<footnote n="1">¹Santa Catarina's State Scientific Police, Brazil</footnote>

<footnote n="2">²Santa Catarina Federal Institute, Florianópolis Campus, Brazil
{luciano.costa, henrique.delabary}@policiacientifica.sc.gov.br, renato.araujo@ifsc.edu.br</footnote>

<footnote>a ID https://orcid.org/0000-0003-1940-5166</footnote>

<footnote>b ID https://orcid.org/0009-0007-9105-4141</footnote>

<footnote>c ID https://orcid.org/0009-0006-5965-334X</footnote>
analyzed, using the official planimetric coordinates
provided in the state legislation as a reference source.

Finally, some recommendations will be presented
for professionals working in environmental forensics,
for public agencies and for Geographic Information
Systems' users in general. The study is structured to
address a key practical challenge: the accurate
specification of boundaries, with a focus on a natural
park. The insights and conclusions drawn may prove
helpful in addressing similar challenges in other areas
of environmental management and geodesy

<figure><img src="image_103.png" /><figcaption>Figure 1: Location map of a Conservation Unit (a legally Protected Area) named Serra do Tabuleiro State Park, located in the eastern region of the State of Santa Catarina (Brazil). Area 84,130 hectares. The Park was created in 1975 (IMA, 2024).</figcaption></figure>

## 2 MATERIAL AND METHODS

### 2.1 Coordinate Reference Systems (CRSs) for Brazil

The Coordinate Reference Systems used in Brazil
were (in historical order): Córrego Alegre, Astro
Datum Chuá, SAD69, and currently SIRGAS2000.
SAD69 stands for South American Datum 69.
SIRGAS2000 stands for “Sistema de Referência
Geocêntrico para as Américas 2000”, or Geocentric
Reference System for the Americas 2000.

Initially, the *Córrego Alegre* coordinate reference
system (CRS) was the first horizontal CRS used in
Brazil, remaining in use until the early 1970s. It was
briefly succeeded by *Astro Datum Chuá*, until the
official adoption of SAD69 in 1977. In 1996, an

updated version, known as SAD69(96), was
introduced to refine the system. Then, in 2005, the
IBGE issued a resolution designating SIRGAS2000
as Brazil's official CRS. A transition period of up to
ten years was established, during which
SIRGAS2000 could be used alongside SAD69. This
transition period concluded in 2015

As it is customary within the geoprocessing
community, Coordinate Reference Systems (CRSs)
are represented by EPSG codes, which are numbers
up to five digits that represent and catalogue
definitions of different CRSs. The acronym EPSG
originated from the now-defunct European Petroleum
Survey Group. In practical terms of EPSG
nomenclature, the following apply for Brazil,
considering only the SAD69 and SIRGAS2000
reference systems, and focusing solely on planimetric
coordinates (UTM):

- **A. EPSG 29192** – SAD69 (Classic Network). According to the system description from the EPSG Geodetic Parameter Dataset website (https://epsg.org), it is possible to verify that, in Brazil, this system was replaced by the SAD69(96) system (code 5532, deprecated). Another source (https://epsg.io) identifies EPSG code 5858 as a replacement for code 29192.

- **B. EPSG 5858** – SAD69(96). The obtained description indicates that this system replaces EPSG 29192 and was, in turn, replaced by EPSG code 31982 (SIRGAS2000).

- **C. EPSG 31982** – SIRGAS2000. This is the EPSG code officially used in Brazil for projects and surveys involving geoprocessing.

Brazil presents an additional complexity, as there
is another reference system not previously mentioned.
In addition to the reference systems with cited EPSG
codes (29192, 5858, and 31982), there is also a
reference system known as “SAD69/96 Doppler
Technique or GPS”.

Brazilian Institute for Geography and Statistics
(IBGE) provides a platform (an application) for
transforming coordinates between different official
geodetic reference systems used in the country. The
application (available as both a desktop and online
version) allows coordinate transformations between
the Córrego Alegre, SAD69, and SIRGAS2000
systems. It is a straightforward and objective tool to
support the geospatial data user community in the
transition to SIRGAS2000. The platform/application
is named ProGRID (IBGE, 2009).
In this work, the three transformations in ProGriD that accept SAD69 as the input reference were used. Thus, in addition to the systems with EPSG codes cited in the previous section (EPSGs 29192, 5858, and 31982), the “SAD69/96 Doppler Technique or GPS” reference system is also a possible input reference for data.

## 2.2 Shapefile Sources for the Park

The polygon that represents the Park's boundaries is available for download on a platform maintained by the management body, which we'll call the Map Library (ML polygon). The information is georeferenced and presented in vector, raster and WMS formats, all designed for use in the UTM SIRGAS 2000 coordinate reference system (EPSG:31982).

The Map Library is widely used by government bodies and society in general. One can download the existing files from the Map Library, from which the layer corresponding to the Park's boundaries was extracted, using the SIRGAS2000 Coordinate Reference System (IMA, 2023).

However, the Park's boundaries are also available on the Management Body's website, catering for those users who don't want to or don't know how to navigate the Map Library, preferring to go directly to the desired shapefile. But we point out that this shapefile unfortunately uses the SAD69 Reference System, although it should represent the same official polygon.

Despite the availability and practicality, unfortunately a practical problem was risen, since there are two official sources (from the Management Body) that provide the Park's delimitation polygons: the Map Library and the website, each with a different Coordinate Reference System (SIRGAS2000 and SAD69).

## 2.3 ML and CB Polygons

The first source (Map Library) shows these features using the Coordinate Reference System cataloged as EPSG:31982, which correctly represents the Planimetric Coordinate System recommended for Brazil: UTM 22J SIRGAS 2000. The delimitation obtained will henceforth be referred to simply as the ML Polygon or simply ML.

The second source, which shows the ‘Current Boundaries’ of the Park, henceforth referred to simply as CB Polygon or CB for short, shows the features using the Coordinate Reference System

cataloged as EPSG:29192, which represents the UTM zone 22S SAD69 Planimetric Coordinate System.

At first, it would be indifferent to use either file, as they should produce the same result. However, when importing the two files (ML and CB) into geoprocessing software, a small difference in positioning was observed between the polygons, which are displaced by distances of around one meters(average difference of 1.2 meters). It should be noted as a very important feature, that the importing of the file with EPSG:29192 can be automatically subjected to an on-the-fly transformation by the geoprocessing software, depending on the settings stipulated by the user.

Naturally, this just-a-few-meters-difference is not acceptable for forensic purposes, and therefore the user will decide to ‘reproject’ the website file (CB polygon, EPSG:29192) to EPSG:31928, which is the CRS recommended for official use in Brazil. Routinely, reprojection is recommended so that the file can be correctly used in geoprocessing software to obtain measurements such as distances and areas. Normally, reprojection does not produce any adverse results.

However, when reprojecting the shapefile obtained from the site file (CB polygon) into EPSG:31928, the processing software QGIS Geographic Information System offers TWO options for data transformation (QGIS, 2023).

Using transformation 1, the same polygon produced by the on-the-fly transformation is obtained (with the same difference of around one meter from the Map Library file ML). The polygon resulting from this transformation will henceforth be called CB1 (Figure 2).

If the user chooses transformation 2, it will be obtained a polygon that is displaced by more than three meters in relation to the Map Library polygon (average difference of 3.6 meters). The polygon resulting from this transformation will henceforth be called CB2 (Figure 2).

In fact, to make an informed choice, users should be aware that there are three distinct "materializations" of SAD69: (1) SAD69 Classic, (2) SAD69/96, and (3) SAD69 based on Doppler Technique or GPS. The user must identify which network the maps and databases belong to, and how they were determined, since these 3 networks (SAD69 Classic Network, SAD69/96 and SAD69 Doppler Technique or GPS) have different distortion patterns. Improper use leads to erroneous results (IBGE, 2009).
<figure><img src="image_104.png" /><figcaption>Figure 2: Three polygons obtained for the Park's boundary, depending on the decisions made by the user. The 3 polygons are slightly displaced geographically, in the order of meters.</figcaption></figure>

Two methods were used by our group in our quest to find the correct polygon among the 3 possible polygons (ML, CB1 and CB2). The *first* method was based on reconstructing the polygon from the planimetric coordinates of the Park's borders, presented in the State Law. The *second* method was based on a comparison between the polygons obtained (ML, CB1 and CB2) and some physical remnants of the Park's boundaries, in the form of some existing fragments of a "Park fence", physically existing in a few remaining regions.

# 3 RESULTS

## 3.1 Reconstruction of the Park's Boundary Polygon from the Coordinates of the Points Reported in the State Law

The Park's polygon was reconstructed based on the planimetric coordinates of the existing points in the State Law (706 points for the polygon). Our group faced some basic problems in the definitions of the points. The first divergence is that point 229 was duplicated and showed different coordinates.

The second divergence is that point 299 was not on the list. It was probably the duplicated point 229 that was typed in wrong. It was then assumed that the entry for the second point 229 was actually point 299.

The third divergence is a series of mistakes. From point 327 to point 340 (14 points in total), each record was duplicated. Fortunately, the coordinates were

also duplicated, except for point 335, where the coordinates of the duplicated points were different (differences of approximately 38.7mE and 52.3mN). In this case, it was assumed that the correct coordinates would be from the first entry of duplicate point 335.

The fourth divergence was point 687, which was left unattributed and therefore assumed not to be part of the polygon.

The fifth problem was that the coordinate reference system for which the coordinates are defined is not fully explained in the State Law text, which only mentions the SA69 system (but as we have already seen, there are three possible SAD69 “flavors”).

The purified set of coordinates was then used as input for reprojection using QGIS software and for reprojection using the ProGriD platform. For each reprojection platform, the three possible alternatives for the input coordinate reference system were considered.

### 3.1.1 Reprojection via QGIS

Considering some of the information in the State Law's text, the user could interpret that this is post-1996 SAD69 data, which would lead to the conclusion that the coordinates would be from the SAD69/96 network (EPSG:5858). By reprojecting the file with this choice of EPSG:5858, a polygon is obtained that coincides with polygon CB1 (Current Boundaries with transformation 1, which is the SAD69(96) transformation).

If the user interprets that the coordinates refer to SAD69 Classic Network points (EPSG:29192), the reprojection may result in the CB1 polygon or the CB2 polygon, depending on the transformation option chosen during the reprojection. The possibilities and results were summarized and schematized in Figure 3.

### 3.1.2 Reprojection via ProGriD

To ensure the transformation with the correct parameters, the planimetric coordinates presented in the State Law were also processed using ProGriD. The platform requires an email address to return the results to, when the input file is very large. Examples of input files on the platform's website are useful for properly formatting the input data (IBGE, 2009). Using ProGriD processing on the same input as before, dully formatted, the following results were obtained (see a summary and illustrative diagram in Figure 3):
a. Assuming SAD 69 Classic Network as the correct input reference, points coinciding with the CB2 polygon are obtained.

b. Assuming SAD 69/96 as the correct input reference, points coinciding with the ML polygon are obtained.

c. Assuming SAD 69/96 Doppler Technique or GPS as the correct input reference, points coinciding with the CB1 polygon are obtained.

<figure><img src="image_105.png" /><figcaption>Figure 3: Representative diagram that summarizes all the alternatives considered and discussed in this work. The quest to reconstruct the Park's polygon by importing and reprojecting the coordinates provided, using QGIS software, leads to the CB1 and CB2 polygons described in the text. By using the transformations available in ProGriD, our group obtained the CB1, CB2 and ML polygons described in the text.</figcaption></figure>

## 3.2 Comparison Between the "Park Fence" and the Park Polygons

A second method that was used to try to define the most suitable polygon was based on the existence of some fragments of a fence that would delimit portions

of the Park. These fragments of the "Park fence" are well known in the surrounding communities and are generally respected as being the local boundaries of the Park, both by the local community and by institutional representatives (Environmental Military Police, Scientific Police and municipal environmental management bodies).

On 13/Dec/2023, several points of the fence were visited and the planimetric coordinates of the base of the fence-posts were taken. In this study, precise planimetric coordinates were obtained using the Real-Time Positioning (RTK) network technique (Figure 4).

<figure><img src="image_106.png" /><figcaption>Figure 4: In December 2023, several remaining fragments of the Park's fence were visited. The precise planimetric coordinates of various points representing the geometry of the fence were recorded.</figcaption></figure>

This technique uses reference stations to generate and transmit corrections to the mobile receiver via the NTRIP protocol. For this survey, the Brazilian Network for Continuous Monitoring of GNSS systems in real time (RBMC-IP), the official IBGE network, was used. The reference station used was located in the campus of the Federal Institute of Santa Catarina (IFSC), in the city of Florianópolis, because it was the closest base. Transmitting the corrections in real time to the mobile receiver in the field allows for a significant reduction in positioning errors, ensuring that coordinates are obtained with high precision, to the nearest centimeter.

The limiting factor for using the networked RTK system is related to the transmission range of the radio waves. In our case, the distance between the base station (IFSC) and the mobile station was no more than 25 kilometers in a straight line, a fully acceptable distance for this positioning method and the topography of the area to be surveyed.

The 20 points sampled (excluding outliers) were imported into QGIS and compared with the polygons mentioned above (ML, CB1 and CB2 polygons). The points obtained do not coincide exactly with any of
the 3 polygons (Figure 5). However, the samples of the fence coordinates show a greater proximity to polygon CB1 (distance μ=0.19 meters σ=0.14 meters). The points fall at greater distances from the ML polygon (distance μ=0.62 meters σ=0.21 meters) and the CB2 polygon (distance μ=2.85 meters σ=1.26 meters).

<figure><img src="image_107.png" /><figcaption>Figure 5: The fragments of the "Park fence" are not exactly aligned with the polygons analyzed, but they are closer to the CB1 polygon. Map produced by the authors, using QGIS 3.30.1-'s-Hertogenbosch. Coordinates are stated in UTM22J SIRGAS2000 (EPSG:31982).</figcaption></figure>

# 4 DISCUSSION

IBGE Resolution No. 1/2005 established SIRGAS2000 as the new CRS for the Brazilian Geodetic System. A transition period of ten years was granted (expired in 2015), allowing SIRGAS2000 to be used concurrently with SAD69The Map Library file ML provides a polygon referenced using the EPSG:31982 UTM SIRGAS 2000 system (which is currently the official coordinate reference system in Brazil).

The 'Current Boundaries' file shows the shapefile referenced using SAD69 Classic Network (EPSG:29192), in non-compliance with IBGE Resolution 1/2005

The "Current Boundaries" file can provide two polygons, depending on the user's interpretation: (1) polygon CB1, assuming that the data is referenced to the SAD69/96 system (as suggested by the dates in the description in the State Law's text); (2) polygon CB2, assuming that the data is referenced to the SAD69 Classic Network system (as suggested by the EPSG:29192 code used in the available shapefile).

When using ProGriD platform to transform the Park's boundaries coordinates (informed in the State Law), the same three polygons are obtained,

depending on the choice between three possible transformations: SAD 69 Classic Network (provides CB2); SAD 69/96 (provides ML); and SAD 69/96 Doppler Technique or GPS (provides CB1).

The polygon in the map library (ML) and the two polygons from the 'Current Boundaries' file (CB1 and CB2) do not coincide, which could lead to problems in the characterization of the Park's border.

Furthermore, the points obtained as a sample for the "Park fence" do not coincide with any of the three polygons, although they are closer to the CB1 polygon.

Although the field GPS data suggest an advantage in using the CB1 coordinate reference system, the challenges in coordinate transformations prevented a clear conclusion regarding the most suitable system for defining the park's boundaries. The discrepancies between the various transformation methods—such as SAD69 Classic Network, SAD69/96, and the Doppler Technique—resulted in inconsistent outcomes. This lack of alignment between the transformations, combined with uncertainties in the input data, complicated the decision-making process. Despite the GPS data favoring CB1, the differences between the transformed polygons and the park's boundary points highlight the difficulty in selecting a definitive coordinate reference system. These issues, along with the lack of consensus in transformation systems, are discussed in greater detail in the alternatives analyzed below.

**Alternative A.** Assuming the data were from SAD69 Classic Network, the CB2 polygon would be obtained using both the QGIS reproject and ProGrid transformation. It is reasonable to assume that the input data is in the SAD69 Classic Network system, since this is stated in the EPSG of the shapefile (EPSG 29192), but the text of the State Law suggests otherwise. Furthermore, the CB2 polygon obtained is the one furthest away from the "Park fence".

**Alternative B.** If the data were from the SAD69/96 system, then the CB1 polygon would be obtained (using the QGIS reproject), but the ML polygon would be obtained (with the ProGriD transformation). It is reasonable to assume that the data is from this reference system, due to the dates declared in the State Law's text. However, this produces two different polygons, with the ML polygon being favored because it was produced by ProGriD, which is the official transformation platform. However, the points collected from the "Park fence" are closer to the CB1 polygon and not the ML polygon.

**Alternative C.** Assuming that the data were originated from the SAD69/96 Doppler Technique or GPS system, then the CB1 polygon would be
obtained (with the appropriate ProGriD transformation). It's reasonable to assume that the data is from this system, because in this case you get the polygon that is closest to the points collected from the "park fence" fragments.

By disregarding alternative A, two viable polygons remain: CB1 and ML. If CB1 were to be established as the legitimate polygon, there would be a number of implications that would make this option costly and not defensible. The first consequence would be the need to reconstruct the Park's polygon to replace the existing one in the Map Library.

When considering the ML polygon as correct, it can be seen that the only way to obtain this polygon from the points in the State Law's text, would be to use ProGriD with the transformation that considers the origin of the data in the SAD69/96 system. If this is the correct polygon, it is important to note that the points collected near the fragments of the Park's fence do not coincide with this polygon, but are a certain distance apart ($μ=0.62$ meters $σ=0.21$ meters).

# 5 CONCLUSION

Under the light of the discussions in this study, our group has produced the following recommendations as a contribution to public authorities and Santa Catarina society.

**Recommendation 1** – Recommendation to the State Government that the Managing Body expressly and unequivocally define which shapefiles (UTM 22S SIRGAS 2000 EPSG:31982) correctly define the Park's polygons.

**Recommendation 2** - Recommendation that users always use shapefiles reprojected to EPSG:31982, refraining from using on-the-fly transformations. Whichever network is selected for importing the points (SAD69 Classic Network or SAD69/96), the user must take appropriate care when "reprojecting" the points to the UTM 22 SIRGAS 2000 reference system (EPSG:31982), as choosing the wrong transformation can lead to positional errors of the order of a few meters.

**Recommendation 3** - Recommendation that users always obtain information or confirmation of the origin of the data when using SAD69 coordinates, so as not to confuse SAD69 Classic Network data, SAD69/96 data and SAD69/96 Doppler Technique or GPS data. When importing into QGIS, pay attention to the deprecated codes. The EPSGs 29192 (SAD69 Classic Network), 5858 (SAD69/96)) and 31982 (SIRGAS 2000) should be used.

**Recommendation 4** – The park's fence does not coincide with the Map Library's original polygon. However, the fence has been respected by the communities and has been used as a marker by the Environmental Military Police and the Scientific Police. The government should clarify the role played by this so-called park fence in the regions where this structure remains.

# REFERENCES

Cain, J. Coordinate Reference Systems (Best Practices for Assignment, Manipulation and Conversion in GIS Systems). ESRI Petroleum GIS Conference: 39. 2013.

IMA. Geoseuc Geographic Information System. Available at http://geoseuc.ima.sc.gov.br/#/, accessed in April/2023.

IMA. Santa Catarina State's Environmental Institute. Instituto do Meio Ambiente de Santa Catarina (IMA). Parque Estadual da Serra do Tabuleiro. Accessed in April/2024: https://www.ima.sc.gov.br/.

EPSG. EPSG Geodetic Parameter Dataset: v11.001. Available at https://epsg.org. Accessed on 01/04/2023.

IBGE. ProGriD - Transformation of Coordinates between Official Reference Systems. Brazilian Institute of Geography and Statistics. 2009. Available at www.ibge.gov.br/geociencias.

QGIS. QGIS Development Team. QGIS Geographic Information System. Accessed on 01/12/2023. Available at https://www.qgis.org.

SANTA CATARINA. State Law 14661/2009. Available at http://leis.alesc.sc.gov.br. Accessed on 01/08/2023.

SANTA CATARINA. State Decree 3446/2010, of 10 August 2010. Available at http://leis.alesc.sc.gov.br. Accessed on 01/08/2023.
# Analysis of Surface Thermal Behavior in Different Local Climate Zones (LCZ): A Case Study in Bragança (Portugal) (2013-2024)

Cátia Rodrigues de Almeida<fnref n="1" />,<fnref n="2" />,<fnref n="3" /> Ⓓª, Artur Gonçalves<fnref n="3" /> Ⓓᵇ and Ana Cláudia Teodoro<fnref n="1" />,<fnref n="2" /> Ⓓᶜ

**Keywords:** Urban Heat Island (UHI), Surface Urban Heat Island (SUHI), Remote Sensing (RS), Albedo, Land Use and Land Cover (LULC), Land Surface Temperature (LST).

**Abstract:** The Urban Heat Island (UHI) effect occurs when temperatures in urban areas are higher than surrounding vegetated areas, especially during the sunset and sunrise. UHI impacts include effects on public health and well-being, changes to the local microclimate, and influence on the local biome. This study evaluates the Land Surface Temperature (LST) and the corresponding Surface Urban Heat Island Intensity (SUHIint) across different Local Climate Zones (LCZs) in Bragança (Portugal) from 2013 to 2024, using images from Landsat 8 and 9 data collected with a portable thermal camera on different surfaces to assess thermal behavior across different scales. The results confirm the existence of the UHI effect in Bragança, where vegetated areas exhibit milder temperatures compared to built areas, especially in summer afternoons. Satellite-derived LST data indicate that the lowest temperature was recorded in an LCZ with vegetation, reaching (-7°C), while the highest minimum temperature was observed in an LCZ with higher density of anthropogenic elements (-3°C). Thermal camera measurements showed that surfaces such as asphalt and exposed soil reached 80 °C in the morning and remained above 60°C in the afternoon. These findings underscore the importance of considering mitigation measures, such as increasing vegetation in urbanized areas or replacing impervious surfaces.

## 1 INTRODUCTION

The viability of human life in each location requires a series of adaptations in Land Use and Land Cover (LULC), presenting a contemporary challenge. With the global population expected to increase from 7.7 billion to 9.7 billion by 2050 (ONU, 2019), it is essential to pay attention to the surface materials incorporated in the construction of houses, railways, commerce, industry, etc. (Oke, 1987; Weng, 2003). Urbanization can contribute to the formation of the Urban Heat Island (UHI) effect, characterized by higher temperatures in urbanized areas compared to surrounding rural areas (Imhoff et al., 2010; Oke, 1982). A specific approach to this phenomenon is the Surface Urban Heat Island (SUHI), which assesses temperature differences based on temporal variability

and surface characteristics (Voogt & Oke, 2003; Weng & Fu, 2014).

SUHI occurs mainly due to the characteristics of anthropogenic elements in the urban landscape, such as asphalt, sidewalks, and streets (Oke, 1988). These materials have a low albedo, meaning they reflect less solar radiation and retain more heat (Coakley, 2003). Figure 1 shows some materials found in an urban context, along with their respective albedo rates. As observed, most anthropogenic materials exhibit low reflectance.

Beyond albedo, the Local Climate Zones (LCZ) methodology was developed to classify urban and rural environments based on the local landscape and physical and thermal characteristics, such as land cover, terrain structure, building density, surface materials. LCZs are divided into various categories:

<footnote n="1">¹Department of Geosciences, Environment and Land Planning, Faculty of Sciences, University of Porto,
Rua Campo Alegre, 687, 4169-007 Porto, Portugal</footnote>

<footnote n="2">²Earth Sciences Institute (ICT), Pole of the FCUP, University of Porto, 4169-007 Porto, Portugal</footnote>

<footnote n="3">³CIMO, LA SusTEC, Instituto Politécnico de Bragança, Campus de Santa Apolónia, 5300 253 Bragança, Portugal
up201600831@fc.up.pt, ajg@ipb.pt, amteodor@fc.up.pt</footnote>

<footnote>ªhttps://orcid.org/0000-0003-4455-9407</footnote>

<footnote>ᵇhttps://orcid.org/0000-0002-4825-6692</footnote>

<footnote>ᶜhttps://orcid.org/0000-0002-8043-6431</footnote>
on one hand, built zones (urban), including blocks of buildings and residential areas with little or no vegetation, and on the other, natural zones (rural), such as forests, parks, water surfaces, etc. (Errea et al., 2023; Stewart & Oke, 2012).

<figure><img src="image_108.png" /><figcaption>Figure 1: The albedo of urban materials (adapted from (EPA - United States Environmental Protection Agency, n.d.)).</figcaption></figure>

Remote Sensing (RS) is a methodology applied to understand UHI using data collected by thermal sensors placed on various platforms. The sensors record the electromagnetic radiation emitted or reflected by surfaces without physical contact with the sampled surface, enabling the calculation of Land Surface Temperature (LST) (Pachêco, 2001; Weng et al., 2004). The satellite mission most commonly used for LST calculation, with freely available data distribution, is Landsat, which has included a thermal sensor since Landsat 4 (operational from July 1982 to June 2001) to Landsat 9 (launched in September 2021) (Almeida et al., 2021; US Geological Survey, n.d.-a; Zhou et al., 2019). In addition to satellite data, portable thermal cameras with in situ data collection can complement the methodology, allowing for a more detailed analysis of surface reflectance.

The objective of this study is to evaluate the UHI effect in Bragança (Portugal) using LST data from Landsat 8 and 9 between 2013 and 2023, along with LST data from a portable thermal camera collected during five field campaigns between 2023 and 2024, divided into summer and winter, to assess the seasonal influence on the results.

The choice of Bragança is based on: i) its mountainous setting, which naturally influences wind circulation and speed, potentially contributing to UHI formation; ii) centralized urban densification with surrounding vegetated areas; iii) its classification into seven distinct LCZs, facilitating the understanding of its LULC; and iv) as a supplement to the existing literature, as Bragança's UHI has already been analyzed using RS images and compared with Air Temperature (Ta) data from a monitoring network

with 23 weather stations distributed across seven LCZs between June 29, 2013, and July 5, 2021. In summer, the correlations were strong and very strong, while in winter, they were very strong in all sensors, confirming similar thermal behavior in both parameters (Almeida et al., 2022, 2023; Menezes, 2017).

Our hypothesis is that the integration of RS data acquired by satellites and thermal cameras can confirm the identification of the UHI effect in Bragança and inspire other researchers to implement similar approaches, contributing to the advancement of knowledge and promoting discussions on the need for sustainable urban planning, including the adoption of high albedo materials and green infrastructure to mitigate the UHI effect (Nations, n.d.). The LST measured by satellites and thermal cameras differs in specificity, with expected errors due to spatial resolution and ground proximity. Landsat 8 and 9 sensors acquire LST at 100m, resampled to 30m, limiting the detection of small transition areas. Additionally, the altitude difference reduces atmospheric influences on ground-level data, as thermal cameras capture at 1.5m, while satellites operate from 705 km.

# 2 METHODOLOGY

## 2.1 Study Area

Located in the northeastern of Portugal (41°48'20" N, 6°45'42" W), the municipality of Bragança has an area of approximately 1,174 km² and, according to the last census conducted in 2021, has a population of 35,341 inhabitants (Bragança, n.d.; INE, n.d.) (Figure 2). It is a region with complex terrain, significant variability in elevation, and heterogeneous LULC, featuring urban covers, agriculture, forested areas, exposed soils, and natural and vegetated areas. Its economy primarily focuses on public sector service companies (Gonçalves et al., 2018).

Bragança is characterized by a Mediterranean climate, with dry summers and mild temperatures (Csb), as defined in the Köppen-Geiger scale (Barceló & Nunes, 2009). Two scales influence the climate: i) regional, with the atmospheric circulation present in transitional zones (between the ocean and continents); and ii) a more detailed scale, influenced by terrain and proximity to water bodies (lakes or rivers). There are two main rivers in and around Bragança: Fervença (which crosses the region) and Sabor, which is located in a peripheral area and has
less local influence (Bragança, n.d.; Gonçalves et al., 2018).

One of the environmental actions implemented in Bragança is the extensive and continuous collection of temperature and relative humidity data through TGP-4500 model sensors, TinyTag, Gemini Data Loggers, Chichester, UK, installed three meters above the ground, enclosed in a protective white box, semi-leveled on light poles, and oriented towards the South.

<figure><img src="image_109.png" /><figcaption>Figure 2. Identification of the 23 Sensors used in this study in Bragança, Portugal: Compact low-rise (CLR) (gray); Compact midrise (CM) (red); Urban green spaces (GAB) (light green); Large low-rise (LLR) (orange); Open midrise (OM) (blue); Rural areas (RCD) (green); Sparsely built (SB) (purple); (adapted from (Almeida et al., 2022)).</figcaption></figure>

The installation criterion for the sensors considered the different LCZs identified in Bragança, ensuring at least three sensors in each class, namely: i) Compact low-rise (CLR): historical city center, including buildings with high density, medium-low height, and presence of embedded bricks/stones (sensors 4 and 6); ii) Compact midrise (CM): modern high-density buildings, with paved surfaces and medium-high height (sensors 3, 7, and 13); iii) Urban green spaces (GAB): predominantly green, covered by trees and low vegetation (sensors 2, 8, 9, and 11); iv) Large low-rise (LLR): industrial and commercial area, with paved parking lots, low or medium density, featuring low and high buildings (sensors 5, 17, and 21); v) Open midrise (OM): streets with isolated or low houses, of medium density (sensors 10, 12, 18, and 22); vi) Sparsely built (SB): transition between rural and urban environments, with scattered houses in forested and agricultural areas (sensors 1, 14, and 15); and vii) Rural areas (RCD): isolated rural spaces surrounding the city (sensors 16, 19, 20, and 23) (Almeida et al., 2022; Gonçalves et al., 2018). The data from these sensors enabled studies on UHI in the

region and surrounding areas, including methodologies with Remote Sensing (Almeida et al., 2022, 2023; Gonçalves et al., 2014, 2018).

## 2.2 Data Used and Processing Conducted

### 2.2.1 Satellite Images

We used Google Earth Engine (GEE) for satellite image processing, a cloud-based platform that processes data remotely, reducing the need for powerful computers and significantly minimizing processing time, which already offers a series of processed products (Hurni et al., 2017; NASA, n.d.-a). We utilized the "LANDSAT/LC08/C02/T1_L2" and "LANDSAT/LC09/C02/T1_L2" collections, with overpass times around 11 AM in Bragança. These satellites operate in a nearly polar and synchronous orbit at a nominal altitude of 705 km, with an orbital inclination of 98.2°, and revisit the same location every eight days, as they operate in alternating orbits, improving temporal resolution for the study areas (US Geological Survey, n.d.-a).

The thermal sensors record data at 100m, which is resampled and made available to users at 30m. In GEE, the collections are provided with atmospherically corrected surface reflectance data, and the LST is derived from the Land Surface Reflectance Code (LaSRC), generated through a single-channel algorithm developed in collaboration between the Rochester Institute of Technology (RIT) and NASA's Jet Propulsion Laboratory (JPL) (NASA, n.d.-a).

We used the code for processing LST (NASA, n.d.-b) and to filter the data for our study area between 2013 and 2024, considering only those with cloud coverage below 20%. Subsequently, we used the quality bands from each collection (QA Pixels) to remove pixels with clouds and cloud shadows.

We utilized the thermal band (B10) from both collections to obtain the LST. The spectral range is from 10.6 to 11.19 µm, and the data are recorded at 100m and resampled to 30m. To convert the LST to Celsius, we applied Equation (1) (NASA, n.d.-a; US Geological Survey, n.d.-b).

$$LST = TIR \times f_s + O - K \quad (1)$$

Where LST is the Land Surface Temperature (°C), TIR is the radiance in the thermal infrared band, f_s is the scale factor of Landsat 8 and 9 (0.00341802), O is the offset of Landsat 8 and 9 (149.0) and K is the conversion constant from K to °C (273.15).
To evaluate the LST by LCZs, we used a .shp file of points in GEE corresponding to the locations of the 23 sensors installed in Bragança to extract the LST values for each of them, and we added a function to export the data to a .csv file. Between 2013 and 2014, we identified 158 images, with the first date being 06/29/2013 and the last being 09/15/2024. Subsequently, we selected results where images provided data for all 23 sensors, resulting in 123 images spanning all seasons. We further filtered the data for summer and winter, obtaining 46 images for summer and 15 for winter, which were used in the processing. The lower number of winter images can be attributed to meteorological conditions, with higher cloud cover and precipitation during winter, aligning with the Köppen-Geiger classification (Barceló & Nunes, 2009).

To process and analyze the results, we calculated the SUHI by applying Equation (2), a methodology already used in other UHI studies (Almeida et al., 2022; Menezes, 2017).

$$SUHI_{\text{Int}} = LST_{LCZ} - \overline{LST_{RCD}} \quad (2)$$

Where SUHI<sub>Int</sub> is the intensity value of the UHI at each point; LST<sub>LCZ</sub> is the LST value extracted from the sensor points; and LST<sub>RCD</sub> refers to the average of the LST values extracted the RCD sensors (16, 19, 20, and 23).

To present the results, we used SPSS software to create boxplots of the SUHI<sub>int</sub> results, separated by LCZ and by seasonality (summer and winter), a methodology already applied in other studies (Almeida et al., 2022; Saher et al., 2021).

### 2.2.2 LST from Portable Thermal Camera

For in situ LST imaging, we used the HT Instruments THT33 portable thermal camera, which has a resolution of 80x80 pixels, a Field of View (FOV) of 21°x21°, a spectral range of 8 to 14 µm, and an Instantaneous Field of View (IFOV) (@1m) of 4.53 mrad. Images were taken at a standard height of 1.5m between the camera and the target to ensure consistent operator-associated error across all samples.

The collections were conducted only on clear, cloud-free days, as these conditions affect measurement accuracy. Consequently, three data collection sessions were conducted: two in winter (02/17/2023 and 02/20/2024) and one in summer (07/27/2023). Since morning data collections had to coincide with satellite overpass times, covering all 23 points with T<sub>a</sub> sensors was not feasible. Instead, ten representative points were selected to capture spatial

diversity and ensure at least one representative for each of the seven LCZs. The selected sensors, their respective LCZs, and the number assigned to each point for identification in the in situ collection are: i) CM: sensor 3, point P08; ii) CLR: sensor 6, point P07; iii) GAB: sensor 9, point P05; iv) SB: sensors 14 and 15, points P09 and P04; v) LLR: sensors 5 and 17, points P06 and P03; vi) OM: sensors 18 and 22, points P02 and P10; and vii) RCD: sensor 23, Point P01.

We assessed the diversity in LULC composition for each selected point and sampled different surfaces within a 10m radius of each sensor. We conducted two measurements at each point: i) between 10 am and 12 pm (within one hour before and after the satellite's overpass time in Bragança); and ii) between 1 pm and 3 pm to analyze temperature variability of the sampled surfaces compared to the first session. For each sampled surface, we took a photograph using a Canon EOS 800D camera to classify it afterward into the following categories: Vegetation; Exposed soil; Road with stones; Sidewalk; Asphalt; Sidewalk with shade; Asphalt with shade; Cobblestone; Iron; Curb and Vegetation with shade. Each surface is assigned a letter to differentiate its LULC for the analysis. For example, five different surfaces were collected at P01, labeled from P01A to P01E.

The HTMercury33 application was used to record, visualize, and export the LST results in .pdf format (User Manual - HTMercury33, n.d.) We created an R routine to extract specific data from the .pdf files to Excel, specifically: S (temperature associated with the fixed central cursor); H (temperature of the hottest point in the image); C (temperature of the coldest point in the image); and Mean (average temperature). Subsequently, we categorized the data by season (summer and winter) and generated boxplots with the maximum and minimum LST.

# 3 RESULTS AND DISCUSSION

## 3.1 SUHI<sub>int</sub>

According to Figure 3, it can be observed that the summer data exhibit greater amplitude and more outliers compared to winter data, which may be associated with a higher availability of solar radiation during this season.

In the summer, most sensors showed negative medians of SUHI<sub>int</sub>, except for sensors 3 (CM), 17 (LLR), and 12 (OM). In these three cases, the surroundings allowed the entry of short-wave
radiation in the early hours of the day from the east (sunset), warming the surface in the early morning. Additionally, the structures of the OM (12) and CM (3) sensors, which contain medium-height buildings, may interfere with heat exchange with the environment due to shadow projection, making the area less prone to heat accumulation (Zheng *et al.*, 2019). In the case of LLR (17), due to the presence of heterogeneous profiles, including industrial and commercial areas and paved parking lots, heat can be absorbed, especially by elements with low albedo (pavements, structures made of cement, etc.) (Zheng *et al.*, 2019). Still in summer, the three sensors that exhibited the lowest medians consist of vegetation elements, as in the case of sensors 9 and 11 (GAB) and 15 (SB), reinforcing that the presence of vegetation aids in regulating SUHI\_int (Dai *et al.*, 2019).

<figure><img src="image_110.png" /><figcaption>Figure 3. Boxplot illustrating the Surface Urban Heat Island Intensity (SUHI\_int) across various Local Climate Zones (LCZ). The zones are categorized as follows: i) Compact Low-Rise (CLR) in gray; ii) Compact Midrise (CM) in red; iii) Urban Green Spaces (GAB) in green; iv) Large Low-Rise (LLR) in orange; v) Open Midrise (OM) in light blue; and vi) Sparsely Built (SB) in purple.</figcaption></figure>

In summer, regarding higher temperatures and disregarding outliers, sensors 17 (LLR), 12, and 18 (OM) exhibited the highest values, reinforcing that the presence of anthropogenic elements may contribute to heat accumulation on the analyzed surfaces (Oke, 1987, 1988; Voogt & Oke, 2003). Finally, the lowest temperature is associated with sensor 9 (GAB), reaffirming the importance of preserving and prioritizing vegetation in spaces to regulate heat (Miles & Esau, 2020).

In winter, the medians were mainly below 0, except for sensors 13 (CM), 2 (GAB), 17 (LLR), and 12 (OM). Most of these sensors are integrated into LCZs with the presence of anthropogenic materials, converging with the values identified in summer, except for sensor 2, which contains vegetation. This

behavior can be explained by the season: in winter, vegetation enters a state of dormancy to protect itself from cooler temperatures, reducing its metabolic activities and, consequently, slowing down heat exchanges. This also limits photosynthesis, as it does in other seasons under similar conditions (Du *et al.*, 2016; French & Inamdar, 2010). There is also leaf shedding in deciduous trees, known as senescence, which may also contribute to this result (Mariën *et al.*, 2019).

As in summer, the highest temperatures are associated with sensors containing anthropogenic elements, as seen in sensors 17 (LLR), 4 and 6 (CLR), and 3 (CM). In contrast, the lowest temperatures occurred in sensors 9 and 11 (GAB), which contain vegetation (Tang *et al.*, 2017).

## 3.2 LST from Thermal Camera

The number of data points available in each measurement session was 220 images (110 in the morning and 110 in the afternoon) on 02/17/2023, 224 on 02/20/2024 (112 in each period), and 186 on 07/27/2023 (93 per period). The differences in sample size are attributed to occasional measurement errors related to the equipment. Figure 4 shows the recorded LST values in summer across different surfaces, with surface identification highlighted in the legend.

In general, summer temperatures are much higher than winter temperatures due to increased electromagnetic energy availability and solar angle, which minimizes shadow projection. Some surfaces reached nearly 80 °C in the morning and above 60 °C in the afternoon, mainly asphalt and exposed soil. In summer, vegetation also retains greater heat in the morning and has cooler temperatures in the afternoon. This can be explained by vegetation's prolonged exposure to electromagnetic energy in the early hours compared to built-up areas, which are shaded by structures and accumulate less heat. In the afternoon, the opposite is observed: vegetation temperatures tend to be cooler due to latent heat exchange with the environment, especially through evapotranspiration, which reduces the heating of these surfaces, in addition to the shading effect provided by tree canopies. In urban areas with low albedo surfaces, such as asphalt and concrete, the radiative balance favors the absorption of short-wavelength radiation, resulting in the surface heating of these areas (Oke, 1988; Yu *et al.*, 2019).

In a more detailed LCZ analysis, at P07 (CLR), we observe that the LULC materials such as asphalt and elements of iron, along with the lack of vegetation, contribute to higher temperatures in summer,
<figure><img src="image_111.png" /><figcaption>Figure 4. Boxplot of Land Surface Temperature (LST) at the ten in situ measurement points using a thermal camera.</figcaption></figure>

particularly in the afternoon, due to the thermal inertia of these materials, which accumulate heat throughout the day. In winter, temperatures are slightly higher than most vegetation classes, especially in the afternoon (Cai *et al.*, 2017). At P08 (CM), asphalt and sidewalk temperatures were higher in summer than in winter, both in the morning and afternoon, which may be associated with the angle of electromagnetic energy. In summer, this energy is more direct, reducing shading from existing buildings in the area (Yu *et al.*, 2019).

At P05 (GAB), it is evident that vegetation created a cooler microclimate, with lower temperatures than densely built areas, especially in summer at both times of the day. The LST values for these locations are slightly lower in winter than in built-up areas, although they represent a less pronounced cooling. This could be associated with the lower radiation levels and the vegetation's dormancy process, which reduces evapotranspiration and impacts natural cooling, as discussed in the UHI\_int results (Rose, 2019).

At points P03 and P06 (LLR), both in summer and winter, the LST showed higher values in paved areas, similar to CM and CLR points, which may be due to the lack of shading in the area. Regarding points P02 and P10 (OM), LST values in summer were more moderate, with warmer asphalt areas and cooler vegetated areas, particularly in the afternoon, potentially helping to regulate local temperature. In winter, due to lower construction density and vegetation presence, temperatures were milder at these points (Oke *et al.*, 2017).

For points P04 and P09 (SB), in both summer and winter, vegetation showed higher values in the morning, possibly due to heat accumulation in the early hours of the day, a similar effect to that observed at point P01 (RCD). The absence of large vertical barriers in open areas facilitates heat exchange with the environment in the night, promoting heat dissipation through convection. However, surface characteristics and solar radiation absorption still significantly influence daytime temperature retention in these areas. For this reason, it is suggested that UHI studies should be conducted after sunset, incorporating additional methodologies such as the use of Unmanned Aerial Vehicle (UAV), which provide greater flexibility regarding the timing of data collection.

Despite scale differences, the results were similar, with higher surface temperatures in built-up areas, reinforcing the need for mitigation measures. Incorporating vegetated surfaces in densely urbanized areas, such as green roofs and substituting impermeable surfaces with alternatives like permeable asphalt, is essential (Cai *et al.*, 2017; Gonçalves *et al.*, 2014; Mullerova & Williams, 2019). It is worth noting that, especially in summer, additional data through field visits are necessary to verify if these initial thermal behavior results are consistent across time.

# 4 CONCLUSIONS

This study used RS data to assess surface thermal behavior across different LCZs and identify the UHI
effect in Bragança. We utilized both satellite data and *in situ* data collected from a thermal camera, comparing the results from both techniques. In terms of spatial resolution and measurement height, the two methods differ: while Landsat is efficient for mapping larger, more homogeneous areas, the use of the thermal camera provides more detailed, localized surface data, being less influenced by atmospheric parameters and more effective at mapping heterogeneous areas. This makes the integration of these techniques particularly valuable for UHI studies.

The SUHI\_int, calculated from satellite data, was higher in summer, especially in LCZs with anthropogenic surfaces, compared to those with vegetation, highlighting the influence of seasonality. *In situ* measurements corroborated these results, showing that impermeable surfaces retained more heat, particularly during summer afternoons.

The analyses suggest that integrating green areas and using high-albedo materials are strategies to mitigate the UHI effect and promote climate resilient cities. Future research could explore long-term UHI dynamics, with additional field data collection and the development of specific mitigation strategies tailored to the different LCZs in Bragança, whose findings could be applied to other urban contexts.

## ACKNOWLEDGEMENTS

This work was supported by national funds through FCT/MCTES (PIDDAC): CIMO, UIDB/00690/2020 (DOI: 10.54499/UIDB/00690/2020) and UIDP/00690/2020 (DOI: 10.54499/UIDP/00690/2020); and SusTEC, LA/P/0007/2020 (DOI: 10.54499/LA/P/0007/2020). The authors are grateful to the Foundation for Science and Technology, I.P., projects UIDB/04683/2020 (https://doi.org/10.54499/UIDB/04683/2020), IDP/04683/2020 (https://doi.org/10.54499/UIDP/04683/2020). Cátia Rodrigues de Almeida was financially supported by Portuguese national funds through FCT (Grant: PRT/BD/153518/2021).

## REFERENCES

Almeida, C. R. de, Furst, L., Gonçalves, A., & Teodoro, A. C. (2022). Remote Sensing Image-Based Analysis of the Urban Heat Island Effect in Bragança, Portugal. *Environments* - MDPI, 9(8). https://doi.org/10.3390/environments9080098

Almeida, C. R. de, Garcia, N., Campos, J. C., Alírio, J., Arenas-Castro, S., Gonçalves, A., Sillero, N., &

Teodoro, A. C. (2023). Time-series analyses of land surface temperature changes with Google Earth Engine in a mountainous region. *Heliyon*, 9(8), e18846. https://doi.org/10.1016/j.heliyon.2023.e18846

Almeida, C. R. de, Teodoro, A. C., & Gonçalves, A. (2021). Study of the Urban Heat Island (UHI) Using Remote Sensing Data/Techniques: A Systematic Review. *Environments*, 8(10), 105. https://doi.org/10.3390/environments8100105

Barceló, A. M., & Nunes, L. F. (2009). Iberian Climate Atlas 1971-2000. In *Atlas Clim{á}tico Ib{é}rico - Iberian Climate Atlas*.

Bragança, M. de. (n.d.). *Estratégia Municipal de Adaptação às Alterações Climáticas*. Retrieved October 1, 2024, from https://www.cm-braganca.pt/cmbraganca2020/uploads/writer_file/document/6050/B_8-\_Publica\_o\_do\_Plano\_Municipal\_de\_Ambiente.pdf

Cai, G., Liu, Y., & Du, M. (2017). Impact of the 2008 Olympic Games on urban thermal environment in Beijing, China from satellite images. *Sustainable Cities and Society*, 32, 212–225. https://doi.org/10.1016/j.scs.2017.03.020

Coakley, J. A. (2003). Reflectance and Albedo, Surface. *Encyclopedia of Atmospheric Sciences*, 1914–1923. https://doi.org/10.1016/b0-12-227090-8/00069-5

Dai, Z., Guldmann, J. M., & Hu, Y. (2019). Thermal impacts of greenery, water, and impervious structures in Beijing's Olympic area: A spatial regression approach. *Ecological Indicators*, 97(July 2018), 77–88. https://doi.org/10.1016/j.ecolind.2018.09.041

Du, H., Song, X., Jiang, H., Kan, Z., Wang, Z., & Cai, Y. (2016). Research on the cooling island effects of water body: A case study of Shanghai, China. *Ecological Indicators*, 67, 31–38. https://doi.org/10.1016/j.ecolind.2016.02.040

EPA - United States Environmental Protection Agency. (n.d.). *Ilhas de calor*. Retrieved October 13, 2024, from http://www.epa.gov/heat-islands

Errea, C. L., Almeida, C. R. de, Gonçalves, A., & Teodoro, A. C. (2023). Remote Sensing Analysis of the Surface Urban Heat Island Effect in Vitoria-Gasteiz, 1985 to 2021. *Remote Sensing*, 15(12), 3110. https://doi.org/10.3390/rs15123110

French, A. N., & Inamdar, A. (2010). Land cover characterization for hydrological modelling using thermal infrared emissivities. *International Journal of Remote Sensing*, 31(14), 3867–3883. https://doi.org/10.1080/01431161.2010.483491

Gonçalves, A., Ornellas, G., Ribeiro, A. C., Maia, F., Rocha, A., & Feliciano, M. (2018). Urban cold and Heat Island in the City of Bragança (Portugal). *Climate*, 6(3), 1–14. https://doi.org/10.3390/cli6030070

Gonçalves, A., Ribeiro, A. C., Feliciano, M., & Maia, F. (2014). *Clima Urbano de la Ciudad de Bragança*.

Hurni, K., Heinimann, A., & Würsch, L. (2017). *Google Earth Engine Image Pre-processing Tool: Background and Methods*. 12. https://www.cde.unibe.ch/e65013/e542846/e707304/e707386/e707390/CDE_Pre-processingTool-UserGuide_eng.pdf