A GIS Procedure to Assess Shoreline Changes over Time Using Multi-temporal Maps: An Analysis of a Sandy Shoreline in Southern Italy over the Last 100 Years

: The aim of the paper is to identify a methodology capable of assessing shore-line changes through a geomatic approach based on the use of GIS (Geographic Information System) software. The paper describes a case study that reports the evolution of a coastline over a period of more than 100 years using medium and large-scale metric maps available in different periods. In fact, the coastlines were obtained from the source maps of the Italian Cadastre (dated 1890), from numerical cartography available on the coastline and acquired in different pe - riod at scales 1:5000 and 1:2000 and, more recently, from the Google Earth Pro platform. To analyse the evolution of the coastline a new procedure has been performed which is based on the use of GIS software, in particular a plugin called DSAS that allows the evaluation of the changes in the coastline and also obtains a statistical analysis of its evolution. The results showed the ease and applicability of the method in determining the evolution of the coastline and the strong erosion of a stretch of coastline with important socio-economic consequences and repercussions was highlighted in the analysed case study.


Introduction
Issues related to environmental changes are now firmly on the agenda, and include aspects such as the reduction of biodiversity and the increasingly unpredictable forces of nature, with dramatic consequences for the entire ecosystem.It is therefore necessary to take concrete measures to protect the environment and regulations with precise indications for planning and conducting activities to control and monitor impacts.
Among the different aspects to be monitored may be, for example, the phenomenon of coastal erosion, especially in areas of particular environmental value, and this has led to a growing interest in the morphological changes of coastal areas [1][2][3].Therefore, the monitoring of these areas through the geomatics approach assumes an important role for the conservation and protection of the environmental heritage and, consequently, has led to an increased demand for accurate knowledge of coastal developments.
The temporal and spatial variability of the coastline represents a particularly open and timely issue.Lira et al. [17] discussed the coastline evolution of low-lying Portuguese sandy coastal area over the last 50 years using digital aerial photographs from the USAF (United States Air Force) 1958 flight and a digital orthophoto from the year 2010.
Yu et al. [18] wrote about the automatic extraction of the coastline using remotely sensed data in an Antarctic environment over several years; the evolution of the coastline represents a clear indicator of the change in the extent and mass balance of ice sheets and shelves.Wang (2019) discussed the evolution of Yellow River delta coastline (China) [19] based Multi-Spectral Scanner (MSS), Thematic Mapper (TM), and Enhanced Thematic Mapper plus (ETM+) images of the study area from 1976 to 2014.
To determine changes in shoreline position, as well as to predict future positioning trends, many studies have used change detection tools connected with a GIS (Geographic Information System) environment [20].In this environment, it is possible to calculate the distance between the furthest shoreline from the baseline and the nearest shoreline for each transect as well as to make an assessment of the annual rate of change or determine other statistical parametres useful for its analysis.Several tools and algorithms, such as the Simple Change Analysis of Retreating and Prograding Systems (SCARPS) [21], BeachTools [22], Digital Shoreline Analysis System (DSAS) [23,24] and Analyzing Moving Boundaries Using R (AMBUR) [25] were developed in GIS environment.
In this research paper, a DSAS tool is used to analyse the evolution of a sandy shoreline and this tool has been used with considerable success in numerous research studies and applications [26][27][28][29][30].For example, Nassar et al. [30] described shoreline change detection along the North Sinai coast in Egypt during the period from 1989 to 2016.In the same way, Das et al. [31] attempted to analyse the potential migration characteristics of Jambudwip Island of the Sundarban biosphere in the Bay of Bengal.
Applications on the Mediterranean coastline have been described in the paper by Kuleli on the Cukurova Delta coast (Turkey) [32], in Khouakhi & Snoussi, 2013 [33] which described the evolution of the bay of Al Hoceima in Moroccan Mediterranean coast and in Armenio et al. [34] in their study on the Gulf of Manfredonia in Southern Italy.
Therefore, the line of research proposed in this paper intends to exploit the efficiency of tools developed in the GIS environment for the management of geospatial data as well as the ability to analyse shoreline changes through specifically dedicated tools.In particular, we want to analyse the evolution of a sandy coastline of an area located in the south of Italy and taking into consideration a broad temporal space, i.e. from 1890 to 2019.
Furthermore, the study aims to assign uncertainty values to the geospatial data available on a territory, with the advantage of effectively analysing changes in the coastline.This methodological approach is also a tool for controlling and monitoring the territory to identify trends and apply coastal risk mitigation interventions and strategies accordingly.

Study Area
The study area is located in the town of Eboli (southern Italy) and concerns a coastline which is approximately 7 km long (Fig. 1a, b); in particular, the area taken into consideration concerns a sandy coastline with widths from 15 to 120 m, with a pine forest behind it and the presence of species of particular value in the Mediterranean area (Fig. 1c, d).
In the sandy areas there are also a dozen establishments with bathing facilities, bars, and restaurants with important socio-economic benefits for the local community.This area is delimited from two rivers, the Tusciano River to the north and the Sele River to the south.From a geological point of view, the deposits in this area constitute a prism of progradational sediments deposited on the Sele River plain front during the Holocene [35].
The choice of this area was determined based on the phenomenon of coastal erosion and the availability of various cartographies from 1890 to 2019 on this territory; in this way, it was possible to obtain a comprehensive long-term overview of this phenomenon.

Geomatics Datasets
The features of the maps taken into account to evaluate the changes in the coastline are shown in Table 1.Images of the cartographies used in the paper can be found in the Appendix (Figs.A.1-A.5).In order to assess the evolution of the coastline over the years, it is necessary to compare several geomatics data in the same reference system.To achieve this aim, the global datum is used in relation to technological progress and data standardisation and studies on the dynamism of the Earth's crust.
Over the years, the Italian territory has been framed within local datum; in particular, the most widely used local datum is the Roma40 datum.The transformation from local datum to global datum takes place through the use of Italian Geographic Military Institute (IGMI) grids, where the correction values of latitude and longitude between the two datum are reported.Since we use plane coordinates, the reference system is the Universal Transverse of Mercator (UTM) projection.Therefore, the geodetic problem consists of the transformation of coordinates from the Gauss-Boaga (Italian local datum) into the UTM33N-ETRF2000 projection.GIS software, by performing a roto-translation with a scale factor (Helmert), allows a transformation between the two datum in automatic way.The advantage of this approach is that no a priori information is required for the seven parameters of the similarity transformation [36].The oldest metric map taken into account is the cadastral cartography, which was realized at the end of 1800s.The Italian cadastral reference system, apart from limited areas, is based on the Cassini-Soldner projection and Bessel 1841 datum; this means that this type of projection assumes various local orientations.In Italy, more than 800 small and large extension origins (axis systems) are present on the territory.This means that several transformations are required to report spatial information in UTM-ETRF2000 plane coordinates; in particular, it is necessary to perform a first transformation from Cassini-Soldner to Gauss-Boaga and then in UTM-ETRF2000.The first transformation can be performed, in turn, into three steps: i) transformations of cadastral coordinates into geographic coordinates with respect to the Bessel ellipsoid, known coordinates of the centre of emanation; ii) transformation of the latter coordinates into geographic coordinates referenced to the Hayford ellipsoid; iii) transformation from geographic coordinates of a generic point (φ, λ) into plane (Gauss-Boaga) coordinates using Hirvonen's formulas [37]: where: E -east coordinate in Gauss-Boaga, N -north coordinate in Gauss-Boaga, FO -false origin in Gauss-Boaga, x -east coordinate in Cassini-Soldner, y -north coordinates in Cassini-Soldner.The coordinates x, y can be calculated by means of the following formula: where λ 0 is the longitude of the central meridian of the fuse, a and b represent the semi-axes of the ellipsoid considered and e′ is the second eccentricity.Considering the Hayford ellipsoid, it is possible to write the numerical value of the following parameters: A 1 = 111 092.08210;A 2 = 16 100.59187;A 3 = 16.96942;A 4 = 0.02226; c = 6 397 376.633; e'2 = 0.0067681702.Lastly, the second transformation (from Gauss-Boaga to UTM-ETRF2000) can be performed through a 7-parameter Helmert transformation.

Georeferencing Maps
In order to make an easy metric comparison between the different geomatics products, it is necessary to build a vector file for each era under consideration.In fact, many of the maps or cartographies are in raster format.For this reason, to obtain the coastline in vector format, it is necessary to perform a georeferencing process.In general, the transformation model that can be adopted is a polynomial function, as reported below [38]:

∑∑
where: x, y -coordinates of source system, n -degree of polynomial, X, Y -coordinates in the target system, a ji -coefficients of the polynomials which are computed using Ground Control Points (GCPs) or points deduced directly from the cartography, were used to improve the intersection of the cartographic grid presented in the maps.
The choice of georeferencing map sheets with a polynomial transformation depends on several factors, such as the potential to model deformation effects due to the state of preservation and intrinsic deformation of materials.
In order to assess the impact of the degree of the polynomial in the georeferencing process, first, second and third order were taken into consideration.

Shoreline 1890 -Cadastral map
The oldest metric map used to determine the coastline is the one on the cadastral map sheets.In Italy, the creation of cadastral maps dates to 1890, as can be seen on the edge of the map.The maps were scanned with a large scanner planner at a resolution of 600 dpi in TIFF format.
Three maps cover the area under investigation: sheet no.45, 55 and 56.The scale of representation is 1:4000 and framed in Cassini-Soldner projection.These maps were geo-referenced in Esri ArcMap software using 70 GCPs and 3-degree polynomial function were used for georeferencing the map sheets.The values of the root mean square (RMS), expressed in terms of pixels, as the polynomial function varies, were reported in Table 2.As shown in Table 2, the different maps showed a different RMS value; this was due to the different state of preservation of the map sheets; this issue is consistent with other works in the literature [39,40].In addition, the 1° transformation shows higher values than the 2° and 3° order transformation; however, it was preferred to use a 2° order transformation because the distortions imposed in some areas of the map are too strong.
Once the map sheets were georeferenced, the 1890 coastline was built in Esri ArcMap software v 10.8, according a polyline.Subsequently, the ESRI shape file was transformed in UTM projection.This task was performed using ConveRgo [41], a piece of open-source software developed by the Italian Regions for the transformation of coordinates in the Italian Geodetic Reference System (Decree 10 November 2011 "Adoption of the National Geodetic Reference System").

Shoreline obtained from 1974 and 1987 (aerial photogrammetry maps)
The first large scale map available on the study area was produced by the "Cassa per le opere straordinarie di pubblico interesse nell'Italia meridionale" (fund for extraordinary works of public interest in Southern Italy) and dated to 1974.The aerial photogrammetric map was drawn in five colours and is framed in the Gauss-Boaga reference system at a scale of 1:5000.
The area of interest is covered by four maps.The coordinates of 18 GCPs (evenly distributed on the map) were obtained by the intersection of cartographic grid.Performing the transformation according the several polynomial function, it was possible to obtain the following RMS error values (Tab.3).In 1987, the municipal administration commissioned the drawing up of a 1:2000 scale map for engineering purposes and for the protection of maritime areas.Taking into account the graphic error, the accuracy of this map is 0.4 m or better 0.2 mm × map scale.
From the cartographic point of view, the Gauss-Boaga projection was chosen for the representation.
The maps were scanned at a resolution of 600 dpi and then, in ArcMap environment, georeferenced in the Gauss-Boaga reference system.RMS error values obtained in this process are reported in Table 4. Subsequently, the shoreline of 1974 and 1987 were transformed in UTM projection using ConveRgo software.

Shoreline 1998 (digital aerial cartography)
In 1998, the municipal administration commissioned the drawing up of a 1:5000 scale map for urban planning of the municipal territory.This latter cartography is already in numerical format but framed in the Gauss-Boaga reference system.Therefore, it was necessary to transform the shoreline into the UTM33N-ETRF2000 cartographic projection by means of ConveRgo software.

Shoreline obtained by Google Earth (GE) Pro in 2019
In order to obtain the shoreline for the year 2019, a polygon from Google Earth (GE) Pro was drawn.GE Pro releases free images with elevated details that may provide some potential for regional land use/cover mapping, especially for those regions with high heterogeneous landscapes [42].In particular, the scene was acquired by satellite in very high resolution (VHR), i.e. satellite able to generate image with a panchromatic resolution lower 1 m; in this way, it was possible to obtain the shoreline with a high degree of accuracy.
In order to obtain a polyline of the 2019 year, it was necessary to carry out the following steps: 1. exporting the kml/kmz file from Google Earth Pro; 2. importing the kml file in GIS environment as a layer file; 3. transform the layer file into a ESRI shape file.

A Statistical Method to Assessing Changes in the Shoreline
Taking into account different shorelines (for a specific year), it is possible to analyse the evolution of the coastline over time.In this paper, the shoreline variation was statistically analysed using the Digital Shoreline Analysis System (DSAS) extension of ArcMap software, which is the main component of Esri's ArcGIS suite.The results of all rate calculations are outputted to a table that can be linked to the transect file by a common attribute field.In this way, further analysis of shoreline changes can be produced.
For each transect, the variance propagation law, for a linear case such as the combination of the different uncertainties ( In the case of the shoreline variation analysis, Equation ( 5) can be explicated as follows: σ are related to the digitalization process.In relation to the terms concerning the scale of representation, the graphical error of the map, expressed by the report, was considered as the variance value [43,44] Since the same coefficient depends exclusively on the scale of representation of the cartography examined, a coefficient 2 ( ) R a equal to 1 was assigned.The RMS value obtained during the georeferencing phase of the individual map was used to evaluate the uncertainty in this process.Since the coastline involves several map sheets, for the calculation of the final variance value, a weighted average of the values of RMS i was performed as a function of the length of the coastline L Ci present on the individual sheet compared to its total length As for the terms concerning the datum change process, the variance values were obtained as shown in Equation ( 9), thus comparing the distances taken in the two different reference systems and then measured on the different maps: For the estimation of the variance values obtained in the digitalization process of the coastlines, the dimensions of the pixel were considered, which are obviously a function of the geometric resolution of each cartography.In fact, during the vectorization phase it is very usual to find oneself in the situation of adjacent pixels that show a clear chromatic change or a less clear segmentation; for this reason, the pixels closest to the ones that identify the coastline have been considered and therefore, the one relative to the size of 3 pixels has been chosen as the variance value.In general, the formula is: where GRM is the geometric resolution of map in post georeferencing process.
In this way, to each coastline is associated a variance value; consequently, it is possible to calculate the EPR (end point ratio) and the EPR unc (uncertainty of end point ratio) through the following formulas: where: NSM -distance between oldest and youngest shorelines, t o-r -time between oldest and most recent shoreline.

( ) ( )
where: uncyA -uncertainty from attribute field of shoreline A, uncyB -uncertainty from attribute field of shoreline B, dateA -date of shoreline A (most recent), dateB -date of shoreline B (oldest).

Evaluation of Coastline Accuracy
For the determination of the 2 CD σ value for the year 1890, a double datum trans- formation was considered: the first one is related to the transformation from the Cassini-Soldner reference system to Gauss-Boaga and the second one is from Gauss Boaga to UTM-ETRF2000.In particular, the uncertainty of the first transformation was carried out through a datum shift relative to the emanation centre (Monte Raione).In other terms, the coordinates of the emanation centre (geographic coordinates referred to the Bessel ellipsoid: φ = 40°40'15".853,λ = 6°08'21".005)calculated through field observations and determined in the Gauss-Boaga system, have been transformed into the Cassini-Soldner reference system.This transformation should lead back to the origin of the Cartesian reference system, i.e.O = (0, 0); in fact this transformation involves a residual of about 0.147 m on the east coordinates and 0.08 m on the north coordinates.
The second value of uncertainty was determined, as for the other years, by comparing the distances measured in the two different reference systems and then deduced from the different cartographies.
Therefore, the uncertainties on the coastline referred to the year 1890, can be summarized as follows: 2 CD σ = 0.188 m from Cassini-Soldner to Gauss-Boaga; 2 CD σ = 0.090 m from Gauss-Boaga to UTM-ETRF2000.
a were assigned a value of 1.In Table 5, the values for the variances are summarized and divided for each epoch.

Workflow for the Statistical Analysis of Shorelines
The DSAS plugin was used to analyse multi-temporal shoreline evolution because is capable of calculating rate-of-change statistics for a time series of shoreline vector data.Indeed, DSAS provides an automated method for establishing measurement locations, performs rate calculations, provides the statistical data necessary to assess the robustness of the rates, and includes a beta model of shoreline forecasting with the option to generate 10-and/or 20-year shoreline horizons and uncertainty bands [45].The first main step for the analysis of the shoreline variation was the creation of a personal geodatabase which is used to manage and store all of the input data; in addition, all the statistical data obtained as outputs are also stored.Moreover, another important aspect for the analysis is that all data are in metre units and UTM projection.In general, the procedure that leads to the result for the determination of the statistical values of the time series data was summarized in the following workflow (Fig. 2).

Fig. 2. Workflow for statistical analysis with DSAS plugin
In order to carry out a statistical analysis of coastline change, it was also necessary to define a positional uncertainty associated with each coastline.In fact, each feature associated to the coastline can come from different cartographic sources, which then lead to measurement uncertainties such as digitalization errors and georeferencing errors.In the present study, the variance values previously calculated for each dataset were inserted in the appropriate section of the plugin, representing the uncertainty value in terms of positioning and measurement.It must be emphasised that in order for the tool to function correctly, the attribute table must be set correctly.

Results and Discussions
The processing of the data conducted in a GIS environment allowed an assessment of the change in the shorelines; in particular, to analyse the accretion or erosion phenomena, four periods were taken into consideration: shoreline changes from 1890 to 1974, from 1974 to 1987, from 1987 to 1998, and from 1998 to 2019.
The analysis of the variation of the shoreline (Fig. 3a) was carried out using the DSAS plugin in order to obtain parametres such as NSM, EPR and EPR unc between the shorelines.In addition, a statistical analysis was performed on the linear regression of the linear regression rate (LRR) of change where all data were computed, regardless of changes in trend or precision.
Subsequently, the results obtained from the LRR parametre were categorized into five different classes; the categorization of the results then allowed the production of a cartography (Fig. 3b) where the erosion or accretion conditions of the coastline can be easily visualized and analysed.In order to evaluate the shoreline change with respect to each individual dataset interval, a further classification was performed based on the obtained values of NSM, i.e., the distance along the transect between two shorelines.Since the value of NSM refers only to the oldest shoreline (in our case 1890) and the most recent shoreline (in our case 2019), the different combinations between the intermediate epochs were performed.
Figure 4 shows the changes of the coastline referred to a specific interval of epochs.
Essentially, the shoreline variation can be characterised into three distinct areas: a northern area characterised by small accretion, a rather stable central area and a southern area characterised by strong erosion with an increasing trend.In addition, in order to evaluate the phenomenon of shoreline retreat and/or advancement, attention was focused on each individual time series.Values referring to the maximum (NSM max ) and mean (NSM mean ) distance between the two shorelines and the mean value (EPR mean ) and uncertainty (EPR unc ) relative to the EPR term were then calculated (Tab.6).
Table 6 shows the erosion trend has been continuous over the years and occurred rather rapidly in the period from 1974 to 1987.Taking into account the erosion and accretion values of the investigated areas, more graphs were built.In particular, Figures 5 and 6 show a boxplot of the EPR mean (end point ratio mean) values classified for the different time series and referred to the erosion and accretion data, respectively.In addition, with reference to the time series, a graph showing the variation over the years of the eroded and accreted areas was built (Fig. 7).In particular, referring to the total historical series evaluated, i.e., the interval 1890-2019, an erosion area of approximately 188,134 m 2 and an accretion area of 113,692 m 2 was detected.The values of shoreline variations are consistent with previous studies conducted in the Mediterranean area and in particular in southern Italy.
In fact, the results of the present research, according to Alberico et al. [46], showed how the phenomenon of coastal erosion in the study area analysed, despite having had an important growth rate in the 1980s or so, has reversed in the last twenty years.
Nevertheless, the phenomenon of the coastline receding inland must be an important element of the evaluation in land management because, according to further studies [47], the study area in the near future could still be transformed, receding inland and modifying the current composition of the dune system.

Conclusions
The paper showed how it is possible to obtain, and in a fast and automatic way, the changes of the coastline that have occurred over the years in a GIS environment.In fact, the rates of change were obtained according to a transect-based approach using the DSAS plugin.One of the advantages of this plugin is that it can evaluate the change of the coastline referred to a historical series of vector data.In the case study, several maps from 1890 until 2019 (georeferenced in a single datum) were considered.By means of these maps, it was possible to obtain the different shorelines, the latter constituting the input vector data of the DSAS plugin.Moreover, to each analysed map, a specific uncertainty value was assigned that takes into account different parametres, such as graphical errors, digitization errors, datum change and georeferencing errors (see Tab. 5); in this way, it was possible to evaluate the measure of the distance between the two shorelines and the related uncertainty on each transect.
The contribution that we wanted to give is to not consider the default parametres of the plugin but to take into account the effective errors associated with the final process of individual shoreline buildings.The attribution of an uncertainty value to each shoreline through a geomatic approach allows us to avoid overestimating uncertainties, especially about the most recent and large-scale maps.
The obtained results lead to a more suitable estimate of the investigated phenomena; however, if on the one hand the analysis conducted does not consider the environmental effects (tides, wind, etc..), on the other hand, these parametres would be insignificant for the time range considered.
From the analysis of shoreline changes, it was demonstrated that in some transects a strong erosion phenomenon is occurred especially in the proximity of the Sele River.In fact, in some transects, the erosion values were more than 200 m.Moreover, analysing the annual rate of change of erosion (Tab.6), it can be discerned that the historical period with the highest erosion rate is the one between 1974 and 1988.
The main origin of the erosion problems is due to the reduction of the solid transport of the Sele River, which has decreased due to the works realized by man along the river (stabilization of upland soils, river water regimentation, etc.) and more generally within the catchment area of the same river.
Considering the high quality of the coastal environment and the tourist industry located in this area, while growth is a benefit for entrepreneurs, the eroded areas not only represent a problem from a socio-economic point of view but also from a naturalistic one; in fact, the erosion of the coast has an impact on both the vegetation (sea lilies and other species of high naturalistic value) and on the pinewood belt with a consequent reduction of the green areas.Indeed, this research intends intended to provide a further contribution to raising awareness among local communities in order to instigate action to mitigate erosion risk.
Finally, from the analysis of the different datasets, it emerged how the cartographic products made available today by platforms such as Google Earth Pro, which are continuously and constantly updated, guarantee that results will be obtained that are comparable to large-scale cartography (e.g.1:2000).Another advantage of the Google Earth Pro platform is the possibility of drawing a series of geometries within the platform itself that can be directly imported into the GIS environment since the maps are framed in a specific reference system.

Fig. 1 .
Fig. 1.Localization of the study area: general framework (a), identification of territory (b), panoramic view (c) and terrestrial photo that shows the erosion (d)

where 2 R a , 2 Rσ 2 G a , 2 Gσ 2 CD a , 2 CDσ 2 D a , 2 D
are related to the scale of representation, are related to the geo- referencing process, are related to the change of datum,

Fig. 3 .
Fig. 3. Shoreline changes in the period 1890-2019: a) categorized by the years; b) categorized by the LRR values

Fig. 5 .
Fig. 5. Boxplot of the EPR mean referring to erosion data

Table 2 .
RMS values obtained on map taking into consideration

Table 3 .
Total RMS error in the map of 1974

Table 4 .
Total RMS error in the map of 1987

Table 5 .
Variance values for different years [m]

Table 6 .
Parametres calculation of the single time series