Using PRISMA Hyperspectral Data for Land Cover Classification with Artificial Intelligence Support

Abstract

Hyperspectral satellite missions, such as PRISMA of the Italian Space Agency (ASI), have opened up new research opportunities. Using PRISMA data in land cover classification has yet to be fully explored, and it is the main focus of this paper. Historically, the main purposes of remote sensing have been to identify land cover types, to detect changes, and to determine the vegetation status of forest canopies or agricultural crops. The ability to achieve these goals can be improved by increasing spectral resolution. At the same time, improved AI algorithms open up new classification possibilities. This paper compares three supervised classification techniques for agricultural crop recognition using PRISMA data: random forest (RF), artificial neural network (ANN), and convolutional neural network (CNN). The study was carried out over an area of 900 km² in the province of Caserta, Italy. The PRISMA HDF5 file, pre-processed by the ASI at the reflectance level (L2d), was converted to GeoTiff using a custom Python script to facilitate its management in Qgis. The Qgis plugin AVHYAS was used for classification tests. The results show that CNN gives better results in terms of overall accuracy (0.973), K coefficient (0.968), and F1 score (0.842).

1. Introduction

The potential for Earth observation (EO) from space for resource management and optimization emerged in the early 1960s and throughout the 20th century, with the advent of digital computers, pattern-recognition technology, and the first artificial satellites. The ability to monitor large areas, while reducing per-unit costs, is one of the main advantages of space-based technologies [1]. Remote sensing (RS) through EO offers the possibility to obtain quantitative information at the pixel level. Having this information helps to understand environmental and socioeconomic trends. Land cover and land use maps can be used for sustainable resource management and the study of climate change phenomena. According to the Eurostat definition, land cover (LC) corresponds to the physical coverage of the land surface, while land use (LU) refers to the socioeconomic use of the land [2]. It is possible to create thematic maps of land cover using different classification methods. Whereas photointerpretation can be ambiguous and subjective, automatic classification is a defined, quantifiable, and repeatable process [3] (p. 359).

Classifying involves assigning a land cover class to all pixels in a digital image. Automatic LULC classification is typically based on the notion that different land cover types have different spectral reflectance behaviors [4]. Classification techniques (e.g., maximum likelihood, decision tree, and neural networks) are then used to define spectral signatures using sample data to discriminate between different classes of selected LULC based on pixel values [5,6]. Image classification, which can be supervised or unsupervised, is the process of extracting semantic Information from raw data, i.e., pixel values, by assigning a class label to each pixel [7,8]. When classifying land cover, supervised approaches offer better performance than unsupervised approaches but require a sufficient number of accurate samples [9].

Machine learning (ML) and deep learning (DL) are two approaches that use artificial intelligence (AI) to classify RS images. Among the classification techniques, support vector machine (SVM), random forest (RF), and maximum likelihood are ML-based algorithms. On the other hand, artificial neural networks (ANNs) can be either ML or DL depending on the layers of the network [10]. All these methods are common supervised classification techniques that are also used for hyperspectral imagery (HSI) [11]. The use of ANN for image classification and change detection has been a well-known technique for many years. Later, the focus shifted to ML algorithms, such as SVM, that could handle large data sets with few training samples, or such as RF, that were easy to use with good accuracy [10]. In the last decade, however, with the rise of DL, there has been a renewed interest in ANNs for their ability to produce good results in image analysis, including land cover classification [10]. These methods are derived from neural networks but have greater computational power. DL methods use deeper layers to extract feature information, particularly in the case of CNN, by considering both spatial and spectral characteristics of images [10]. DL models provide promising results for object recognition or classification of hyperspectral data, by better handling images with high spatial and spectral resolution [11,12].

Multispectral data are a key resource in remote sensing because they provide more than one spectral measurement per pixel [13] (p. 4). Most multispectral satellites record information in the visible and near-infrared regions of the electromagnetic spectrum, in a number of bands ranging from three to six or more, improving the ability to distinguish man-made surfaces, vegetation, clear and turbid water, rocks, and soil [14]. However, bands in multispectral sensors are not contiguous across the spectrum and often have bandwidths of 100–200 nm. As a result, unlike hyperspectral sensors, they do not have sufficient spectral resolution to directly identify materials with diagnostic characteristics [15] (p. 306). Hyperspectral data cover a wide range, from visible (0.4–0.7 nm) to shortwave infrared (SWIR) (2.4 nm) and are useful for detailed land cover legends. The ability to distinguish features with similar spectral signatures is enhanced by the availability of multiple bands [16]. In addition to improving the ability to discriminate between similar objects, hyperspectral data make it possible to perform advanced studies, for example, predicting the type and amount of crops traits by estimating grassland biochemical parameters [17,18].

HSI from aircraft has been in development since the 1980s, but the first hyperspectral space EO missions for civil and scientific purposes have only been available since the early 21st century [19]. The main providers of space-based hyperspectral data over the last few decades have been Earth Observer-1 (NMP/EO-1) [20], launched under NASA’s New Millennium Program in 1999, with the Hyperion spectrometer, and PROBA (Project for On-Board Autonomy), launched in 2001 with the Compact High-Resolution Imaging Spectrometer (CHRIS) developed by the European Space Agency [21]. Among the latest EO hyperspectral missions are PRISMA (Hyperspectral Precursor of the Application Mission) and EnMAP (Environmental Mapping and Analysis Program). EnMAP is a German imaging spectroscopy mission that was launched in April 2022 and recently completed its commissioning phase [22,23]. PRISMA is developed by the ASI and has been in orbit since March 2019, with commissioning completed in January 2020 [24]. The PRISMA mission is expected to contribute to the advancement of environmental RS by providing hyperspectral data for various applications, such as monitoring of agricultural crops, forest resources, and inland and coastal waters; mapping of natural resources, soil properties, and soil degradation; climate change studies; and environmental research [25].

Working with HSI raises issues that are well known in the scientific community, related to the size of the data. The very high dimensionality of hyperspectral data, due to the large amount of information recorded in different bands, has some disadvantages: high storage costs, redundancy, and degraded performance [26]. The dimensionality reduction is addressed using techniques based on band selection or feature extraction, as we will see below. However, the reduction in the size of the data must preserve the most relevant information it contains. To predict the separability of two classes of materials, the statistical distance between two spectral bands must be measured, and one of the most common methods is the Bhattacharyya distance [1]. The other issue is related to the notion that the number of samples that are used to train a classifier has an impact on its accuracy. Therefore, in order to maintain accuracy, the number of training pixels per class needs to increase as the dimensionality of the data increases, according to the curse of dimensionality or the Hughes phenomenon (Hughes 1968) [13,27].

The literature on the use of PRISMA data, particularly in the field of land cover classification, is still limited. Most of the research shows the potential of PRISMA data for specific purposes, such as forest conservation with wildfire fuel mapping [28] or fire detection [29], geological applications [30], cryospheric applications [31], urban surface detection [32], and mapping methane point emissions [33]. There are also interesting studies in the agricultural field dealing with specific crop or vegetation type discrimination [34,35,36,37]. However, the possibility of using PRISMA data to distinguish LULC classes has not been fully investigated.

This research aims to evaluate the possibility of classifying permanent agricultural crops (i.e., orchards, fruit trees, olive groves, etc.) in a highly fragmented agricultural area using PRISMA data. The purpose of this study is to discriminate entities with very similar spectral signatures using HSI. From this point of view, high spectral resolution becomes an advantage. However, the handling of hyperspectral images may not be an easy task due to the high dimensionality of the data [26,27]. Therefore, as explained in the following sections, two techniques were tested to achieve data dimensionality reduction. For the purposes of the research, it was decided to compare three different AI-based classifiers by evaluating the accuracy of their results. In particular, among the consolidate methods for supervised classification tasks, RF, ANN, and CNN were chosen. This choice was made to test the PRISMA data using algorithms with known performance, and which are expected to produce different results [38]. In addition, it was found useful to report the processing time for data of different dimensionality, used as input in each classifier.

2. Materials and Methods

The main steps of the proposed methodology are shown in Figure 1. The starting point of the procedure is the selection of the study area. Then, on the one hand, sample data are collected, both in the field and through photointerpretation, for the construction of a training data set, which is also crucial for the selection of the LC classes. On the other hand, the PRISMA satellite cube downloaded from the ASI portal was processed, including dimensionality reduction. Finally, three types of classifiers (RF, ANN, and CNN) were used for the supervised classification of the hyperspectral cube. The classification maps obtained were subjected to an accuracy evaluation phase, which allowed the selection of the best result.

2.1. Study Area

The study area is in Italy, between the provinces of Naples and Caserta, and covers an area of approximately 900 km² north of the city of Naples, as shown in Figure 2.

The area is largely coincident with what was once Campania Felix and is still known as “Terra di lavoro” [39] (pp. 16–17). This is a large flat agricultural area, on the left side of the Volturno River, partially crossed by a network of drainage canals known as Regi Lagni. It is an area characterized by a high degree of fragmentation. The complexity of agricultural uses in this area can be derived from the Campania Agricultural Land Use Map (CUAS) data, as shown in

Table 1. CUAS surfaces classification. Source: Campania Region 2009 [40].

CUAS Class	Hectares	%
Urbanized environment and artificial surfaces	25,373	28.30%
Orchards and small fruit trees	15,931	17.77%
Spring–summer crops—vegetable crops	12,594	14.05%
Spring–summer crops—industrial crops	8632	9.63%
Natural pasture and upland grassland	4401	4.91%
Vernal fall crops—cereals	3703	4.13%
Deciduous forests	3354	3.74%
Spring–summer crops—grain crops	2801	3.12%
Permanent pasture, meadow	2588	2.89%
Complex cropping and plots	2512	2.80%
Olive groves	2245	2.50%
Protected crops—horticultural and fruit crops	1240	1.38%
Temporary crops associated with permanent crops	1180	1.32%
Areas with sparse vegetation	655	0.73%
Grasslands	649	0.72%
Water	508	0.57%
Fall green crops—tubers	397	0.44%
Shrubs and bushes	220	0.25%
Vineyards	161	0.18%
Pastures unused or of uncertain use	138	0.15%
Pasture	127	0.14%
Poplar groves, willow groves, other deciduous trees	92	0.10%
Bare rocks and outcrops	82	0.09%
Citrus groves	44	0.05%
Coniferous forests	16	0.02%
Areas of natural colonization	11	0.01%
Other permanent crops or orchards	6	0.01%
Areas degraded by fire and other events	3	0.00%
Artificially recolonized areas (reforestation)	2	0.00%
Mixed broadleaf and coniferous forests	1	0.00%
Tot.	89,666	100%

[40]. The CUAS is an official map, accredited at the conventional scale of 1:50,000. It has been used in this research because it allows, although with unproven accuracy, the surveying of much smaller minimum mappable units than the Corine Project, which has a minimum mapping unit of 25 ha [41]. In 2009, the reference year of the CUAS, the area was mainly occupied by arable land (31%), artificial areas (28%), and orchards (18%). This study area is particularly suitable for testing the classification of complex mosaics due to its extreme fragmentation.

Table 2. CLC surfaces classification. Source: Corine project 2012 [41].

CLC Class	Code_12	Hectares	%
Non-irrigated arable land	211	20,982	23%
Complex cultivation patterns	242	16,687	19%
Continuous urban fabric	111	14,062	16%
Fruit trees and berry plantations	222	10,641	12%
Discontinuous urban fabric	112	6351	7%
Industrial or commercial units	121	3877	4%
Broad-leaved forest	311	3215	4%
Permanently irrigated land	212	3038	3%
Natural grasslands	321	2399	3%
Olive groves	223	1987	2%
Land principally occupied by agriculture […]	243	1441	2%
Annual crops associated with permanent crops	241	1416	2%
Transitional woodland–shrub	324	954	1%
Sclerophyllous vegetation	323	553	1%
Mineral extraction sites	131	431	0%
Pastures	231	411	0%
Road and rail networks and associated land	122	368	0%
Airports	124	225	0%
Water courses	511	195	0%
Dump sites	132	157	0%
Green urban areas	141	141	0%
Construction sites	133	48	0%
Burnt areas	334	43	0%
Sport and leisure facilities	142	26	0%
Coniferous forest	312	18	0%
Tot.		89,666	100%

shows the data obtained from the Corine Land Cover Classification (CLC) [41] for the year 2012. The largest percentage of the study area (30%) is occupied by anthropic areas (code_12 “111” to “142”). Agricultural areas (code_12 “211” to “112”) occupy 23%, while areas with permanent tree crops (code_12 “212” to “223”) occupy about 17% of the area.

Based on the study area configuration, 18 classes were selected for ground truth collection and image classification. As shown in

Table 3. The following table shows the correspondence between the level III classification of the CLC project and the classes used in this study.

Code_12	CLC Class	Matching Class	ID
111	Continuous urban fabric	Impervious	1
112	Discontinuous urban fabric
121	Industrial or commercial units
122	Road and rail networks and associated land
124	Airports
131	Mineral extraction sites	Bare soil	2
132	Dump sites	Impervious	1
133	Construction sites
141	Green urban areas *
142	Sport and leisure facilities
211	Non-irrigated arable land	Bare soil—crops or low vegetation **	2–3
212	Permanently irrigated land
221	Vineyards	Vineyards	18
222	Fruit trees and berry plantations	Citrus grove	9
		Apricot orchard	10
		Cherry orchard	11
		Persimmon orchard	12
		Hazelnut orchard	13
		Walnut orchard	14
		Peach orchard	15
		Poplar orchard	16
223	Olive groves	Olive grove	17
231	Pastures	Bare soil—crops or low vegetation **	2–3
241	Annual crops associated with permanent crops
242	Complex cultivation patterns
243	Land principally occupied by agriculture […]
311	Broad-leaved forest	Woods and forest	4
312	Coniferous forest
321	Natural grasslands	Bare soil—crops or low vegetation **	2–3
323	Sclerophyllous vegetation	Woods and forest	4
324	Transitional woodland–shrub
334	Burnt areas	Bare soil	2
511	Water courses	Water	7
-	-	Clouds	5
-	-	Shadows	6

* Considering pixel spatial resolution in terms of 30 × 30 m of PRISMA image, it was preferred to assign this class to the class containing impervious surfaces. ** Depending on the case, there could be correspondence for class ID 2 or 3, as this study was not based on multitemporal analysis. Some CLC classes (i.e., Code_12: 211;212;231;241;242;243;321) may correspond to the class “ID 2—Bare soil” or alternatively to the class “ID 3—Crops or low vegetation”, depending on the presence or absence of vegetation cover at the time of PRISMA image acquisition (e.g., because some areas have been tilled with agricultural machinery). The authors did not combine the two classes when classifying, but the two classes need to be aggregated to compare them with CLC classes in this table.

, the level of detail was decreased or increased starting from the subdivision corresponding to Corine project level III. The level of detail was reduced for classes that were not of interest for the case study if these classes had similar characteristics. Conversely, the level of detail was increased for categories of relevance to the case study, particularly in the case of fruit trees.

2.2. PRISMA Hyperspectral Image

According to the technical specification document [42], the PRISMA satellite orbits at an altitude of 615 km and can deliver 223 images per day. The PRISMA spectrometer is integrated with a medium-resolution panchromatic camera and can provide images with high spectral resolution: the radiation reflected from the Earth’s surface is recorded in the following spectral ranges: VNIR (400–1010 nm) and SWIR (920–2505 nm). The hyperspectral cube consists of 240 frames with a bandwidth less than or equal to 12 nm, with a SWAT of 30 km at nadir and a GSD of about 30 m in the hyperspectral and 5 m in the panchromatic. The PRISMA products, available at http://prisma.asi.it/, are provided at different levels [43], in particular:

• Level 0 (L0) refers to raw data that include ancillary data from satellite instruments;
• Level 1 (L1) are top-of-atmosphere (TOA) radiance images, radiometrically corrected and calibrated in physical units;
• Level 2b (L2b), obtained by L1 product, spectral radiance is transformed into geophysical parameters, images are geolocated and geocoded at the bottom of the atmosphere;
• Level 2c (L2c) are geolocated at bottom-of-atmosphere reflectance images; Angstrom correction is applied;
• Level 2d (L2d) product is a geocoded version of L2c; images are orthorectified by the ASI using a DEM.

The study area is bounded by 30 × 30 square kilometers corresponding to the scene “PRS_L2D_STD_20191130100450_20191130100454_0001” acquired by the PRISMA satellite on 30 November 2019 and shown in Figure 3.

For this study, it was decided to use the L2d level of the ASI data provided with reflectance values, geocoded and orthorectified. “SRTM90_CGIAR_4.1” is the digital elevation model (DEM) used by the ASI to orthorectify the L2d product of the scene. The STRM acronym stands for “Shuttle Radar Topography Mission” and refers to the joint mission of NASA and the German and Italian space agencies that started in 2000, with the goal of obtaining a high-quality global digital elevation model [44]. Version 4.1 of the SRTM DEM is distributed by the Consortium for Spatial Information (CGIAR-CSI) [45]. It is important to note that global DEMs are generally less accurate than local DEMs and may contain errors, the main sources of which are the post-processing or generation techniques used [46]. The area chosen for this study is largely flat, so negligible errors in altitude can be accepted. However, it was observed that the geocoded data were also not very accurate in spatial distribution (+/−300 m in latitude and +/−150 m in longitude). In addition, functionality to handle the HDF5 format of PRISMA data was not implemented in all mainstream software until a few years ago. For these reasons, a Python script (Prisma Tool) was developed [47]. The Prisma Tool allows data to be converted from HDF5 to the more common GeoTiff format, making the data easier to use in GIS software. At the same time, it allows for better georeferencing of images through the input of the correct vertex coordinates.

Dimensionality Reduction

As mentioned above, hyperspectral data have great potential for distinguishing features with similar spectral signatures but also some disadvantages related to data dimensionality. Dimensionality reduction techniques can generally be divided into two types: those based on feature extraction and those based on band selection. Those based on feature extraction, such as principal component analysis (PCA), which selects a subset of features with high variance from the original data, involve a loss of information. This is why feature selection is the preferred dimensionality reduction technique in HSI, where a subset of bands is selected without loss of their physical meaning [48]. Several techniques have been proposed in the literature to solve the band selection challenge [49].

In this case study, starting with the original 240-band cube (Cube-O, where “O” stands for original), two different methods were used.

• The first is based on an approach like supervised band selection techniques that allow us to choose which bands to keep. This step was aimed at eliminating bands corresponding to atmospheric windows with low or zero radiance values, rather than greatly reducing the size of the data. Random spectral signatures associated with different types of surfaces were compared as shown in Figure 4. In accordance with previous approaches [48,50], the bands with the lowest signal-to-noise ratio and atmospheric water vapor absorption were selected for removal. Through this process, a new hyperspectral cube of 156 bands (Cube-R1, where “R1” indicates the first reduction) was derived from the original 240 bands. The bands removed were those related to the ranges 1290–1490 nm, 1700–2050 nm, and 2350–2495 nm.

Figure 4. Random spectral signature of the original PRISMA Cube-O; the red rectangles correspond to the area in which the bands with the lowest values have been removed—authors’ elaboration using Qgis with the EnMAP-Box plugin described in the following section.

Figure 4. Random spectral signature of the original PRISMA Cube-O; the red rectangles correspond to the area in which the bands with the lowest values have been removed—authors’ elaboration using Qgis with the EnMAP-Box plugin described in the following section.

• The second method is based on the application of the statistical method of PCA [51] (pp. 115–128), which is widely used in the scientific community [49] to reduce the cube dimension, removing bands that are strongly correlated to others. The original hyperspectral cube was used as input data. At the end of the processing, 20 principal components were selected, preserving a percentage of variance in the data of 98.7% (Cube-R2, where “R2” indicates the second reduction). The AVHYAS plugin described in the following section was used for PCA.

It is worth noting that the Cube-R2 obtained through the PCA has kept a high percentage of the variance that was computed on Cube-O. PCA, which is very useful to distinguish between classes that are very different from each other, could mask the information recorded in the original data about small differences between similar classes. For the classification task, the authors decided to alternatively test the two images (Cube-R1 and Cube-R2) obtained using the models described above to evaluate this issue.

2.3. Software Used in Research

Qgis has been used as an open-source GIS software for the purpose of this study. This software is an official project of the Open-Source Geospatial Foundation (OSGeo) [52] and is one of the most widely used open-source GIS software worldwide due to its ability to easily implement a variety of functions through plugins. Qgis offers several plugins for supervised classification of satellite images.

However, the ability to classify hyperspectral data, especially using DL-based approaches such as neural networks, is still limited. A solution is offered by the Advanced Hyperspectral Data Analysis (AVHYAS) plugin [53]. The AVHYAS plugin has been developed by the Space Applications Centre (SAC) of the Indian Space Research Organization. It is a HSI processing and analysis plugin based on Python-3. The plugin extends the existing Qgis GDAL functionality. It uses several modules such as TensorFlow for DL and Scikit-Learn for machine learning. The purpose of the plugin is to exploit the visualization capabilities of Qgis by integrating them with the most advanced techniques available in the literature for hyperspectral data management and classification. The ability to experiment with different supervised classification techniques based on ML and DL within the same application was another reason for choosing AVHYAS.

Another very useful plugin was EnMAP-Box, developed by the Humboldt University of Berlin and the University of Greifswald, to handle data from the German hyperspectral satellite EnMAP [54]. This plugin constituted the starting point for the development of the AVHYAS plugin mentioned above. It combines the visualization capabilities of Qgis with the ability to manage and process hyperspectral data, as in the previous case. Some of its main features are listed below. The plugin further provides (i) the ability to easily visualize HSI, both as a two-dimensional image and as a hyperspectral cube; (ii) the exploration of high-dimensional data and hyperspectral libraries; and (iii) the collection of spectral signatures in spectral libraries. It also allows data to be processed using some of the most popular algorithms in the scientific community (i.e., support vector machines or RF-based raster classification, regression, and clustering approaches from the scikit-learn library).

For field sample data collection, the QField mobile application [55] has been used; it is based on Qgis and was developed specifically for field work. The main goal is to help users perform the tasks required in the field, such as data collection and storage.

2.4. Sample Data

In order to apply supervised classification techniques, it was necessary to collect a dataset of labeled data, to be used both for training the classifiers and testing the results. In addition, as mentioned above, the number of training samples affects the accuracy of a classifier. Using a large number of bands requires a large number of observations to maintain classification accuracy [20]. Since it was chosen to maintain a large number of bands as the input data, a consistent number of ground truths was collected. Specifically, three sets of labeled data were used. Each set was collected at different times and in different ways.

(A) The first set was collected on 23 October 2019. The coordinates of the collected points, corresponding to the centroids of the agricultural plots, were obtained using differential GPS and then entered into a CSV file. The CSV file, which also included attributes related to the type of LC, was then imported into GIS software. Using the coordinates in the attribute table, the points were automatically georeferenced. Thus, a vector datum of the points was obtained, which in turn made it possible to derive polygon vector data, considering the perimeter of the cultivated areas. Photointerpretation and geotagged images collected during the survey made this step possible. During the first survey, 160 control points were surveyed and classified simultaneously as the satellite image was acquired. Of these, 106 involved agricultural parcels, which in turn were divided into 43 points for permanent crops and 63 for temporary crops. Figure 5 shows the preservation of the canopy for almost all the permanent crops and trees classified in this study at the time of the first survey. This study is not based on a multitemporal analysis of the study area. The classification was made on the basis of a November image. During this period, the canopy of the trees may have had a low number of leaves. However, the leaves of the plants selected for this study were not yet in the phase of leaf abscission, as can be seen in the photos below.

(B) Later, research was suspended due to the COVID-19 pandemic. However, it was necessary to increase the number of labeled data to have more confidence in the statistical validation of the results. It was assumed that ground truths could still be integrated for the arboriculture sector after three years of image acquisition. In the survey, great care was taken in the estimation of the probable date of planting from the tree characteristics. In this way, only plantings that were older than the date of the satellite image was included in the survey. The second survey took place in April 2023 and allowed the detection of 150 additional ground truths. This brings the total number of GCPs related to permanent cultures to 180. In this case, the QField app was used during the survey to collect point vector data directly in the field, using the AGPS of a smartphone. In an attempt to open up to the world of citizen science and crowdsensing data collection, this type of approach was preferred over the traditional one. The accuracy of AGPS provided by common smartphones was considered adequate for this study, given the 30 m spatial resolution of the PRISMA image. Using the QField app, the collected points could be entered directly into the GIS environment in their correct location [55]. This is possible through the activation of a satellite base map and smartphone localization in the app during the survey. The georeferenced vector point datum was then used to derive a polygon datum using the same methodology as the first labeled dataset (A).

(C) Photointerpretation also allowed for the addition of new classes such as clouds, shadows, and water. It also increased the area for some of the classes already present, such as impervious surfaces, bare soil, and temporary agricultural surfaces. Photointerpretation was performed using the results of different false-color compositions of the PRISMA image. A total of six different false-color band compositions were used. The main objective was to improve the ability to distinguish between bare soil and temporary or permanent vegetation by photointerpretation. The selection of the bands is based on the combination method of the raster layer styling panel of the Qgis EnMAP-Box plugin [54]. Table 4 shows the RGB bands used for each composition, and Figure 3 shows the resulting maps.

Table 4. This table shows the bands used in false-color composition to obtain maps 1 to 6, shown in Figure 6.

		Red Channel		Green Channel		Blue Channel
	MAP	Band Number	nm	Band Number	nm	Band Number	nm
Natural color	1	34	0669.5	22	0562.5	13	0492.5
False color	2	34	0669.5	50	0833.4	22	0562.5
Color infrared	3	50	0833.4	34	0669.5	22	0562.5
Agriculture	4	128	1616.5	50	0833.4	13	0492.5
Healthy vegetation	5	50	0833.4	128	1616.5	13	0492.5
Recent harvest	6	193	2198.8	50	0833.4	22	0562.5

Table 4. This table shows the bands used in false-color composition to obtain maps 1 to 6, shown in Figure 6.

		Red Channel		Green Channel		Blue Channel
	MAP	Band Number	nm	Band Number	nm	Band Number	nm
Natural color	1	34	0669.5	22	0562.5	13	0492.5
False color	2	34	0669.5	50	0833.4	22	0562.5
Color infrared	3	50	0833.4	34	0669.5	22	0562.5
Agriculture	4	128	1616.5	50	0833.4	13	0492.5
Healthy vegetation	5	50	0833.4	128	1616.5	13	0492.5
Recent harvest	6	193	2198.8	50	0833.4	22	0562.5

Figure 6. Maps resulting from false-color combination: 1—natural color; 2—false color; 3—color infrared; 4—agriculture; 5—healthy vegetation; 6—recent harvest.

Figure 6. Maps resulting from false-color combination: 1—natural color; 2—false color; 3—color infrared; 4—agriculture; 5—healthy vegetation; 6—recent harvest.

Table 5. Ha surfaces of labeled data distinguished by class and group: A—ground truths surveyed in 2019; B—ground truths surveyed in 2023; C—labeled data obtained by photointerpretation.

ID	Class	A	B	C	Tot.
1	Impervious	8.23	0.00	203.85	212.08
2	Bare soil	26.35	8.75	86.81	121.90
3	Temporary crops	19.83	1.29	47.29	68.41
4	Woods and forest	6.49	0.82	42.60	49.91
5	Clouds	0.00	0.00	176.35	176.35
6	Shadows	0.00	0.00	217.53	217.53
7	Waters	0.52	0.00	17.90	18.42
8	Greenhouses	1.97	0.00	27.65	29.62
9	Citrus groves	0.00	2.88	0.00	2.88
10	Apricot orchard	0.73	9.28	0.00	10.00
11	Cherry orchard	0.00	4.34	0.00	4.34
12	Persimmon orchard	4.53	8.82	0.00	13.35
13	Hazelnut orchard	0.00	18.14	0.00	18.14
14	Walnut orchard	5.99	12.61	0.00	18.60
15	Peach orchard	4.87	1.44	0.00	6.32
16	Poplar orchard	2.16	2.73	0.00	4.89
17	Olive groves	11.29	5.40	0.00	16.70
18	Vineyard	1.11	3.68	0.00	4.79
Tot.		94.07	80.18	819.98	994.23

shows the dimensions in hectares of the surfaces in the sample dataset by class and group.

2.5. Supervised Classification Techniques

The primary purpose of this study is to evaluate the possibility of using PRISMA data to classify land cover classes, in particular different types of orchards, as mentioned in the introduction. For this purpose, three of the most common [11,38] classification techniques were investigated: RF, ANN, and CNN. Classification tests were performed in Qgis 3.14.1 software using the AVHYAS plugin [53]. The different techniques were tested by classifying the same study area. The images used as input were both the one corresponding to the Cube-R1 with 156 bands and the one corresponding to the Cube-R2 obtained by the PCA. The same number of labeled data was used in each test run. In the case of neural networks, the training data set was randomly divided into training (70%) and test (30%) sites by AVHYAS prior to each test run. For random forests, k-fold cross-validation with 3 iterations was used for validation. The AVHYAS plugin [53] provides some of the most common networks in the literature [38] through a dedicated DL classification module. Reflecting the evolution of the literature, 13 different neural network architectures are available within this module [53]. The default setting was used in the selection of some parameters for the proposed study, since the aim was not to define the optimal scenario in the use of this tool. In addition to different input data, the number of trees for the RF algorithm and the number of epochs for the neural networks differed between tests.

2.5.1. Random Forest

The RF algorithm developed by Breiman [56] in 2001 is an ML-based approach. Classification is based on individual decision trees working as an ensemble. A random subset of labeled samples is used to train each tree. The randomization of the training data reduces the probability of error as the number of trees increases. Each tree expresses only one vote per instance, and the final classification is given by the majority of votes from all trees [56,57]. The choice of the number of trees affects the classification result. Considering similar works, a suitable number to balance accuracy and time spent is around 500 [58,59]. However, in this study, six different tests with 100, 500, and 1000 trees were performed, alternately using both Cube-R1 and Cube-R2 as input, to test the performance of the used data.

2.5.2. Artificial Neural Network

Some tests were carried out with a classification approach based on the ANN, using the model of the basic network proposed by the plugin. In this case, the “Baseline (Fully Connected NN)” model was chosen. This is a sequential model of a fully connected network with 3 hidden layers, where each neuron is fully connected to all those in the previous layer. The learning rate has been left at the default value of 0.0001. This value, considered small, is slower in that it requires more epochs for training, but it is generally more suitable for dealing with complex problems [60]. A total of 6 trials were run. The number of epochs chosen from time to time was 10, 50, and 100, where Cube-R1 was used as input in 3 cases, and Cube-R2 was used in the remaining 3 cases. Some tests were performed using an ANN-based classification approach using the DL module available in AVHYAS [53].

2.5.3. Convolutional Neural Network

The utility of using CNNs for HSI classification has been widely described in the literature [38,59,61]. Thanks to specific convolutional kernels, CNN can use both spatial and spectral information to classify images. In the present work, among the models available in the AVHYAS plugin, the “Hu (1-D CNN)” was used. This model is based on the one proposed by Hu [62]. The architecture of the network is similar to the one shown in Figure 7. It is a one-dimensional network and considers only spectral signatures by using a convolutional layer and a fully connected layer. The default value for the learning rate, which was kept in this case, is 0.01. Six trials were run with the number of epochs 10, 50, and 100. Cube-R1 and Cube-R2 were used alternately as input data.

2.6. Accuracy Assessment

In case studies involving the use of remotely sensed data to produce land cover maps, assessment of the accuracy of the resolutions is essential. For quantitative assessment of the accuracy of the results, an error matrix can be generated, which can be used for a number of descriptive and analytical statistical analyses. The most common terms are overall accuracy (OA), kappa coefficient (K) (Cohen 1960), user’s accuracy (UA), and producer’s accuracy (PA) [63]. OA is the percentage of correctly classified pixels and is given by the ratio of the number of correctly classified pixels to the total number of test pixels used. Cohen’s K is a measure of the agreement between the predicted values and the actual values. This metric measures how much better the classification that was made is in comparison to a classification made by chance. UA estimates the commission error and indicates the probability that a pixel predicted to be in a certain class actually belongs to that class. PA estimates the omission error and indicates the probability that a pixel belonging to a given class has been correctly classified. Another statistical indicator used is the F1 score, which is the harmonic mean between precision and recall for each class, varying from 0 (worst case) to 1 (best case) [64].

The error matrix automatically generated by the AVHYAS plugin in the final report was used to evaluate the results in this case study. The classification results were evaluated globally in terms of OA and K. Results for individual classes were evaluated in terms of UA, PA, and F1 score. In the case of random forest, cross-validation was used, so the validation was carried out on the total number of labeled pixels. In the case of neural networks, 30% of the randomly selected pixels for the test phase were used to validate the results.

3. Results

This section presents the results of the classification techniques used. For each one of the three methods of classification used (RF, ANN, and CNN), the results are very different in terms of OA, K, and F1 score indices and also for PA, UA, and F1 evaluated for each class.

Table 6. RF, ANN, and CNN comparison.

	Input Data	OA	K	F1	h
RF	Cube-R2 (20 PCA)	0.89	0.87	0.60	0.1
ANN	Cube-R1 (156 bands)	0.96	0.96	0.77	67.2
CNN	Cube-R1 (156 bands)	0.97	0.97	0.84	2.4

shows a comparison of the best results of the methods used for classification in terms of global accuracy (OA, K, and F1) for all the classes and running time (h).

Figure 8 shows a comparison graph between the RF, ANN, and CNN performances in terms of the F1 score obtained for the classes of permanent agricultural cover that are mainly considered in this study (ID from 9 to 18).

The following subsections report the values of the best results obtained in the different tests using RF, ANN, and CNN. In the first table of each subsection, in addition to OA and K, the data used as input and the processing time are also reported. The second table in each paragraph shows the results for each of the classes considered and the area in Ha obtained.

3.1. Random Forest Results

Table 7. RF classification results, overall evaluation.

Input	n. Trees	OA	K	F1	Hours
Cube-R2	500	0.887	0.867	0.603	0.052

shows the overall values for the best result of the RF algorithm obtained with the PCA-derived Cube-R2. Figure 9 shows the resulting classification map, while a detail is shown in Figure 10.

Table 8. RF classification results, per class evaluation.

ID	Class	UA [%]	PA [%]	F1 [%]	Ha
1	Impervious	93%	88%	91%	37,203
2	Bare soil	92%	96%	94%	4922
3	Temporary crops	83%	80%	81%	16,671
4	Woods and forest	84%	81%	83%	5444
5	Clouds	100%	100%	100%	2698
6	Shadows	99%	99%	99%	53,888
7	Waters	72%	99%	84%	225
8	Greenhouses	91%	88%	90%	1002
9	Citrus groves	27%	69%	39%	103
10	Apricot orchard	31%	28%	30%	2690
11	Cherry orchard	6%	16%	9%	123
12	Persimmon orchard	40%	38%	39%	2402
13	Hazelnut orchard	72%	62%	67%	377
14	Walnut orchard	40%	33%	36%	3539
15	Peach orchard	29%	38%	33%	342
16	Poplar orchard	18%	71%	28%	213
17	Olive groves	59%	60%	59%	2566
18	Vineyard	14%	78%	24%	78

shows the values by class for the best result obtained with the RF algorithm.

3.2. Artificial Neural Network Results

Table 9. ANN classification results, overall evaluation.

Input	n. Epochs	OA	K	F1	Hours
Cube-R1	100	0.963	0.956	0.766	67.166

shows the overall values for the best result obtained with ANN, using Cube-R1 with 156 bands as input. OA and K index values are very high but not homogeneous, because they are influenced by the results obtained in the cases of the first eight classes. Figure 11 shows the resulting classification map, while a detail is shown in Figure 12.

Table 10. ANN classification results, per class evaluation.

ID	Class	UA [%]	PA [%]	F1 [%]	Ha
1	Impervious	100%	98%	99%	34,668
2	Bare soil	99%	100%	99%	4628
3	Temporary crops	94%	96%	95%	13,118
4	Woods and forest	97%	99%	98%	2500
5	Clouds	100%	100%	100%	3247
6	Shadows	99%	100%	100%	51,221
7	Waters	98%	100%	99%	664
8	Greenhouses	96%	100%	98%	1275
9	Citrus groves	60%	55%	57%	1224
10	Apricot orchard	91%	45%	60%	6589
11	Cherry orchard	29%	50%	36%	297
12	Persimmon orchard	77%	83%	80%	3375
13	Hazelnut orchard	98%	79%	88%	1703
14	Walnut orchard	71%	79%	75%	3296
15	Peach orchard	64%	67%	65%	2209
16	Poplar orchard	35%	55%	43%	862
17	Olive groves	87%	87%	87%	2852
18	Vineyard	0%	0%	0%	540

shows the values by class for the best ANN score.

3.3. Convolutional Neural Network Results

Table 11. CNN classification results, overall evaluation.

Input	n. Epochs	OA	K	F1	Hours
Cube-R1	100	0.973	0.968	0.842	2.399

shows the global values for the best result obtained with CNN, using Cube-R1 with 156 bands. Figure 13 shows the resulting classification map, while a detail is shown in Figure 14.

Table 12. CNN classification results, per class evaluation.

ID	Class	UA [%]	PA [%]	F1 [%]	Ha
1	Impervious	99%	100%	99%	30,850
2	Bare soil	99%	100%	100%	5429
3	Temporary crops	99%	95%	97%	13,956
4	Woods and forest	98%	96%	97%	4290
5	Clouds	100%	100%	100%	3050
6	Shadows	100%	100%	100%	53,105
7	Waters	100%	100%	100%	570
8	Greenhouses	99%	98%	98%	1733
9	Citrus groves	80%	100%	89%	1217
10	Apricot orchard	53%	71%	61%	3533
11	Cherry orchard	36%	50%	42%	622
12	Persimmon orchard	77%	77%	77%	2571
13	Hazelnut orchard	93%	93%	93%	718
14	Walnut orchard	84%	79%	81%	3778
15	Peach orchard	68%	75%	71%	1778
16	Poplar orchard	76%	57%	65%	1643
17	Olive groves	96%	85%	90%	3396
18	Vineyard	40%	86%	55%	2027

shows the values by class for the best result obtained with CNN.

4. Discussion

Land cover classifications based on multispectral data allowed us to distinguish more general land cover classes. The use of hyperspectral data has advanced the identification of detailed maps of LULC. However, as noted in the introduction, most of the research focuses on temporary crops [34,61,65] or a few tree species [57,66]. The results show that a detailed legend for land cover maps can be obtained by supervised classification using HSI. In particular, it is possible to extend the classification to the fourth level, starting from the third level of CLC. This study shows how the PRISMA data can be used to distinguish between many types of permanent crops. In this case study, CLC class “222—Fruit and berry plantations” was divided into eight subclasses. The results show that in the case of CNN, using Cube-R1, the F1 score is higher than 0.7 for the following classes: hazelnut orchard—0.93; olive groves—0.90; citrus groves—0.88; walnut orchard—0.81; persimmon orchard—0.77; peach orchard—0.71. Based on these results, it can be said that the main objective of this study has been demonstrated. Supervised classification using HSI provides an opportunity to obtain a detailed legend for land cover maps. Furthermore, to understand the actual applicability of the data in the application domain, the worst-case scenario was assumed (non-multitemporal analysis with a single autumn image), but even in this case, the results were very interesting. It is also true, as the photos in Figure 5 show, that the climatic conditions of the study area were conducive to the maintenance of the canopy over a longer period of time.

Among the different techniques used for the classification, the best choice seems to be the use of CNN. In terms of OA, K, PA, UA, and F1, this technique gave the best global result and per-class evaluation.

However, the classification with the RF algorithm gave interesting global results (OA—0.887; K—0.867; F1—0.603), using Cube-R2, but not comparable to those of the neural networks. It is an excellent alternative for land cover mapping with aggregated classes due to its very fast computation times. Indeed, excellent results were obtained with F1 scores above 80% for impervious, bare soil, temporary crops/low vegetation, and woods and forests. With respect to the classes of interest (ID from 9 to 18), the best results concern the classes hazelnut orchard (67%) and olive grove (59%).

Classification with ANN produced excellent overall results (OA—0.963; K—0.956; F1—0.766) using Cube-R1. The F1 score, which is higher than in the RF case, for some of the classes of interest (ID from 9 to 18), such as cherry and poplar groves, is less than 60% and 0% in the case of vineyards. It may be possible to obtain better results for these classes with an image from a different season. But in the case of olive, hazelnut, and persimmon groves, results higher than 80% in terms of F1 seem to be very favorable. Among the disadvantages, it is necessary to highlight the onerous cost of calculation in terms of hours spent.

In the case of convolutional neural networks, using Cube-R1 as in the previous case, the OA and global K coefficients are very high (OA—0.973; K—0.968; F1—0.842). This is because they are influenced by the very high level of accuracy in the case of aggregated classes such as impervious. However, in this case also, UA, PA, and F1 values are higher compared to the previous techniques. The convolutional network seems to guarantee better classification results, especially for the following classes: hazelnut, olive, persimmon, and walnut. On the contrary, for cherry, vineyard, and apricot, the classification is less accurate. For the classes of interest (ID from 9 to 18), the classification output shows percentages higher than 70% in terms of the F1 score in 6 out of 10 cases (hazelnut orchard—0.93; olive grove—0.90; citrus grove—0.88; walnut orchard—0.81; persimmon orchard—0.77; peach orchard—0.71).

Convolutional technology produces an overall F1 score that is almost 10% better than ANN and 25% better than RF. When evaluating the results obtained in the classes of interest (ID from 9 to 18), CNN guarantees almost 30% better F1 results compared to RF in all classes. Compared to ANN, there are improvements of more than 20% in the vineyard, citrus, and poplar classes; they are equivalent for the apricot and persimmon classes; while for all the others, there are slightly better results with CNN.

Based on the results obtained, it has been proven that the convolutional network guarantees excellent results with PRISMA hyperspectral data, especially with the cube resulting from the band selection (Cube-R1) and with a number of epochs equal to 100.

5. Conclusions

This study shows that the HSI data, by achieving very high classification results in terms of F1 scores for some classes, can discriminate between different types of permanent crops (ID from 9 to 18). As shown before, the results in terms of F1 score are higher than 0.70 in 6 out of 10 cases (hazelnut orchard—0.93; olive grove—0.90; citrus grove—0.88; walnut orchard—0.81; persimmon orchard—0.77; peach orchard—0.71). This shows that satellite imagery can be used for level IV classification and confirms studies based on PRISMA data to distinguish crops [34,67], vegetation [36], or forest types [35].

The most promising technology for this kind of application is based on neural network methods, especially CNNs (CNN F1 results are almost 30% better compared to RF for permanent crops). Indeed, various studies have demonstrated the potential of CNNs for hyperspectral data classification compared to other techniques [38,59,61,62]. The best results were obtained with CNN using Cube-R1 (156 bands). This confirms that strong dimensionality reduction is not necessary when using DL-based methods [37].

When dealing with high dimensionality of data, methods based on ML such as RF have limitations [68]. This aspect necessitates information reduction, which inevitably has an impact on the ability to discriminate classes with similar characteristics. However, as shown in [58], better results can be obtained with focused band selection.

The processing times may not be relevant for a case like this, where only a single data set is processed. However, this aspect could become crucial for multitemporal analyses or when experiments are carried out with different sample data. Although the best results were obtained with CNN, this method is more time consuming than RF. In fact, the main disadvantages of this method are the computational time and computer power consumption [59,69]. However, it is important to note that the algorithms have been used on open-source tools, so most of the processing/analysis is single-threaded. The time required could be significantly reduced as open-source software evolves to a multithreaded perspective.

In terms of reducing the high dimensionality of the data, the best results were obtained using Cube-R1 derived by band selection. Eliminating bands in the three regions of the spectrum where the transmittance is low is always useful [48,49,50]. Data obtained with this method preserve more than 150 bands of the original 240. It guarantees good results when used as input for neural networks, but not for RF. Feature extraction using the PCA method is essential in this case. Starting from the original Cube-O, PCA can significantly reduce data dimensions while preserving a high percentage of variance. Therefore, it is useful for discriminating between very different classes, but as stated in the previous paragraphs, it tends to overlook the few differences between similar classes.

As mentioned in the introduction, there is to date a lack of studies using PRISMA for LULC classification. The reasons are mainly due to the difficulties in accessing and managing the data. However, future improvements are expected in both data availability, partly due to new hyperspectral missions such as EnMAP [23], and data quality, using higher-accuracy DEMs for orthorectification. This would improve the reflectance accuracy, and thus even more accurate classifications are likely. It is expected that the interest in HSI will increase soon, especially because of the possibilities and advantages in the field of natural and agricultural land classification and monitoring.

The analysis of multi-seasonal imagery will certainly be a possible future application of PRISMA hyperspectral data for land cover mapping. The most promising way seems to be the use of convolutional networks. Future tests will certainly favor CNN by including different models for the network architecture. A chance to improve the results could certainly be achieved by making use of 3D CNNs that take advantage of the combination of spatial and spectral information [70]. Also, a different dimensionality reduction method with more specific band selection may be used. This would allow for better results with ML algorithms that are less time-consuming, such as RF, or less sensitive to the Hughes phenomenon, such as SVM [71].