Identification and analysis of long-term changes in river deltas and riparian zones using time-series multispectral data

Today’s remote sensing data and technologies offer the capability to effectively monitor diverse and challenging environments around the world, such as coastal river and riparian zones. Coastal riparian zones and river deltas usually demonstrate extreme coastline changes in terms of the extent of water coverage of inland territories due to flood events, low and high tides, the climate, specific environmental characteristics, etc. In this paper, we exploit freely available multispectral time series data for previous decades, utilizing Landsat missions in order to develop an open-source-based image processing pipeline for the extraction of the actual yearly average coastline status of riparian river delta areas. The latter present significant temporal coastline changes between years, semesters, and months. Average mean maps are generated and then compared to several temporal levels in order to distinguish long-term significant changes and ecosystem threats. Additionally, a custom long short-term memory (LSTM) neural network is deployed to forecast the evolution of the coastline by exploiting the average value for each pixel across all available images as a training sample and producing a forecast output for the next period. The network achieves accuracy scores of 89.77% over 'non-water' depicting pixels and 84.26% over 'water' depicting pixels, for regions that present frequent changes between land and water coverage over time. The predicted map presents high statistical agreement with the respective average map generated in two different validation approaches, with kappa coefficients of 85.9% and 91.4%, respectively.


Introduction
Coastlines are the interfaces between land, ocean and inland waters, and are strongly influenced by human activity (Winarso et al. 2001). Additionally, coastal and riparian zones host a rich variety of biodiversity and ecosystem types and are therefore valuable from environmental, social and economic (Feyen et al. 2020) aspects. The European Commission (Boteler 2014), the United Nations (UN), and many other political and environmental organizations have issued a series of guidelines and practices regarding the protection of these areas, such as the European Green Deal, the Directive for Maritime Spatial Planning, the UN's Sustainable Development Goals (SDGs) 13 and 14, as well as the newly released Coastal Zones dataset from the Copernicus Land Monitoring Service. Specifically, the UN's SDG 14 highlights the importance of accurate shoreline information for coastal conservation.

3
Since a large and growing portion of the world population live in flat coastal areas which are environmentally heavily stressed, the assessment of spatiotemporal coastline changes is of major importance. To this end, numerous relevant studies have aimed to address different environmental risks in coastal zones. Chatziparaskeva et al. (2022) efficiently analyzed the concentrations of microplastics on marine and coastal waters in the Mediterranean region, whereas many other studies (Tsai et al. 2022;Wu et al. 2022) have monitored the coastline changes of estuary environments by using multispectral time-series datasets. More specifically, coastline monitoring assessment studies have been carried out through different approaches (remote sensing, photogrammetry, GIS) across the wider Mediterranean region (Travers et al. 2010;Nicholls and Hoozemans 1996) as well as in specific countries and areas, e.g., within Spain (Llorens et al. 2018), Greece (Skilodimou et al. 2021), Italy (Fabris 2019) and Cyprus (Andreou et al. 2017).
In parallel, the plethora of earth observation (EO) assets-data and technological tools-provide a unique opportunity to establish novel methodological frameworks for the development of efficient and sustainable monitoring solutions for environmentally sensitive regions. In many cases, even data with medium spatial resolution, such as data derived from Copernicus, the ESA, and Landsat, are deemed adequate (Bresciani et al. 2021). Through the years, many different solutions have been proposed for mapping flooding, water extraction and severe coastline changes due to soil erosion, a summary of which was documented by Anastopoulos and Nikolakopoulos (2021). Such solutions usually make use of the unique spectral characteristics of water surfaces through specific indices such as mNDWI and NDVI, or even use classification approaches (Makantasis et al. 2018;Sivanpillai et al. 2010;Sheng et al. 2008;Huang et al. 2014). Furthermore, Abdelhady et al. (2022) successfully deployed an automated shoreline extraction pipeline from VHR multispectral imagery, in which lidar-derived data were utilized for validation purposes. The technological advances of the past decades, such as the development of AI algorithms for satellite EO data, the deployment of high-end computing architectures, effective data-handling procedures (Munawar et al. 2022), as well as the introduction of additional non-satellite approaches such as UAV monitoring and lidar captures (Rovithis et al. 2017), have further expanded the potential use cases to coastline environments. The continuous progress made in new processing techniques and technological approaches alongside the availability of more open-source solutions is allowing the parallel monitoring of several case studies or even climate phenomena all round the world (Donchyts et al. 2016;Ruheili et al. 2021). Supervised and unsupervised classifications (Taveneau et al. 2021), automatic change detection techniques, and machine and deep learning (DL) applications using sophisticated artificial neural network (Makantasis et al. 2018) architectures, among other methods, are constantly being improved using more sophisticated mathematical models or fusion approaches (Zhu et al. 2021). To this end, algorithmic solutions specifically focusing on the application of DL approaches to flood/water extraction case studies are being tested and fine-tuned (Sadiq et al. 2022).
All of the above emphasize the significant risks to which coastal zones are exposed as well as the need for sustainable management of these areas. The importance of defining a methodological framework that utilizes modern tools and technologies and can identify significant changes and risks to which coastal zones are subject in real or near-real time is highlighted. Such a methodology would also assist the competent authorities and stakeholders in making coastal areas more adaptable and resilient (Bertoldi 2018).
In this study, a simplified open-source methodology is proposed. This methodology utilizes multispectral Landsat datasets and (i) periodically determines the tide-independent status of a river delta coastline and (ii) forecasts the future coastline based on previous observations and trends. This solution, which was developed via Python open-source code, exploits widely used multispectral data and indices and offers potential flexible integration with external operational products and platforms. The present study does not aspire to provide benchmark classification results in terms of overall score and algorithmic fine-tuning. On the other hand, it does provide a tangible methodological and analytical processing pipeline that is easily comprehensible and can be adapted to similar fragile environmental ecosystems.

Case study
The study area for the EPIPELAGIC project is depicted in Fig. 1 and corresponds to an area of 762 km 2 . It concerns the coastal area of the Thermaikos Gulf and includes part of the city of Thessaloniki and the Axios Delta National Park (Axios National Park-Loudia-Aliakmon). This park is of high ecological importance, and it is protected by the Ramsar Convention as a Wetland of International Importance and is in the network of Natura 2000 sites. Considering its environmental importance, it has been also highlighted in several other scientific publications (e.g., Karageorgis et al. 2001;Nikolaidis et al. 2009).
However, in this work, we first evaluated our approach using a smaller area of interest (AOI) that presents similar characteristics and faces the same climate change hazards, due to the large extent of the EPIPELAGIC AOI. In this paper, all processes refer to the wider region of the Evros River delta in NE Greece, (Fig. 2), which is considered a wetland area of unique natural beauty and of high environmental importance. This AOI was selected based on its characteristics-a coastal area near a river delta combining a variety of land uses (urban, industrial, commercial) and a variety of economic activities.
The Evros River delta AOI consists of swampy areas with low vegetation and intense seasonal changes as well as temporary floods based on meteorological phenomena. The lack of a clear coastline and the extent of the seasonal changes make the study of coastal changes particularly difficult; for this reason, it is considered an ideal case study for this multitemporal study.

Data acquisition
A continuous flow of satellite data and information is required to be able to monitor the coastal environment for the full period of the past 20 years. An additional target of this study is to make the proposed methodology as interoperable as possible in other AOIs by utilizing (i) a commonly used set of multispectral data and (ii) a series of homogeneous time-series data for the last 20 years. To this end, data from three different Landsat missions were selected: Landsat 5TM, 7 and 8.
Multispectral datasets from Landsat sensors have been evaluated in hundreds of studies over the years (Hemati et al. 2021), even for the wider Aegean Sea (Ekercin 2007).
Overall, for the data acquisition process, a selection criterion of < 5% cloud cover was set, and a quality control check was performed in order to discard no-data stripes and/or corrupted scenes. 159 Level 1 Collection 1 Landsat datasets were selected for further processing, as seen in Table 1.

Datacube creation and processes
The workflow of the developed methodology for the creation of mean maps is presented in Fig. 3. All implemented technical steps carried out and reported in this section  were executed through a series of Python-developed functions in a unified script. Consequently, Python libraries such as gdal, openCV, rasterio, sklearn, pytorch, and pandas were widely utilized in this study. As typically implemented in most pixel-based multispectral data handling processes, we did not treat each spectral band separately; we proceeded with the creation of spectral datacubes for each scene. The initial step prior to atmospheric correction was to crop all datasets with a basic Python function over the selected AOI. An unsupervised classification using a k-means algorithm was executed for each generated datacube to produce the classification maps.

Atmospheric correction
As mentioned above, the selected datasets consisted of Level 1 Collection 1 data, which represent top-of-atmosphere values but are considered ideal for any time-series analysis in general. An appropriate workflow was developed at the start of the datacube-generation Python script in order to apply the necessary atmospheric corrections (Fig. 4).
The coefficients given in the metadata for each scene (.MTL) were used for this correction. Specifically, this process was divided into two parts: (i) the conversion of the original values of the image into radiation values (DN to radiance) and (ii) the conversion of radiation values into reflectance values (radiance to reflectance) (US Geological Survey 1998).
For the first step, the "gain and bias" method described by Green Edmund et al. (2000) was implemented as follows: where L is the radiance value in W/(m 2 ·srad·μm), DN is the 8-bit digital number radiance value as captured from the sensor, with values ranging from 0 to 255, and G and B are the coefficients "channel gain" and "channel offset," available from each scene's metadata file.
The final surface reflectance value was then calculated through where ρ is the reflectance index, d 2 is the squared Earth-Sun distance in astronomical units, ESUN is the mean value of the out-of-atmosphere radiance, and SZ is the vertical scene angle.

Index selection
A crucial data-handling step was to enhance the spectral capacity through the implementation of datacube composites by adding selected indices on top of the existing spectral bands, thus acknowledging the fact that coastline identification is a challenging task. Some widely adopted selected indices for the extraction of water bodies are NDWI, MNDWI, and NDVI:

NDVI ∶= (NIR − RED)∕(NIR + RED)
NDWI ∶= (Green − NIR)∕(Green + NIR) In many cases, the correct use of such indices (McFeeters et al. 1996, Xu et al. 2006) alongside the effective configuration of the threshold value can produce significantly good results through the exploitation of the characteristics of the spectral reflectance of water. Specifically, in the study of Kyriou and Nikolakopoulos (2015), satisfactory results were achieved using the NDVI, NDWI and MNDWI indices in a flood mapping analysis of the Evros River case study using Landsat data.
Similarly, in this study, the spectral bands in the visible and the infrared parts of the spectrum were isolated from each scene (Table 2), while the aforementioned three indices/bands were calculated for the production of each datacube composite (Fig. 5).

Unsupervised classification
The k-means clustering algorithm is a well-established tool for performing unsupervised classifications of remotely sensed geospatial and satellite multispectral data, as it allows the automatic separation of a cloud of points into subgroups/ clusters (Tou et al. 1974;Lv et al. 2010). In this study, we chose to apply the k-means algorithm to the atmospherically corrected datacubes for the identification of temporal changes in the coastline and water bodies through the years.
Each scene was classified by the k-means algorithm into two different classes for their separation into water bodies and land. These classifications were compiled sequentially per year and then every 5 years, so the average of these classification results was calculated, and consequently the probability of matching each pixel to water or land was produced as well. It should be noted here that a key step in this processing chain was the creation of a cloud mask to remove atmospheric water vapor-using the BQA channel of the Landsat receivers-in order to remove incorrectly sorted pixels. Pixels corresponding to clouds or cloud shadows mNDWI ∶= (Green − Middleinfrared)∕(Green + Middleinfrared) ✓ Band 10 -Band 11 - were not removed from the results of the classification. The band quality assessment band (BQA band) was used to remove these pixels and create a mask that gives these pixels a value of 0. The binary result of the classification for each scene was a 1516 × 1108 raster GeoTIFF image in which the pixels corresponding to water had a value of 2, the pixels corresponding to land had a value of 1, and the pixels corresponding to clouds and cloud shadows had a value of 0 (Fig. 6). The exact selection of the correct threshold value for binary mask creation was carried out after considering the respective literature and other studies (Pan et al. 2020), and was finalized as illustrated in Fig. 6, below: The final product was a set of classification maps of the area: a map for each year (20 maps) and a map for every 5 years (4 maps) for the period 2000-2020. Additionally, maps were produced for the identification of important changes during every consecutive 5-year period. At the final stage of our processing pipeline, all output maps were georeferenced into the Hellenic Geodetic Reference System 1987 (HGRS87).

Forecasting
Complementary to the aforementioned processing pipeline of EO data time-series analysis, a module was implemented for predicting the status of the water bodies in our AOI at the next monitoring period. To this end, an artificial neural network (ANN)-specifically, a long short-term memory (LSTM) network (Hochreiter et al. 1997)-was used. In recent years, deep learning (DL) architectures, including LSTMs, have been extensively used in change detection and forecasting applications based on geospatial data (Papadomanolaki et al. 2021;Wang et al. 2019;Campos-Taberner et al. 2020).
The main goal was to train a model that can predict the future state (no-water/water) of each individual pixel based on the time series produced by the previous EO analysis for the years 2000-2019 (see "Datacube creation and processes"). To achieve this goal, the following steps were carried out: a) The mean value and variance of water occurrence observations for each pixel were estimated for the entire time span (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018)(2019) in order to classify the pixels into three classes: "almost always land," "variable," and "almost always water." More specifically, the variances of the available water occurrences for each pixel were computed and then classified as "variable" if their variance was higher than thres, where and minvar and maxvar are the minimum and maximum variance values, respectively. Nonvariable pixels were then directly assigned to the "almost always land" and "almost always water" classes, based on the average water occurrence for the specific pixel. b) All available water occurrence observations were grouped within fixed time intervals (every year between 2000 and 2020) and the average water occurrence value was estimated. Taking into account that pixels that are "mostly land" or "mostly water" may exhibit different behaviors, we further segmented the "variable" class into two subclasses, namely "mostly land" and "mostly water," based on the mean total water occurrence of each pixel (less or more than 50%, respectively). Consequently, the prediction of the future states of the pixels in each subclass was performed by using a simple twocell LSTM model for each subclass. More specifically, as presented in Fig. 7, the first LSTM cell had one input (the observed water occurrence for a specific year) and N = 64 hidden feature outputs, while the second LSTM cell had N = 64 inputs and N = 64 outputs. Finally, a third linear layer (64 × 1) was used to provide the final output, i.e., the observed water occurrence for the next year (monitoring period). Each model was trained using the set of yearly estimated water occurrences for the years 2000-2020, i.e., a time series with 21 measurements. For the Evros delta case study, the "mostly land" subgroup contained 10,166 pixels and the "mostly water" subclass contained 6098 pixels. The first 80% of the pixels in each time series were used for training, while the remaining 20% were reserved for validation, i.e., for predicting the value in 2020 based on the previous monitoring periods (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018)(2019). For both subclasses, thres = minvar + (maxvar − minvar) × 0.2 Fig. 7 Diagram of the deployed LSTM network the training and validation errors for predicting the last available water occurrence value converged rapidly (the mean square error for both validation sets was approx. 0.055 after 5 epochs). When the yearly water occurrences for the "variable" class were quantized to binary values using a threshold of 0.5, the corresponding validation accuracies after 10 epochs were 89.77% for "mostly land" pixels and 84.26% for "mostly water" pixels.

Results
In this section, all of the results derived from our study and the aforementioned methodology are presented, namely: (i) yearly average maps of the status of a coastal area during a calendar year, as extracted from a series of multitemporal EO data; (ii) change detection maps for different monitoring periods; and (iii) the prediction map for the AOI, based on previous observations and a pre-trained LSTM neural network.

Average maps
In the present study, a big volume of EO data is handled and several outputs can be extracted at several intermediate implementation stages. Based on the workflow presented in previous sections ("Data acquisition" and "Datacube creation and processes"), the main outputs produced from the time-series datacubes are the yearly average probability maps. All of the generated average maps are presented in Fig. 8.
Based on the above, and by examining the results through visual interpretation, it is already clear that changes appear over the long term close to the river basin buffer zone and the general coastline. Of course, many of those changes are strongly influenced by the overall river management by local authorities in both Greece and Turkey. Such an example is highlighted in Fig. 9 for the eastern riverbank, where the cultivation of rice fields can dramatically change the probability that this area is water covered.

Temporal changes
Accordingly, a set of 5-year composite maps were generated and are presented in Fig. 10.
The methodology applied is similar to the thresholding process described in the "Forecasting" section; essentially, the areas presenting high variance and extreme temporal changes are categorized into two separate classes: low probability of water/flood and high probability of water/flood. Again, many of the regions categorized as "high variance" in each 5-year period are either parcels in which flood-dependent cultivation is performed or changes along the coastline and the river basin that are not easily distinguishable and are implemented over the years.
Consequently, in this analysis, temporal change maps for this period were created, which revealed the more pressured areas with (i) receding water bodies or (ii) new water bodies. These change maps are demonstrated accordingly in Fig. 11. Finally, a change map for the reference years of 2000 and 2020 was also created; it revealed the gross changes that occurred through this period (Fig. 12).
It is clear that the following features are efficiently identified by this change map:

Forecasting result
The final output of this study concerns the forecasting image derived from the methodology described in the "Forecasting" section and the module implemented to predict the status of the water bodies at each monitoring period using an LSTM neural network. The training of the model considered all average map outputs for the period 2000-2019, while the result is a map showing predictions for the year 2020. Below, both the average map produced and the LSTM prediction are presented accordingly.

Forecasting validation
The aforementioned results (Fig. 13) are compared for the purpose of internally validating the LSTM prediction. These two instances seem identical (no major differences), which is expected since most of the AOI is permanently covered by either water or land. Hence, an extracted set of pixels was defined to make this validation process both more accurate and more challenging. Specifically, two different validation tests are hereby presented. The first is applied to the high-variance regions of the year 2020, which present frequent changes between land and water coverage over time. The second is a buffer zone 250 m around the coastline. Firstly, a raster mask of the selected pixels is created, upon which the validation/ comparison process is performed. The kappa coefficient and its respective confusion matrix between the yearly average map and the LSTM prediction for 2020 are calculated for both tests. The two selected validation masks are presented in Fig. 14.
For the high-variance regions, the two instances (classification and LSTM prediction) present good agreement, even though this selection of pixels covers some cultivation parcels and some areas that periodically flooded due to anthropogenic actions that cannot be easily predicted. This test achieved a kappa coefficient score of 85.9%, and its corresponding results (percentage and pixels) are presented in the confusion matrices of Table 3.
For the 250 m coastal buffer zone, the two instances (classification and LSTM prediction) present similar results, with an overall kappa coefficient score of 91.4%. Again, the corresponding results (percentage and pixels) of this test are presented in the confusion matrices of Table 4.

Conclusions
In this study, a novel end-to-end methodological framework for the effective monitoring of coastal and river riparian zones, such as river deltas, was proposed. The results indicate the workability of this end-to-end tool, even when applied to challenging case studies; however, it does not aspire to present benchmark data or algorithms in terms of accuracy/classification scores.
This framework is easily adaptable and can be integrated into existing solutions, as it exploits open access data and methodologies. Multispectral time-series data for the past  1 3 20 years were analyzed using simplified Python automated and semi-automated image processing workflows in order to assess a series of change detection maps of the coastline of Evros River in Greece. Additionally, an LSTM network was deployed to perform forecast analysis based on the extracted classification maps. The comparison indicated high statistical agreement, with kappa coefficient scores of 85.9% and 91.4% for the two different validation procedures, respectively. Despite the fact that the workflow provides promising results as a proof of concept, it still has limitations in terms of (i) EO data usage, (ii) the absence of in-situ sensor measurements for a more thorough and accurate validation procedure, and (iii) correlation with observed and documented environmental hazards in the region.
Future work for the authors following this study includes (i) the application of this methodology to wider and more challenging case studies, such as the EPIPELAGIC project case study, (ii) the integration of additional multitemporal information derived from in-situ sensors and measurements for verification/validation purposes, (iii) the use of additional EO multispectral data such as Copernicus Sentinel 2 and the recently launched Landsat 9, and (iv) the investigation of different network architectures and the use of ground truth information as an additional validation procedure.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.