Time series clustering using trend, seasonal and autoregressive components to identify maximum temperature patterns in the Iberian Peninsula

Time series (TS) clustering is a crucial area of data mining that can be used to identify interesting patterns. This study introduces a novel approach to obtain clusters of TS by representing them with feature vectors that define the trend, seasonality and noise components of each series in order to identify areas of the Iberian Peninsula (IP) that follow the same pattern of change in regards to maximum temperature during 1931–2009. This representation allows for dimensionality reduction, and is obtained based on singular spectrum analysis decomposition in a sequential manner, which is a well-developed methodology of TS analysis and forecasting with applications ranging from the decomposition and filtering of nonparametric TS to parameter estimation and forecasting. In this approach, the trend, seasonality and residual components of each TS corresponding to a specific area in the Iberian region are extracted using the proposed SSA methodology. Afterwards, the feature vectors of the TS are obtained by modelling the extracted components and estimating their parameters. Finally, a clustering algorithm is applied to group the TS into clusters, which are defined according to the centroids. This methodology is applied to a climate database with reasonable results that align with the defined characteristics, enabling a spatial exploration of the IP. The results identified three differentiated zones that can be used to describe how the maximum temperature varied: in the northern and central zones, an increase in temperature was noted over time, whereas in the southern zone, a slight decrease was noted. Moreover, different seasonal variations were observed across the zones.


Introduction
Climate change is an increasingly noticeable global problem that has a significant impact on society as a whole and on ecosystems. As a result, research on climate change has become increasingly important, especially studies relating to temperature, and has been expanding. As stated in the report of the Sixth Assessment of the Intergovernmental Panel on Climate Change (IPCC), the global mean surface temperature (GMST) increased by 1.1 °C between 2001-2021 and 1850-1900, accelerating its rate after the 1970s (IPCC, 2021). Moreover, in the very near future , the GMST may warm up as much as 0.25 °C per decade, according to the climate model predictions by Samset et al. (2020) and Tebaldi et al. (2021).
Although these numbers may seem low, the changes and effects are really remarkable, as global warming manifests prolonged droughts, heat waves and forest fires. For instance, over the last fifteen years (2003,2010,2015,2018), Europe has experienced extreme heat waves (Kuglitsch et al. 2010;Russo et al. 2015;Molina et al. 2020). According to Calheiros et al. (2021), weather and climate conditions such as high temperatures, moderate annual precipitation and prolonged dry spells are the cause of a large number of forest fires and burnt areas around the globe, particularly in areas of the Iberian Peninsula (IP). According to climate predictions, these conditions are expected to continue for the foreseeable future, with some projections predicting that more intense, prolonged, and frequent extreme heat events will occur in Europe in the twenty first century, with a higher impact on the IP and Mediterranean regions (IPCC 2014;King and Karoly 2017;Dosio et al. 2018;Vicedo-Cabrera et al. 2018).
Analysis of the temperature changes experienced by the IP, as well as the projections that have presented concerning this issue, will be a more manageable process if studies that were conducted to observe how these changes or variations in temperature have occurred in the IP are included. This will certainly be the case if these changes are defined by zones or sub-regions, since the temperature in the IP presents spatial variations strongly influenced by distance from the sea and complex orography, promoting a marked climate gradient from north to south (Lorenzo and Alvarez 2022). On the other hand, considering analyses that take into account temperature extremes will be more beneficial, since most studies in the climate field have focused on mean climate trend analysis, which does not tell us about unusual changes (Gebremichael et al. 2022).
The extreme temperature changes experienced in a geographical area or subregion can be understood by obtaining time series (TS) clusters of these temperatures and defining them in points or areas distributed over a geographical area. Since TS clustering is used to identify interesting patterns in TS data sets, finding TS clusters can be valuable in different domains, such as responding to anomalies or novelties, detecting discordance, integrating applications in dynamic change recognition in TS, and predicting, recommending, and discovering patterns (Aghabozorgi et al. 2015). Different studies on TS clustering (Warren Liao 2005;Rani and Sikka 2012;Aghabozorgi et al. 2015;Ergüner Özkoç 2021) agree that there are three main categories or approaches to TS clustering, depending on whether one is working directly with raw data, indirectly with features or characteristics extracted from the raw data, or indirectly with models built from raw data. Some studies of different domains have used TS clustering to define patterns or find matching TS (Keogh and Smyth 1997;Huhtala et al. 1999;Wang and Wang 2000;Möller-Levet et al. 2003;Huiting et al. 2006;Guo et al. 2008;Lee et al. 2010Lee et al. , 2018Li and Wei 2020;Shi et al. 2021).
This study is framed within TS clustering based on the approach of extracting features from data, and proposes a procedure to cluster TS by trend, seasonality, and main autocorrelations to ensure that patterns of change in maximum temperature (TMAX) can be identified for each zone in the IP during the period of 1931-2009. The novelty of this methodology is due to the decomposition of TS using singular spectrum analysis (SSA), a well-developed TS analysis and forecasting methodology whose applications have a wide scope ranging from nonparametric TS decomposition and filtering to parameter estimation and forecasting (Golyandina and Korobeynikov 2014). First, in this decomposition process, three components associated with the trend, seasonality and residual of the initial TS are reconstructed, allowing the parameters that describe these components to be extracted. Second, the representation of each TS is obtained from a feature vector generated on the basis of the calculated parameters, which allows the TS to be clustered using unsupervised learning algorithms, such as k-means (Hartigan and Wong 1979), k-medoids (Park and Jun 2009), hierarchical agglomerative (HA; Lukasová 1979) and Kohonen self-organising maps (SOM; Kohonen and Oja 1996), which are well-known and representative conventional algorithms that use the Euclidean distance. Finally, in the experiment on a climatic database, after comparing the clusters obtained with the different methods, a hybrid approach that combines HA and k-means, called hkmeans, was selected as a clustering algorithm to identify TS that are similar and follow a pattern. A set of 1776 points from a grid of 25 × 25 km 2 elaborated through spatial interpolation kriging by the "Servicio de Desarrollos Climatológicos" of the Meteorological Spanish State Agency was used. This grid includes points distributed throughout Spain, Portugal and the closest areas of the Atlantic Ocean and the Mediterranean Sea; for each point, a monthly TS of TMAX from January 1931 to 2009 is considered. The way in which the TS are represented here allows the clustering process to be optimised, since one of the advantages of implementing the SSA decomposition is the ability to eliminate the noise of the series (one of its main applications), together with the study of the spectral profile. In addition, costs and speed are improved with the reduction of the dimensionality. Therefore, SSA is used not only to decompose a series into a number of components to determine its spectral profile or to estimate its parameters, but also as a basis for recognizing characteristics and patterns of TS ensembles. The clustering results of this study made it possible to identify three differentiated zones: zone 1 is situated in the north of the IP, in areas with the lowest TMAX and a higher proportion of increase compared to the other areas; zone 2 is located more to the south, in areas with the highest TMAX, as well as a slight decline over an extended period; and zone 3 is stationed more towards the centre, in areas with intermediate TMAX, showing an increase over time. In addition, it was observed that the identified zones show different seasonal variations.

3
The remainder of this document is organised as follows. Section 2 describes the TS decomposition method, which uses SSA in a sequential manner. Section 3 proposes the new method for defining the trend, seasonality and autoregression patterns of TS, including extracting TS trend, seasonality and residual components and clustering and defining patterns. Section 4 presents the results of the method. Section 5 presents the main conclusions of the paper.

Sequential SSA decomposition method
In the following section, the theory of the sequential SSA decomposition method for extracting TS components is briefly presented.
SSA is a technique that is known for its application in TS analysis and prediction and has recently been used to analyse digital images and other objects that are not necessarily flat or rectangular and may contain gaps. SSA is a very unique and attractive methodology for solving a wide range of problems in different areas related to TS and digital images, as it naturally combines parametric and nonparametric techniques (Golyandina et al. 2018).
This technique is based on the singular value decomposition (SVD) of a specific matrix obtained from a TS and aims to decompose original TS into a small number of interpretable components, such as a trend that is smooth and slowly varying, with oscillatory components that are periodic, pure quasiperiodic or amplitude-modulated, and noise without any pattern or structure (Golyandina et al. 2001;Zhigljavsky 2011;Xiao et al. 2014;Golyandina and Korobeynikov 2014).
As stated in Golyandina et al. (2001), SSA consists of four steps: embedding, SVD, grouping and diagonal averaging. In some references, steps 1 and 2 of the generic SSA scheme are combined in the decomposition stage, whereas steps 3 and 4 are combined in the reconstruction stage. Several additive components of the original TS are obtained through SSA. In the following, the SSA method is presented formally.
Input: = t 1 , t 2 , … , t N , the initial TS, which is a one-dimensional N-order TS.

Step 1: Embedding
The so-called trajectory matrix is obtained through the equation X = T( ) , where T is a linear map that transforms the TS into a matrix of order L × K , and where L is an integer that is called the window length, 1 < L < N , and K = N − L + 1.
The set of all possible path matrices can be denoted as Hankel matrices, M (H) L,K , where all elements along the diagonal are equal. If N and L are fixed, then there is a biunivocal correspondence between the path matrices and the TS.
The trajectory matrix X constructed from lagged vectors, which are generated from the TS , can be represented in the following way:

2.2
Step 2: Decomposition of X into the sum of the rank 1 matrices The result obtained in this step is the following decomposition: where U i ∈ R L and V i ∈ R K are vectors, such that ‖U i ‖ = 1 and ‖V i ‖ = 1 for all i and i denotes nonnegative numbers.
If such a decomposition is performed using conventional SVD, the corresponding SSA method is "Basic SSA." In addition, the SVD of the matrix X is calculated via the eigenvalues and eigenvectors of the matrix S = XX T of size L × L . Here, 1 , … , d denotes the eigenvalues of the matrix S listed in decreasing order of magnitude , U d denotes the orthonormal system of the eigenvectors of the matrix S corresponding to these eigenvalues, considering that d = L . Since are elementary matrices of rank 1. Thus, the SVD of the trajectory matrix can be written as the following: is called an SVD eigenvector of order i and consists of the singular value (equal to √ ii ), an eigenvector U i (the left singular vector) and a factor vector V i (the right singular vector).

Step 3: Grouping
The input of this step consists of expansion (2) and specification of how to cluster the components of (2). The index set {1, 2,…,d} must be segmented into m disjoint subsets. I 1 , The resulting matrix X I corresponding to group I is defined as the following: Thus, if a partition is specified in m disjoint subsets of the index set {1, 2, … , d} , then, by expansion (2), the result of the grouping step leads to the following decomposition: The above procedure for choosing the sets I 1 , I 2 , … , I m is called the eigentriple grouping procedure. The grouping of expansion (2), where I j = {j} , is called elementary. (1) 1 3

Step 4: Reconstruction
In this step, each matrix X Ik from lumped decomposition (5) is transferred into the form of the input object , which is a TS of length N. To do this, each matrix X Ik is hankelised and, by means of one-to-one correspondence between Hankel matrices and TS, is transformed into a new series of length N . Thus, applying diagonal averaging to X Ik produces a reconstructed series ̃ k of order N (for more details, see Sect. 1.1.2.6 of Golyandina et al. 2018).
Consequently, the resulting decomposition of the input object is the following: If the grouping is elementary, the reconstructed objects ̃ k in (6) are called elementary components.
The SSA parameters, i.e., the length of the window L and the way in which X ik matrices are grouped, are very important for the outcome of the decomposition and depend on the properties of the initial TS and the objective of the analysis. One aspect that helps in choosing these parameters is the notion of separability. The separability of two TS, (1) N and (2) N , means the possibility of extracting (1) N from an observed sum (1) N + (2) N . According to Golyandina et al. (2001), SSA can approximately separate signals, noise, sinusoidal waves with different frequencies, trends and seasonality, etc.
It is possible to obtain recommendations for the choice of window length based on the (approximate) separability conditions. For example, the value of L should be large enough ( L ≈ N∕2) , and if extracting an existing periodic component of a TS with one or several known periods is needed, it is advisable to choose a window length that is proportional to the highest period. In Golyandina (2010), the choice of SSA parameters is discussed.
SSA can be performed sequentially, which is recommended when the TS structure is complex (Golyandina et al. 2012). Sequential SSA consists of two stages: in the first stage, the extraction of the TS trend with a small L is performed, and in the second stage, the periodic components of the residue are detected and extracted with L ≈ N∕2.

Trend, seasonality and autoregression SSA-based TS pattern identification
In TS data mining, feature extraction is one of the dimensionality reduction procedures. Features extracted from series concisely represent the relevant features of each TS as a finite set of inputs for a clustering algorithm that can discern the similarities and differences of TS (Wang and Hyndman 2006). In this paper, features from complete series, rather than subsequences, are extracted to look for complete TS that have similar patterns (i.e., similar trend, stationarity and autoregression patterns).

Pattern identification algorithm
The algorithm used to identify trend, seasonality and autoregression patterns in TS that is proposed in this study can be summarised in the following steps: 1 Perform sequential SSA to extract the three components of the initial series that are associated with trend, seasonality and residual. 2 Model the extracted series in such a way that their associated characteristics can be extracted.
2.1 Trend component: from trend = − t , estimate and using a linear regression, where trend is the extracted trend series and t is time. 2.2 Seasonality component: from seasonal = c 1 sin(2 t∕w) + c 2 cos(2 t∕w) , estimate c 1 and c 2 using linear regression, where w is the period, seasonal is the extracted seasonal component and t is time. 2.3 Residual component: from residual , obtain an AR(T) and calculate the autocorrelations 1 , 2 , 3 and w , where w is the period and residual is the extracted residual component. For each initial series, a feature vector is constructed by considering the estimated parameters.
3 Use a conventional clustering algorithm to obtain a similar TS. 4 Average the initial series of each group to obtain their representative patterns based on the defined characteristics.
The new algorithm is explained step by step below.

Extracting the trend, seasonality and residual TS
The series for trend and seasonality and those associated with the residual of the original TS are extracted using a sequential SSA-based decomposition approach from the TS trend, seasonality and autoregression pattern identification algorithm proposed in this paper. In sequential SSA decomposition, TS are decomposed into a small number of components (i.e., trend, seasonality and residual) in two stages (Golyandina et al. 2012). When the trend shape is complex, it is impossible to completely decompose the TS at once, resulting in the decomposition process being performed sequentially (Golyandina and Korobeynikov 2014). Thus, sequential SSA is applied to each initial TS to ensure that, in the first stage, the component associated with the trend is extracted by choosing an L value that is as small as possible, (i.e., divisible by the identified period w), to make the TS containing a periodic component smooth. In the second stage of sequential SSA, by considering the maximum L value ( L ≈ N∕2 and it is divisible by w ) for greater separability, the seasonality of the residual generated in the first stage can be extracted, generating a new TS for the residual part.
As noted in Sect. 2, the four steps of the SSA are highly connected, with each step immediately linked to the previous one, making each of them fundamental to the representation process proposed here, which starts with the extraction of the series of trend , seasonal and residual (all of order N ) described in this section. Thus, a correct identification of w and a correct definition of L are essential.
This process generates three TS for each initial TS, allowing them to be modelled and the parameters of interest to be calculated.

Constructing a feature vector for each TS
After modelling the extracted TS according to steps in 3.1.2 of the proposed TS trend, seasonality and autoregression pattern identification algorithm and to the estimated parameters, a feature vector for each TS is constructed.
Following the decomposition of TS i by means of sequential SSA to obtain the reconstruction of the trend, seasonal and residual components, which are modelled based on steps 3.1.2.1. to 3.1.2.3. of the proposed algorithm, the feature vector or representation of i , denoted as

Obtaining similar TS
A variety of methods have been developed to assess whether two TS are similar or in the same group. In this paper, the Euclidean distance is employed as a similarity criterion using the unsupervised learning algorithm hkmeans to group feature vectors that translate from the trend, stationarity and residual TS.
Given a dataset of TS with n points 1 , 2 , … , n , where each i is represented by a feature vector C i = i , i , c i1 , c i2 , i1 , i2 , i3 , is , the unsupervised hkmeans partitioning process P = P 1 , P 2 , … , P k , where the C i vectors are clustered according to the Euclidean distance as a similarity measure and are called characteristic-based or representation-based TS clustering. P i is called a cluster, where D = ⋃ k i=1 P i and P i ∩ P j = ∅, ∀i ≠ j. In the clustering process, each reconstructed C i vector has been normalised by centring to a mean of 0 and scaling to a standard deviation of 1 (Bro and Smilde 2003). Afterwards, a similarity matrix is obtained by taking into account these vectors and considering the Euclidean distance as the similarity measure. Finally, the similarity matrix is used as input to the clustering algorithms. Consequently, the clustering process is optimised, noise is eliminated, and computational costs are reduced.

Representative patterns of each group
After clustering with a cluster algorithm, the original TS of each group is utilised to calculate a new TS gi by averaging the TS of each group g i , with each gi serving as a representative prototype of group g i . Subsequently, steps in 3.1.2 of the proposed algorithms are applied to define the trend, seasonality and autoregression patterns.

Data preparation
A set of 1776 points from a grid in an IP of 25 × 25 km 2 elaborated through spatial interpolation kriging conducted by the "Servicio de Desarrollos Climatológicos" of the Meteorological Spanish State Agency was used. According to Fig. 1, this grid includes points distributed in Spain, Portugal, Southern France, North Africa and the closest areas of the Atlantic Ocean and the Mediterranean Sea, and considers a monthly TS of TMAX from January 1931 to December 2009 for each point, with 948 observations each.

TS analysis
As shown in Fig. 2, the TS of the monthly mean TMAX at the points located in such a grid of the IP are clearly very seasonal, with temperature peaks occurring in the summer months and temperature drops occurring in the winter months. Moreover, a variable trend is observed throughout the period, as shown in Fig. 6.

Decomposition and reconstruction of SSA
Sequential SSA was applied as a data preprocessing tool due to its capacity for component separability, which facilitates dimensionality reduction, TS representation and result interpretability. Figure 6 shows that the trend shape of the TS studied is complex, indicating that the decomposition is performed sequentially. First, the trend is extracted. For such a changing trend shape, its extraction is similar to smoothing (Golyandina and Korobeynikov 2014), starting with choosing a window length L = 12 to smooth the TS containing a periodic component, as would be done in the moving average procedure. The window length must be the minimum possible length and must be divisible by the period, which in this case is 12, as seen in the periodogram in Fig. 3, suggesting that the seasonality consists of sine waves with the indicated period. Since it is a monthly series, the horizontal axis of this graph varies from 0 to 6, and the peak at point 0 is related to a long-term trend, whereas the peaks at points 1 and 2 are related to the annual and semi-annual cycles, respectively. Generally, the same was noted in all TS in this study.
After performing the first decomposition, it is confirmed that the first eigentriple corresponds to the trend, whereas the other eigentriples contain high-frequency components, meaning that they are not related to the trend. This can be seen in Figs. 4 and 5, which show the shape of the six main eigenvectors and the result of the reconstruction carried out by each of the six eigentriples. Additionally, it can be seen that the principal eigenvector has practically constant coordinates; thus, it corresponds to pure smoothing by the Bartlett filter (Golyandina et al. 2012). Notably, the first reconstructed principal eigenvector in Fig. 5 produces the same exact reconstructed trend shown in Fig. 6.
After extracting the trend, the seasonality is extracted from the residual generated in the first decomposition stage. Figure 7 shows the periodogram of this residual, demonstrating that there is seasonality consisting of sinusoidal waves with periods of 12 and 6, and that it is not possible to observe the peak at the point 0 shown in Fig. 3 since the trend was removed. Since the length of the TS is N = 948 , to guarantee better separability, the window length L = 468 is taken as the maximum window length, ensuring that L ≤ N∕2 and L is divisible by 12.
To properly identify the searched sinusoidal waves, eigenvalue plots, eigenvector scatter plots and the correlation matrix of the elementary components ( W − matrix ) are used. Figure 8 shows several steps produced by approximately equal eigenvalues, each of which generates a pair of eigenvectors corresponding to a sine wave. This is confirmed in Fig. 9, which shows 2 almost regular polygons ET1-2 and ET3-4 corresponding to periods of 12 and 6, which occur due to seasonality and are clearly explained by the periodogram in Fig. 7.
The components are organised according to the order of the periodogram values at these frequencies (Golyandina and Korobeynikov 2014). The correlation matrix in Fig. 10 demonstrates that the indicated component pairs show high within-pair correlations, but between them, the correlation is almost zero. Some pairs of eigentriples that satisfy the same properties are observed, but they are referred to as noise because identifying the period to which they correspond is not interpretable for monthly data.
Thus, from the previous identification, the seasonality is extracted by reconstructing the ET1-4 clustering result. As shown in Fig. 11, the original TS corresponds to the residual of the 1st stage, TS F1 is the extracted seasonality and TS Residuals is the final noise or residual generated in the 2nd stage. The noise residuals obtained are heterogeneous (Golyandina and Korobeynikov 2014).
The resulting sequential SSA decomposition is shown in Fig. 12.

Parameter estimation for the representation
Considering the trend, seasonality and noise TS extracted through sequential SSA, the parameters that correspond to each one are estimated, as indicated in steps 3.1.2 of the proposed algorithm. The parameters estimated by modelling each component of the initial TS are shown in Table 1. The procedure is applied to the distinct TS considered in the dataset, and a feature vector is constructed for each TS.

TS clustering
When clustering methods are applied to any dataset (with random structures or not), the data is divided into groups. Thus, the clustering tendency of the new dataset obtained with the feature vectors of each TS is evaluated. The Hopkins statistic test evaluated the spatial randomness of the data by measuring the probability that the dataset was generated by uniform data distribution, enabling it to be used to evaluate the clustering tendency of a dataset (Cross and Jain 1982;Banerjee and Davé 2004). If the value of such a statistic is close to 0.5, the data is considered random, but if it is close to 1.0, it indicates that the dataset contains very well-defined clustered data. For the dataset, a Hopkins statistic value of 0.8550996 was obtained, which indicates that the dataset has natural clusters.
Since the dataset shows a clustering tendency, this study proceeds with clustering. Initially, four clustering algorithms, i.e., k-means, k-medoids, HA and SOM, are compared to the dataset. Internal validation measures of clustering are used.
Internal clustering measures the goodness of a clustering structure without taking into account external information that is not present in the data (Kremer et al. 2011;Song and Zhang 2008;Brun et al. 2007). This approach is usually based on compactness and separation criteria because clustering aims to arrange similar objects  (Kraus et al. 2011;Zhao and Karypis 2002). Three such internal validation indices are selected: the silhouette index (Rousseeuw 1987), the Dunn index (Dunn 1974) and the connectivity index (Liu et al. 2013). These indices suggest that among the applied algorithms, the best algorithm for clustering data into two groups is HA. However, although there are no objective criteria for choosing the number of clusters, when the popular silhouette (Rousseeuw 1987), elbow (Thorndike 1953) and GAP (Tibshirani et al. 2001) methods are applied, they estimate that the optimal number of clusters is 3. Looking at the information extracted from the selected internal validation indices with respect to three clusters, the methods to consider include HA and k-means clustering, as noted in Table 2. The connectivity index corresponds to the extent to which elements are placed in the same cluster as their nearest neighbours in the data space and should be minimised. In contrast, the Dunn index should be maximised, since, according to its calculation, a dataset has compact and well-separated clusters if the diameter of the clusters is small and the distance between clusters is large. The silhouette index should also be maximised because it measures how well an observation is clustered, in addition to estimating the average distance between clusters.
Considering what is suggested by these indices, in order to improve the clustering results, a hybrid method, called hierarchical k-means (hkmeans), was selected to cluster the TS according to their extracted features. The hybrid hkmeans method combines the two algorithms suggested by the indices in three steps. In the first step, it calculates the hierarchical clustering results and cuts the tree into k clusters. In the second step, it calculates the centre or mean of each group. In the third step, it calculates the k-means using the set of cluster centres defined in the previous step as the centre of the initial clusters (Kassambara 2017;Lee et al. 2010). Although k-means is one of the most popular clustering algorithms, it has some limitations; for example, the number of clusters needs to be specified in advance, the initial centroids need to be selected randomly, and the final clustering solution is sensitive to that initial random selection of centroids. However, for hkmeans, the selection of the initial centroids for k-means is defined using the hierarchical clustering result, which reduces these limitations and improves the final results of k-means clustering.
Here, hkmeans is applied to the set formed by the feature vectors of the TS. In the first step, an HA algorithm with Ward's linkage method is considered. The final clusters are shown in Fig. 13. To retrieve information on the groups obtained, the centroids of each group are accessed, which are the averages of each characteristic used to define the characteristic vectors that represent each cluster. Thus, Cluster 1 is identified as containing points in IP with lower temperatures than the points in other clusters, although the point values of this cluster vary greatly over time, increasing over time. The seasonal component describing these TS typically corresponds to a cyclical variation, and is described by c 1 = −4.62810 and c 2 = −5.20363 . In addition, the main autocorrelation of the residual component is negative. On the other hand, Cluster 2 contains points in IP where the average TMAX reaches the highest values, but the values decrease slightly over time, with a seasonal component whose amplitude of cyclical variation is defined by c 1 = −4.96816 and c 2 = −5.13317 on average. In addition, the main autocorrelation of the residual component is positive and has a higher absolute value than that of the other groups. Cluster 3 is characterised by points in IP where the average TMAX is at an intermediate level with respect to the temperatures of the other clusters and experiences a positive change over time. The seasonal component generally shows a cyclical variation described by c 1 = −6.07224 and c 2 = −8.18860 , and its main autocorrelation in the residual component, which is negative, is the lowest absolute value compared to the other clusters. Figure 14 illustrates the distribution of the points in IP according to the clusters obtained (the grid is composed of longitudes and latitudes in UTM coordinates).
As can be seen in the map in Fig. 14, Cluster 1 includes areas of northern Spain, northern Portugal and southern France, as well as the Mediterranean region. In Cluster 2, most of the points are distributed in the Mediterranean Ocean, with the remaining points in southern Spain, northern Africa and the Atlantic Ocean. In Cluster 3, the points are primarily distributed in the Spanish and Portuguese territories.
When analysing the resulting patterns using the series obtained from the centroids of each cluster, the three zones in the IP are clearly differentiated: in the north and central zones, an increase in temperature was noted over time, whereas in the south, a slight decrease was observed. The north zone of the IP, where the areas with the lowest TMAX are found, experienced a 0.2034 °C increase in its TMAX per decade between 1931 and 2009, whereas the central zone showed an increase of 0.135 °C per decade. In contrast, the south zone, where the areas with the highest TMAX are located, showed a slight decline.
Moreover, different seasonal variations were noted for each zone: the north zone shows its largest variations in winter months, whereas the central zone shows its variations in spring and autumn months. The south zone does not show any marked differences in monthly variation.

Conclusion
This paper proposes a novel method for clustering TS data that considers their trend, seasonality, and residual components. This approach involves representing TS data as feature vectors that are constructed by extracting the trend, seasonality and noise components using SSA decomposition. The results demonstrate reasonable groupings based on the defined features.
The proposed procedure can be applied to discover patterns in TS datasets, extract valuable information, and perform exploratory analysis on large TS datasets to support modelling efforts. The experiments allowed for spatial exploration and description of the variations of TMAX in the IP from 1931 to 2009 based on different zones defined by their trend and monthly variation.
Furthermore, this method could be used to test TS datasets with varying lengths or seasonal periods, as it is not restricted to TS data with uniform characteristics. Future research could explore the applicability of this method in multivariate TS analysis, as SSA can also be utilised for decomposition of such types of series.
material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.