Mapping of Sicilian Pocket Beaches Land Use/Land Cover with Sentinel-2 Imagery: A Case Study of Messina Province

Abstract

Pocket beaches (PBs) are among the most attractive tourist sites and economic development contributors in coastal areas; however, they are negatively impacted by the combined effects of climate change and anthropogenic activities. Generally, research on PBs is conducted from the beach towards offshore. Studies on the land use/land cover (LULC) of PBs are limited and currently lacking. Such studies deserve more investigation due to the importance of LULC in PBs’ functioning. In this study, supervised classification methods were investigated for LULC mapping of the PBs located in the province of Messina. Sentinel-2B satellite images were analyzed using maximum likelihood (MaL), minimum distance (MiD), mahalanobis distance (MaD) and spectral angle mapper (SAM) classification methods. The study was conducted mainly in order to determine which classification method would be adequate for small scale Sentinel-2 imagery analysis and provide accurate results for the LULC mapping of PBs. In addition, an occurrence-based filter algorithm in conjunction with OpenStreetMap data and Google Earth imagery was used to extract linear features within 500 m of the inland buffer zone of the PBs. The results demonstrate that information on the biophysical parameters, namely surface cover fractions, of the coastal area can be obtained by conducting LULC mapping on Sentinel-2 images.

1. Introduction

Pocket beaches (PBs) are small beaches limited by natural rocky headlands or artificial structures [1,2,3]. These headlands prevent longshore sediment transport and bestow natural protection to the PBs. Depending on the level of their protection, PBs are often considered as closed systems. These enclosed beaches are the most attractive tourist sites and economic development contributors in coastal areas [4,5]; however, the combined effects of climate change and anthropogenic activities put pressure on the coastal landscapes [6] of these ecological niches, consequently impacting their economic value. Within a buffer zone of 500 m from the shoreline, samples representing the main land cover classes of the PBs were selected and extracted from Sentinel-2 imagery. The area of study was chosen in order to evaluate if the adopted approach would be an adequate study model for other different PBs. For each PB, a number of polygon shapefiles for a class of LULC were manually digitized. We separated the samples for the training of the classification model from the samples used for validation and accuracy assessment of the classification results. The validation samples were selected by adopting photointerpretation [7]; Sicilian Corine Land Cover, Sentinel-2B, and Google Earth imagery were acquired for different dates. The part of the image attributed to water also has minimal percentages because the sea was previously cut out; however, for the PBs located near cliffs, a small part of the sea water was taken into consideration.

Generally, research on PBs is conducted from the beach towards offshore areas [8,9,10]. Most research has focused on planview geometry and morphological characteristic analysis of PBs [11,12] or morphohydrodynamics investigations [10,13]. Studies on the land use of PBs [14,15,16] are limited and currently lacking. However, such studies deserve to be considered and conducted due to the importance of LULC in PB function.

Earth observation satellite imaging is the most cost-effective source of data for managing, mapping, and monitoring both offshore [17] and onshore [18,19] coastal areas [20]. It is useful in monitoring anthropogenic coastal activities, such as forestry and agriculture [21], construction [22], mining [23], transport and navigation [24], tourism and recreation [6], fishing and aquaculture [25,26,27], and shipping operations [28]. A significant amount of the world’s population occupies or otherwise visits coastal areas [29], and their activities are increasing and putting pressure on the local environment. It is essential to conduct regular monitoring to identify and evaluate the impact of anthropogenic activities in these areas. Globally, many studies have been carried out in coastal areas using freely available earth observation data [6,30,31], although these studies were conducted on a large scale.

This study is a part of the products of the Pocket Beach Management and Remote Surveillance System (BESS) projects, cofinanced by the European Regional Development Fund, within the Interreg Italy–Malta 2014–2020, and coordinated by the University of Messina in collaboration with the Ministry of Gozo (MG), Euro-Mediterranean Center on Insular Coastal Dynamics—University of Malta (UM), University of Palermo (UNIPA), and the National Institute of Geophysics and Volcanology (INGV).

In this study, supervised classification methods were investigated for LULC mapping of the 24 PBs located in the province of Messina. Sentinel-2B satellite images were analyzed using maximum likelihood (MaL) [32], minimum distance (MiD) [33], mahalanobis distance (MaD) [34] and spectral angle mapper (SAM) [35] classification methods. The study was conducted mainly in order to determine which classification method would be most suitable for small scale Sentinel-2 imagery analysis, aiming for the extraction of information on biophysical parameters, namely surface cover fractions, by conducting LULC mapping, and to provide accuracy results for LULC mapping in the 500 m buffer zone of the PBs.

2. Materials and Methods

2.1. Area of Study

Within the BESS Project, we studied 132 PBs located on the main and minor islands of Sicily, Malta, and Gozo. We collected information about sediments (grain size, color, and composition) and coastal evolution using drone and satellite images. The shape parameters [11,12] of the bay were manually digitized based on 25 cm spatial resolution orthophotos obtained from the 2013 virtual hub dataset (Web Map Service) [36], which was gathered via bathymetry by direct survey or indirect analysis, and a 500 m land use buffer derived from satellite images. Geographic information system (GIS) analysis was applied for data processing and mapping to provide both management and maintenance systems based on the integrated knowledge of these environmental niches [37].

The PBs are closed environments that, considering their degree of confinement, can be considered relicts: a witness of an exceptional storm or a tsunami; the transition of sediments coming from an overflow; the continuous feeding of streams; or the sudden fall of a cliff. The interest of our project was to include an analysis of the potential impact of territorial pressure on the evaluation. In a territory with an estimated density of 192.3 habitants per km² [38] and with significant tourism (14–15 million tourists in Sicily in the pre-COVID-19 era [39]) frequenting the beaches, both “standard” or pocket, the ecological status greatly suffered as a result. For this reason, it is useful to perform an LULC analysis by taking advantage of recent technological development in remote sensing by using the available means—satellite images and commercial software algorithms. To define the methodology, the 24 PBs in the Province of Messina (Figure 1) were chosen because it has the longest coastal perimeter, about 412.65 km (more than ¼ of the whole considered territory), and its coastline is exposed to the waves of two different seas: the Ionian from the east and the Tyrrhenian from the north. Furthermore, the coastal area shows a different state of urbanization because of its variability of both morphology and geodiversity. Alluvial plains derived from the dismantlement of the back chain are characterized by different grades of metamorphic outcropping. From 2001 to 2019 [40], the area of study saw slight demographic changes (

Table 1. Population and area of coastal municipalities of the province of Messina affiliated with the pocket beaches (PBs) studied in the BESS project.

Municipalities	Population in 2001	Population in 2006	Population in 2010	Population in 2013	Population in 2016	Population in 2019	Area (km2)
Gioiosa Marea	7229	7198	7209	7198	7014	6880	26.48
Taormina	10,778	11,026	11,076	11,050	10,909	10,844	13.13
Milazzo	32,083	32,590	32,601	31,882	31,473	30,860	24.7
Lipari	10,556	10,894	11,386	12,500	12,753	12,836	89.71
Malfa	847	872	943	996	956	987	8.73
Patti	13,108	13,391	13,611	13,420	13,347	13,066	50.07

), the population decreasing on the main island and increasing on the Aeolian Islands. During the summer, the population everywhere increases by approximately 30% (based on media information and personnel observations; no punctual scientific data), with the highest peaks on the minor islands.

Along the coastline facing the Ionian Sea, there are 3 PBs in two different locations: 1 PB formed by 2 beaches in a double-curved structure due to the central presence of Isola Bella in the bay of Taormina (from the name of the homonymous municipality) and 1 PB in the bay of Mazzarò. These beaches are limited by the railway line, hotels, restaurants, and beaches of the most important tourist destination of Sicily.

Along the Tyrrhenian coast, 21 PBs were selected: 8 PBs on the main island and 13 on the Aeolian Archipelago. Three beaches, grouped into 1 PB, are located around the Milazzo Peninsula, exposed to the east, north, and west waves, respectively. The whole peninsula is part of a marine reserve, where fishing and bathing activities are strictly regulated. Three beaches, grouped into 1 PB, are located along the northwest side of Cape Tindari and part of a natural protected area; only the westernmost PB is accessible by a very bumpy white road, and the others are accessible only by the sea. Due to this, they maintain a remarkably high degree of naturalness. The two PBs of Capo Calavà are divided by a short promontory—the east one is a very small and wild beach and very difficult to reach, whereas the western one is the longest PB in our study. The longest PB is a wide beach with the presence of the railway line situated above sea level and few activities in the most internal part of the beach itself.

The Aeolian Archipelago contains seven islands: the island of Salina, which has three municipalities—2 PBs from our study are located on the coast of the municipality of Malfa; the other six islands are under the jurisdiction of the municipality of Lipari, which contains 2 PBs along the Volcano Island, 4 PBs along the Lipari Island, and 1 PB along the coast of the Panarea and Filicudi islands, respectively. For each beach, a BESS Project code (SIC01ME01, SIC01ME02 …) was assigned with letters and numbers: “SIC” to distinguish the Sicilian ones from the Malta and Gozo beaches. The first “01 …” represents the number of one of the 132 PBs studied along the whole Sicilian coastline, ME represents the province of Messina, and the second “01” represents the number of one of the 24 PBs studied in this province.

2.2. Data

We analyzed four bands of multispectral Sentinel-2B satellite images at a spatial resolution of 10 m for the land-use mapping of 24 PBs in the coastal municipalities of the province of Messina. Sentinel-2B was placed in orbit on 7 March 2017, together with its twin Sentinel-2A, launched on 23 June 2015; they provide Earth observation data with five days of temporal resolution. Their high spatial resolution of 10, 20, and 60 m allows crucial information to be obtained on LULC in coastal areas. For each PB in Sicily and the Maltese islands, such information can be obtained every 5 days, allowing regular monitoring of the anthropogenic activities on the enclosed beaches. The Sentinel-2B satellite imagery used in this research (

Table 2. Sentinel-2 imagery analyzed for land use on PBs of the province of Messina. Three images acquired on 20 October 2019 at 05:50:29 UTC by Sentinel-2B satellite were analyzed for land use/land cover mapping of the PBs located in six municipalities of the province of Messina. These images cover three tiling grids for each municipality with cloud cover percentage varying from 0.43 to 1.88. Four bands (blue, green, red, and near infrared) with a spatial resolution of 10 m were used.

Municipalities	Sentinel-2 Tiling Grid ID	Cloud Cover (%)
Gioiosa Marea, Lipari, Patti, Malfa	33SVC	0.43
Taormina	33SWB	1.65
Milazzo	33SWC	1.88

) was downloaded for free from the SciHub (https://scihub.copernicus.eu/dhus/#/home) website accessed on 28 May 2020. We downloaded the S2MSI2A product type data—these are level-2 atmospheric corrected satellite images available for the general public.

In addition to Sentinel-2B images, very high-resolution Google Earth pro images acquired on 7 July 2019, 1 August 2019, and 23–24 May 2020 and OpenStreetMap (OSM) data were used to conduct photo interpretation and manual digitization.

OSM data are free and open access volunteered geographic information data and are often used in LULC mapping as an additional source of information for mapping and visualization [41,42]. It is advised to use OSM data in conjunction with recognized data, such as sentinel, because OSM still lacks a formal standard.

2.3. Methodology for Data Processing

In order to contribute to the current effort aiming for the protection of Mediterranean coastal ecosystems, in general, and their PBs, within BESS Project, the latter has been created to take advantage of the latest land surveillance technology and to implement a management system for specific environmental niches in order to find solutions to manage both Maltese and Sicilian PBs.

The satellite image methodological approach adopted in this study is presented in Figure 2. The flow chart shows the steps followed to extract LULC classes and their accuracy assessment. The methodology can be divided into two steps. The first one concerns data download and preparation; the second step is LULC mapping. As we downloaded atmospherically corrected S2MSI2A products, the first step consisted of elaborating on the layer stacking by combining 4 bands into 1 multispectral image and the selection of training samples for each class (

Table 3. Description of classes used for anthropogenic activity mapping on PBs of the province of Messina.

Classes	Description
Shadows	Missing data
Built-up	Residential buildings and asphalt surfaces
Scattered vegetation	Sparse vegetation with 50% grassland covering surface
Vegetation	Vegetated arable lands or heterogeneous agricultural areas
Dense vegetation	Permanent crops, plantation trees, natural or semi-natural forest
Bare land	Nonvegetated agricultural and nonagricultural areas
Bare land (rocks)	Rocky outcroppings
Water	Seawater and swimming pools

) available for each PB.

Both images were processed in the ENvironment for Visualizing Images (ENVI) software, version 5.5 [43]. In addition to LULC, the main linear structures, e.g., streets, were extracted using occurrence-based filtering tools available in ENVI software version 5.5 [44]. Accuracy assessments of LULC classification results were conducted based on validation training samples. The photo interpretation process was carried out by selecting the training sample, whereas the main streets were manually digitized based on OSM data. The results, in percentage, of LULC at 10 m of spatial resolution within 0.5 km inland for each PB revealing indirect anthropogenic activities are presented.

The parameters used to evaluate the quality of LULC results are: (1) the producer’s accuracy, which indicates the degree of confidence for correctly classified areas based on the validation of the training samples of a given class; (2) the user’s accuracy, which is the percentage of correctly classified areas within all samples classified; (3) the overall accuracy (OA), which is the degree of confidence for the correctly classified areas based on all validation training samples; and (4) the Kappa coefficient, which is the agreement of the classification results with the validation training samples.

2.4. Training and Validation Samples Collection

Training and validation samples are the most important input data required for supervised classification of satellite imagery. They determine the number and the names of classes for output results, they allow evaluation of the quality and reliability of classification results, and, last but not least, they are the backbone of any supervised classification work. It is important to ensure that these training and validation samples are sufficient class representatives and cover the whole area of study. In this study, training and validation samples, respectively representing 30% and 45% of the studied buffer zone, were manually selected and digitized on Sentinel-2 images, based on photointerpretation [7] by combining the same satellite images with very high-resolution Google Earth imagery and Sicilian Corine Land Cover version 2011 mapped at 1:10,000 scales [45,46].

Training and validation samples were collected by selecting and digitizing portions of areas of each LULC class on the Sentinel-2 imagery. The process was performed by comparing the class that we intend to create on Sentinel-2 image with Sicilian Corine Land Cover and Google Earth imagery. As some classes for Sicilian Corine Land Cover would have changed, we preferred to divide them into vegetation, dense vegetation, urban vegetation, and scattered vegetation. Such classes are easily identifiable on Sentinel-2 in its false color near infrared composite image as well as on very high-resolution Google Earth imagery. In addition, we sampled data related to bare land, bare land (rocks), built-up, built up (urban), water, and quarry, which are also noticeable on Sentinel-2 in its false color near infrared composite image. To ensure the quality and accuracy of the output results, the training and validation segments were independently created.

2.5. Class Pair Separability (Jeffries–Matusita or Transformed Divergence)

In supervised classification, the training data are used to characterize the classes and to train a predictive model, which is then used to predict land cover in the rest of the image [47]. In order to evaluate the heterogeneity of LULC classes, the separability between two classes was computed. The separability or Jeffries–Matusita distance is an effective metric which evaluates the training samples [48] and measures the quality of how well two classes are distinguished. The separability values range between 0 and 2, with high values indicating better separability, which leads to the precise classification results, while lower values indicate worse separability. For each class pair, we used the value of separability of 1.52 for training and validation samples (Figure 2).

2.6. Brief Description of the Proposed Image Classification Algorithms

The analysis of Sentinel-2 imagery for LULC mapping of PBs was conducted by using 4 classification algorithms.

2.6.1. Maximum Likelihood Classification

The maximum likelihood classifier is a supervised classification method available in different geographic information system software, such as ArcGIS, ENVI, and Erdas Imagine. In our study, we used the maximum likelihood algorithm available in ENVI version 5.5. Maximum likelihood is commonly used for LULC [49] mapping using Sentinel-2 satellite image analysis [50]. During satellite image processing, each pixel is assigned to the class with the highest likelihood. This assignment depends on a threshold value set by the user—in our case, the value threshold was lower, and thus some pixels were unclassified. More details about the maximum likelihood classifier are presented by Ahmad and Quegan [51].

2.6.2. Minimum Distance Classification

The minimum distance classifier available in the ENVI software is a supervised classification algorithm which uses the mean of a region of interest’s shapefiles for each class and calculates the Euclidean distance from each unknown pixel to the mean region of interest of each class [52,53]. In our study, no distance threshold value was specified, which allowed us to perform the classification of the pixels of input data to the nearest class. The classification was conducted in ENVI software version 5.5.

2.6.3. Mahalanobis Distance

The Mahalanobis distance is a direction-sensitive distance classifier. It uses statistics for each class of input data. Mahalanobis distance is similar to the maximum likelihood classification, but it assumes that all class covariances are equal, and therefore, it is a faster method [52]. During the processing in the ENVI software, no specific distance threshold value was used, and this allowed the algorithm to classify all the pixels to the closest training data.

2.6.4. Spectral Angle Mapper

The spectral angle mapper is a physically based spectral classification that uses an n- (number of spectral bands) dimensional angle to match pixels to training data. It is a supervised classification machine learning method that is used to determine the spectral similarity between two spectra (a given spectra image and reference spectra) by calculating the angle between the spectra and treating them as regions of interest in a dimensional space equal to the number of bands [54,55]. The pixels are classified to the classes of the smallest angle between the training mean region of interest and each pixel region of interest in n-dimensional space. The classification was performed in ENVI software version 5.5 with reflectance data calibrated to the shapefile of the area of interest, which allowed us to perform a classification relatively insensitive to illumination and albedo effects.

2.7. Occurrence-Based Filtering

We applied texture filtering to the Sentinel-2B imagery using an occurrence-based filter available in ENVI software version 5.5. For more than three decades, texture has been of great interest to remote sensing analysts [56]. For each pixel, texture analysis generates spatial variability measurements, which produce data that is easily classifiable at the pixel level [55]. This information helps to identify linear structures composed of contiguous pixels with the same tone values. The first-order occurrence of measured textures, namely mean and variance, were computed. This allowed for the distinction of streets and roads from other LULC classes by measuring the distribution of gray tones around the mean (Figure 3). The measure of mean and the variance were computed with 3 × 3 kernels as follows (Equations (1) and (2)) [57]:

Mean:

$$M = {\Sigma }_i{\Sigma }_{j}X\left(i,j\right)$$

(1)

$$\mathrm{Variance}:\hspace{1em}\hspace{1em}\hspace{1em} V={\Sigma }_i{\Sigma }_j\left(i-M\right)²X\left(i,j\right)\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}$$

(2)

where M is the mean, V is the variance, i is gray tone, and X is a matrix of gray tone of the image.

The results of the occurrence-based filter were then compared with OSM data to identify and digitize the main street in a 500 m buffer zone from the PB. The process followed a photo-interpretation process which allowed us to classify and manually digitize the streets.

3. Results

The results obtained show that Sentinel-2 imagery provides useful LULC information that could also be used to evaluate and better understand the impact of anthropogenic activities in small coastal ecosystems.

In our study, we focused on these 24 PBs (

Table 4. Descriptions of 24 PBs of the province of Messina.

Name of the Beaches	Municipality/Island	Codes of the Beaches	Coordinates (Decimal Degrees)
			Latitudes	Longitudes
South Isola Bella	Taormina	SIC01ME01	37.850035	15.297601
North Isola Bella	Taormina	SIC01ME02	37.852593	15.300546
Mazzarò	Taormina	SIC02ME03	37.855315	15.301613
East Milazzo	Milazzo	SIC03ME04	38.263882	15.243056
North Milazzo	Milazzo	SIC03ME05	38.269745	15.237303
West Milazzo	Milazzo	SIC03ME06	38.264215	15.236534
East Tindari	Patti	SIC04ME07	38.151473	15.040707
Central Tindari	Patti	SIC04ME08	38.151205	15.037139
West Tindari	Patti	SIC04ME09	38.150371	15.03108
East Capo Calavà	Gioiosa Marea	SIC05ME10	38.190424	14.920996
West Capo Calavà	Gioiosa Marea	SIC06ME11	38.188102	14.910634
Punta dell’Asino	Lipari/Vulcano	SIC07ME12	38.370453	14.998676
Punta Bandiera	Lipari/Vulcano	SIC08ME13	38.373471	15.003052
La Forbice	Lipari/Lipari	SIC09ME14	38.450786	14.960705
Pignataro di Fuori	Lipari/Lipari	SIC10ME15	38.477989	14.972756
Sabbie Bianche	Lipari/Lipari	SIC11ME16	38.498898	14.961744
Lido Blu	Lipari/Lipari	SIC12ME17	38.501816	14.962953
Punta Scario	Malfa/Salina	SIC13ME18	38.582392	14.834378
Pollara	Malfa/Salina	SIC14ME19	38.580203	14.807179
Zimmaro	Lipari/Panarea	SIC15ME20	38.628785	15.065786
West Preistorico	Lipari/Panarea	SIC16ME21	38.625762	15.061491
Central Preistorico	Lipari/Panarea	SIC16ME22	38.62564	15.063061
East Preistorico	Lipari/Panarea	SIC16ME23	38.625176	15.063807
Le Punte	Lipari/Filicudi	SIC17ME24	38.555222	14.58378

) because they seemed to be a good representation of Mediterranean PBs located in different coastal situations; otherwise, similar results have also been obtained for the other 132 PBs investigated within the BESS Project.

The results obtained from the multispectral Sentinel-2B satellite image analysis are presented in maps and tables. The maps represent the location of LULC classes for the main Island of Sicily (Figure 4) and for the minor Islands of the Aeolian Archipelago (Figure 5), whereas the table contains statistical analysis and accuracy assessment results (

Table 5. Overall land use/land cover classification results and analysis obtained for the PBs located on the coast of Messina province.

Codes of the PBs	Classes	Area		Producer’s Accuracy (%)	User’s Accuracy (%)	Overall Classification Accuracy (%)	Kappa Coefficient
		m2	(%)
SIC01ME01, SIC01ME02, SIC02ME03	Shadows	53,700	6.05	100	27.78	96.98	0.93
	Built-up	426,100	48	100	98.62
	Vegetation	407,800	45.94	95.70	100
SIC03ME04, SIC03ME05, SIC03ME06	Shadows	61,500	6.27	100	76.47	89.07	0.85
	Built-up	119,600	12.21	86.36	96.94
	Vegetation	399,700	40.8	96.35	84.62
	Dense vegetation	189,400	19.33	81.48	97.06
	Bare land	209,300	21.36	95.12	75.00
SIC04ME07, SIC04ME08, SIC04ME09	Shadows	59,600	5.45	56.1	96.72	87.06	0.79
	Scattered vegetation	496,200	45.42	91.01	93.08
	Dense vegetation	256,600	23.49	99.32	82.21
	Bare land	184,500	16.89	68.25	89.84
	Bare land (rocks)	95,400	8.73	83.56	46.56
SIC05ME10, SIC06ME11	Shadows	94,100	7.37	98.08	100	89.63	0.86
	Built-up	428,800	33.62	92.24	99.55
	Scattered vegetation	214,200	16.79	74.71	65.33
	Vegetation	322,000	25.24	88.10	81.70
	Dense vegetation	112,700	8.83	89.39	100
	Bare land	103,500	8.11	100	98.54
SIC09ME14	Built-up	167,100	36.64	92.83	89.80	92.63	0.89
	Vegetation	111,600	24.47	95.73	94.58
	Bare land	173,500	38.04	94.32	93.79
	Water	3800	0.83	67.65	100
SIC10ME15	Vegetation	333,800	71.93	100	83.55	94.81	0.87
	Bare land	130,200	28.06	92.96	100
SIC11ME16, SIC12ME17	Built-up	7088	26.77	97.74	98.30	98.76	0.97
	Scattered vegetation	424,400	53.44	99.15	99.09
	Vegetation	157,100	19.78	100	97.17
SIC07ME12, SIC08ME1	Built-up	82,400	8.58	82.83	93.18	96.08	0.93
	Vegetation	421,900	43.93	97.51	99.49
	Bare land	448,400	46.69	100	94.67
	Water	7600	0.79	66.67	100
SIC15ME20, SIC16ME21, SIC16ME22	Built-up	83,500	12.44	98.42	95.04	92.80	0.88
	Vegetation	350,200	52.21	98.31	84.30
	Bare land	237,000	35.33	86.44	99.14
SIC17ME24	Shadows	19,800	3.47	98.63	100	90.45	0.82
	Built-up	99,000	17.35	94.59	100
	Scattered vegetation	338,600	59.35	87.84	98.26
	Dense vegetation	113,100	19.82	95.06	70
SIC13ME18	Built-up	226,800	51.53	96.99	99.56	97.76	0.95
	Vegetation	176,900	40.19	100	92
	Water	36,400	8.27	97.97	99.32
SIC14ME19	Shadows	153,400	23.92	97.40	100	94.89	0.92
	Built-up	158,500	24.71	92.48	72.32
	Vegetation	223,400	34.84	99.08	98.19
	Bare land	38,700	6.03	63.43	96.52
	Water	67,200	10.48	100	95.63

). The results cover 24 PBs located in six coastal municipalities of the Messina province, including classification maps and statistical parameters, and the quality and extent of classes. We noticed that in the 500 m buffer zone of the PBs, information on the biophysical parameters, namely surface cover fractions, can be obtained by conducting LULC mapping on Sentinel-2 images. For most of the PBs, a large area is occupied by vegetation, built-up areas, and bare land. For each PB, more detailed descriptions are presented in the Supplementary Material (Figures S1–S12 and Tables S1–S12). Information on the biophysical parameters represented by the surface cover fraction indicates the main anthropogenic activity conducted in nearby each PB.

4. Discussion

Currently, open available geospatial data are being used for land use/land cover (LULC) change analysis of coastal areas [6,30,31,32]. In this study, we conducted a LULC mapping analysis of the PBs located on the coast of Messina province. The area of study was chosen in order to evaluate if the adopted approach would be an adequate study model for other different PBs. Our research methodology evaluates four classification methods with the aim of determining which one could be used for an accurate LULC mapping of the PBs. Within a buffer zone of 500 m from the shoreline, samples representing the main land cover classes of the PBs were selected and extracted from Sentinel-2 imagery. They were used during image classification procedures and eventually for extracting information on the biophysical parameters, namely surface cover fractions. For most of the PBs, a large area is occupied by vegetation, built-up areas, and bare land. The identification of surface cover fractions indicates the anthropogenic activity conducted nearby each PB.

4.1. Effect of Class Combination and Topography on Classification Accuracy and Precision

An evaluation of the effect of class combination and topography on the classification accuracy and precision was performed. For each PB, the class pair separability was computed. Separability values indicate how LULC classes can be distinguished on satellite images. The values of separability of two combined classes vary between 0 and 2, with higher values indicating the best separability of classes leading to more precise classification results. The area of study presents topographic variations resulting in shadows near the cliffs or hill slopes. An interesting issue is the presence of the shadows. Some images may have surfaces that are not immediately detectable, which are classified as shadows. In most of the images analyzed, the percentage of shadow is less than 8%; in six images, it is practically absent; and in the last image, it reaches almost 24%. When the images are taken at the same time, the issue of the shadows could be resolved by integrating freely available multitemporal satellite images with moderate spatial resolution, such as those from the Landsat program. In any case, the separability of different LULC features was guaranteed, and it was taken into account in the application of classification methods. These shadows reduce the separability of LULC classes, eventually reducing the precision of the results. In this study, we used the value of separability of 1.52 for training and validation samples.

Figure 6 shows the difference between the performances of the four different classification methods evaluated. For each method, overall classification accuracy and kappa coefficient are presented. We noticed that MaL performs better than MiD, MaD, and SAM classifiers; however, for the classification of only two classes, SAM performs better than other classifiers. We noticed that the range of overall classification accuracy is predicted by the number of classes as well as the values of the class pair separability. For a small number of classes and a high value of separability, the MaL classifier demonstrates the best ability for mapping. However, for only two LULC classes, the SAM classifier performs better with high overall classification accuracy and kappa coefficient values, due to the smallest angle between the classes and spectral bands.

4.2. Limitation and Suitability of Conducting Supervised Classification Procedures

Phiri et al. [58] reviewed the contribution of Sentinel-2 to land cover/use classification and noticed that Sentinel-2 has a positive impact on land cover/use monitoring, specifically in the monitoring of crops, forests, urban areas, and water resources; Xu et al. [59] used Sentinel-2 multispectral imagery to extract a high-precision urban impervious surface. In our case study, we noticed that for LULC mapping activities conducted on a small area, the ability to differentiate all land cover objects is limited by spatial resolution and the presence of missing data pixels in mountainous landscapes. The results of the current study demonstrate that information on the biophysical parameters, namely surface cover fractions, of the coastal area can be obtained by conducting LULC mapping on Sentinel-2 images. The proposed approach contributes to coastal research studies, especially research on the LULC of PBs. We demonstrate that the supervised classification procedures rely mainly on three components, namely classifier algorithms, training and validation samples, and dataset imagery. While recent improvements have been made to machine learning classifier algorithms and several datasets of imagery are available, more effort should be made for the collection and the preparation of training and validation samples.

5. Conclusions

In this study, satellite images of Sentinel-2B were processed for LULC mapping of the PBs located on the coast of Messina province. The area of study was chosen in order to evaluate if the adopted approach would be an adequate study model for other different PBs. The obtained maps of the LULC allowed us to determine the location of different anthropogenic activities conducted near the PBs. They are essentially dominated by agriculture and forestry. The construction activities dominating built-up areas and urban areas represent the second most prevalent human activities. Considering the applied methodology, information on the biophysical parameters, namely surface cover fractions, can be retrieved from Sentinel-2 images. Indeed, due to the global coverage of Sentinel-2, its images can be used for the LULC mapping of other PBs. Such information will effectively improve the monitoring of coastal areas.