Geostatistics and Structure from Motion Techniques for Coastal Pollution Assessment along the Policoro Coast (Southern Italy)

: The coastal system, with its physical chemical and biological components, is the place where the lithosphere, atmosphere and hydrosphere intersect and interact, and in which human activity has a considerable impact on the balance of the whole ecosystem. The aim of this study is to assess the health of the coastal environment in a natural protected area in Policoro (Italy). Sand samples have been collected and analysed to detect the presence of environmental contaminants and pollutants (heavy metals). A photogrammetric survey was carried out using an unmanned aerial vehicle (UAV), and an updated orthophoto of the area was obtained using the Structure from Motion (SfM) processing technique. Geostatistical techniques have been used to assess the distribution of the contaminants in the study area. information on the spatial distribution of the contaminants investigated in the sampling area. Statistical analysis was performed using the statistical tool of Microsoft O ﬃ ce Excel 2019.


Introduction
Environmental pollution is defined as the presence of contaminants, in any chemical and physical state, in the natural environment [1]. There are two main types of pollution: Natural pollution, caused by natural processes, e.g., CO 2 and sulfur emissions caused by volcanic eruptions; and anthropic pollution caused by human activities. The last one is due to industrial processes, incinerator emissions, thermoelectric power plants, oil refineries, iron and plastic production processes, chemical fertilizers, radioactive waste, etc. The global and constant growth of these activities over the last century, in conjunction with a demographical increase, has contributed to the exposure of new areas, both inland and coastal, to pollution risk. Modeling environmental data is a fundamental step to understanding phenomena characterized by variability [2]. An environmental variable is a quantitative or descriptive measure of the natural environment [3]. In order to understand the environmental status of a site, sampling has to be carried out in the study area. In this study, heavy metals were chosen as analytes because they play an important role in the ecosystem. They don't degrade and have a high toxicity, both for humans and for all living species, because of the bioaccumulation process. This physical-chemical process allows metal to bind to cellular structures, hindering the performance of the organism's vital function. [4][5][6][7]. The collected samples are useful for obtaining measures which are representative of the studied phenomenon. Geostatistics is a class of statistics used to analyse values associated with spatial phenomena, useful for describing spatial patterns and interpolating values. The theory of Geostatistics, first elaborated by Matheron [8] and then disseminated by Isaaks and Srivastava [9] Webster and Oliver [10], Goovaerts [11] allowed the development of a theoretical system for the study of a soil variability, spatial interpolation and the choice of the optimal sampling scheme. At present, Geostatistics is widely used to quantify the spatial distribution of environmental variables. Starting from a dataset of samples, Geostatistics allows the prediction of the value of a study parameter in unknown positions and offers the possibility of converting punctual and discrete data into areal and continuous data (gridding). The geostatistical technique that allows this transformation is the spatial interpolation, which can be global or local and exact or approximated. In the global interpolation, all the data domain is used for the estimation of a value in a given point, and the variation of a single input value affects the entire output surface; in the local interpolation, only values which are proximal to a given point, are used to estimate the value in that point [12]. In this case, the variation of a single input value only affects the output in the local region. Spatial interpolation is defined as exact when the interpolated surface passes through the measurement points, while it is approximate when it does not pass through the measurement points, that is, when margins of uncertainty are expected on the sampled data. A spatial interpolation method is the deterministic one, which uses a deterministic function to predict value at an unknown point. The most important and most used deterministic methods are the Inverse Distance Weighted (IDW) and the Spline. In the IDW method [13,14], the influence of input data is inversely proportional to the distance between the input datapoint and the point at which the value has to be predicted. Moreover, the output surfaces are quite realistic. The Spline method [15] performs an estimate of the values through a mathematical function (applied to the proximal data) that minimizes the curvature by creating a very smooth surface. This methodology is widely used when dealing with parameters that have a low variability (concentration of pollutants, temperatures, piezometric surfaces). This study is intended as a useful preliminary study to test this methodology before applying it on a larger area. The proposed area of study has been chosen because of the presence of a drainage canal, which might convey polluting chemicals derived from production activities in the area.

Production Activities in the Study Area
The production activities found on the site mostly concern zootechnical activities and agriculture. Within the protected area, the livestock consistency is very high (about 1200 cattle and 1300 sheep and goats are present). The bovine breeding is of an intensive type and intensive agriculture is carried out both within the protected area and in the immediate vicinity. Arboreal (orchards) and herbaceous crops (fodder) are widespread. These activities have a strong influence on the conservation of priority species and habitats present on the site, as they contribute significantly to the alteration of the delicate balances that regulate the ecosystems. In particular, reference is made to the following issues: • Supplement of nutrients: The types of breeding and agriculture described above significantly increase the contribution of organic material which, because of the washout, reaches the drainage canals, causing a considerable increase in the concentration of nitrates in the soil; • Use of pesticides: The massive use of pesticides affects the quality of the water in the drainage channels. Particularly in summer, as a result of the excessive evaporation, these compounds can reach high concentrations, causing severe damage to insectivorous species after bioaccumulation.

Flora and Fauna
The presence of relict populations of species linked to an extinct hygrophilous environment and typical species of consolidated sand dunes is significant. Twenty-seven species of community interest fauna have been identified: Among these, 21 species are included in the Council Directive 79/409 EEC on the protection of wild birds, while two species of mammals, two species of reptiles and two species of invertebrates are included in the Council Directive 92/43 EEC on the conservation of natural habitats and of wild fauna and flora. The list of birds includes many migratory species, whose presence in the area confirms the importance of the area as a stop-over site during migrations. Moreover, along the neritic zone there is a consistent transit of pelagic migratory species, observed in large numbers comprising up to 500 individuals, during the post-reproductive migration. The small temporary puddles that form along the shoreline as a result of storm surges are also used by several bird migratory species. Water birds are also important: They colonize the reedbeds at the mouth of the Sinni river and the hygrophilous formations located along the drainage canals. These biotopes are very important for the wintering of some migratory species, which colonize the reedbeds near the canals. [16].

Geological and Geomorphological Setting
The area covered by this study is located in a Site of Community Importance and in a Special Protected Area (SIC-ZPS IT9220055), named "Bosco Pantano di Policoro and Ionic Coast Foce Sinni" on the Ionian coast of the Basilicata Region ( Figure 1) in a section bordered to the south by the mouth of the river Sinni and to the north by the mouth of the Agri river.
Geosciences 2020, 10, x FOR PEER REVIEW 3 of 11 small temporary puddles that form along the shoreline as a result of storm surges are also used by several bird migratory species. Water birds are also important: they colonize the reedbeds at the mouth of the Sinni river and the hygrophilous formations located along the drainage canals. These biotopes are very important for the wintering of some migratory species, which colonize the reedbeds near the canals. [16].

Geological and Geomorphological Setting
The area covered by this study is located in a Site of Community Importance and in a Special Protected Area (SIC-ZPS IT9220055), named "Bosco Pantano di Policoro and Ionic Coast Foce Sinni" on the Ionian coast of the Basilicata Region ( Figure 1) in a section bordered to the south by the mouth of the river Sinni and to the north by the mouth of the Agri river. The study area is oriented in a SW-NE direction along the Ionian coast of Basilicata, in the western sector of the Gulf of Taranto, proximal to the area between the coastal plain of Metaponto and the alluvial areas of the main Lucanian rivers. From a geological point of view (Figure 2), the presence of Upper Pleistocene marine terraced deposits [17][18][19], located at the intersection of the front of the chain and the Ionian coast of Lucania, is one of the main features. Holocene deposits (sandy, sandy-gravelly and sandy-silty deposits) form the coastal plain. They are located in the area between the most recent terraced marine deposits and the innermost dune belt, and consist of a sedimentary succession of sands, clayey sands and gravels, with thickness variable from 20 to 40 m, deposited on the local basement represented by Argille subappennine (Upper Pliocene-Lower Pleistocene). The coastal dunes consist of fine sand, varying in color from gray to yellowish, and form a strip of cordons, parallel to the current shoreline, towards Marina di Ginosa, where the dunes reach heights of up to 13 m. The stretch of coast (Idrovora of Policoro-Ginosa Marina) in which the study area is located, is characterized by a beach composed of coarse and dark sand towards the Policoro idrovora, growing gradually finer and clearer towards Ginosa Marina. The dune belt gradually increases in extent and height towards Ginosa Marina, bordered towards the inland by dark, clayey-silty deposits. Sedimentological, mineraolgical and anemometric studies, conducted in an area that includes the one in question, identify a prevalent direction of sediment dispersion towards NE [20,21]. The study area is oriented in a SW-NE direction along the Ionian coast of Basilicata, in the western sector of the Gulf of Taranto, proximal to the area between the coastal plain of Metaponto and the alluvial areas of the main Lucanian rivers. From a geological point of view (Figure 2), the presence of Upper Pleistocene marine terraced deposits [17][18][19], located at the intersection of the front of the chain and the Ionian coast of Lucania, is one of the main features. Holocene deposits (sandy, sandy-gravelly and sandy-silty deposits) form the coastal plain. They are located in the area between the most recent terraced marine deposits and the innermost dune belt, and consist of a sedimentary succession of sands, clayey sands and gravels, with thickness variable from 20 to 40 m, deposited on the local basement represented by Argille subappennine (Upper Pliocene-Lower Pleistocene). The coastal dunes consist of fine sand, varying in color from gray to yellowish, and form a strip of cordons, parallel to the current shoreline, towards Marina di Ginosa, where the dunes reach heights of up to 13 m. The stretch of coast (Idrovora of Policoro-Ginosa Marina) in which the study area is located, is characterized by a beach composed of coarse and dark sand towards the Policoro idrovora, growing gradually finer and clearer towards Ginosa Marina. The dune belt gradually increases in extent and height towards Ginosa Marina, bordered towards the inland by dark, clayey-silty deposits. Sedimentological, mineraolgical and anemometric studies, conducted in an area that includes the one in question, identify a prevalent direction of sediment dispersion towards NE [20,21].

Materials and Methods
In order to update the cartographic and geomorphological context, aerial photogrammetric surveys were carried out using an unmanned aerial vehicle (UAV) Phantom 4 Pro, equipped with a 20 MP camera, CMOS 1″ sensor, FOV 84° 8.8/24 mm (35 mm eq.) f/2,8-f/11, ISO 100-3200, flying at a 45 m quote, with a corresponding Ground Sampling Distance (GSD) of 2 cm. The flights were carried out on 16 July 2019, between 12:30 and 13:00, in a range of time without the influence of high or low tide, in order to obtain a reliable and updated morpho-topographic model of the area, unaffected by tidal phenomena.

Materials and Methods
In order to update the cartographic and geomorphological context, aerial photogrammetric surveys were carried out using an unmanned aerial vehicle (UAV) Phantom 4 Pro, equipped with a 20 MP camera, CMOS 1" sensor, FOV 84 • 8.8/24 mm (35 mm eq.) f/2,8-f/11, ISO 100-3200, flying at a 45 m quote, with a corresponding Ground Sampling Distance (GSD) of 2 cm. The flights were carried out on 16 July 2019, between 12:30 and 13:00, in a range of time without the influence of high or low tide, in order to obtain a reliable and updated morpho-topographic model of the area, unaffected by tidal phenomena.
The entire survey of the study area has produced 395 images that were processed in a Linux environment using MicMac, an open source software which was developed by the French Geographic Institute (IGN). Its modular structure allowed the precise regulation of all the obligatory and optional parameters to process the dataset with the Structure from Motion (SfM) technique [22]. The SfM is a technique that allows the obtaining of 3D models using a dataset of 2D images, generating ortophotos and models based on point clouds, in which each point has spatial and RGB information. The workflow used is summarized in Figure 3. The entire survey of the study area has produced 395 images that were processed in a Linux environment using MicMac, an open source software which was developed by the French Geographic Institute (IGN). Its modular structure allowed the precise regulation of all the obligatory and optional parameters to process the dataset with the Structure from Motion (SfM) technique [22]. The SfM is a technique that allows the obtaining of 3D models using a dataset of 2D images, generating ortophotos and models based on point clouds, in which each point has spatial and RGB information. The workflow used is summarized in Figure 3. . Structure from Motion processing workflow.The aerial photogrammetry was useful to identify morpho-topographic features that are not clearly visible from the ground level. Moreover, it was a time-efficient method to produce cartographic outputs usable in GIS. The most important output obtained was the georeferenced orthophoto which was imported in the QGIS open-source software for the environmental characterization phase. This methodology allows a choice of sampling points using an updated map of the study area, which is affected by frequent coastline and morphological variations due to storm surges. Two main strategies can be followed to choose the location of the sampling points:.

•
Convenience sampling: the choice is based on an examination of previous studies. This choice is preferred in complex sites when the historical information available allows for the prediction of the location of the most vulnerable areas and the most probable sources of contamination; • Systematic sampling: the choice of the sampling points is made on the basis of a random or statistical criterion, based on a predefined or random grid. This choice is preferred when the size of the area, or the lack of historical information on the site, make it difficult to obtain a preliminary environmental characterization.
Since there are no historical data available in the study area, the strategy of systematic location of the sampling points was chosen. A total grid, with lateral dimensions of 150 × 60 m (9000 m 2 ), and eight internal mesh having side dimensions of 37.5 × 30 m (1125 m 2 ), was defined ( Figure 4). In each mesh, the sampling point choice was also based on the presence of environmental matrices (channel water and marine water) not subject to this study (Table 1, Figure 5). The aerial photogrammetry was useful to identify morpho-topographic features that are not clearly visible from the ground level. Moreover, it was a time-efficient method to produce cartographic outputs usable in GIS. The most important output obtained was the georeferenced orthophoto which was imported in the QGIS open-source software for the environmental characterization phase. This methodology allows a choice of sampling points using an updated map of the study area, which is affected by frequent coastline and morphological variations due to storm surges. Two main strategies can be followed to choose the location of the sampling points: • Convenience sampling: The choice is based on an examination of previous studies. This choice is preferred in complex sites when the historical information available allows for the prediction of the location of the most vulnerable areas and the most probable sources of contamination; • Systematic sampling: The choice of the sampling points is made on the basis of a random or statistical criterion, based on a predefined or random grid. This choice is preferred when the size of the area, or the lack of historical information on the site, make it difficult to obtain a preliminary environmental characterization.
Since there are no historical data available in the study area, the strategy of systematic location of the sampling points was chosen. A total grid, with lateral dimensions of 150 × 60 m (9000 m 2 ), and eight internal mesh having side dimensions of 37.5 × 30 m (1125 m 2 ), was defined ( Figure 4). In each mesh, the sampling point choice was also based on the presence of environmental matrices (channel water and marine water) not subject to this study (Table 1, Figure 5). Geosciences 2020, 10, x FOR PEER REVIEW 6 of 11    The surface sand samples were taken on 23 July 2019 from 11:00 to 11:20 in a time range not influenced by the daily tide excursion. The samples were taken using a stainless-steel sampler which did not affect the physical-chemical characteristics of the sand; each sample was stored and transported in a sterile LDPE sample bag with an ID label.
The methodology used for chemical analyses is compliant with the Italian provision of law D.M. 13th September, 1999, n • 248, 21st October 1999. In order to define Cr VI, the methodology IRSA-CNR Qd. 64 vol. 3 n. 16/1986 wasapplied. About 500 g of sand were taken from each sampling point, and stored in the refrigerator at a temperature of 4 • C as they were transported to the Serveco S.r.l. analysis laboratory. Considering the human activities in the study site and the ecosystemical importance of the protected coastal area, the analysis focused on the concentration of these heavy chemical elements: Sb, As, Be, Cd, Co, Cr, Cr VI, Hg, Ni, Cu, Pb, Se, Sn, Tl, V, Zn. In the data spatial exploration phase, decided to two types of interpolation (IDW and Spline) were performed in order to obtain information on the spatial distribution of the contaminants investigated in the sampling area. Statistical analysis was performed using the statistical tool of Microsoft Office Excel 2019.     According to the correlation matrix between the analytes (Table 2), a strong correlation can be observed between some analytes (Sb-Cr, Sb-Ni, Sb-V, Sb-Zn, As-Co, As-Cu, Co-Cu, Cr-Ni, Cr-Zn, Cr-V, Ni-V, Ni-Zn, V-Zn). This correlation highlights a strong tendency to assume the same spatial pattern. The significance analysis of the correlations between the analytes was carried out considering 14 degrees of freedom of the system and obtained the following significance levels: 0.01 bilateral hypothesis (Sb-Ni, Sb-V, Sb-Zn, Ni-V, Ni-Zn, V-Zn), 0.05 bilateral hypothesis (Sb-Cr, Cr-Zn, Cr- Ni ± 11; Cu ± 1.9; V ± Zn ± 6.2).

Results and discussion
Geosciences 2020, 10, 28 8 of 11 The metal Sb appears to be present in greater quantities in samples C1 and C2, where it reaches~7 mg/kg, while it assumes lower values in the other samples, up to a minimum value of~1 mg/kg. The highest values of the as metal are found in the C7 and C8 samples while they are lower than 10 mg/kg in the other six samples.
The metal Co results have concentrations higher than 5 mg/kg in all the samples, and in particular reaches a maximum of 12-13 mg/kg in samples C1 and C7. Cr has concentrations with high spatial variability and reaches very high values. The points where the highest values are reached are in samples C1 and C2, where they exceed 500 mg/kg.
The Ni metal assumes a spatial variability similar to that of Cr, and also, in this case, high values are reached considering the type of metal. Values exceed 60 mg/kg in all points with maximum values close to 150 mg/kg in sample C1. Cu has values in the range 5-14 mg/kg with points C3, C4 and C6 less affected by this metal.
The concentration of the metal V is in the range of 9-27 mg/kg and reaches the maximum values at points C1 and C2, while the minimum values are found at points C3 and C4. Zn metal concentrations are more constant at all points and values between 25-52 mg/kg are measured.

Sb
As The correlation matrix between the sampling points (Table 3), which represents the total concentration of all heavy metals, shows a very high correlation index among all the sampling points. Only C5 seems to be almost unrelated to C1 and C2 (0.1-0.05) and this could be explained by its being the furthest point from the ones where the highest contamination was found, meaning it is not affected by mechanisms of natural or anthropogenic dispersion, due to the presence of the drainage channel.
Spline and IDW interpolation algorithms were performed. The area actually interpolated for both cases was 5500 m 2 because of the presence of environmental matrices, which are not the subject of this study, and because of the interpolation limits. Figure 7 shows the application of the two interpolation algorithms to the study area. In particular, for each sampling point, the concentrations of all the heavy metals analyzed were summed, and the results were interpolated in a GIS environment, in order to obtain a map of the area most affected by the contamination. The figure shows how the two interpolation algorithms are consistent in showing the presence of the most polluted area between points C1 and C2 (in red), while, moving towards the North and East, there is a gradual decrease in the contamination. Spline and IDW interpolation algorithms were performed. The area actually interpolated for both cases was 5500 m 2 because of the presence of environmental matrices, which are not the subject of this study, and because of the interpolation limits. Figure 7 shows the application of the two interpolation algorithms to the study area. In particular, for each sampling point, the concentrations of all the heavy metals analyzed were summed, and the results were interpolated in a GIS environment, in order to obtain a map of the area most affected by the contamination. The figure shows how the two interpolation algorithms are consistent in showing the presence of the most polluted area between points C1 and C2 (in red), while, moving towards the North and East, there is a gradual decrease in the contamination.  Table 4 describes the two algorithms interpolation errors, related to the concentrations detected on site. In the first column, the IDs of the sampling points are reported; in the second column, the total concentrations of all the heavy metals for each sampling point; in the third and fourth column, the concentrations calculated respectively with the algorithm are reported, IDW and Spline; in the fifth and sixth columns, the interpolation errors of the two algorithms used were calculated and reported. These errors were calculated by the difference between the real values and the values calculated by the two interpolations for each sampling point.   Table 4 describes the two algorithms interpolation errors, related to the concentrations detected on site. In the first column, the IDs of the sampling points are reported; in the second column, the total concentrations of all the heavy metals for each sampling point; in the third and fourth column, the concentrations calculated respectively with the algorithm are reported, IDW and Spline; in the fifth and sixth columns, the interpolation errors of the two algorithms used were calculated and reported. These errors were calculated by the difference between the real values and the values calculated by the two interpolations for each sampling point. In the IDW method, the average error of interpolation, compared to the real concentration measured, was −0.19 mg/kg, while in the Spline methodology the error was 0.19 mg/kg ( Table 4). The standard deviation calculated for the Spline methodology was 1.24, while it was 0.17 for the IDW method.
Both interpolations show a more evident red area, in which samples C1 and C2 are located, which is the most contaminated by heavy metals. The contamination decreases moving away from this area, reaching lower concentrations in the area subtended by samples C3, C4, C5, C6. Using the interpolation techniques described for the mapping of real concentrations in the investigated area, it is possible to calculate the potential concentrations of the analytes at points located within the same area. These "expected" concentrations can be deduced by using GIS software to create iso-concentration curves of the studied analytes.

Conclusions
The processes that influence coastal environment dynamics and the strong temporal variability that characterize them were the starting points for the development of this work. The aim of this study was to assess the health of the coastal environment in a natural protected area in Policoro (Italy). A photogrammetric survey of the study area was carried out using a UAV and processing the images using the SfM technique, obtaining an updated orthophoto in order to make an accurate choice of environmental sampling points. Eight samples have been collected in the area and analysed in the laboratory. Geostatistical techniques, supported by GIS techniques, have been applied to estimate the spatial distribution trend of heavy metals in the area, by interpolating punctual and discrete measurements derived from the sampling. This work highlighted some critical environmental issues in this coastal area, which could be due to human activities in the area.
In fact, the concentration of heavy metals is higher in proximity to the drainage channel (Figure 7), suggesting that the pollutants have been transported by the channel itself, while it's lower in proximity to the shoreline, probably because of the longshore current dispersion action.
Testing this methodology in a small area was a useful preliminary study, and the results can be considered as a starting point, before applying the SfM and geostatistical techniques in a wider area to evaluate the applicability of this method to assess the dispersion of pollutants at a different scale. Considering the environmental and ecological importance of the site, these first results could be very interesting for further studies.