A Comparison between Standard and Functional Clustering Methodologies: Application to Agricultural Fields for Yield Pattern Assessment

Abstract

The recognition of spatial patterns within agricultural fields, presenting similar yield potential areas, stable through time, is very important for optimizing agricultural practices. This study proposes the evaluation of different clustering methodologies applied to multispectral satellite time series for retrieving temporally stable (constant) patterns in agricultural fields, related to within-field yield spatial distribution. The ability of different clustering procedures for the recognition and mapping of constant patterns in fields of cereal crops was assessed. Crop vigor patterns, considered to be related to soils characteristics, and possibly indicative of yield potential, were derived by applying the different clustering algorithms to time series of Landsat images acquired on 94 agricultural fields near Rome (Italy). Two different approaches were applied and validated using Landsat 7 and 8 archived imagery. The first approach automatically extracts and calculates for each field of interest (FOI) the Normalized Difference Vegetation Index (NDVI), then exploits the standard K-means clustering algorithm to derive constant patterns at the field level. The second approach applies novel clustering procedures directly to spectral reflectance time series, in particular: (1) standard K-means; (2) functional K-means; (3) multivariate functional principal components clustering analysis; (4) hierarchical clustering. The different approaches were validated through cluster accuracy estimates on a reference set of FOIs for which yield maps were available for some years. Results show that multivariate functional principal components clustering, with an a priori determination of the optimal number of classes for each FOI, provides a better accuracy than those of standard clustering algorithms. The proposed novel functional clustering methodologies are effective and efficient for constant pattern retrieval and can be used for a sustainable management of agricultural fields, depending on farming systems and environmental conditions in different regions.

1. Introduction

Automatic detection of constant patterns in soil and crop properties over extensive spatial and temporal scales is important for a range of agricultural, ecosystem science, and natural resource management applications in order to plan interventions and to predict scenarios [1]. Satellite optical remote sensing nowadays provides suitable data for gathering information in a spatially- and temporally-continuous manner [2].

Multi-year time series of remotely-sensed images of vegetation and soils usually include a directional trend, seasonal regular fluctuations and additional irregular (random) short-term fluctuations due to local factors [2]. Depending on the research focus, it can be either the seasonal, short-term, or trend component which is of special interest—or all three simultaneously [3]. However, separating changes taking place due to soil characteristics and fertility and crop phenological cycles, inter-annual climatic variability, human activities (fertilization, irrigation, etc.), and long-term trends is challenging unless a data set of sufficient duration and frequency exists. In the last decade, some space agencies (e.g., NASA and USGS with Landsat 8–30 m of spatial resolution or ESA with Sentinels—10 to 30 m of spatial resolution) have adopted free access policies to mission imagery leading to a diffuse access to decametric optical data. In this way, they also triggered a special interest in long-term monitoring and in the development of creative, novel techniques to create robust multi-temporal products.

Satellite-based farm monitoring systems at the field scale (i.e., suitable for precision agriculture), which are rarely available (see e.g., [4,5]), are aimed at the optimization of agricultural practices for increasing crop production, reducing costs, increasing farm profitability, and allowing a more rational and sustainable use of fertilizers and agro-chemicals [5]. One of the most commonly practiced methodology for taking into account field spatial variability in precision agriculture is the identification of crop management zones as relatively homogenous sub-units of a field that can each be managed with a different, but uniform, agronomic practice [5,6]. In this application framework, mapping constant patterns in agricultural fields by exploiting satellite image archives has been investigated by several authors [7,8,9,10] and applied for different international projects (e.g., the crop yield monitoring system within the European Commission MARS project [11] and Global Agriculture Monitoring—GLAM [12]).

Several authors have also investigated uni- and multivariate geostatistical methods for the identification of uniform management zones, or sub-field units of similar yield potential, by using data fusion techniques [13,14] combining different data layers. Such layers typically include data acquired using ground-based proximal sensors, as well as yield maps and soil-sampling data, so the feasibility of such methodologies will depend on the availability of this kind of information, which is normally quite rare even amongst the most advanced and intensively managed farming systems.

Since cropping systems are very different in their within-fields spatial distribution and phenology, as well in their management, depending on farming structures and environmental conditions in different regions, an image-based approach relying only on freely available remote sensing data is appealing, in order to find a more flexible and general solution. In fact, freely-available data from routine acquisitions, with a medium revisit time (4–16 days) from spaceborne sensors, are considered to be more flexible as compared to commercial satellite acquisitions, which are typically on-demand and do not guarantee the availability of long time series. This is considered as a requirement for the type of multi-temporal clustering analysis proposed in this paper.

Cluster analysis relies on algorithms used for grouping pixels according to measured or perceived intrinsic characteristics or similarity, and not using a priori labelling. The absence of categorical information distinguishes data clustering (unsupervised classification) from supervised classification or discriminant analysis [15,16,17,18]. The aim of clustering is to group similar objects, starting only from data, thus, it is a trial approach to determine similar patterns (K groups or clusters) and to assign each data point to a cluster. Since clustering is an unsupervised procedure, the performance of clustering algorithms is usually assessed by validity indices that combine a measure of intra-cluster variability (to be minimized) and a measure of inter-cluster separation (to be maximized). Such validity indices could also be used to identify the optimal number of clusters for a given dataset. A general intrinsic problem of this type of unsupervised algorithms lies in the fact that the data may contain clusters of different shapes and sizes whose final meaning is generally unknown, for example the type of soil/vegetation cover [19,20].

Clustering and, in particular, K-means, which was proposed over 50 years ago, has been widely used in the last decades and a large number of clustering algorithms have been published. This proves both the effort in developing a versatile clustering algorithm and the ill-posed problem of clustering, i.e., to find the correct number of clusters [19,20,21]. As regards the application to remote sensing, Richards et al. [22] introduced a weighted method for clustering the individual units of geologic maps obtained from remote sensing images (Magellan Spacecraft radar images of the planet Venus) to provide planetary geologists a method to numerically test the consistency of a mapping with the entire multidimensional dataset of the study area. Ma et al. [23] in their experimental results on three different remote sensing images (Landsat TM and Quickbird satellite imagery, and AVIRIS airborne imagery) demonstrated that the use of a spectral-spatial clustering with a local weight parameter determination method can provide higher clustering accuracy in comparison to other clustering methods. In the study of [24], remote sensing vegetation indices were extracted from MODIS satellite data over the period 2000 to 2015 and then used to identify homogeneous hydrological watersheds in Iran. A fuzzy method (fuzzy c-means) was used to cluster the watersheds based on the extracted indices; results indicated three homogeneous regions, which are consistent with the variations of topography and climate of the study area and that will likely need similar management considerations and measures.

In recent years, a very popular alternative to the traditional K-means algorithm has been given by spectral clustering [25,26]. In this approach, the clustering problem is revisited as a partitioning problem for the graph representing the data points connected by their pairwise similarities. Then, the Laplacian graph is built and its first K eigenvectors computed. Conventional K-means on the rows of the eigenvector matrix finally gives the clustering. Several algorithms for spectral clustering have been proposed, according to different choices for the similarity distance and different normalizations for the Laplacian graph. All of them outperform the traditional K-means algorithm, allowing for non-convex clusters. Moreover, their computational complexity depends on the sparsity of the similarity matrix rather than on the size of the dataset. Recently, several clustering algorithms [27,28,29] have been proposed, improving upon the spectral clustering approach, and have been tailored for remote sensing data. Among these more sophisticated clustering methods, Rodriguez and Laio [27] proposed the Fast Search and Find of Density Peaks Clustering algorithm; Elhamifar and Vidal [28] applied the Sparse Manifold Clustering algorithm; Murphy and Maggioni [29] describes the application of the Diffusion Learning algorithm to high-dimensional remote sensing data. A broad survey of methods for determining K appears in [30,31], who have also proved in their study theoretical guarantees for the number of clusters K for a wide variety of data. Despite its many successful implementations, spectral clustering is not particularly suited to time series of multispectral images. The application of spectral clustering to such data is not so common in the literature: as is thoroughly discussed in [32], the risk is to disregard the inherent time structure of data while clustering. Moreover, the computational cost of this method is remarkably higher: without specific assumptions heavily constraining the sparsity of the similarity graph, the clustering procedure becomes unfeasible even for medium-sized (20,000 pixels) fields.

A sensibly different approach to clustering, though related to the K-means algorithm, is to shift to a functional data paradigm. Functional data analysis, indeed, is a very active topic in statistics since many applied problems are more naturally modeled in functional terms, even though only a limited number of observations could be available. Such functions represent the trajectories of a stochastic process and the clustering is performed on the whole curves, or on a suitable parameterization of them (e.g., basis function expansions, functional principal components). A general survey of the related literature can be found in [33]. Functional clustering has been successfully applied in different research fields and received particular attention from statisticians in the last decade ([34,35] among others). Indeed, the recent literature [36,37,38,39,40] describes several applications to environmental, climatological, as well as remote sensing data, so that exploring the performance of functional clustering algorithms on multivariate data is one of the purposes of the present study. Specifically, we applied the following clustering procedures to time series of satellite imagery: (1) standard K-means applied to multispectral reflectance data (MK); (2) functional K-means (MFK), where the clustering is performed on the smoothed time series after a suitable basis functions expansion, as in [33]; (3) multivariate functional principal components clustering analysis (FPCA), based on the multivariate functional PCA proposed by [34]; and (4) hierarchical clustering, based on the multilevel FPCA (MFPCA) proposed by [35,41,42], to analyze random functions observed for a large number of units—in our case the images pixels—at multiple subunits—in our case the different wavelengths.

The different approaches presented in this manuscript have been successfully applied to real data in several research fields, such as online auctions, spectrometry, cerebral vascular geometry, gene expression profiles, just to cite a few (see the thorough review by [34] and references therein). For environmental applications, Nguyen et al. [43] proposes a novel unsupervised clustering method, specifically developed for segmenting dense terrestrial laser scanning data in coastal marsh environments. The authors implemented unsupervised clustering with the well-known K-means algorithm by applying an optimization to determine the “k” clusters. The developed clustering method provided promising accuracy results for vegetated land cover segmentation in coastal marsh environments. In [44], remote sensing imagery and mobile phone positioning data were integrated to retrieve urban functional zones by hierarchical clustering obtaining good results.

The main objective of this study is to assess the capability of different clustering procedures in mapping constant patterns, at the field level, using as input multispectral satellite data, and, in particular, by the use of multivariate functional clustering for remote sensing time series data that was, in our knowledge, not applied before for agricultural studies.

Since the clusters obtained from long satellite time series can be considered as the representation of different constant potential within the field [1], thus, not affected by seasonal fluctuation, they should match the crop yield patterns for each crop season. Consequently, crop yield maps collected from 2007 to 2017 were used to compare the capability of each applied clustering procedure to detect useful constant patterns of crop vigor within agricultural fields. The polygon kriging method [45] was chosen to compare the differences between the detected patterns, thus evaluating if a significant difference between the yields of each class existed.

2. Materials and Methods

2.1. Study Area and Landsat Dataset

The study area is located in the Maccarese S.p.A. farm (lat. 41°52′18″N, long. 12°14′05″E, alt. 8 m above sea level) near Rome (Central Italy) (Figure 1), within a completely flat agricultural area, originating from land reclamation works carried out in the 1920s. The soil is classified as Cutanic Luvisol [46], with sandy clay loam texture, becoming more clayey towards the northeast of the site. The main crops cultivated in this area are maize, durum wheat, wheat, triticale, barley, carrots and broad beans. These type of crops indicate that there are two main active crop growth seasons in this area, which were considered for this study, from February to April for winter-spring crops, such as wheat, barley, triticale, rapeseed, and broad beans, and from June to August for summer crops, such as maize.

Within the farm, 94 fields of interest (FOI) between 10 and 20 ha were chosen for this study. Figure 1 shows the 94 FOI (in blue color) overlaid on the 08/07/2015 Landsat 8 imagery (R: 655 nm; G: 561 nm; B: 483 nm).

The optical dataset acquired on the study area in Italy is composed of 37 Landsat Enhanced Thematic Mapper Plus (ETM+) (with no data gaps on the FOIs) and 13 Landsat 8 Operational Land Imager (OLI) (GSD 30 m) clear sky (no cloud contamination) images covering all the 94 FOIs, acquired from 2008 to 2015. All the images were spectrally and spatially resampled retaining four spectral bands (blue, green, red, and NIR bands) with a spatial resolution of 30 m and an approximate scene size of 13 km north-south by 11 km east–west. This dataset was used as input in the processing chain for the different clustering procedures, which were applied to produce the constant pattern maps at the field level. It is worth mentioning that crop yield maps collected from 2007 to 2017 were used to compare the capability of the different clustering procedures to detect constant patterns of crop vigor within the selected FOIs. This is because the clusters obtained from satellite time series can be assumed as the representation of different constant potential within each FOI, thus not affected by seasonal fluctuation, and which should match the crop yield patterns for different crop seasons.

2.2. Landsat Time Series Processing Chain and Standard K-Means Clustering

The methodology, developed in IDL code [47] for processing Landsat time series, consists of three main steps (Figure 2).

The first step involves the setup of the Landsat time series (T₁–Tn, i.e., 2008–2015) for further processing. In detail, each image is first atmospherically corrected by using the FLAASH tool based on the MODTRAN4 transfer code implemented in the ENVI image analysis software version 5.0 (Harris Geospatial Solutions, Inc.; Broomfield, CO 80021, United States of America) [47,48]. Next, to ensure an accurate co-registration between the images and the FOI in vector format, each image was spatially resampled (bilinear method) to a higher spatial resolution corresponding to a ground sampling distance (GDS) of 5 m/pixel. The resampling step of the data is just a workaround to assure the correct positioning of the final map with respect to the actual FOI borders.

The second step is devoted to the extraction of the time series (T₁–Tm, where m indicates a subset of the n images pertaining to the time series) for each FOI (i), to be clustered in the third step, taking into account that two different datasets, one pertaining to soil and the other one to vegetation, should be extracted. In fact, only images of the active vegetation crop season were considered in this study, for the calculation of vegetation indices and the recognition of related biomass patterns within the same FOI. To this aim, the Normalized Difference Vegetation Index (NDVI), one of the most commonly used vegetation indices for estimating plant “greenness” or photosynthetic activity [49], was calculated for each FOI throughout time:

$$NDVI=\frac{NIR-RED}{NIR+RED}$$

(1)

where RED means the Landsat sensors’ red band (630–690 nm), while the NIR band is between 760 and 900 nm. The value of this index ranges from −1 to 1. The common range for green vegetation is 0.2 to 0.8 [46]. The selection of a FOI (at a certain date) as pertinent to the soil or vegetation time series was accomplished by calculating the percentage of pixels exceeding a selected threshold on NDVI, and setting another threshold on this percentage. This is done for each FOI over the whole Landsat time series. In particular, for the study area, in order to be sure to identify FOIs occupied by green vegetation, we set a rather high NDVI threshold, i.e., ≥0.6. When NDVI was higher than 0.6 for at least 70% of the pixels of a FOI, at a given acquisition time, the FOI was considered as covered by green vegetation. Conversely, when the NDVI was <0.2 for at least 70% of the pixels, the FOI was considered as occupied by bare soil (or senescent vegetation). This procedure allows to identify, for each image, three types of FOI: (i) green vegetation FOI; (ii) soil FOI; and (iii) mixed FOI. The latter two were excluded from further processing in the present study. It is important to remark that this procedure does not take into account the type of crop (e.g., crop rotation) occurring in the FOI, so that the index trend could include variability related to different FOIs use.

The third processing step consists in the application, to the “Vegetation Fields” layer, of specific indices, suitable to describe the variability at the field scale of the vegetation optical properties (e.g., NDVI, simple ratio index; [49,50]). The standard NDVI was used for this methodology. Furthermore, before applying the clustering procedure, on the NDVI layers, a Z-score spatial normalization [2,51] was used to improve the comparison of spatial patterns occurring within a given FOI, of images acquired on different dates/years. The following equation was used:

$$Z-scor{e}_{i,j}(t)=\frac{NDV{I}_{i,j}(t)-{\mu }_j(t)}{{\sigma }_j(t)}$$

(2)

where NDVIi, _j(t) is the value of the NDVI value of the pixel i, field (FOI) j at date t; µ_j(t) is the mean of the NDVI values of the FOI j at date t; σ_j(t) is the standard deviation of the NDVI values of the FOI j at date t. By using this normalization, the pixels of each FOI will have zero mean and standard deviation equal to one on each considered date, instead of having NDVI values that depend on the phenological period in which the image was acquired. This should improve the robustness of the identified patterns.

The last step of the procedure consisted of deriving, from the normalized NDVI “Vegetation Field” time series dataset, constant patterns at the field level, by using the standard K-means algorithm, a simple, but effective, commonly-used non-hierarchical clustering algorithm [21].

The K-means algorithm, implemented in IDL programming language (Harris Geospatial Solutions, Inc.; Broomfield, CO 80021, United States of America) [47], calculates initial class means evenly distributed in the data space, and then iteratively clusters the pixels into the nearest class by using a minimum distance method. In each iteration, class means are recalculated and pixels reclassified with respect to the new means, i.e., the within-cluster variation is used as a measure to form homogenous clusters. All the pixels are assigned to the nearest class, unless a standard deviation or distance threshold is specified, in which case some pixels may be left unclassified if they do not meet the selection criteria.

In this study, for the standard K-means classifier applied for this clustering procedure, we fixed a priori (as requested by the algorithm) the number of clusters at three for each FOI. The choice is reasonable taking into account the limited FOIs’ extents (10–20 ha) and their within-field variability. To this aim, for each FOI the average normalized NDVI (temporal series trend) class and the classes above and below the FOI average were calculated.

2.3. Clustering Methods for Multilevel Functional Data

A different clustering approach was also considered for comparison (right part of Figure 2). It uses as the input dataset the atmospherically-corrected and resampled Landsat images (T₁–T_n) obtained in the first step of the processing chain described in Section 2.2 (Figure 2). Specifically, for each date in the Landsat time series, the reflectance images corresponding to the four spectral bands centered at 0.482 µm, 0.561 µm, 0.654 µm, and 0.864 µm were extracted. Then, different clustering procedures (described in the following sections) were applied to this multispectral dataset. As is well known, the choice of the optimal number of classes is a crucial issue. Here, we decided to evaluate, for any classification method with a fixed number of classes, the silhouette index [52,53,54], a measure of how the objects are well matched to their own cluster and poorly matched to neighboring clusters. Then, we choose the classification corresponding to the optimum value of the silhouette index for each FOI.

We developed and tested the following functional clustering procedures: (1) standard K-means on multispectral reflectance data; (2) functional K-means; (3) multivariate functional principal components clustering; and (4) hierarchical clustering.

2.3.1. Silhouette Index

Several cluster quality measures, or cluster validity indices, have been proposed in the literature [15]. External indices compare cluster labels to externally supplied class labels (also known as gold standard data or ground truth), whereas internal ones are based solely on cluster similarity. In general, internal measures reward compactness (small distances within clusters) and separation (large distances between clusters). They mainly differ in the way these two aspects are weighted and related, and how sensitive they are to outliers. We chose to adopt the very popular silhouette index because, as thoroughly discussed in [16] where 30 cluster validity indices are compared, “there is a group of about 10 indexes [sic] that seems to be recommendable and Silhouette, Davies-Bouldin and Calinski-Harabasz are in the top of this group”. Moreover, as clearly stated in [55] the silhouette index is less prone to outliers than other internal indices, such as Davies-Bouldin index. The average silhouette width provides an evaluation of clustering validity, and might be used to select an ‘appropriate’ number of clusters [55]. This index is a normalized summation-type index: the compactness, or cohesion, is measured based on the distance between all the points in the same cluster while the separation is measured based on the nearest neighbor distance. For a generic pixel i, we define a(i) as its average distance from the other pixels in its own cluster and b(i) as the minimum among the average distances from the pixels in the other clusters. Then, the silhouette index of pixel i is given by:

$$s(i)=\frac{b(i)-a(i)}{max\left\{a\right(i),b(i)\}},$$

(3)

Thus, a s(i) value close to 1 indicates that pixel i is appropriately clustered. The average of the s(i) over all data of a cluster is a measure of how tightly grouped all the data in the cluster are, and the average of s(i) over all data of the entire dataset is a measure of how appropriately the data have been clustered.

2.3.2. Multivariate K-Means (MK)

First of all, we implemented a standard K-means algorithm on the multispectral reflectance time series Y_i(t) (with i = 1, 4 in our case), i.e., multivariate data, to check if the algorithm (MK) is capable of distinguishing between bare soil and vegetation, starting from the whole time series of reflectance data, and also to obtain a baseline to assess the improvements due to the introduction of the functional clustering techniques, described in the following subsections. Then, while the clustering methodology is still the traditional K-means, the differences with the algorithm described in Section 2.2 lie in the use of the whole time series, without eliminating images representing non-vegetated periods, and in the representation of each image pixel by four reflectance values at different wavelengths instead of a compact vegetation index, such as the NDVI.

2.3.3. Functional K-Means (MFK)

According to the functional clustering approach, for each pixel the multispectral reflectance time series Y_i(t) (with i = 1, 4 in our case) is considered as a discrete sampling of a multidimensional time trajectory, the latter being a realization of a stochastic process Y(t). Such functional data can be represented in finite dimensional space by approximation of the functions in a finite basis (generally, by using B-splines) or by dimension reduction techniques (functional PCA), also in a multilevel/hierarchical framework.

Inspired by the work of [35] we perform a B-spline smoothing of the multivariate data Y_i(t) as provided by the R package fda [53], resulting in multivariate smoothed functions Ŷ_i(t) defined on the whole time interval; to cluster these smoothed functions, the simplest choice is to obtain discrete samples of the Ŷ_i(t) by evaluating them on an equispaced time grid and then to apply the K-means algorithm to these discrete samples. Several distances can be chosen (e.g., Euclidean, L2, H1), to assess the role of the curves and of their first derivative in distinguishing different clusters. We adopt the L2 distance [53]:

$$d(Y^{(1)},Y^{(2)})=\sqrt{\sum _i{\int }_T(Y_i^{(1)}(t)-Y_i^{(2)}(t){)}^{2}dt}$$

(4)

where ${\hat{Y}}^{(1)},{\hat{Y}}^{(2)}$ are the multivariate smoothed time series for pixel (1) and (2), respectively, with the subscript i denoting their components; the integral is to be evaluated on the whole time interval T or on its discretization on the equispaced grid mentioned above.

2.3.4. Multivariate Functional Clustering (FPCA)

We also implemented a multivariate FPCA, based on the approach proposed by [42]. Classical multivariate FPCA concatenates the multiple functions into one, to perform univariate FPCA. This method should perform well on the considered dataset, since the components of multivariate functional data are measured in the same units and have similar variation. However, normalization can also be considered in the FPCA in order to account for differences in degrees of variability and in quantity among the components of the multivariate random functions. Specifically, we truncate the Karhunen Loeve expansion of the multidimensional smoothed functions Ŷ_i(t) at the first q terms and cluster on the principal component scores. Thus, FPCA reduces high-dimensional functional data into low-dimensional principal component score vectors, while retaining most information: the number of components kept, q, is such that the cumulative percentage of variance explained is greater than 0.9. Then, clustering is applied to the estimated scores.

2.3.5. Hierarchical Clustering (MFPCA)

Hierarchical clustering is a well-acknowledged methodology consisting of a series of partitions, obtained proceeding by successive fusions of subjects (agglomerative clustering) or successive separations (divisive clustering). In such a broad field, the multilevel functional PCA (MFPCA, [45]) was specifically designed for multilevel functional data, i.e., data with a natural hierarchy, for example multiple recordings of the same subject. Since the covariance of such data has a multilevel structure, this method is able to decompose the total functional variation into between-subject and within-subject variation via the functional analysis of variance. Then it can capture dominant modes of variations and reduce dimensions by conducting FPCA at both levels. Specifically, we implemented the hierarchical clustering proposed in [41] assuming two levels in our data.

We state that level 1 clustering refers to the average behavior across all dimensions, whereas level 2 clustering is related to dimension-specific behaviors. At level 1, two units are in the same cluster if their unit-specific means are similar in shape; while at level 2, two units are in the same cluster if their corresponding deviations from the unit-specific means are dynamically similar or they move together over time. Again, MFPCA reduces high-dimensional functional data into low dimensional principal component score vectors: at each level, only a few components are kept and clustering is applied to the estimated scores at this level. Obviously, all these procedures depend on the fixed number of K clusters. In this case study we chose K = 4, by constraining the procedure to distinguish bare soil (one cluster) from vegetation and, within vegetated FOIs, to discriminate, as for the standard K-means algorithm, three clusters representing different levels of crop vigor (mean class, above the mean, and below the mean).

2.4. Clustering Accuracy Assessment

The clustering accuracy assessment is not a trivial task, for example internal cluster validity indices could not be used—as they are part of the clustering procedure in most of the approaches, and external validity indices could not be applicable too. In fact, external validity indices ask for ground truth, which, in this circumstance, could not be found. The reason lies in the fact that clustering results, though very useful for crop management, do not represent an observable or measurable entity with an evidence—such as land cover or leaf chlorophyll content—it is a notional category, theoretically indicative of long-term yield potential/trend pattern. This is the reason why, taking into account this hypothetical relation, yield maps were used as an indirect representation of the desired output. This has to be considered an indirect information since: (i) yield is related not only to crop limiting factors, but also to agro-practices; and (ii) yield is affected by seasonality and the maps are available just for few years, making it difficult to generalize such data to long-term potential.

Therefore, the reliability of the clustering methods, for the assessment of within-field patterns of crop productivity, was tested in this work by using crop yield maps for fourteen FOIs of the Maccarese S.p.A. farm. In four FOIs, the yield maps were available for two years (FOI 10, 13, 24 and 54), whereas in all the others yield maps were only available for one year, resulting in eighteen validation test cases. Wheat yield data/maps between 2007 and 2017 were collected by a combine harvester equipped with a differential global position system (DGPS) and a grain flow sensor, providing spatially referenced yield data [56,57,58]. The raw yield data were refined by removing yield values obtained at the edge of the FOI, or in case the harvester cutting width was incorrectly set, or for those values collected when the combine harvester was on a low speed or without differential GPS data.

The refined yield data were used to evaluate if there are significant differences in yield between areas classified by the clustering methods. For this comparison we used the polygon kriging method [42], which allows to take into account the spatial correlation in the FOI, hence obtaining for each cluster (an irregularly shaped polygon) a spatially weighted average estimate (the expected value) and the standard deviation of the yield [49,50]. The polygon kriging method is similar to block kriging, but the former is applied to an irregular polygon, discretized into regular cells, instead of a square block. The average covariance function for each polygon is obtained by the computation of a weighted discrete summation of the point covariance function, where the weights are based on the proportion of the cell inside the polygon [45]. The differences between clusters (polygons) were determined by observing the expected values and their 99% confidence intervals (lower and upper limit) for each polygon.

We compared the different clustering methods observing the optimum number of clusters for each FOI assuming that this optimum corresponds to the maximum number of clusters detected among all the clustering methods, which are also statistically different in terms of crop yield.

3. Results

3.1. Clustering Results

3.1.1. Standard K-Means Clustering Results

The results obtained by applying the standard K-means clustering algorithm using as input the concatenated images extracted from the Landsat images time series (i.e., the vegetation field (i) dataset (T₁–T_m) in Figure 2), are presented in Figure 3a. For comparison, the results of the hierarchical MFPCA (Figure 3b), as applied on the spectral reflectance time series, is also shown. The choice is justified by the fact both methods rely on a fixed number of clusters for each FOI to be defined in advance: three clusters for the standard K-means and four clusters for the MFPCA.

The results obtained for the biomass clustering for each investigated FOI do not show temporally-structured results, since the K-means algorithm is not able to identify specific temporal behaviors, but just similar spatial behavior for each FOI across Landsat time series. It is evident from Figure 3 that although in some fields there is a general agreement between the results of the two clustering algorithms, in the identification of areas with crop vigor higher (or lower) than the average; in many cases the results are rather different.

3.1.2. Functional Clustering Results

All the clustering procedures from Section 2.3 have been applied to the complete Landsat spectral reflectance time series, as described in Section 2.1, after the atmospheric correction and spatial resampling of the images.

The clustering maps for each procedure have been produced and the MK, MFK, and FPCA functional clustering maps are presented in Figure 4, while maps obtained from MFPCA are presented in Figure 3b. In these figures, in order to simplify the comparison of the results, we choose to adopt a color scale ranging from red to green. Such a color scale was determined by post-clustering evaluation of the NDVI mean value of each cluster and a subsequent ordering of the clusters from the lowest (red) to the highest (green) NDVI. It should be stressed that NDVI was not used within the clustering algorithms, but just as a post-clustering for labelling purposes. For each procedure, the optimal number of clusters, according to the related mean silhouette index, is generally very low: for two thirds of the considered FOIs this number is two or three (specifically, 60 FOIs for MK, 65 FOIs for MFK and 68 FOIs for FPCA). Only for a few specific FOIs all the procedures suggest a higher number of clusters: for 11 fields the estimated optimal number of clusters is six, or higher. In general, FOIs borders are not identified as a single cluster. However, some spatial structures in the vegetation, such as stripes parallel to the FOI borders, are often highlighted and detected.

3.2. Ability of Clustering Methods to Capture Actual Yield Patterns

The differences between each cluster (polygon) of the validation FOIs were statistically evaluated taking into account the expected yield values obtained applying the polygon kriging on the 14 validation FOIs. The number of classes was set in advance for K-means (3 classes) and MFPCA (four classes). The classes identified by the two methods were significantly different in 16 validation tests for K-means and in 15 cases for MFPCA (

Table 1. The number of expected classes (E) detected by the five clustering methods (see text for the abbreviations) and the corresponding number of significant different classes (S) obtained from polygon kriging yield comparisons over the whole validation set. The last column sows the optimum number of classes (OPT), assumed as the maximum of S for each row.

FOI	Area (ha)	Year	K-Means		MK		MFK		FPCA		MFPCA		OPT
			E	S	E	S	E	S	E	S	E	S
10	9.0	2007	3	2	5	4	4	4	4	3	4	4	4
10	9.0	2010	3	3	5	4	4	3	4	4	4	4	4
13	15.6	2007	3	3	2	2	2	2	2	2	4	4	3
13	15.6	2010	3	3	2	2	2	2	2	2	4	4	3
24	15.6	2010	3	3	5	5	6	6	3	3	4	3	6
24	15.6	2016	3	3	5	4	6	5	3	3	4	4	6
25	2.8	2009	3	2	2	2	3	2	3	3	4	2	3
30	7.8	2010	3	3	2	2	2	2	2	2	4	3	3
36	9.6	2009	3	3	5	4	5	4	5	4	4	4	4
39	16.9	2009	3	3	4	3	8	7	7	6	4	4	7
54	20.1	2009	3	2	2	2	2	2	2	2	4	4	2
54	20.1	2010	3	3	2	2	2	2	2	2	4	4	3
27	15.3	2016	3	3	4	4	4	4	2	2	4	4	4
64	55.6	2017	3	3	8	8	5	5	2	2	4	3	8
77	41.1	2017	3	3	5	5	6	5	2	2	4	4	5
78	52.6	2017	3	3	2	2	2	2	2	2	4	4	4
83	45.4	2017	3	3	3	3	2	2	2	2	4	4	4
94	32.8	2017	3	3	5	4	3	3	3	3	4	4	4

). The other three clustering methods (FPCA, MK, and MFK) did not provide a constant number of classes. The number of detected and significantly different classes was the same for thirteen validation tests, both using MK and MFK methods, whereas using the FPCA method the number of significantly different classes matched the number of detected classes 16 times (Table 1).

We remark that “optimal” (OPT in Table 1) means here the maximum number of statistically different clusters that can be identified across all the considered methods. This is assessed by considering as optimal the highest number of clusters, but only after statistical verification that these clusters correspond to significant differences in yield maps. As Table 1 shows, in the fourteen validation FOIs, this number (last column in Table 1, OPT) ranges from two to eight, so that the methods that do not fix the clusters number in advance (i.e., FPCA, MK, and MFK) are clearly to be preferred. Indeed, while it is reasonable to generally assume three or four clusters, the data can suggest different scenarios on different FOIs.

Examining the validation results, the lowest class in terms of crop vigor is always localized close to the FOI edges when using the K-means method (Figure 3), except for some FOIs in the southern part of the Maccarese farm. For the other clustering approaches, where the number of classes was not set in advance, the low class (crop vigor less than field average) is not confined to the FOIs’ boundaries, and there are many FOIs where only two classes were detected (average field crop vigor and crop vigor higher than the field average). Consequently, the “low” class is not present (Figure 4). The MFK method, in comparison to the FPCA, seems to provide a spatial distribution more in agreement with the K-means, even though the results provided by MFK are more detailed, i.e., this method provided the best results in terms of “optimal” number of statistically different clusters retrieved. In fact, the number of statistically different clusters for the MFK method is often very close to the optimal number (Table 1).

As an example, to give an agronomical insight into the meaning of the clusters, we showed the comparison between recorded yield maps and the clusters detected by the five methods in six FOIs (Figure 5).

The yield map of the FOI 24 shows lower productivity in the northeastern, southwestern, and central parts of the FOI (Figure 5). The K-means method reflects the spatial distribution of the wheat yield, especially at the FOI boundaries, however, both MFK and MK methods identified the highest number of significantly different classes (Figure 6). In particular, MFK allowed to identify six clusters consistent with the yield map and statistically different among them (Figure 5 and Figure 6). The brown (cluster 2) and red (cluster 3) areas, respectively showed an expected value of 4.2 and 3.7 t·ha⁻¹, and these values are consistent with measured yields values within these clusters (ranging from 1.9 to 5.3 t·ha⁻¹). The yellow cluster is characterized by an expected value of 4.6 t·ha⁻¹, a value very close to the measured mean value within the FOI (4.8 t·ha⁻¹). The three green areas, in the northern and southern parts of the FOI, actually detected regions having the highest measured yield values. The other two clustering methods showed different problems: FPCA detected only three classes that do not describe all the spatial variability of the FOI, especially in the southern part of the FOI, while the clustering obtained by MFPCA did not match up well with the yield map.

Although the recorded yield values are very low in the FOI 30, this FOI is clearly divided into two main parts: high productivity in the northeast part, where the average measured yield value is 4.5 t·ha⁻¹, and low productivity in the southwest, showing an average yield value of 1.9 t·ha⁻¹ (Figure 5). The FPCA, MK, and MFK methods provided very similar results: only two clusters that are consistent with the two abovementioned yield patterns. The two classes detected by MFK showed an expected value of 2 t·ha⁻¹ for the yellow cluster and 4.6 t·ha⁻¹ for the green one, while for the FPCA and MK both the expected values of the two clusters are slightly higher than those estimated using MFK. The K-means identified a third cluster at the edges, but although all the clusters are significantly different, the expected yield values of the medium and high class are very similar: 3.3 and 3.4 t·ha⁻¹, respectively. The four classes detected by MFPCA seem inadequate, in particular, the two lowest classes are not significantly different between them and both have a very small dimension (Figure 5).

In the FOI 36, low yield zones are localized in the southwest and northeast (yield: 0.88–2.0 t·ha⁻¹), while in the central area higher yield values were registered (3.2–6.8 t·ha⁻¹). This spatial distribution is well described by the MFK, FPCA, and K-means clusters. However, the two low-productivity clusters of MFK and FPCA are not statistically different. The MFPCA method was able to distinguish four classes showing a reasonable visual match with the yield map, even if the lowest three classes are very close among them.

The optimal number of clusters in the FOI 27 is four; both MK and MFK methods provided the best correspondence between yield patterns and clusters (Figure 5). However, the clusters detected by K-means method also show a clear link with measured yield values, while the two other methods (FPCA and MFPCA) were not able to display the within-field variability. Observing the yield map, a low yield area can be clearly observed at the northern edge of the field that corresponds to cluster 1 of the MK and K-means methods. This cluster is adjacent to an area showing higher yield values, which matches cluster 3 of MK, MFK, and K-means methods, where the expected value for the three methods ranges between 4.3 (MFK) and 5.0 (MK and K-means) t·ha⁻¹. The southern part of the field is more productive, especially in the southwest section, which corresponds to cluster 2 for MK (6.1 t·ha⁻¹) and cluster 1 for MFK (6.2 t·ha⁻¹).

FOI 64 is a very large field, composed of many parcels, thus, the optimal number of clusters is quite high and it was detected by the MK method. The other clustering approaches underestimate the field variability, in particular FPCA, which provided only two clusters and MFPCA, which identified four classes, but only two have a relevant extension (Figure 5). The yield map of FOI 64 shows a narrow belt in the central area having low yield values, which were included in cluster 4 of MK (expected value 4.7 t·ha⁻¹). In the southeast part of FOI 64 the highest yield values can be observed (average value of 11.5 t·ha⁻¹), which correspond to clusters 5 and 6 of MK, having, respectively, an expected value of 13.4 and 10.7 t·ha⁻¹.

FOI 77 consists of two large parcels, which show different productivity levels: high yield values in the southwest part (mean yield value of 7.8 t·ha⁻¹) and low and medium values in the northeast part of the FOI (mean yield value of 4.7 t·ha⁻¹). The optimal number of clusters of this FOI is five, which was obtained by the MK and MFK methods. However, clusters 1, 4, and 5 obtained by the MK method have very similar expected values, which range between 6.3 and 6.5 t·ha⁻¹. Cluster 2 of MK and cluster 4 of MFK match the area with low productivity in the northeast part of the field, where the measured yield values range between 1 and 5 t·ha⁻¹. In this low productivity cluster the expected yield values are 3.1 for MK and 3.3 t·ha⁻¹ for MFK. The other clustering methods provided few clusters, which do not match the spatial variability displayed by the yield map.

4. Discussion

The present study illustrates the potential of the proposed clustering methodologies for applications in the area of precision land management. In this context, the identification of the most suitable methodologies for capturing the stable patterns underlying the spatial variability observed in soil and vegetation variables is a very active research area. Most of the approaches that have been proposed rely on complex multivariate geostatistical techniques of data fusion of ground-based data, e.g., from geophysical [56,59] and proximal sensors [13,14], which are not always available with remote sensing data. Therefore, there is a great interest in the development of methodologies relying only on remote sensing data, which could be rapidly developed into high-level spatial products delivered to land managers. However, in order to be trusted by users, the validity of the methodologies proposed for the identification of stable, clear, and manageable zones within a field needs to be assessed. This is always rather difficult, since the observed spatial variability is superimposed onto a generally even higher temporal variability. Additionally, the validity of the identification of uniform zones has to be considered in relation to the objectives of the user, e.g., choosing classes of uniform management.

In the literature there are several sophisticated clustering methods based on spectral clustering algorithms that are used for different applications such as those proposed by [27,28,29]. However, in the present work, we chose, as an objective assessment variable of the goodness of the compared clustering methodologies, the ability to identify from multispectral satellite time series statistically significant yield differences among the clusters. Now, the application of spectral clustering or diffusion learning algorithms to remote satellite time series is quite rare in the literature. The main reason, in our opinion, is that these methods are not well suited for handling sequential information, so that the risk is to disregard the time structure of data while clustering. Moreover, the computational cost of algorithms such as spectral clustering is sensibly higher than both the standard K-means and functional clustering. Our approach assumes a functional model for the data at hand (a time series of multispectral images). Such functions represent the trajectories of a stochastic process and the clustering is performed on the whole curves, or on a suitable parameterization of them. Then, the observations are supposed to be a finite sample of points extracted from these curves, which belong to an infinite dimensional space. To efficiently represent such curves, we adopt the coefficients of their basis expansion, as in the MFK method (where the dimension depends on the chosen number of basis functions), or also their principal components scores for the FPCA method. Both these two last methods reduce the dimension of the clustering problem: as a comparison, in the MK method, the dimension of the data is given by the length of the time series times the number of considered wavelengths. That dimension for our dataset is about 200.

Clustering methodologies results obtained for the chosen agricultural area in Italy using Landsat multispectral satellite time series as input can be summarized and discussed as follows: The standard K-means clustering, applied to NDVI Landsat time series in the first procedure, proposed in Section 2.2, generally seems to retrieve a cluster dimension and distribution that are in accordance with the FOI’s dimension and crop yield within-field spatial distribution of the study area. Although the K-means classes and the yield maps often match, setting the number of classes to three could be restrictive in large fields (e.g., see FOI 36 in Figure 5) and excessive in small fields (e.g., see FOI 30 in Figure 5). On the other hand, multivariate K-means (MK) is theoretically very similar to the standard K-means, but it was applied to the spectral reflectance Landsat time series and, therefore, not only related to the difference between near-infrared and red light (i.e., NDVI). Moreover, this method has the potential to exploit the information yielded by the whole multispectral signal (Landsat bands). In our experiments, the clustering differences between standard K-means and MK are mainly due to the freedom in the choice of the number of clusters. On the other hand, the MFK and FPCA clustering methodologies take a great advantage from the functional approach: data are smoothed and their dimension reduced before clustering. Consequently, these procedures are weakly affected by irregularities in the time series. MFPCA seemed to be a promising approach, because of the hierarchical nature of the information to be retrieved, i.e., a large spectral difference between bare soil and vegetation and a smaller difference among the different crop yield in the FOIs. However, the procedure is still to be refined and it is quite sensitive to the irregular temporal resolution of the time series of images; the clustering results on the considered dataset are definitely not satisfying with respect to the yield maps.

In general, as clearly reported in the Results section, the clustering approaches that do not require setting the number of cluster in advance (i.e., MK, MFK, and FPCA), have proved to be more suitable to describe the within-field variability as compared to the other methods (standard K-means and MFPCA). It has to be noted that the silhouette index, applied here to choose the optimal number of clusters, while being undoubtedly a reasonable proposal, has not always provided good results. Specifically, when applied to the FPCA method, it often suggested just two clusters, even in very large FOIs. The individuation of a priori criteria for determining the optimal number of clusters remains then, in our opinion, an open issue.

Another point to be stated is that functional clustering methods have shown to be able to detect constant patterns in time series data and to identify statistically different clusters in the reported experiments. Among these methods, MFK has shown the best performance, both in terms of statistical significance of the different clusters detected in the analyzed FOIs and in terms of clustering accuracy as evaluated by consistency with the yield maps used for validation. Another advantage of MFK clustering methodology is its ease of implementation and low computational cost with respect to other approaches mentioned in the Introduction section.

5. Conclusions

In this study, the capability of different clustering procedures in mapping constant patterns, at the field level, using as input multispectral satellite data, and, in particular, by the use of multivariate functional clustering for remote sensing time series data, has been assessed. Archived Landsat 7 and 8 time series have been considered to this aim. Crop yield maps collected from 2007 to 2017 have been used to assess the capability of each applied clustering procedure to detect constant patterns related to within-field crop vigor.

The main conclusion of our study is that the MK, MFK, and FPCA clustering approaches, which do not require setting the number of clusters in advance, are more suitable to describe the within-field variability as compared to the other methods (standard K-means and MFPCA).

In particular, functional clustering methods demonstrated: (i) as concerns the temporal dimension, to be less prone to the seasonal noise, and to be able to capture the long term signal, which is the information of interest; and (ii) as regards the spatial dimension, to be suitable to detect constant patterns in satellite time series identifying the optimal number (OPT values in Table 1) of clusters, with a good match with the yield maps’ patterns.

The MFPCA is an exception, as it appears to be quite sensitive to the uneven temporal distribution of satellite acquisitions, and since it asks for an a priori definition of the number of classes.

In fact, the possibility of adaptively choosing the number of clusters proved to be crucial. This is grounded in the fact that the differences in the granularity of field variability, as observed in the real parcels, are great, even in neighboring ones. Nevertheless, the selection of a priori criteria for determining the optimal number of clusters (in this study the silhouette index) remains, in our opinion, an open issue to be addressed in the future.

In conclusion, on the basis of the obtained results, the usefulness of multivariate functional clustering for the analysis of remote sensing time series for agricultural studies has been proved.

Among the compared approaches, MFK has shown to outperform the other methods, both in terms of statistical significance of the identified clusters—within the single FOIs, and of clustering accuracy—as evaluated by consistency with the yield maps used for validation. Moreover, MFK, with respect to other approaches, offers faster implementation and lower computational costs.

Future research will be advantaged by the free availability of a large amount of satellite time series data, such as Sentinel 2A,B (European Space Agency) multispectral data, which have operated simultaneously since 2017 (with a revisit time of five days and spatial resolution of 10–20 m).

In view of the results obtained by using Landsat bands for the MFK algorithm, instead of a spectral index involving two bands in the VNIR as for the standard K-means algorithm, current multispectral (Sentinel and Landsat 8) and future hyperspectral satellite sensors (e.g., PRISMA, GF-5, EnMAP, ECOSTRESS, VENµS, and Sentinel 10; [60]) will provide a further improvement in clustering methods and performance. The availability of new data with increased spatial, temporal and spectral resolution will certainly require further experiments on these clustering procedures to more robustly assess their performance.