1. Introduction
Coastal areas are defined as those where an interface or a transition between land and sea exists. These areas constitute important and sensitive environments with significant complexity, which are, however, difficult to define due to their diversity in form and dynamics as well as to their spatial boundaries [1,2].
Although coastal areas account for 2% of the planet’s surface, they contain almost 10% of the Earth’s population [3]. Coastal erosion is the wasting of land along the shoreline due to diverse natural or human causes, such as wave action, wave and tidal currents, high winds, earthquakes, landslides, dam construction, drainage, and other human-induced causes of erosion. Coastal urbanization pressure, rising sea levels, and the increasing frequency of storms and cyclones are also maximizing coastal erosion and putting people’s lives and properties at risk. Bird (1996) estimated that 70% of sandy shorelines are eroding [4]. However, this research had no robust estimation of shoreline change rates because of the lack of constant and long-term monitoring. Thus, recent studies suggest a lower rate of constant erosion for the world’s sandy beaches, in the order of 24% [5].
The usefulness of remote sensing imagery has been practically proven in coastal erosion detection and mapping for many decades. Over the years, remote sensing imagery has been transforming as a result of the evolution of aerial platforms from airplanes and satellites to unmanned aerial vehicles (UAVs). The first attempts to map the coastal zone with aerial photography can be traced back to the 1920s, just after the termination of the First World War, during which the method of photogrammetry was first developed. From the outset, the advantages of aerial photos become obvious: they are easily interpreted, they cover large areas, and they offer the possibility to view and map in 3D and to perform very accurate photogrammetric measurements. In one characteristic example, diachronic aerial images, covering a period of forty years, were photogrammetrically processed for monitoring the shoreline evolution at the southern New Jersey coast [6]. Additionally, aerial photographs were used to examine the role of the presence of Posidonia oceanica on shoreline morphology changes in Sardinia, Italy and Corsica, France [7,8].
The next big step was the launch of the first satellites during the early 1970s. The launch of the Landsat series, despite its low spatial resolution, opened new horizons for coastal area mapping due to the existence of multispectral bands, the repeatability of the scenes captured, and the large area coverage. Landsat Thematic Mapper and Enhanced Thematic Mapper images from different years, covering a period of 25 years, were processed in order to map the shoreline changes in Egypt [9]. Landsat imagery covering a 30-year period has been used for the continuous inspection of coastline changes in western Florida [10]. Similar datasets have been interpreted for detecting changes of the shoreline in India [11], China [12], and Namibia [13]. Furthermore, Landsat data, along with Google Earth images and aerial photographs, were processed in order to perform a diachronic survey of the evolution of the coastline in Lower Casamance and southern Gambia, Africa [14]. The development of very high resolution (VHR) satellite sensors (i.e., Ikonos and Quickbird) around 2000 rebooted the use of remote sensing data in coastal monitoring. The new sensors present special characteristics, such as four to eight multispectral bands, programmability in image acquisition globally, high repeatability (between 24 to 72 h), and storage of the data in digital format. An automatic procedure for shoreline extraction from Ikonos images using a mean shift segmentation algorithm was proposed by Reference [15]. More recently, Sekovski et al. [16] processed a WorldView-2 multispectal image to portray the coastline by combining supervised and unsupervised image classification methods. In the last five years, the explosion of the commercial UAV market has opened new avenues in coastal monitoring, as described in Reference [17]. UAVs present some obvious advantages for coastal surveying. These advantages include the best spatial resolution among all the other remote sensing surveys, a lower cost compared to all other ground or remote sensing surveys, potential to transmit the acquired photos to the remote controller, the ability to repeat acquisition in case something goes wrong, and very low security risks as the drones are unmanned.
As mentioned in the previous paragraph, different types of remote sensing data have been used in coastal environment applications. High resolution satellite data and air photos were used to monitor the coastline’s changes or to measure volume changes in coastal areas [18,19,20,21], sea surface temperatures [22,23], sea ice coverage [24,25], and coastal change processes [26,27]. A review paper summarizing the remote processing methods performed in shoreline delineation has been presented [28]. Recently, the use of UAVs for coastal monitoring was investigated in many studies [29,30,31,32,33,34,35,36,37]. In one of these studies, UAV digital photogrammetry was applied in order to develop point clouds and 3D models of the coastal cliffs in East Sussex, monthly for a one-year period, in order to map their recession [38]. UAV images and structure from motion (SfM) photogrammetry was also tested at Dongshan Island, China to map coastal changes [39].
In order to map the diachronic evolution of a coastline, there is a need to use multitemporal and multisensor data with diverse spatial resolution. In such a study, air photos and the Landsat and Spot data series covering a period of 45 years were combined in order to map the shoreline evolution in Mexico [40]. Diverse multidisciplinary survey techniques, such as terrestrial laser scanning, terrestrial photogrammetry, and structure from motion photogrammetry using UAV data, have been tried in order to evaluate their potential in cliff erosion monitoring [41]. Aerial photos and very high-resolution satellite data, such as that from Ikonos, Quickbird, Worldview-2, and Pleiades, covering a 14-year period, have been used for estimating/defining coastal erosion in southern Thailand [42]. Topographic maps, Google Earth images, and Landsat Thematic Mapper data were combined to map shoreline changes along the mangrove ecosystem of East Indonesia [43]. A study combining aerial photographs from three different dates (1960, 1973, and 1988) and Quickbird imagery from 2014 focused on the changes in coastline location in the Bay of Jijel (eastern Algeria) [44]. Aerial photographs and a variety of VHR satellite images (Quickbird, Worldview-1, and Worldview-2) were analyzed to map the shoreline changes from 1943 to 2012 in Papua New Guinea [45]. Historical maps, analogue aerial photos, and airborne Light Detection and Ranging (LiDAR) data were compared in order to evaluate the shoreline change for the period 1881 to 2015 [46]. Historical aerial photos and VHR satellite data were processed for shoreline mapping on the Marshal Islands [47,48]. In the context of air photos and satellite images in combination with a geographical information system (GIS), a Digital Shoreline Analysis System (DSAS) and field surveys were utilized for the assessment of coastal erosion and the kinematics of the coastline in order to understand how natural and anthropogenic factors affect the diachronic coastline evolution [49].
In this study, we use a 73-year inventory period motivated by the strong accretion and erosion rates in the southwestern coastal area of Lefkada (Figure 1). In Figure 1, the scale remains at the same position at different dates, indicating the decrease or increase of the beach width. From 1945 to 1972, the erosion rate was quite small; between 1972 and 2008, the erosion increased; and from 2008 to 2016, the phenomenon is reversed and a high accretion rate is observed. Multi-date analogue aerial photos, digital aerial photos, declassified satellite imagery provided by United States Geological Survey (USGS), Pleiades satellite data, and UAV images were combined for the accurate mapping of the western Lefkada coastline during the last 73 years. This is the first time that such a range of data is being analyzed and compared simultaneously. As far as we know, studies on Pleiades triplet data focused mainly on the vertical accuracy of the derived digital surface model and the horizontal accuracy of the produced orthomosaics [50,51,52,53,54]. Only in one case study [42] were Pleiades data combined with other VHR satellite images in order to map shoreline changes. In the current paper, Pleiades triplet data are processed and compared to air photos and UAV data. Another novel characteristic of the current study is the photogrammetric processing of declassified satellite imagery for coastline evolution. Although more than 24 years have passed since the liberalization of these datasets, there are no studies describing their photogrammetric processing. Furthermore, aerial, satellite, UAV photogrammetry, and geographical information systems are handled together for the remote sensing data processing and the shoreline erosion mapping. In addition to the different remote sensing data used in this paper, the geologic environment is equally fascinating as most of the accumulated sediments in the coast come from active landslides. Thus, the specific work tries to correlate the coastline evolution with the earthquakes through time.
The remainder of the current manuscript is divided into five sections: In the next section, the geological and geomorphological status of Lefkada Island is described. Section 3 presents the materials used and the performed methods. In Section 4 the results are stated, and Section 5 discusses and interprets the results. The final section presents the conclusions.
2. Study Area
The northwestern-most Hellenic arc, where the Ionian Sea islands are located, is regarded as a site of complex continent–continent to continent–ocean convergent plate margins [55]. The collision zone in the north is separated from the continent–ocean convergence margin in the south by the Cephalonia Transform Fault zone (CTFZ). The CTFZ (Figure 2a) is a major dextral strike-slip structure [56], considered as a highly active seismotectonic structure hosting most of the earthquakes in Cephalonia and Lefkada Islands (Figure 2a and Table 1). The fault is separated in a southern segment that strikes SW–NE offshore to the west of Cephalonia Island. The northern segment of the CTFZ strikes SSW–NNE offshore to the west of Lefkada Island. The total length of the CFTZ is about 60 km [57]. Most known large earthquakes in the study area (Table 1) are the cluster of events during the 1953 Ionian earthquakes that seriously damaged the Cephalonia, Zakynthos, Ithaki, and Lefkada Islands [58,59]. Several large earthquakes have also damaged Lefkada Island in the historical period as well as in the instrumental era of seismicity [58]. Many historical earthquakes have also been documented with the largest of them occurring from 1577 up to 1869 AD at magnitudes ranging from 5.9 to 6.5 [60]. The 2003 Mw6.3 [61] and the 2015 Mw6.4 strong earthquakes [62] are regarded as key earthquakes for understanding the seismotectonics of the northern segment of the CTFZ. Most of these earthquakes are large enough and located close enough to Lefkada to trigger significant landslides in the study area [63,64,65]. Seismicity in the area shows temporal gaps of between 12 and 30 years (Table 1).
In addition to the CTFZ, the so-called Athani Fault is mapped onshore in the study area controlling the western Lefkada coast [66,67]. The Athani Fault is defined for almost 20 km, crossing lengthwise the west coast of Lefkada Island (Figure 2b). The Athani Fault represents a dextral strike slip fault that juxtaposes Paxos unit limestone against Pleistocene to Holocene deposits. In addition to the Athani Fault, other smaller faults, being almost parallel with the CTFZ and Athani Faults with similar kinematics, define steep slopes hosting landslides during earthquake activity [65]. The November 2015 earthquake of Lefkada, 12 years after the 6.3 earthquake of 2003, affected Lefkada’s coastal cliffs. The landslides and the coastal area, on which this work focuses, are at obvious structural proximity with two major faults in the area, namely, the Athani Fault and CTFZ. Thus, the western part of Lefkada Island is an ideal study site characterized by high and frequent rates of seismicity (Table 1) and crossed by active faults.
3. Materials and Methods
3.1. Materials
In the current study, different datasets of airborne, satellite, and UAV platforms were combined and processed using diverse photogrammetric techniques. Acquisition dates, number of data, spatial resolution, and other characteristics of the datasets are presented in Table 2.
The first historical dataset is an orthophoto mosaic of 1945 created for the needs of the Greek Cadastre. The specific digital orthomosaic was created with photogrammetric methods from analogue aerial photographs acquired in 1945. Having a pixel size of 1 m, it covers the whole Greek territory. The specific dataset was developed by the National Greek Cadastre and Mapping Agency. We did not perform any further process on it.
The second dataset includes 3 sets of declassified satellite imagery obtained in the late 1960s–1970s by American military missions and declassified for the first time after 1995. The specific declassified photographs were acquired by U.S. military satellites and give globally important information of the planet surface [70]. This imagery was, for many years, classified as top secret, and it was delivered for free use in 1995. The specific archive contains more than 990,000 photographs. The photos were acquired between 1959 and 1980. The images present a variety of scales and quality. Very often, the scenes are cloudy. The pixel size of the photos varies and ranges between 6 and 30 feet (2–9 m). The film is scanned at 7 microns, and the final total size of an image can surpass 1.3 GB [71]. The first image of the specific dataset was collected during the CORONA KH-4B mission in May 1972 (Table 2). The specific satellite collected the photographs with a telescopic camera, and the film returned to the earth via recovery capsules. While the first Corona missions carried a single panoramic camera or a single frame camera, the next satellites (KH-4, KH-4A, and KH-4B) had double panoramic cameras on board, each looking at 15 degrees forward or backward from the satellite orbit. The pixel size of the image was around 6 feet, but according to previous studies, the KH-4B missions provided the best spatial resolution (1.83 m at nadir [72,73]. The second and the third images of the specific dataset were collected during the CORONA KH-9 program missions. The KH-9 program was active for 7 years (1973 to 1980). The specific mission collected images according to a predefined grid in order to eliminate the image distortion. The images were acquired with overlap for stereoscopic analysis. The KH-9 system produced 9 × 18 inch photos with a pixel size of 20 feet. During the 7-year mission, KH-0 collected almost 29,000 images in 12 accomplished space journeys. The second image of the declassified satellite imagery dataset was acquired in June 1974 during the mission 1208-5. The third image was acquired in October 1980 during the mission 1216-5. The spatial resolution of both images is better than 20 feet, and they were both scanned at 7 microns.
The third dataset comprises 8 analogue aerial photographs from 1992, at 1/15,000 scale accessed through the Hellenic Military Geographical Service (HMGS). With 60% along the track overlap, these photos developed 7 stereopairs and were selected as they combine the best spatial resolution and are cloud free.
The fourth dataset comprises eight analogue aerial photographs from 2000, at the same scale (1/15,000), and from the same source (HMGS).
During the period 2007–2009, digital air photos covering the whole of Greece were acquired by the National Greek Cadastre and Mapping Agency. After photogrammetric processing, an orthomosaic with a pixel size of 50 cm and a digital surface model with a spatial resolution of 5 m were produced. Both covered the whole country, and they are used as basemaps in many studies as they present the highest horizontal and vertical accuracies. The orthophoto and the Digital Surface Model of National Greek Cadastre and Mapping Agency covering Lefkada Island belong to the fifth dataset used in this study.
The sixth dataset consists of 2 scenes of Pleiades satellite imagery. The Pleiades system comprises two satellites in the same orbit but with 180° offset. The first satellite was launched in 2011, while the second satellite was launched 1 year later. Each satellite simultaneously collects 1 panchromatic and 4 multispectral bands. The panchromatic band has a pixel of 0.7 m, while the 4 multispectral bands have a pixel size of 2.8 m [74]. Lefkada Island is covered by two scenes (triplets). The first scene was acquired on 26 September 2016, and the second was acquired on 13 October 2016. Each scene (triplet) was composed of 3 images (tri-stereo), and both scenes were totally free of clouds.
A total of 462 UAV photographs comprise the seventh and final dataset. These were acquired by the University of Patras team during field work in November 2018. A photogrammetric flight was performed in the Egremni beach area. The flight’s altitude was 60 m above ground level (agl). The along-the-track overlap of the photos was 90%, while the respective across-the-track overlap was 75%. Table 3 presents all the flight details.
3.2. Methods
The 7 datasets noted in the previous paragraph were processed using different photogrammetric and computer vision techniques discussed below. We used two software packages, ERDAS IMAGINE Leica photogrammetry Suite (LPS) and Agisoft Photoscan Professional, for the processing of the various remote sensing data. The Greek Cadastre orthophoto and Digital Surface Model (DSM) of 2008 were used as a basemap and reference for all the other datasets. According to the Greek Cadastre specifications, the orthomosaics have a horizontal (planimetric) accuracy of 0.71 m while the vertical accuracy of the respective DSM is better than 2 m.
3.2.1. Analogue Air Photo Data Processing
In LPS (2014 release), a block was created and the analogue air photos were imported into it (Figure 3). The internal triangulation of the air photos was calculated using the fiducial points on each photo. Furthermore, there is a need for external Ground Control Points (GCPs) to orthorectify the imagery (Figure 3). There are many possible error sources in an air photo block, such as lens distortion, film distortion, and atmospheric refraction [75]. Each of these may possibly decrease the accuracy of aerial triangulation results and the accuracy of the final orthomosaic. The major parameters which increase or decrease the accuracy of the triangulation and the accuracy of the final orthomosaic are the following: the position of the air photos (block geometry), the allocation, the quantity and the accuracy of the points that are used (control or tie points), and the existence of any other random error. The whole procedure is described in detail in a previous study [75] and cannot be repeated in the current paper. As the block geometry and the existence of random errors are predefined, the main duty of the user is to select with high precision the GCPs and the tie points.
In the first LPS block of analogue air photos of 2000, we used 71 GCPs, while in the second LPS block of 1992 air photos, we used 54 GCPs. The coordinates of GCPs were obtained from the Greek Cadastre orthomosaic, while the respective altitude was obtained from the respective DSM. The GCPs had a very good distribution around the broader area with an emphasis given on the area across western Lefkada beach (Figure 3). Then, LPS was used to calculate the root mean square error (RMSE) values for the entire block, as seen in Table 4. The RMSEs calculated for the two blocks are assumed acceptable in order to proceed to the orthorectification of the photos and the creation of the final orthomosaics, one for the May 1992 air photos and one for the September 2000 air photos.
3.2.2. Declassified Satellite Imagery Processing
The 3 declassified satellite scenes were processed in 3 different blocks in LPS, one for each date. In a previous study [76], LPS was used for the analyses of CORONA KH-4 imagery. That study demonstrated the huge contribution of these data process for archaeological studies in the Near East. Stereo models, digital elevation models (DEMs), and orthomosaics have been successfully produced utilizing a simplified frame-model in LPS for orthorectification of the CORONA KH-4 images. The GCPs used in that study [76] were selected from freely available datasets. LPS was also used for the orthorectification of CORONA KH-4 data in another study [77]. Ikonos Geo product imagery and Shuttle Radar Topography Mission DEMs were used for the orthorectification of the declassified image. CORONA imagery has been processed in ERDAS IMAGINE OrthoBASE Pro in order to produce a DEM [78]. Sohn et al. [73], among others, used the same software to assess the quality of CORONA imagery derived, and Reference [79] proposed a methodology for CORONA image processing based on ERDAS IMAGINE toolboxes. CORONA KH-9 imagery was processed with satisfactory results in the past [71,80].
In the current study, the 3 declassified scenes were orthorectified using LPS (2014 release). Around 50 points were selected from the Greek Cadastre orthophoto to serve as the GCPs for each declassified scene. Particularly for the 1972 imagery, 58 ground control points allocated all over the scene were used (Figure 4). The declassified imagery of 1972 was orthorectified with a pixel size of 2 m, while the respective images of 1974 and 1980 were orthorectified with a pixel size of 4 m.
3.2.3. Pleiades Data Processing
Pleiades tri-stereo images were processed using also LPS. The exact procedure for Pleiades triplets processing is described in Reference [54]. Initially, all 6 panchromatic images that group the two triplets were inspected for radiometric unconformities during the image acquisition as proposed in Reference [54]. Commonly, radiometric anomalies show up only in one of the two or three images (stereo pair or triplet). They are caused by the existence of surfaces of highly reflective objects (e.g., metal roofs and water surfaces). A different acquisition angle may produce sun glint in one of the images and not in others [50]. None of these effects was noted in the present set of Pleiades images. The coordinates of 180 GCPs were obtained from the orthomosaic provided by the National Cadastre and Mapping Agency S.A., while the elevation was retrieved by the respective DSM.
3.2.4. UAV Data Processing
UAV images were imported into the Agisoft PhotoScan software. The software combined computer vision techniques and structure from motion (SfM) photogrammetry to achieve direct georeferencing or bundle adjustment with ground control points (GCPs) [75,81,82] or simple similarity transformation over the whole block without GCPs. Alike to classic photogrammetry from air photos or satellite stereopairs, SfM photogrammetry uses overlapping images taken from different points of view. The major variation between the two approaches is that SfM defines the internal camera characteristics, the position of the camera, and the orientation of each image in an automatic way without the need for any prior knowledge or grid determination [83]. All the necessary calculations for the internal orientation are estimated automatically by a repeat procedure called “bundle adjustment”. Bundle adjustment reflects the automatic localization of common characteristics in a group of overlapping photos [75]. The geometry of the whole scene is built as more overlapping images are processed and more mutual objects are detected and related. The requirement for many overlapping photos in order to cover an entire area of interest generated the procedure name of structure from motion photogrammetry or photogrammetry produced from a moving camera [82].
During the flight campaigns, artificial targets were spread out in the broader area and measured with a differential Global Navigation Satellite System (GNSS) receiver (Leica GS08). These artificial targets were easily detectable in the UAV images and used as ground control points. The use of artificial targets as ground control points is described in detail in a previous study [17]. The allocation of these targets is presented in Figure 5. As described in more detail in another case [75], the GNSS sensor was receiving corrections—through the GSM network—from the Greek Hellenic Positioning System (HEPOS). The receiver collects signals from GPS (L1, L2, and L2C frequencies), GLONASS (L1 and L2 frequencies), and Satellite-Based Augmentation Systems like Wide Area Augmentation System (WAAS), European Geostationary Navigation Overlay Service (EGNOS), etc. All the GCPs’ measurements have a horizontal accuracy better than 1.3 cm and a vertical accuracy better than 2 cm.
3.3. Historical Shoreline Analysis
To quantify shoreline changes, we separated the 73-year period from 1945 to 2018 into 3 subperiods. These were 1945–1972, 1972–2000, and 2000–2016. This time separation enabled the detailed analysis of the erosion and accretion in rates and in absolute lengths highlighted along the shoreline trace. For the shoreline analysis, two different approaches were used. Firstly, the shoreline was digitized from all the datasets. Then, two different flowcharts—one fully automatic and one semiautomatic—were followed for the shoreline evolution mapping, as described in the next paragraphs. The decision to use two different flowcharts was based on a previous study [84] which proved that different techniques can provide very correlated results for smooth shorelines but less correlated results for irregular shorelines.
3.3.1. Digital Shoreline Analysis System (DSAS)
A software package called the Digital Shoreline Analysis System (DSAS) was used in order to measure the shoreline rates-of-change by comparing vectorized shorelines from diverse dates [85]. DSAS is an add-on tool of ArcGIS software, developed by the USGS in conjunction with the TPMC Environmental Services. The DSAS software creates transect lines perpendicular to the coastline with reference to a specific baseline set by the user. The user also determines the spacing of the transect lines across the shoreline [40]. The shift of the shoreline seaward or landward, with reference to the baseline, is characterized as accumulation or erosion at each transect, respectively, and the statistical values are considered as positive for accumulation and negative for erosion. The software calculates diverse statistical values in order to measure the shoreline position change [85]. The main values are the end point rate (EPR), which measures the rate of the coastline change between two successive shorelines, and the net shoreline movement (NSM), which calculates the total distance between the successive shorelines. The specific method for coastline change detection has been used in the past with different types of remote sensing data. For example, DSAS was used with multi-date satellite images from the Indian sensors (IRS P6 and LISS-III) in order to extract the shoreline change on the western coasts of India [86]. The same software was used with diverse Landsat Thematic Mapper and Enhanced Thematic Mapper images during a 25-year period to derive the shoreline changes in Egypt [7]. In another example, diachronic aerial photographs and Quick-bird satellite data were used in combination with DSAS [44] in order to detect and measure the shoreline movements at the Bay of Jijel (eastern Algeria). Multitemporal and satellite dataset have been interpreted with DSAS for the change detection in 112 km of shoreline in India [11]. Aerial photographs and a variety of VHR satellite images (Quickbird, Worldview-1, and Worldview-2) were analyzed to map the shoreline changes from 1943 to 2012 in Papua New Guinea [45] with DSAS. Multi-date Landsat images were also processed with the same software [12]. Historical aerial photos and VHR satellite data were processed for the shoreline mapping in the Marshal Islands using DSAS [47,48].
In the current study, DSAS utilizes the digitized shorelines in reference to a baseline projected to the same reference system (Greek Grid). The shorelines cover the following years: 1945, 1972, 1974, 1980, 1992, 2000, 2008, 2016, and 2018. For each shoreline vector, the software requires that the date is predefined in the format year/month/day. The software creates transects which intersect the multi-date shorelines at specific points, which are used for the calculations of the EPRs. The EPR calculation is presented in Figure 6, while the general flowchart of DSAS is presented in Figure 7.
3.3.2. Semiautomatic Processing in GIS Environment
For the shorelines of Egremni and Gialos in SW Lefkada Island, we also applied semiautomatic processing in a GIS environment for extracting shoreline alterations through time. The methodological procedure was based on the creation of 102 transect lines spaced every 100 m. Their contact point with the studied shorelines identifies the coastal point viewed in a distance diagram. The baseline, which was the road from Porto Katsiki to Athani and from Athani to Komilio that can be identified in images more recent than 1972, was projected to images of all different ages. In Figure 8, the digitized road from the 1972 orthophoto is overlaid in the orthomosaics of 2000 and 2016. The digitized roads confirm the accuracy of the derived orthomosaics of different ages because their spatial adjustment in all the orthomosaics can be identified easily and they also act as testimonies of the used methodology (Figure 8). Thus, the semiautomatic procedure used in this paper represents an evaluation tool for the DSAS methodology and attaches great importance in the study to the shorelines of Egremni and Gialos of Lefkada Island in terms of accretion and erosion.
4. Results
The images of 1945, 1972, 1974, 1980, 1992, 2000, 2008, and 2016 and a UAV flight campaign during 2018 were used to analyze the rate of change at two significant and popular tourist beaches, Egremni and Gialos, in which both include similar sandy and rocky parts (Figure 2b). These shorelines are regarded as significant for the geomorphological evolution of the west part of the island because a series of factors seem to have a common impact on their evolution. These factors are the wave action; the seismicity; and, for the years after 2000, the human factor. The 73 years of inventory of both beaches revealed remarkable changes over time across the shoreline in the previously referred periods. In these periods, we mapped high erosion and/or accretion rates and absolute length changes. We also defined shoreline trends from erosion to accretion and vice versa (Figure 9). Analytically, our results over discrete periods are summarized in Figure 9. In total, the shoreline mapping results are presented in three figures (Figure 9, Figure 10 and Figure 11) and in one table where the maximum erosion and accretion rates are highlighted in each time period (Table 5). In Figure 10 and Figure 11, maximum erosion (up to 149 m) and maximum accretion (up to 44 m) can be observed. Over the entire coast, from the Porto Katsiki headland to the northernmost end of the studied shoreline (Gialos headland), erosion affects 6300 m and accretion affects 3700 m, indicating that erosion has prevailed over accretion during the last 73 years.
For the period of 1945–1972, erosion was dominant in six areas where the erosion annual rate surpasses 0.50 m and three large areas where accretion surpasses 0.50 m per year (Figure 9a). Overall, 4800 m of the coast is under erosion while accretion extends to 5200 m, suggesting that the area of erosion is almost equal to the area of accretion. In the south (Egremni beach), accretion in the period 1945–1972 shows a maximum of 51 m while maximum erosion is 54 m (Figure 10a). In the Egremni beach, the accretion is active in 2800 m, or 60% of the total length of the beach, with erosion prevailing in 1900 m, or 40% of the beach length. For the same period in the Gialos beach, accretion is active in 2400 m, or 45% over the beach length, and erosion in 2900 m, or 55% of the beach length (Figure 11a). Maximum accretion in Gialos beach is as much as 22 m and maximum erosion is 30 m.
During the period 1972–2000, we defined the change from erosion into accretion in comparison to the 1945–1972 period. The new pattern depicts erosion into areas of accretion and vice versa (Figure 9b). At the southern end of the shoreline, erosion prevailed in both periods but now has an impact on a wider area (Figure 10b). At Egremni beach during this period, erosion is remarkably high, up to 70 m in the north and up to 38 m in the south (Figure 10b). Overall, the accretion is active in 1400 m, or 30% of the total beach length, while erosion is prevailing in 3300 m, or 70% of the beach length (Figure 10b). In the Gialos beach, accretion is active in 2800 m, or 53% of the beach length, and erosion in 2500 m, or 47% of the beach length (Figure 10b). Maximum accretion in Gialos beach is as much as 62 m, and maximum erosion reaches 37 m (Figure 11b).
Over the 2000–2016 period, there was again a change from an erosion to an accretion trend in comparison with the previous period (1972–2000). The new pattern once more shows erosion moves into areas of accretion and vice versa (Figure 9c). The southernmost end of the shoreline remains constantly under erosion in this period (Figure 9c). At the Egremni beach during this period, erosion is remarkably high, up to 42 m in the north and up to 90 m in the south (Figure 10c). Maximum accretion is 38 m at the north end of the Egremni beach (Figure 10c). Overall, the accretion is active in 2000 m, or 42% of the total length of the beach, while erosion is prevailing in 2700 m, or 58% of the beach length (Figure 10c). At the Gialos beach, the accretion affected an area 3400 m long, or 64% of the beach length, and erosion affected 1900 m, or 36% of the beach length, in the period 2000–2016 (Figure 11c). Maximum accretion in Gialos beach is as much as 22 m, and maximum erosion is 20 m. In summary, over the 73 years, both beaches were under erosion, with the highest erosion observed at the southernmost end of the Egremni beach near the Porto Katsiki headland (149 m) (Figure 10d), and the highest accretion was observed at the southern end of Gialos beach near Egremni headland (44 m) (Figure 11d).
During the 73 years, there are periods related to intense seismicity and periods with low or absence of seismicity (Table 1). Significant seismicity is recorded at the beginning of the first period, during 1948 and 1953. In the second period, the biggest earthquake occurred in the middle of the 1972–2000 period. The most recent period of seismicity was recorded in 2003 and 2015. Based on these data, our mapping of shorelines corresponds to significant, low, and very high seismicity. Based on our unpublished mapping data relating to the 2003 and 2015 earthquakes, significant landslides and mass wasting across the shoreline were triggered due to these two earthquakes.
5. Discussion
The phenomenon of shoreline changes is quite complex, affected by oceanographic and meteorological conditions and geological parameters. In general, the study area is influenced by west to northwest prevailing winds that rarely blow above a maximum speed of 6 Bf [87]. The tide wave action is characterized by low amplitude as in the rest of Mediterranean Sea, and the average offshore wave height is almost 0.79 m [88]. However, in this study, as oceanographic and meteorological data and studies are missing for Lefkada Island, we have to focus exclusively on geological parameters such as seismicity, active deformation, and mass wasting across the west cliff of Lefkada in order to explain the shoreline changes. Those changes were mapped using multi-date and multisensor remote sensing data combining, for the first time, aerial photos, declassified satellite imagery, Pleiades data, and UAV data. The photogrammetric processing of these data and the very low RMSE (Table 4) resulted in the excellent georeferencing of the diachronic data (Figure 8). The extensive field work provided the necessary number of ground control points for the accurate georeferencing of the UAV data in order to combine them with the other remote sensing data. This is fully in accordance with previous studies [89,90] that noted that the registration of UAV imagery with any other remote sensing coarser product depends on the use of ground control points that succeed in highly accurate georeferencing. However, the requirements of this procedure are by far more demanding in time and ground support [89].
The study area includes two sweeping beaches, located at southwestern Lefkada Island, separated by a headland developed between the Gialos and the Egremni beaches (Egremni headland). To the south, the study area ends at the Porto Katsiki headland and, to the north, at Gialos headland. The beach to the north of the Egremni headland is characterized by maximum accretion of 44 m (Figure 11d), while north of the Porto Katsiki headland, erosion dominates, reaching up to 149 m in the time period from 1945 to 2016 (Figure 10d). Since erosion or accretion is constantly sustained through the 73-year period in these two locations, we consider that the beach topography and oceanographic conditions are the main mechanisms controlling erosion or accretion at these specific sites. In contrast to these sites, all other shoreline lengths show remarkable changes from erosion to accretion and vice versa. Trying to highlight these changes in a better way, we considered seismicity as a major factor controlling sediment equilibrium and active deformation through earthquake induced landslides along the Egremni–Gialos beaches.
Indeed, the study of different time imageries confirms that the landslides in the overhanging cliff are directly related to earthquake occurrence and affect the forest development or retreat (Figure 12). Figure 12 presents the diachronic evolution of one of the landslides in the overhanging cliff in Egremni beach. It is obvious that from 1945 to 2018 that the extent of the landslide was almost doubled, offering waste material to the Egremni beach. The usual material that covers Lefkada beaches is mostly coarse sand with a small percentage of pebbles [87]. However, both Egremni and Gialos Beaches present a significant percentage of medium to large pebbles and granules (Figure 12) attributed to mass wasting phenomena evolving in both beaches. Thus, our multidate analysis indicates that, primarily, the seismic events of 1948, 2003, and 2015 [69] intercalated in two out of three time periods that represented crucial parameters for cliff erosion and sediment accumulation along the studied shoreline. However, Reference [69] indicated that the analyzed shoreline during the most recent seismicity has subsided a few centimeters and Ganas et al. 2016 [91] suggested extensive coastal road failures related to the outbreak of numerous earthquake-induced landslides. Thus, the changes of the two periods 1945–1972 and 2000–2016 are regarded as the result of the strong earthquakes of 1948, 1953, 2003, and 2015 (Figure 2a and Table 1). Their immediate impact is the prevalence of shoreline accretion, especially in the sandy part of both beaches (Figure 10a,c and Figure 11a,c). In contrast, during the 30-year period of seismic quiescence (from 1953–1983) and especially in the 1972–2000 period, only one earthquake affected the area (Figure 2a and Table 1); hence, the main process is the shoreline erosion (Figure 10b and Figure 11b). On the contrary, the two other seismic quiescence periods appear incapable of producing significant changes on the accretion/erosion equilibrium.
It is worth mentioning that both the fully automatic DSAS method and the semiautomatic method gave the same results (Table 5). The specific result is in accordance with the results derived in a previous study for a smooth shoreline in the Tyrrhenian Sea, Italy [84]. Both methodologies can be applied to vectorized coastlines independently from their originality (raster basemap used for digitizing the coastline). DSAS was, in the past, used with shorelines digitized from aerial photographs and VHR satellite images [44,45,47,48] or medium resolution satellite data [7,11,12,86]. The main concept in all these studies is that the raster data should have a similar spatial resolution and very accurate georeferencing. Those two prerequisites have been fully accomplished in the current study, as presented in Table 2 and Table 4.
Particularly important for understanding the phenomenon of erosion and/or accretion is the last period of monitoring with UAV in the north part of the Egremni beach. During this monitoring, we mapped extensive earthquake-induced landslides and prevalence of erosion. These highly accurate UAV data indicate beyond any doubt that, for a better understanding of the mass wasting after each earthquake, further remote sensing data acquisition and analysis is needed in order to gain knowledge about the exact time that a strong earthquake effect stops its impact as a sediment feeder for the affected shoreline. The basic advantages of the UAV SfM photogrammetry (low cost, flexibility, and high accuracy) proved the feasibility of the specific technology in coastal change mapping. The same conclusion was derived in a recent study [92] using UAV to monitor wind- and wave-driven morphological changes on a beach-dune at Truc-Vert in southwest France.
6. Conclusions
The photogrammetric processing of diverse remote sensing data (air photographs, satellite images, and drone-acquired footage) and further analysis in a GIS environment indicates that the southwest shoreline of Lefkada Island is under dynamic equilibrium. This equilibrium is strongly controlled by geological parameters such as the subsidence of the studied shoreline during co-seismic deformation and mass wasting. The following observations can be noted:
Headlands appear to control accretion and/or erosion along the sweeping beaches of SW Lefkada Island.
Periods of strong seismicity predate shoreline accretion, while periods of relatively long seismic quiescence, in the order of 30 years, are related to erosion.
Periods of strong earthquakes give rise to the change of the former prevailed condition. Areas of erosion change to accretion areas and vice versa.
Mass wasting is dramatically related to the first years after earthquakes, and sediment dispersal signifies the long-term phenomenon of sediment drift along and across the shoreline.
Both famous beaches will maintain their sediment budget, and an analogous evolution will be sustained as long as strong earthquakes occur in the future.
Furthermore, for the first time, it was demonstrated that photogrammetric processing is appropriate for declassified satellite imagery provided by USGS. In order to map the diachronic evolution of a coastline, diverse remote sensing data can be combined. UAVs can be successfully combined with other remote sensing data and have proven to be a very cheap, accurate, and flexible method for coastline monitoring.
Author Contributions
Conceptualization, Konstantinos Nikolakopoulos; methodology, Konstantinos Nikolakopoulos and Vasiliki Zygouri; software, Konstantinos Nikolakopoulos, Aggeliki Kyriou, Vasiliki Zygouri, and Dionysios Apostolopoulos; validation, Konstantinos Nikolakopoulos, Aggeliki Kyriou, and Vasiliki Zygouri; field work, Konstantinos Nikolakopoulos, Ioannis Koukouvelas, Vasiliki Zygouri, and Aggeliki Kyriou; data curation, Aggeliki Kyriou, Vasiliki Zygouri, and Dionysios Apostolopoulos; writing—original draft preparation, Konstantinos Nikolakopoulos, Ioannis Koukouvelas, and Vasiliki Zygouri; writing—review and editing, Konstantinos Nikolakopoulos, Ioannis Koukouvelas, and Vasiliki Zygouri; project administration, Konstantinos Nikolakopoulos and Ioannis Koukouvelas.
Funding
This research received no external funding.
Acknowledgments
The authors would like to thank the “EnCeladus” Hellenic Supersite (
Conflicts of Interest
The authors declare no conflict of interest.
Figures and Tables
Enlarge this image.Figure 1. Diachronic images showing a rock outcrop in the central part of the Egremni beach: Note the characteristic erosion and accretion of the shoreline over the years. The scale bar is located at the same point in all the images. (a) Air photo mosaic of 1945. (b) Declassified satellite imagery of 1972. (c) Air photo mosaic of 1992. (d) Air photo mosaic of 2000. (e) Cadastre mosaic of 2008. (f) Unmanned aerial vehicle (UAV) mosaic of 2018.
Enlarge this image.Figure 2. (a) The main structural features that affect the seismotectonics of the Ionian Sea: The stars provide the epicenters of the seismic events presented in Table 1. (b). A simplified display of the main geological characteristics of Lefkada island: The black box delineates the study area.
Enlarge this image.Figure 3. Selection of ground control points during the orthorectification process of the air photos of 2000: (a) Cadastre orthomosaic displaying the study area. (b) One of the air photos of 2000 displaying the study area. (c) Enlargement of Figure 3a. (d). Enlargement of Figure 3b. (e). Focus on one ground control point on the Cadastre orthomosaic. (f) Focus on the respective ground control point on air photo of 2000.
Enlarge this image.Figure 4. Allocation of the ground control points on the declassified satellite imagery.
Enlarge this image.Figure 5. Allocation of the ground control points measured with differential Global Navigation Satellite System (GNSS) receiver on the UAV orthomosaic.
Enlarge this image.Figure 6. An example of the end point rate (ERP) calculation using 1972 and 2000 shorelines at the study area.
Enlarge this image.Figure 7. Flowchart of the algorithm of Digital Shoreline Analysis System (DSAS) v4.3.
Enlarge this image.Figure 8. (a) The study area where the DSAS and semiautomatic processing in geographical information system (GIS) were performed: Orthomosaics from (b) 1972, (c) 2000, and (d) 2016. In red color, a digitized road from the 1972 orthophoto is overlaid in all mosaics.
Enlarge this image.Figure 9. The processing results of DSAS software for the study area: Red colors represent erosion rate per year, while green colors represent accretion rate per year. (a) Six areas where the erosion annual rate overpasses 0.50 m and three large areas where accretion overpasses 0.50 m per year for the period of 1945–1972. (b) Erosion into areas of accretion and vice versa for the period of 1972–2000. (c) One more change from the erosion to accretion trend in comparison with the previous period (1972–2000). (d) Areas where erosion and accretion take place over the period of 1945–2016. (e) Highest erosion overpasses 1.50 m per year observed in the southernmost end of the Egremni beach for the period of 2016–2018. (f) Study area in Lefkada Island.
Enlarge this image.Figure 10. The processing results of the semiautomatic method for Egremni beach: Red colored bars represent total erosion in each transect for the specific time period, while green colored bars represent total accretion for the same time period. The blue triangle states the Egremni rock outcrop shown in Figure 1. The yellow line represents the erosion/accretion rate along the Egremni beach for each time period. N and S are abbreviations for north and south compass points. Each diagram represents a specific time period: (a) 1945–1972 period, (b) 1972–2000 period, (c) 2000–2016 period, (d) 1945–2016 period, and (e) 2016–2018 period.
Enlarge this image.Figure 11. The processing results of the semiautomatic method for Gialos beach: Red colored bars represent total erosion in each transect for the specific time period, while green colored bars represent total accretion for the same time period. The blue triangle states the Gialos reference in Figure 9. The yellow line represents the erosion/accretion rate along the Gialos beach for each time period. N and S are abbreviations for north and south compass points. Each diagram represents a specific time period: (a) 1945–1972 period, (b) 1972–2000 period, (c) 2000–2016 period, and (d) 1945–2016 period.
Enlarge this image.Figure 12. The diachronic evolution of one of the landslides in the overhanging cliff at Egremni beach: From 1945 to 2018, the extent of the landslide almost doubled, offering waste material to Egremni beach. (a) The landslide extent as mapped from the 1945 orthophoto. (b) Multicolor polygons represent the landslide extent at different dates. The background image is from the 2008 Cadastre orthophoto. (c) The landslide extent as mapped from the 2018 UAV orthophoto.
The most significant earthquakes affecting Lefkada Island during the studied time period. For more information on the earthquakes see References [62,64,68,69].
| Date | Longitude | Latitude | Depth | Mw |
|---|---|---|---|---|
| 22 April 1948 | 38.68 | 20.60 | Unknown | 6.5 |
| 30 June 1948 | 38.80 | 20.60 | Unknown | 6.4 |
| 12 August 1953 | 38.30 | 20.80 | Unknown | 7.0 |
|
|
||||
| 17 January 1983 | 37.97 | 20.25 | 9 km | 6.5 |
|
|
||||
| 14 August 2003 | 38.79 | 20.56 | 12 km | 6.3 |
|
|
||||
| 17 November 2015 | 38.67 | 20.60 | 11 km | 6.5 |
The datasets used in this study, their source, and their spatial resolution.
| Date | Data Type | Source | Reference System | Number of Photos | Spatial Resolution |
|---|---|---|---|---|---|
| 1945 | Air photo mosaic | National Greek Cadastre and Mapping Agency | Greek Grid | Orthomosaic | 1 |
| 27/05/1972 | Declassified satellite imagery | United States Geological Survey | No reference system | 1 | 2 m |
| 05/06/1974 | Declassified satellite imagery | United States Geological Survey | No reference system | 1 | 4 m |
| 05/10/1980 | Declassified satellite imagery | United States Geological Survey | No reference system | 1 | 4 m |
| 11/09/1992 | Analogue air photos | Hellenic Military Geographical Service | No reference system | 8 | 1 m |
| 30/05/2000 | Analogue air photos | Hellenic Military Geographical Service | No reference system | 8 | 1 m |
| 2008 | Digital colored air photos | National Greek Cadastre and Mapping Agency | Greek Grid | Orthomosaic | 0.5 m |
| 2016 | Pleiades Imagery | Enceladus Greek Supersite | World Geodetic System 84 | 2 triplets | 0.5 m |
| 3/11/2018 | Unmanned Aerial Vehicle Imagery | University of Patras | World Geodetic System ‘84 | 462 | 0.1 m |
Characteristics of the UAV flight campaign and pixel size of the produced orthophoto and Digital Surface Model.
|
|
|
|
|
|
|
|
|---|---|---|---|---|---|---|
| 60 | 90 | 75 | 462 | 4.2 cm | 8.4 cm | 3000 × 4000 |
The number of ground control points used in each Leica Photogrammetry Suite block, the respective number of tie points, and the aerial triangulation root mean square errors.
| LPS Block | No of Images | No of Ground Control Points | No of Tie Points | RMSE (m) |
|---|---|---|---|---|
| Declassified 1972 | 1 | 58 | 0 | 0.0044 |
| Declassified 1974 | 1 | 45 | 0 | 0.0727 |
| Declassified 1980 | 1 | 38 | 0 | 0.1668 |
| Air photo 1992 | 8 | 54 | 83 | 1.9753 |
| Air photo 2000 | 8 | 71 | 92 | 1.8602 |
| Pleiades 2016 | 6 | 180 | 224 | 0.1263 |
Maximum and minimum rates of EPR along the study area.
| Maximum/Minimum Rates of EPR (DSAS Method) | Maximum/Minimum Rates of Accretion/Erosion (Semiautomatic Processing Method) | |||||
|---|---|---|---|---|---|---|
| Period | Maximum Rate (m/yr) | Minimum Rate (m/yr) | Egremni Beach | Gialos Beach | ||
| Maximum Rate (m/yr) | Minimum Rate (m/yr) | Maximum Rate (m/yr) | Minimum Rate (m/yr) | |||
| 1945–1972 | 1.92 | −1.35 | 1.89 | −1.26 | 0.81 | −1.11 |
| 1972–2000 | 2.45 | −2.53 | 1.07 | −2.10 | 2.21 | −1.32 |
| 2000–2016 | 2.24 | −5.73 | 2.38 | −5.62 | 1.44 | −1.25 |
| 1945–2016 | 0.55 | −1.53 | 0.21 | −2.10 | 0.62 | −0.17 |
© 2019 by the authors.
Details
; Zygouri, Vasiliki; Apostolopoulos, Dionysios