Keywords: Unsupervised classification, k -means, land cover, multivariate normal density, spectral imaging data

Land cover classification of remote sensing data is a fundamental tool to study changes in the environment such as deforestation or wildfires. A current challenge is to quantify land cover changes with real-time, large-scale data from modern hyperor multispectral sensors. A range of methods are available for this task, several of them being based on the k -means classification method which is efficient when classes of land cover are well separated. Here a new algorithm, called probabilistic k -means, is presented to solve some of the limitations of the standard k -means. It is shown that the new algorithm performs better than the standard k -means when the data are noisy. If the number of land cover classes is unknown, an entropybased criterion can be used to select the best number of classes. The proposed new algorithm is implemented in a combination of R and C computer codes which is particularly efficient with large data sets: a whole image with more than 3 million pixels and covering more than 10,000 km 2 can be analysed in a few minutes. Four applications with hyperspectral and multispectral data are presented. For the data sets with ground truth data, the overall accuracy of the probabilistic k -means was substantially improved compared to the standard k -means. One of these data sets includes more than 120 million pixels, demonstrating the scalability of the proposed approach. These developments open new perspectives for the large scale analysis of remote sensing data. All computer code are available in an open-source package called sentinel.

Introduction

Monitoring environmental changes has become critical for many issues related to sustainable development. Deforestation, wildfires, and other land use changes have profound impacts on human activities, our environment, and biodiversity [START_REF] Pettorelli | Using the satellite-derived NDVI to assess ecological responses to environmental change[END_REF][START_REF] Newbold | Has land use pushed terrestrial biodiversity beyond the planetary boundary? A global assessment[END_REF][START_REF] Betts | Global forest loss disproportionately erodes biodiversity in intact landscapes[END_REF]. For instance, there is some evidence suggesting that pathogen outbreaks are linked to changes in land cover and particularly deforestation [START_REF] Jones | Zoonosis emergence linked to agricultural intensification and environmental change[END_REF][START_REF] Morand | Changing landscapes of Southeast Asia and rodent-borne diseases: decreased diversity but increased transmission risks[END_REF].

There are two main approaches to track on-going environmental changes: either by monitoring and measuring land uses directly in the field, or by remote sensing with satellites, aircrafts, or other airborne devices (e.g., unmanned aerial vehicles).

Although the second approach has some limitations compared to the first one, it has some definite advantages that cannot be matched by field data. In particular, satellites can cover the whole surface of the Earth with a frequency of a few days or weeks [START_REF] Li | A global analysis of Sentinel-2A, Sentinel-2B and Landsat-8 data revisit intervals and implications for terrestrial monitoring[END_REF][START_REF] Wulder | Land cover 2.0[END_REF]. Furthermore, the most recent satellites are equipped with high-resolution sensors which are able to record a wide range of information such as reflectance at different wavelengths, altitude, or temperature [START_REF] Fu | Remote sensing satellites for Digital Earth[END_REF]. During the last decade, there has been a remarkable increase in the resolution of these sensors. To illustrate this progress, the University of Twente maintains a database listing 334 satellites (some being out of service) and 396 sensors with a number of bands ranging between 1 and 16,921 (https://webapps.itc.utwente.nl/sensor/; accessed 2021-08-31). Among these sensors, 43 (11%) are indicated to have a resolution of one meter or less (until 1.25 cm), and 90 (23%) others are listed with a resolution between 1 m and 10 m. Spectral imaging sensors record electromagnetic waves and provide data in two broad categories: hyperspectral imaging (HSI) where reflectance is recorded for several hundreds of narrow bands (typically a few nanometres wide), and multispectral imaging (MSI) where reflectance is recorded for a few bands (usually less than 20) each with a width of few tens or hundreds nanometres. Both HSI and MSI usually record wavelengths beyond visible light (e.g., ultraviolet, infrared). During the last decade, a range of open-access remote sensing data have been made available (e.g., https://developers.google.com/earth-engine/; see also Guo et al., 2020). As an example of these developments, the Sentinel program is made of seven satellites currently in orbit around the Earth (https://sentinel.esa.int/). Two of them, Sentinel -2A and Sentinel -2B, are equipped with an MSI sensor which records reflectance in thirteen bands from ultraviolet (UV) to infrared (IR) including three bands in visible light [START_REF] Gascon | Copernicus Sentinel-2A calibration and products validation status[END_REF]. The Sentinel program stands apart from other similar programs because the data are available publicly in near-real time through the Copernicus datahub (https://scihub.copernicus.eu/). Each satellite covers the same location every two weeks, so the same location is potentially covered every week giving the opportunity to monitor environmental and land use changes at relatively high temporal resolution [START_REF] Li | A global analysis of Sentinel-2A, Sentinel-2B and Landsat-8 data revisit intervals and implications for terrestrial monitoring[END_REF].

One of the applications of spectral imaging data is to infer land cover and land use. Two types of approaches are used for this objective. In supervised classification methods, there is a reference sample with known land cover which is used to "train" the classification procedure in a first step, and the sample with unknown land cover is then classified in a second step. In unsupervised classification, there is no reference sample: classes or groups are defined following different criteria (see [START_REF] Wulder | Land cover 2.0[END_REF], for a recent review). As discussed below, both approaches have their respective advantages. For unsupervised classification, the k -means algorithm has been widely used in various contexts (see next section).

The objective of the present paper is to present a new method, called the probabilistic k -means, to analyse large-scale, spectral imaging data. The most important original feature of this method is to take into account variance heterogeneity among groups. Furthermore, a specific aim was to perform analysis of images with several millions of pixels in reasonable times. For instance, an image (or product) of Sentinel -2 (about 10,000 km 2 ) at a 10-m resolution has more than 120 million pixels.

The method is available in a computer package called sentinel (which also includes functions to query, manage, and download data from the Copernicus datahub). Before detailing the proposed methodological development, the next section presents a review of the recent literature on the applications of the k -means method to the analysis of remote sensing data. Four applications are then presented with two HSI data sets and two MSI data sets. The discussion gives further comparisons with previous contributions on unsupervised classification. The perspectives of current and future developments are also discussed.

Literature review

This review focuses on developments and applications of the k -means method in remote sensing data analysis published during the last ten years. Where possible, the sizes of the imaging data and the software used have been noted.

Several papers attempted to develop methods aimed to improve the properties of the k -means method. [START_REF] Galluccio | Graph based k-means clustering[END_REF] developed a method which assumes there are modes (areas of highest densities of observations) in the distribution of reflectance. These modes are found in the multivariate density space using the link lengths of a minimum spanning tree. Basically, the goal of their method is to initialise the centres of the k -means algorithm. They applied it to image data from Paris (512 × 521 pixels, 7 bands) and from Mars (300 × 120 pixels, 256 bands).

Another study found that the standard k -means algorithm usually performs poorly on HSI data [START_REF] Zhang | Comparative analysis of bat genomes provides insight into the evolution of flight and immunity[END_REF]. These last authors define the pure neighbourhood index (PNI) to perform neighbourhood-constrained k -means which adds steps to the iterations of the standard k -means with a weight function defined with the PNI. They applied this method to a 200 × 200 pixels image with 80 bands. [START_REF] Haut | Cloud implementation of the K-means algorithm for hyperspectral image analysis[END_REF] used the MapReduce computational framework to analyse two images from Indian Pines (145 × 145 pixels and 2678 × 614 pixels, both with 220 bands).

They programmed their analyses with Apache Spark for distributed computing and Python Scikit for the k -means. However, they did not assess the effect of different numbers of groups. [START_REF] He | Enhanced land use/cover classification using support vector machines and fuzzy k-means clustering algorithms[END_REF] showed that support vector machine (SVM), a supervised classification method, performs very well even with a small training data set. On the other hand, fuzzy k -means (FKM) was found to have a reduced usefulness with large data sets. These authors proposed a fusion of the two methods where the entropy is used to find the appropriate number of groups (see below for details about the use of entropy). They applied their method on two SPOT6 images (1982 × 1630 pixels and 2113 × 2151 pixels) each with six reference classes. Their analyses were implemented in ENVI and IDL (ver. 4.8). [START_REF] Zhang | Hierarchical feature learning with dropout k -means for hyperspectral image classification[END_REF] used an object-based approach defining a hierarchy from the pixels up to the scene. Their analyses used a combination of principal component analysis (PCA) on HSI images, k -means with drop-out, and SVM. The code was implemented in LIBSVM. They applied their approach to the Indian Pines data (145 × 145 pixels, 220 bands) and the University of Pavia data (610 × 340 pixels, 103 bands). They concluded that the drop-out k -means improves efficiency of the standard k -means with a small computational burden. They also demonstrated that the spatial information contained in the neighbourhood of pixels is useful, although their results did not relate this improvement with the identification of physical objects on the ground. Similarly, in another study [START_REF] Kavzoglu | An experimental comparison of multi-resolution segmentation, SLIC and K-means clustering for object-based classification of VHR imagery[END_REF] used k -means to perform image segmentation in a framework of object-based image analysis. They applied their approach to an image with 5000×3700 pixels and 8 bands. They found that k -means generally performs well for image segmentation using different specific algorithms. They implemented their computations with ENVI and MATLAB.

Image matching and indexing are also applications of k -means. [START_REF] Cao | Towards information-theoretic K-means clustering for image indexing[END_REF] used k -means to perform image indexing based on the Kullback-Leibler discrepancy.

They provided code in C++ and Matlab. [START_REF] Sedaghat | Very high resolution image matching based on local features and k -means clustering[END_REF] performed image matching using k -means in a second step to classify images into groups. They used MATLAB to implement their method.

Several papers used k -means to perform fine-scale spatial structure analyses. [START_REF] Kuo | Leaf segmentation based on k -means algorithm to obtain leaf angle distribution using terrestrial LiDAR[END_REF] analysed canopy structure by quantifying leaf angle distribution using a combination of k -means and an octree data structure: they analysed point cloud data (PCD), a kind of LiDAR (light detection and ranging) data which can reconstruct 3-D structures. The PCD were first split into octree subspaces so that each single octree unit contained no more than 1500 points. Each octree unit was then analysed with a standard k -means. Direct observations led these authors to infer that a leaf used between 500 and 1500 points, which helped them find the number of groups in the k -means analyses. [START_REF] Reza | Rice yield estimation based on K-means clustering with graph-cut segmentation using low-altitude UAV images[END_REF] used graph-cut and k -means to identify rice grains and estimate their sizes: they first applied k -means on the red-green-blue (RGB) image data after converting them to the Lab colour space, and then used a graph-cut algorithm to identify the rice grains. The best value of number of groups in the k -means analyses was found with the histogram method [START_REF] Kanthana | Automatic grayscale classification using histogram clustering for active contour models[END_REF]. The analysed images had 600 × 400 pixels. [START_REF] Wang | Realization of remote sensing image segmentation based on K-means clustering[END_REF] used k -means for image segmentation to identify roads from satellite images: the image data were converted from the RGB space into the the Lab colour space and then analysed with k -means fixing the number of groups to three (no information on image size was given). Some authors used k -means to quantify temporal changes from several images.

Kesikoglu et al. (2013) combined PCA with a fuzzy version of k -means called c-means

to analyse temporal changes from image differencing, so there were effectively only two groups in their c-means analyses. [START_REF] Lv | Novel land cover change detection method based on k-means clustering and adaptive majority voting using bitemporal remote sensing images[END_REF] used k -means with adaptive majority voting (AMV) to quantify change magnitude image (CMI). Their method starts from a "central" pixel, and builds a region around it. In a second step, a k -means analysis is done in the region with two groups (changed vs. unchanged pixels). In a third step, the region is extended with the AMV algorithm. They analysed four images ranging in size from 412 × 300 to 950 × 1250 pixels.

Overall, k -means is a widely used method in image and remote sensing data analysis; it is often used in combination with other data analysis methods (e.g., PCA). A remarkable diversity of approaches have been developed during the past decade most of them with different objectives. The sizes of the data are generally moderate, and very little open-source software has been contributed by these studies.

Methods

Data

Remote sensing data are usually arranged in a rectangular raster with variables associated with each pixel of the raster. These variables may be univariate (i.e., a single value is associated to each pixel) or multivariate. In this paper, we consider a multivariate setting where these variables are the values of reflectance measured in different wavelengths (the bands). In the present study we do not consider the spatial arrangement of the pixels in their respective rasters, so that the pixels are assumed to be independent. Therefore, the data under consideration below are denoted as X with the values of reflectance arranged in a matrix with n rows and p columns, where n is the number of pixels of the raster (i.e., the product of the number of rows by the number of columns of the raster), and p is the number of bands of the image. Measures of reflectance are usually more or less noisy [START_REF] Chavez | An improved dark-object subtraction technique for atmospheric scattering correction of multispectral data[END_REF][START_REF] Zhang | Characterization of Sentinel-2A and Landsat-8 top of atmosphere, surface, and nadir BRDF adjusted reflectance and NDVI differences[END_REF]. The exact values measured by the sensor depend on land cover and also on several factors such as the satellite or aircraft position, the time of the day, the atmospheric conditions, and so on.

Probabilistic k -means

The k -means method is a widely used, unsupervised classification procedure [START_REF] Hastie | The elements of statistical learning. Data mining, inference, and prediction[END_REF]. It requires specification of the number of groups (or clusters), denoted as K here, then the algorithm proceeds by assigning observations to a group depending on the distance to the group means [START_REF] Lloyd | Least squares quantization in PCM[END_REF]. If the values of these means are unknown (which is the most common case), some initial values are chosen randomly, the observations are assigned as explained above, the group centres are recalculated, and the whole procedure is repeated until group assignments are stable. The method works with multivariate data using a multivariate distance such as the Euclidean distance.

Standard k -means algorithms work well when within-group variances are homogeneous so that group assignments using distances are likely to be valid. However, when variances are heterogeneous, this is likely to result in misclassification of observations. Figure 1 shows a small simulated example with two groups each with 200 observations drawn randomly from two normal distributions with means 0 and 6, and standard-deviations (SD) 2 and 0.1, both respectively for each group. Even though the two means are very different, the large SD of the first group is likely to result in mixing of observations from both groups, and thus a k -means-based classification may be in error for these observations. The standard k -means indeed resulted in 15 misclassified observations in this case.

A solution to this problem is to rely on a probabilistic approach when classifying observations in the different groups. In the above simple simulated case, it is straightforward to apply this approach: after running a standard k -means classification, the means and SDs of both groups are calculated, then the probability densities are calculated for all observations using parameters of both groups: each observation is reclassified to the group for which it has the highest density. This can be represented graphically with a classification limit where the inferred density curves intersect (Fig. 1B). Note, on the other hand, that the limit for the standard k -means is defined by the equidistant point between the two group means. The reclassification procedure can be repeated until the overall classification is stable. In this simple case, a single iteration is enough and results finally in four observations misclassified.

This approach can be generalised to multivariate data using the densities of multivariate distributions. However, this requires estimation of a number of parameters which is likely to grow substantially with the number of variables. For instance, a multivariate normal distribution with p variables has 2p + p(p -1)/2 parameters: p means, p SDs, and p(p -1)/2 covariances. Therefore, the number of parameters is proportional to p 2 . A way to avoid having to estimate too many parameters when p increases is to first perform a PCA on the matrix X. PCA is usually used to perform dimension reduction in order to obtain a number of variables smaller than p that maximise the overall variation in X. Another property of PCA is that these principal components (PCs) are orthogonal: in other words, the coordinates of the observations (here the pixels of the image) on these PCs have zero covariances. We denote the matrix of these PC-based coordinates as Z. From a geometrical point of view, a PCA resulting in p PCs is a global rotation of the axes defined by the original p variables with the constraint that the covariances of the PCs are equal to zero. Therefore, this considerably simplifies the calculations of multivariate normal densities since it is now needed to estimate only 2p parameters (p means and p SDs) for each group of the classification.

Another crucial difference with PCA as commonly used in data analysis is that it is important here to not scale the original variables (i.e., divide them by their respective SD) before performing the PCA. If one of the variables has a large variance compared to the others, then it will contribute overwhelmingly to the PCA and will pull the overall variation in the data compared to the patterns from the covariances. This is the reason why variable scaling is usually recommended before running a PCA (e.g., [START_REF] Venables | Modern applied statistics with S (fourth edition)[END_REF]. However, the present goal is to discriminate groups with the calculated PCs where the overall variance is actually the consequence of the existence of these groups. So, in order to not erase this overall variance, the variables should not be scaled.

A word should be said about the choice of the form of the density distribution.

The present work assumes that the rows of Z (not X) follow a multivariate normal distribution. Furthermore, it is assumed that this distribution is non-homogeneous: its parameters (means and SDs) vary among the K groups and are assumed to be homogeneous within each group. In practice these assumptions may not be valid and other distributions may reflect more accurately the distribution of reflectance within each group. However, at the moment there is no theoretical or empirical justification for one distribution rather than another. Furthermore, the crucial point here is to assess variation among groups of different land cover and the normal distribution with its two parameters (mean and SD) may be flexible enough to accommodate such variation.

Selecting the number of groups

Likelihood-based information criteria

With unsupervised classification, there are two possible situations: the number of groups may be known a priori, or this number must be inferred from the data.

In the second situation, a parametric, probabilistic approach makes possible to use standard statistical tools such as Akaike's information criterion (AIC, [START_REF] Akaike | Information theory and an extension of the maximum likelihood principle[END_REF] which requires to compute the likelihood of the data. We must take care that group assignment is uncertain and has to be considered explicitly when calculating the likelihood function. Thus, we have to calculate the probability for the ith row of Z (z i ) using the parameters (means and SDs) estimated for group j multiplied by the probability that pixel i belongs to group j. These products are then summed over all K groups for each pixel. Finally, the log-likelihood is the sum of the log-transformed probabilities over all n pixels:

L = n i=1 ln   K j=1 fj × ξ(z i |μ j , σj )   , (1)
where ξ if the multivariate normal density function, fj if the estimated proportion of pixels in group j, and μj and σj are the estimated parameters for group j. There are thus 2pK + K -1 parameters estimated from the data: mean and SD for each column of Z and for each group, and K -1 proportions (since K j=1 f j = 1). We may now calculate the AIC:

AIC = -2L + 2(2pK + K -1). ( 2
)
The value of K resulting in the smallest value of AIC is to be preferred. Another criterion which can be used is the Bayesian information criterion (BIC) defined by [START_REF] Schwarz | Estimating the dimension of a model[END_REF]:

BIC = -2L + (2pK + K -1) × ln n. (3)
A simulation study is presented in the Supplementary Information which shows that both criteria are not robust to non-normality of the data. In particular, if the observations follow a uniform distribution and there is no heterogeneity (i.e., all observations are generated from the same distribution, thus K = 1), then both AIC and BIC will select a model with K > 1. Furthermore, in this situation the values of AIC and of BIC tend to decrease continuously when K is increased (see Supplementary Information for details).

Informational entropy

Another procedure for selecting the value of K is based on the principle of entropy [START_REF] Burrough | High-resolution landform classification using fuzzy k -means[END_REF]. This approach can be applied if there is a measure of uncertainty in the assignment of observations to the groups: in that case each observation is given a value of membership to each group with the constraint:

K j=1 m ij = 1, (4)
where m ij is the membership value of observation i for group j. The entropy, H, for a given value of K, is then calculated with:

H = - 1 n n i=1 K j=1 m ij × ln m ij . (5)
The value of K resulting in the smallest value of H gives the best description of the data. Membership values have been defined in the context of FKM [START_REF] Burrough | High-resolution landform classification using fuzzy k -means[END_REF][START_REF] He | Enhanced land use/cover classification using support vector machines and fuzzy k-means clustering algorithms[END_REF], but they can be adapted in a straightforward way to the probabilistic k -means developed here using the probability densities. Furthermore, the computation of the multivariate normal densities on a log-scale (see next section) makes possible to calculate H even when densities can reach very low values (see Supplementary Information).

Computational details and implementation

The overall workflow is summarised on Figure 2. The whole procedure was implemented in code written in the R and in C computer languages. The PCA was performed by singular value decomposition (SVD) which is faster and numerically more stable than the usual eigendecomposition [START_REF] Venables | Modern applied statistics with S (fourth edition)[END_REF]. With HSI data, it was observed that a relatively substantial number of PCs had nearly zero variance so that keeping all PCs made the computations much slower for no benefit: the number of PCs selected was set to keep at least 99% of the overall variance.

For MSI data, all PCs are kept. The coordinates on the p PCs are first analysed with a standard k -means using [START_REF] Hartigan | Algorithm AS 136: a K-means clustering algorithm[END_REF] algorithm which is particularly efficient and fast. The means and SDs are calculated for the p PCs and each group. The multivariate normal densities are calculated on a logarithmic scale which avoids numerical underflows and considerably simplifies the calculations (the overall densities are calculated with sums instead of products if full densities were used). Furthermore, the mathematical expression is factorised to avoid repeating redundant computations (e.g., the terms -ln( √ 2πσ j ) were computed once for all observations). These factorisations result in running times around 2.5 times faster than using the internal log-density function. Finally, the densities are evaluated separately for each pixel and only its classification is stored, avoiding to store all densities which would require an array of npK real values (amounting to 4.1 GB of memory with n = 3.3 × 10 6 , p = 13, and K = 12). Furthermore, this makes the overall memory requirement independent of the value of K. The running times are predicted to be proportional to n, p, and K (i.e., O(npK)). It was evaluated that a single iteration of the algorithm takes ≈ K 5 sec on a standard laptop with n = 3,348,900 and p = 13. On the other hand, the number of iterations required to reach convergence depends on the data: analyses of data sets with strong structure converge quickly (typically less than 10 iterations with K = 2), whereas if there is no structure convergence takes longer to reach.

The probabilistic reclassification is iterated until convergence. Furthermore, two stopping criteria have been defined: the maximum number of iterations can be fixed by the user (e.g., 200); or the procedure can be stopped when less than a fixed proportion of pixels are reclassified (e.g., if this proportion is zero, then iterations are stopped when no pixel is reclassified). This probabilistic k -means has been coded in a C routine called from R.

All code is available in a package named sentinel distributed on GitHub (https:

//github.com/emmanuelparadis/sentinel). Some code is also provided to use the standard k -means method in order to ease comparisons with the present method.

This package includes code to query the SciHub repository where the Sentinel data are stored.

Applications

Five data sets were analysed (Table 1). They are described in details below.

Hyperspectral Data: Pavia University and Okavango

Two hyperspectral data sets were considered. The Pavia University and Okavango data are two HSI data sets that have been preprocessed (http://www.ehu.eus/ ccwintco/index.php/Hyperspectral_Remote_Sensing_Scenes; accessed 2021-07-07). Both data sets are associated with reference data defined as "ground truth" with 9 and 14 classes of land cover, respectively (Tables S1-S2). Two analyses were performed with both data sets. First, the ground truth data only were analysed with the standard and the probabilistic k -means with K set equal to the known number of classes of land cover for each data set. The classification performance of each method was quantified with the overall accuracy as defined by [START_REF] Olofsson | Good practices for estimating area and assessing accuracy of land change[END_REF]. Because both k -means algorithms are unsupervised, the reference and inferred land cover values were matched with the Hungarian algorithm, a method which aims to maximise the values on the diagonal of a matrix, as implemented in the package RcppHungarian [START_REF] Silverman | RcppHungarian: solves minimum cost bipartite matching problems[END_REF]; the diagonal values of the matrix output were used to calculate the accuracy. Second, the complete data set was analysed with the probabilistic k -means using increasing values of K: the value of H was computed for each value of K and the final maps were drawn for both standard and probabilistic k -means using the value of K giving the smallest value of H.

Southern France

An image data taken by the satellite SPOT6 was analysed. The image was taken on 2019-06-27 above the South of France and had no cloud cover. The area of the image was estimated to be 3207 km 2 . A preliminary analysis of the CORINE land cover database over this area found that it is covered by 34 distinct land classes (as defined by the CORINE database). Out of these 34 classes, 17 were represented by less than 0.5% of the area, whereas 14 classes were represented by at least 1% of the area (Tables S3). The data were analysed with the probabilistic k -means with increasing values of K: the value of H was computed for each value of K and the final maps were drawn for both standard and probabilistic k -means using the value of K giving the smallest value of H.

Eastern Thailand

One area was selected in Thailand extending from N 14°27 56 to N 13°27 the period 2016-2018 identified twelve main land cover classes (Tables S4).

4 Results

Pavia University

The overall accuracies were 0.55, and 0.62 for the standard and probabilistic kmeans, respectively. The maps drawn with the ground truth data only show that some areas are not correctly identified with both methods (Fig. 3). However, some areas look more homogeneous with the probabilistic than with the standard k -means which is consistent with with the better overall accuracy of the former.

The smallest value of H was observed for K = 2, and then the value of H increased when the value of K increased as well (Fig. 4). However, two local minima of H were observed for K = 6 and K = 9. Figure 5 shows the overall maps inferred from the k -means methods assuming K = 9.

Okavango

The overall accuracies were 0.62, and 0.79 for the standard and probabilistic kmeans, respectively. The ground truth data were more dispersed than for the previous data so that the improvement accuracy of the probabilistic k -means was not so clearly visible (Fig. 6).

A pattern similar to the previous data was observed with respect to the relation between H and K: the smallest value of H was observed for K = 2, and then the former increased when the latter increased (Fig. 4). Local minima were observed for K = 6, K = 9, and K = 11. Figure 7 shows the overall maps inferred from the k -means methods assuming K = 11.

Southern France

The value of entropy was very low for K = 2 (H = 0.004) meaning that classification between the two groups was practically almost perfect (Fig. 8). However, a very substantial portion of the study area was covered by water which could readily explain this result. Indeed, as for the previous data sets, the value of H increased when K increased. However, a local minimum was observed for K = 15. The maps inferred from both k -means methods assuming K = 15 show some interesting differences (Fig. 9). Particularly, the coastal lagoons which were found to be covered with different classes by the standard k -means were all grouped in the same class by the probabilistic k -means (Fig. 9).

Eastern Thailand

The values of H varied with respect to K in the same way than for the previous data sets: the smallest value was observed for K = 2 and local minima were observed for K = 12 at the finest resolution (10 m) and for K = 15 at the coarsest resolution (60 m; Fig. 10). These two values of K were selected to infer the maps at their respective resolutions . Overall, the maps inferred with probabilistic kmeans show better delimitation of the fields compared to the results obtained with the standard k -means, particularly for the paddy fields on the west of the study area.

Computational efficiency

With 3,348,900 pixels and 13 variables, each iteration with K = 12 took around 2.4 sec. Therefore, 200 iterations (the default limit set in the code of sentinel) took 8 min. Furthermore, it was observed that in all cases, with either real or simulated data, there was convergence to a stable classification with no further reassignment.

In all cases reassignment was around 0.001% of the pixels after 200 iterations.

Discussion

The present work has contributed a new k -means method which appears as an improvement compared to currently available implementations with respect to three points: better accuracy, possibility to identify the number of groups, and ability to handle and analyse very large data sets. Each point is discussed below.

With both HSI data sets for which ground truth data were available, the method proposed here showed better overall accuracy compared to the standard k -means.

The improvement was particularly substantial for the Okavango data set. Although no ground data were available for the two MSI data sets analysed here, the maps suggest improved classification with the probabilistic method compared to the standard k -means. These results clearly suggest that the proposed probabilistic method has improved performance compared to the standard k -means for land cover classification of spectral imaging data. The fact that the assumption of homogeneous variance is relaxed in this method is certainly an important factor to explain this improvement (see further below).

The present results emphasise the importance of selecting the number of groups, K. The above review suggests that this issue did not receive a lot of attention in the recent literature. Although the information criteria presented above seem good candidates to select the best value of K in a probabilistic framework, this was not conclusive (see Supplementary Information). Clearly, the lack of statistical robustness shown by this approach is problematic and needs to be investigated further. On the other hand, the entropy-base criterion, H, proposed by Burrough et al. (2000, and previous references therein) appears a good alternative. However, some care must be taken when using it. It was observed that the smallest value of H was always obtained with K = 2 groups. This could make sense considering that spectral imaging data often show a strong discrimination between two broad classes of land cover (e.g., land vs. water, urban vs. vegetation), so that it is expected that classification with K = 2 yields essentially very good results so that all values of membership, m ij , are either zero or one. On the other hand, in all applications the values of H showed a local minimum for more realistic values of K. Therefore, it is suggested here that the entropy-based criterion is useful provided it is used within a range of realistic values of K (i.e., avoiding too small values).

The main feature of the approach adopted in this paper is to relax the assumption of homogeneous variance which underlies the standard k -means algorithm. The assumption of homogeneous variance is an important feature of the ISODATA method [START_REF] Ball | ISODATA, a novel method of data analysis and pattern classification[END_REF] which is derived from the standard k -means. [START_REF] Memarsadeghi | A fast implementation of the isodata clustering algorithm[END_REF] when implementing the ISODATA assumed that 'the clusters are wellseparated, that is, the probability that a point belonging to one cluster is closer to the centre of another cluster than to its own cluster centre is negligible.' As illustrated above, if the variances are homogeneous this is likely to result in misclassifications.

Interestingly, [START_REF] Memarsadeghi | A fast implementation of the isodata clustering algorithm[END_REF] made no parametric assumption on the distribution of the data within groups (or clusters). Indeed, if the groups are well separated and their variances small enough, there is no need to make any such assumption and the standard k -means algorithms are very likely to perform very well.

The DBSCAN method [START_REF] Ester | A density-based algorithm for discovering clusters in large spatial databases with noise[END_REF][START_REF] Li | DBSCAN clustering algorithms for non-uniform density data and its application in urban rail passenger aggregation distribution[END_REF] is another unsupervised clustering method related to k -means which takes into account the spatial configuration of the data as well as noise. However, the DBSCAN, although closely related to probabilistic k -means, is more complex and current implementations have been explored only with limited data sizes, typically with a few ten thousand observations [START_REF] Hahsler | dbscan: fast density-based clustering with R[END_REF].

There has been substantial research on applying the k -means method to the analysis of remote sensing data (e.g., [START_REF] Lv | Parallel K-means clustering of remote sensing images based on MapReduce[END_REF][START_REF] Pascucci | A Comparison between standard and functional clustering methodologies: application to agricultural fields for yield pattern assessment[END_REF] above review of the recent literature). Besides these applications to remote sensing, an approach has recently been developed to take noisy data into account in the context of quantum computing [START_REF] Kerenidis | q-means: a quantum algorithm for unsupervised machine learning[END_REF][START_REF] Khan | K-means clustering on noisy intermediate scale quantum computers[END_REF]: these proposals can be compared to the method proposed here in the sense that they aim to deal with overlapping clusters; however, they treat this issue quite differently. [START_REF] Ma | Spectral-spatial clustering with a local weight parameter determination method for remote sensing imagery[END_REF] proposed an elaborate method named spectral clustering which seems to outperform other classification methods. However, spectral clustering appears to be a computationally costly method and seems unfeasible even with a few ten thousands pixels [START_REF] Pascucci | A Comparison between standard and functional clustering methodologies: application to agricultural fields for yield pattern assessment[END_REF]. On the other hand, the method proposed in this paper is economical in terms of computations as it only requires evaluation of densities for each pixel and each group. For instance, Rodriguez and Laio (2014) developed a clustering method based on densities but also requiring calculation of distances among observations. More recently, [START_REF] Liu | L p -norm probabilistic K-means clustering via nonlinear programming[END_REF] proposed a method, with a name similar to that presented here, which is based on a probabilistic modelling of the standard k -means which is solved by numerical optimisation (see also [START_REF] Li | Probabilistic k-means clustering via nonlinear programming[END_REF]. However, similarly to the standard k -means, and by contrast to the present method, they assumed homogeneous variance among groups. These authors implemented their method in MATLAB and presented several applications with artificial and real data sets of modest sizes (several thousands observations or less)

with known number of groups. [START_REF] Richards | Weighted model-based clustering for remote sensing image analysis[END_REF] present a method that shares some similarities with the present one as it is based on the multivariate normal distribution. However, these authors proposed to maximise the log-likelihood function by expectation-maximisation [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm (with discussion)[END_REF]. Besides, they weight the contributions of the pixels to the likelihood function with respect to their spatial contiguity, which was not considered in the present work (but see the perspectives below).

In terms of running times, the probabilistic k -means has an attractive feature.

As reported above, the analysis of a complete Sentinel -2 data set covering more than 12,000 km 2 at a resolution of 60 m with 13 bands takes around 8 min. Additionally, three analyses with several million pixels (including one with more than 10 8 pixels) are reported showing how the method presented here is scalable to very large data sets.

In addition to the developments on k -means, the package sentinel presented in this paper adds to the software tools available for the analysis of Sentinel -2 data. [START_REF] Ranghetti | sen2r": an R toolbox for automatically downloading and preprocessing Sentinel-2 satellite data[END_REF] presented another package written in R, sen2r, to handle and manage Sentinel -2 data. By comparison, sentinel makes it possible to search, download, and manage products and data from all Sentinel satellites. Besides, the R environment makes it possible to read the different file formats used in Sentinel products thanks to the packages rgdal [START_REF] Bivand | rgdal: bindings for the 'Geospatial' Data Abstraction Library[END_REF] and ncdf4 [START_REF] Pierce | ncdf4: interface to Unidata netCDF (version 4 or earlier) format data files[END_REF].

These integrated tools have the potential to contribute to a software environment for time-series analysis of remote sensing data [START_REF] Gray | Consistent classification of image time series with automatic adaptive signature generalization[END_REF][START_REF] Cai | Enhancing MODIS land cover product with a spatial-temporal modeling algorithm[END_REF][START_REF] Gómez | Optical remotely sensed time series data for land cover classification: a review[END_REF].

The present paper aims at developing and implementing a fast unsupervised classification method to analyse multispectral data. The ultimate goal of this work is to be able to analyse large scale remote sensing data to infer changes in forest cover over large areas (see e.g., [START_REF] Hermosilla | Disturbance-informed annual land cover classification maps of Canada's forested ecosystems for a 29-Year Landsat time series[END_REF]Paradis, 2020b,a). The approach presented in this paper offers several perspectives of future development in several directions. Although k -means is basically an unsupervised method, it could be extended into a supervised method by defining known groups and evaluating the a priori distribution of reflectance. This poses some difficulties since it is difficult to find reference sites with relevant information to use as 'training' data. Another direction which is currently pursued by the author is to include spatial contiguity in the model. [START_REF] Richards | Weighted model-based clustering for remote sensing image analysis[END_REF] used this information to calculate weights in their likelihood function. Currently, an approach using edge detection with Prewitt's and Sobel's operators [START_REF] Wang | Evaluating edge detection through boundary detection[END_REF]) is under study.

Figure 1 :

 1 Figure 1: (A) Two normal distributions with mean and standard-deviation 0 and 2 (red) and 6 and 0.1 (blue). (B) Two hundred observations, shown under the x-axis, were simulated from each distribution in A. After a standard k -means classification, 15 observations were misclassified. After a probabilistic k -means, four observations remained misclassified. The curves show the densities inferred from the observations.

Figure 2 :

 2 Figure 2: Workflow of the probabilistic k -means for the analysis of spectral imaging data.

Figure 3 :

 3 Figure 3: Maps of the Pavia University data set considering only the reference data.

Figure 4 :

 4 Figure 4: Values of entropy (H) with different numbers of groups (K) for the Pavia University and Okavango data sets.

Figure 5 :

 5 Figure 5: Maps of the Pavia University data set.

Figure 6 :

 6 Figure 6: Maps of the Okavango data set considering only the reference data.

Figure 7 :

 7 Figure 7: Maps of the Okavango data set.

Figure 8 :

 8 Figure 8: Values of entropy (H) with different numbers of groups (K) for the Southern France data set.

Figure 9 :

 9 Figure 9: Maps of the Southern France data set.

Figure 10 :

 10 Figure 10: Values of entropy (H) with different numbers of groups (K) for the Eastern Thailand data set.

Figure 11 :

 11 Figure 11: Maps of Eastern Thailand data set at the 10 m resolution. Scales are UTM-based (in km).

Figure 12 :

 12 Figure 12: Maps of Eastern Thailand data set at the 60 m resolution. Scales are UTM-based (in km).

Table 1 )

 1 was analysed with the same procedure than for the Southern France data. With Sentinel -2 data, four bands are available at a resolution of 10 m, six at 20 m, and three at 60 m. Two data sets were built from this image: us-

	49 , and

ing the highest resolution bands (10 m, 4 bands) and using all bands aggregating the highest resolution bands at 60 m (13 bands). Similarly to the Southern France data, there was no ground truth data for this data set. An analysis of land cover data from the European Spatial Agency Climate Change Initiative (ESA/CCI; http://maps.elie.ucl.ac.be/CCI/viewer/index.php; accessed 2019-11-27) for

Table 1 :

 1 Main features of the data analysed in this study.

	Site	Number of pixels	Resolution (m) Bands
	Pavia University	207,400 (610 × 340)	1.3	103
	Okavango	377,856 (1476 × 256)	30	145
	Montpellier	89,161,101 (9799 × 9099)	6	4
	Eastern Thailand 120,560,400 (10980 × 10980)	10	4
	Eastern Thailand	3,348,900 (1830 × 1830)	60	13

Acknowledgements

The author is grateful to two reviewers for their constructive comments. The calculations used for data analyses benefited from the ISEM computing cluster platform. This is publication ISEM 2021-332.

Code availability

https://github.com/emmanuelparadis/sentinel.

Declaration of data availability

The data that support the findings of this study are openly available from http:// www.ehu.eus/ccwintco/index.php/Hyperspectral_Remote_Sensing_Scenes (Pavia University and Okavango), from https://sso.theia-land.fr/ (Southern France), and from European Space Agency's Copernicus SciHub at https://scihub.copernicus. eu/ (Eastern Thailand).

Declaration of interest statement

The author declares no conflict of interest.