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

: 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 classiﬁcation 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) classiﬁcation methods. The study was conducted mainly in order to determine which classiﬁcation 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 ﬁlter 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.


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 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 2 [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 Land 2021, 10, 678 3 of 20 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 1 4 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), 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.
Land 2021, 10, x FOR PEER REVIEW 3 of 20 habitants per km 2 [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), 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.

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) 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. 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 (%) 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.

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) 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. Both images were processed in the ENvironment for Visualizing Images (ENVI) sof ware, version 5.5 [43]. In addition to LULC, the main linear structures, e.g., streets, wer extracted using occurrence-based filtering tools available in ENVI software version 5 [44]. Accuracy assessments of LULC classification results were conducted based on val dation training samples. The photo interpretation process was carried out by selecting th training sample, whereas the main streets were manually digitized based on OSM dat The results, in percentage, of LULC at 10 m of spatial resolution within 0.5 km inland fo each PB revealing indirect anthropogenic activities are presented.
The parameters used to evaluate the quality of LULC results are: (1) the producer accuracy, which indicates the degree of confidence for correctly classified areas based o the validation of the training samples of a given class; (2) the user's accuracy, which is th percentage of correctly classified areas within all samples classified; (3) the overall accu racy (OA), which is the degree of confidence for the correctly classified areas based on a validation training samples; and (4) the Kappa coefficient, which is the agreement of th classification results with the validation training samples. 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.

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.

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).

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.

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].

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.

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.

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.

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)

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.

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: Variance: where M is the mean, V is the variance, i is gray tone, and X is a matrix of gray tone of the image.
(a) (b) (c) (d) 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.

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. Mean: Variance : 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.

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) 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). 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.

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.

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.

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

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 ; OA and kappa coefficient for the Patti-Tindari's PBs East (SIC04ME07), Center (SIC04ME08), and West (SIC04ME09) (c); OA and kappa coefficient for the two PBs of East Capo Calavà (SIC05ME10) and West Capo Calavà (SIC06ME11) on the coast of Gioiosa Marea (d); OA and kappa coefficient for the PB of Lipari-La Forbice (SIC09ME14) (e); OA and kappa coefficient for the PB of Lipari-Pignataro di Fuori (SIC10ME15) (f); OA and kappa coefficient for the two PBs of Lipari-Sabbie Bianche (SIC11ME16) and Lipari-Lido Blu (SIC12ME17) located on the east coast of Lipari (g); OA and kappa coefficient for the two PBs of Volcano Island, Punta dell'Asino (SIC07ME12) and Punta Bandiera (SIC08ME13) (h); OA and kappa coefficient for four PBs of Zimmaro (SIC15ME20), West Preistorico (SIC16ME21), and Central Preistorico (SIC16ME22) of Panarea Island (i); OA and kappa coefficient for one PB of Filicudi-Le Punte (SIC17ME24) (j); OA and kappa coefficient for one PB of Punta Scario (SIC13ME18) located on the coast of Malfa on Salina Island (k); and OA and kappa coefficient for one PB of Pollara (SIC14ME19) located on the coast of Malfa (l).

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.