Regional Spatial Mean of Ionospheric Irregularities Based on K-Means Clustering of ROTI Maps

Abstract

In this paper, we investigate and propose the application of an unsupervised machine learning clustering method to characterize the spatial and temporal distribution of ionospheric plasma irregularities over the Western African equatorial region. The ordinary Kriging algorithm was used to interpolate the rate of change of the total electron content (TEC) index (ROTI) over gridded 0.5° by 0.5° latitude and longitude regional maps in order to simulate the level of ionospheric plasma irregularities in a quasi-real-time scenario. K-means was used to obtain a spatial mean index through an optimal stratification of regional post-processed ROTI maps. The results obtained could be adapted by appropriate K-means algorithms to a real-time scenario, as has been performed for other applications. This method could allow us to monitor plasma irregularities in real time over the African region and, therefore, lead to the possibility of mitigating their effects on satellite-based location systems in the said region.

1. Introduction

Ionospheric electron density irregularities during the nighttime are crucial components of equatorial ionospheric variability studies. Research has revealed that the ionosphere over the Equatorial Ionization Anomaly (EIA) region is the most complex region with the highest variability, and its morphology is much more different from other latitudes. The effect of nighttime ionospheric plasma irregularities renders space-based systems almost unusable in these equatorial latitudes [1,2]. Ionospheric irregularities’ occurrence, as shown by different techniques, indicates a marked seasonal (with peaks during equinoctial months) and longitudinal variability. The Western low-latitude longitudinal sector of Africa, considered in this paper, has a higher occurrence of irregularities [3]. It has to be noted that the higher occurrence of irregularities during equinoctial months at low latitudes corresponds to the seasonal period of increased development of the EIA [4].

The rate of change of TEC index (ROTI), defined by [5], has been identified to be an important parameter to characterize the occurrence and magnitude of ionospheric plasma irregularities all over the globe. This index has been found to correspond to the ionospheric scintillation, a common radio signal phenomenon of the EIA region [6]. Chandra et al. in [7] used the ROT and ROTI to identify the signals severely affected by ionospheric gradients. Abe et al [8] analyzed the signature of plasma irregularities in an SBAS system’s performance in Central–West African equatorial and low-latitude regions through the ROTI. Among other results, they found that irregularities greatly contributed to the nighttime degradation of the SBAS system’s performance in the African equatorial and low-latitude regions, and 62.5% of the nighttime degradation could be associated with ionospheric irregularities over the region. The same authors demonstrated in [9] that the ROTI is a good parameter to estimate ionospheric irregularities during the storm-time period in the low-latitude African region.

However, effective evaluation of the ROTI over a region is challenging due to the scarcity of ground-based data for estimating ionospheric plasma irregularities. In some regions, the ROTI can only be obtained through the interpolation of sparse measurements. Although various interpolation algorithms are available globally, the accuracy of the estimated points tends to decrease as the distance between the pierce points increases. Therefore, the choice of interpolation algorithms and the quality of input data significantly impact the reliability of ROTI estimation. Predicting ionospheric parameters from ground-based GNSS measurements with enough accuracy is a key challenge in African ionospheric studies, particularly due to the irregularly spaced receiver locations. The Kriging interpolation technique appears to offer a solution to the estimation problem, based solely on the assumptions and knowledge of the variogram, resulting in a low residual error variance that guarantees confidence bounds [10]. Sayin et al [11] describes Kriging as an algorithm suitable for the spatial interpolation methodology in environmental geoscience. It estimates the variance and the spatial correlation structure of the surface from the measurements. Many researchers have applied the Kriging algorithm to estimate global and regional maps of ionospheric parameters (e.g., [12,13,14,15]). Regional ionospheric TEC and ROTI maps have been evaluated and observed to provide a reasonable corrective measure for ionospheric delays [16,17,18,19] and monitor the occurrence and severity of large-scale irregularities [20,21].

Machine learning is a branch of artificial intelligence that considers that algorithms can identify patterns and make decisions from datasets with little human intervention. Unsupervised machine learning is the process of finding hidden patterns in datasets. The starting point of the unsupervised machine learning process is the division of data points into groups of similar values called “clusters” that are different from the data of other groups. Representing data by a few clusters leads to data simplification.

The K-means algorithm [22,23] is an unsupervised machine learning technique that starts with an initial group of randomly selected centroids, which are used as the beginning points for every cluster, and then performs iterative calculations to optimize the positions of the centroids. It stops the process of creating and optimizing clusters when either: (a) the centroids stabilize or (b) the defined number of iterations is achieved.

K-means has been extensively used in geophysics in areas like weather and climate predictions [24,25,26,27,28,29] and seismology [30,31]. K-means has also been used to model the daily variations in TEC by Pongracic et al [32].

A reliable and accurate index to monitor ionospheric irregularities is a major challenge in the performance improvement of space-based systems in regions where irregularities are frequent. Therefore, this paper focuses on the provision of a spatial mean based on a technique for the optimal stratification of regional ROTI maps using post-processed data. Appropriate K-means-based algorithms could be used to test real-time data acquisition, monitoring facilities, and mapping techniques, as has been performed in other applications [33,34,35]. Achieving this could lead to possible mitigation effects of ionospheric plasma irregularities on satellite-based systems.

This paper is organized as follows: Section 2 covers the description of the data along with the ROTI, Kriging, and K-means techniques. Section 3 presents and discusses the application of these techniques for ROTI mapping and the proposed sampling method to obtain a spatial mean in a couple of cases. We finish by drawing a series of conclusions in Section 4.

2. Materials and Methods

2.1. Data and ROTI Calculation

To achieve the set objective, GNSS ground-based data were obtained from 22 stations located within the West African equatorial region through a series of public networks like the International GNSS Services, IGS, AFREF, NIGNET, and SONEL. The geographical descriptions of the stations are given in Figure 1 and

Table 1. GNSS stations used in this study with their coordinates and networks.

#Station	Station Code	Country	Geog. Lat. (°)	Geog. Lon. (°)	Network
1	BJAB	Benin	7.11	2.00	GPS-AFREF
2	BJCO	Benin	6.23	2.23	GPS-AFREF
3	BJKA	Benin	11.07	2.55	GPS-AFREF
4	BJNA	Benin	10.15	1.22	GPS-AFREF
5	BJNI	Benin	9.57	3.12	GPS-AFREF
6	BJPA	Benin	9.21	2.37	GPS-AFREF
7	BJSA	Benin	7.56	1.60	GPS-AFREF
8	OUAG	Burkina Faso	12.35	−1.95	GPS-IGS
9	YKRO	Cote D’Ivoire	6.87	−5.24	GPS-IGS
10	NKLG	Gabon	0.35	9.67	GPS-IGS
11	ACRA	Ghana	5.6	−0.20	GPS-IGS
12	CGGN	Nigeria	10.12	9.12	GPS-IGS
13	OSGF	Nigeria	6.92	11.18	NIG-NET GPS
14	ULAG	Nigeria	6.52	3.39	NIG-NET GPS
15	FUTA	Nigeria	7.2	5.3	NIG-NET GPS
16	RUST	Nigeria	4.8	6.98	NIG-NET GPS
17	UNEC	Nigeria	6.42	7.51	NIG-NET GPS
18	FUTY	Nigeria	9.35	12.5	NIG-NET GPS
19	BKFP	Nigeria	12.47	4.23	NIG-NET GPS
20	ABUZ	Nigeria	11.15	7.65	NIG-NET GPS
21	FG07	Sao Tome	0.35	6.74	SONEL GPS
22	DAKR	Senegal	14.75	−17.49	GPS-IGS

.

The ROTI was estimated over these stations in order to have a proxy for the plasma irregularities encountered by the GPS signals passing through the region. This parameter was quantified in each visible GPS satellite by taking the standard deviation of the ROT (time derivative of the TEC) at every 5-minute interval (Equation (1)).

$$\mathrm{R}\mathrm{O}\mathrm{T}\mathrm{I}=\sqrt{<{\mathrm{R}\mathrm{O}\mathrm{T}}^2>-{<\mathrm{R}\mathrm{O}\mathrm{T}>}^2}$$

(1)

where $\mathrm{R}\mathrm{O}\mathrm{T}={\displaystyle \frac{{\mathrm{T}\mathrm{E}\mathrm{C}}_{\mathrm{k}}^{\mathrm{i}}-{\mathrm{T}\mathrm{E}\mathrm{C}}_{\mathrm{k}-1}^{\mathrm{i}}}{{\mathrm{t}}_{\mathrm{k}}-{\mathrm{t}}_{\mathrm{k}-1}}}$ is the rate of change of TEC computed at every 30 s interval and then converted into TECu/min; i is the visible satellite; k is the epoch; and TEC is the total electron content. The ROT and ROTI were measured in TECu/min.

The ionospheric plasma irregularities were studied considering disturbed (2 October 2013; Am = 75) and quiet geomagnetic conditions (28 October 2013; Am = 4). The case of the moderate storm on 2 October 2013 was selected since it produced an enhancement of the irregularities in the region. The solar wind plasma speed, interplanetary magnetic field along the z-axis (IMF-Bz), and disturbance storm time (Dst) indices are shown in Figure 2 as the storm time parametric indices. The increment in the solar wind plasma speed on 2 October 2013 (DOY 275) from about 350 km/s to 525 km/s impacted the interplanetary magnetic field, and it fluctuated to around −9 nT. The Dst underwent a negative excursion to around −72. A more in-depth analysis of this storm can be found in [36,37].

Figure 3 shows an example of the ROT and ROTI variations throughout the storm day (2 October 2013, day of year (DOY) 275) for two West African stations: FUTY (left) and NKLG (right). The fluctuations evident after sunset reveal the presence of ionospheric irregularities and were observed by all visible satellites at both stations.

2.2. ROTI Map Estimation

To evaluate the spatial distribution of the ionospheric irregularities over the West African sector, the ROTI was mapped over the sector using the ordinary Kriging technique at a grid resolution of 0.5° by 0.5° longitude and latitude. We need to take into account that the critical challenge of ionosphere studies in the African sector is to predict the ionospheric parameters from the observed ground-based GNSS measurements at irregularly spaced locations. Therefore, the interpolation was a necessary step to perform the group division of the maps based on their similarity in the ROTI (clustering) in order to obtain a spatial mean of the irregularities’ occurrence in the region.

Blanch [17] reported that Kriging provides a solution to the problem of estimation based on the assumptions and knowledge of a variogram. In this work, we used the ordinary Kriging model as described in Blanch et al.; Rodríguez-Bouza et al; and Walter et al. [18,19,38] to estimate the ROTI at a nearby ionosphere pierce point (IPP) for any given point (Equation (2)).

$${\mathrm{R}\mathrm{O}\mathrm{T}\mathrm{I}}_{\mathrm{e}\mathrm{s}\mathrm{t}}={\sum }_{\mathrm{i}=1}^{\mathrm{n}}{\mathrm{w}}_{\mathrm{i}}{\mathrm{R}\mathrm{O}\mathrm{T}\mathrm{I}}_{\mathrm{I}\mathrm{P}\mathrm{P}\left(\mathrm{i}\right)}$$

(2)

where ROTI_est is the estimated ROTI at a given point; w is the weighting matrix, which is the coefficient that describes the planar fit trend; ROTI_IPP is the ROTI at the IPP; and i is the satellite in view.

$$\mathrm{w}=\left[\mathrm{W}-\mathrm{W}\mathrm{G}{\left({\mathrm{G}}^{\mathrm{T}}\mathrm{W}\mathrm{G}\right)}^{-1}{\mathrm{G}}^{\mathrm{T}}\mathrm{W}\right]\mathrm{c}+\mathrm{W}\mathrm{G}{\left({\mathrm{G}}^{\mathrm{T}}\mathrm{W}\mathrm{G}\right)}^{-1}$$

(3)

where W is the reciprocal of the covariance of all the measurements; G is the quantity that describes the geometry of the measurements; and c is the covariance of the vector and scalar field, which contains the characteristics of the planar trend and the decorrelation between the neighboring measurements. The scalar field determines the accuracy of the algorithm in terms of the multivariate Gaussian distribution, residual error analysis, and correlation distance between the neighborhood measurements.

$$\mathrm{W}={\sum }_{\mathrm{i}=1}^{\mathrm{n}}{({\mathsf{\sigma }}_{\mathrm{I}\mathrm{P}\mathrm{P}\left(\mathrm{i}\right)}^2+{\mathsf{\sigma }}_{\mathrm{d}\mathrm{e}\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{r}}^2)}^{-1}$$

(4)

where ${\mathsf{\sigma }}_{\mathrm{I}\mathrm{P}\mathrm{P}\left(\mathrm{i}\right)}^2$is the variance of the ROTI at each IPP, and ${\mathsf{\sigma }}_{\mathrm{d}\mathrm{e}\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{r}}^2$ is the decorrelation factor, which was estimated to be 35 cm using the planar fit for the WAAS at the middle latitude [38,39]. This value could not hold in the equatorial region, where the ionosphere is driven by many factors and the correlation with geomagnetic activity is not linear, unlike at middle latitudes [40]. In order to estimate the decorrelation factor, the approach in [38] and [41,42] was followed, with the relation found in [43]:

$${\mathsf{\sigma }}_{\mathrm{d}\mathrm{e}\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{r}}^2(el,{\mathrm{d}}_{\mathrm{i},\mathrm{j}})=\mathrm{k}$$

(5)

where el is the elevation angle, and d_i,j is the distance between the ionospheric measurements. The constant factor k that relates the ionospheric delay values among the measurements was averaged to 70 cm.

The variance was defined as:

$${{\mathsf{\sigma }}^2}_{\mathrm{I}\mathrm{P}\mathrm{P}}={\displaystyle \frac{{\sum }_{\mathrm{i}=1}^{\mathrm{N}}{\left({\mathrm{R}\mathrm{O}\mathrm{T}\mathrm{I}}_{\mathrm{i}{\mathrm{p}\mathrm{p}}_{\mathrm{i}}}-{\mathsf{\mu }}_{\mathrm{R}\mathrm{O}\mathrm{T}\mathrm{I}}\right)}^2}{\mathrm{N}}}$$

(6)

where ${{\mathsf{\sigma }}^2}_{\mathrm{I}\mathrm{P}\mathrm{P}}$shows the accuracy of the ROTI map over the study area, ${\mathsf{\mu }}_{\mathrm{R}\mathrm{O}\mathrm{T}\mathrm{I}}$ is the mean of the ROTI over the study area, N is the IPP number, and ${\mathrm{R}\mathrm{O}\mathrm{T}\mathrm{I}}_{\mathrm{i}{\mathrm{p}\mathrm{p}}_{\mathrm{i}}}$ is the ith element in the set.

In addition, we set the maximum search radius distance and the minimum number of IPPs to predict a value at a 0.5° × 0.5° grid resolution to 600 km and 8, respectively, against the 2100 km and 10 used generally at middle latitudes [44,45].

Figure 4 shows an example of the ROTI at an IPP over the West African region at a given epoch (2 October 2013 2000 UT) on the left side and the resultant interpolated ROTI map after applying the Kriging technique on the right. It can be noticed that the IPPs around the DAKR station in Figure 4 were not considered in the interpolation. For the reasons explained above, the regional maps of the ROTI were further reduced to the densest area of IPPs.

2.3. Optimization Sample Technique

The majority of ionosphere variable indices are determined through regular grids. However, in cases where the area lacks a regular pattern of sample locations or measurements, this solution might be too restrictive. Since the maps of the ROTI in the region are irregularly shaped in each epoch, we applied an optimization sample method based on an irregular grid. This approach is widely used in the survey and monitoring of natural resources and presents some advantages, one of them being to avoid the border effects due to the increase in Kriging variance. Moreover, an irregular pattern was expected to lead to a more precise spatial mean [46]. This process consisted of splitting the study area into clusters or cells as compact as possible using an unsupervised machine learning technique, as described in the Introduction.

For this task, we used a specialized R library, spcosa [47], to implement spatial coverage sampling and, thus, obtain the spatial means of the mapped ROTI in the study region. The spcosa package made use of a well-established method for geographical partitioning and cluster analysis, K-means, to find compact clusters in a multivariate property space. The clusters were represented by their centroids, and the method aimed to minimize an objective function (Equation (7)), which is the mean shortest squared distance (MSSD):

$$\mathrm{M}\mathrm{S}\mathrm{S}\mathrm{D}={\displaystyle \frac{1}{\mathrm{N}}}\sum {\mathrm{m}\mathrm{i}\mathrm{n}}_{\mathrm{j}}\left({\mathrm{D}}_{\mathrm{i}\mathrm{j}}^2\right)$$

(7)

where N is the number of raster cells, and (D²_ij) is the minimum squared distance between the ith cell and all the cluster centroids.

K-means proceeds in an iterative way, starting with an initial solution, which can be a random configuration of cluster centroids. Then, 2 steps are repeatedly alternated: one for the reallocation of raster cells from one cluster to another, and another for the recalculation of the coordinates of the centroids. The iteration process stops when the MSSD function cannot be lowered.

3. Results and Discussion

3.1. Regional ROTI Maps under Different Geomagnetic Conditions

In this section, we present some results of the performance of the application of ordinary Kriging geospatial interpolation over the West African sector. The nighttime values of the ROTI are used to showcase the validity of the interpolation technique and monitor the distribution of the electron density over the region under study during two different geomagnetic condition events. Figure 5 and Figure 6 show examples of the time regime evolution of the plasma irregularities over the West African region using hourly regional ROTI maps (left column) and their variance (right column) during two geomagnetic conditions on 2 and 28 October 2013 with different signatures. The variance maps indicate the performance of the interpolation technique over the region. The complete set of ROTI maps every 5 min for this case is provided as a movie in the Supplementary Materials.

Figure 5 presents the results of the ordinary Kriging interpolation of the instantaneous ROTI values in the West African regional maps for 2 October 2013 (DOY 275). The figure shows the non-uniform distribution of the ROTI over the study area on the left, indicating the presence of nighttime plasma irregularities in yellow–red from 1900 UT, while the magnitude of the irregularities was stronger during the hours of 2000–2200 UT, as indicated by the higher ROTI values. In general, as seen from the variance maps in Figure 5 (right), it is observed that the prediction values corresponded well to the measured values in nearly all the interpolated areas considered. There was an exception in some parts at the border near the sea, which were far away from the measurements, where high values of variance were observed, indicated by the red regions.

In the same vein, Figure 6 shows the results of another test case that shows the performance of ordinary Kriging interpolation on the regional map. The ROTI data obtained during the quietest geomagnetic conditions on 28 October (DOY 301) 2013 were used. The results show the presence of plasma irregularities from the southern crest of the anomaly that spread westward. The pattern of the irregularities insulating the eastward side of the region was validated with the variance maps (Figure 6, right column). Also, in this case, low variance values near the IPPs where interpolation was performed can be seen. The histogram in Figure 7, corresponding to the variance values of all the ROTI maps every 15 min for day 301 exhibits mean, standard deviation, and 99% and 75% percentile values of 0.2, 0.63, 0.53, and 0.65, respectively. This confirms that the Kriging interpolation algorithm’s modeled values may be accurate compared with experimental values.

3.2. Optimum Sampling Applied to Test Cases

When the K-means clustering algorithm is applied, there exists a trade-off between the number of clusters (k-value) and computational efficiency. Increasing the k-value enhances the capacity to discern underlying structures within the data yet, concurrently, it escalates the computational complexity and time. In contrast, diminishing the k-value streamlines the computation but may result in the omission of significant patterns. Therefore, striking an optimal balance between the number of clusters and computational efficiency is crucial when selecting the k-value. Even though various methods exist to find an optimal number of clusters, in our application of the K-means method, we decided to apply a visual inspection of the resulting clusters for different k-values based on the data structures and the appropriate cluster size for a region with a radius that varies from 500 to 600 km. In our test cases, we determined, after some trials, that stratification into 30 strata could represent a good approximation of the structures of the irregularities in the studied region in order to use their centroids as a regional index. It has to be taken into account that the number of clusters determined the computation time of the algorithm.

Figure 8 (right) shows the spatial coverage sampling for the test case on 2 October 2013 (DOY 275) at 20.20 UT associated with the ROTI map (left). The centroids of the clusters found after the stratification were taken as the sample points and the corresponding ROTI value for each one was obtained.

After the stratification and sampling, the locations and values of the centroids were reported in

Table 2. Coordinates of the cells in Figure 5 with ROTI values at centroids.

#Cluster	Lon.	Lat.	Roti
1	1.0406142	5.4264093	0.98796
2	9.0656952	5.4496616	0.89953
3	0.7574102	0.7424731	0.24533
4	13.9909615	12.8552702	0.32832
5	1.84456	15.9886426	0.28087
6	9.4122007	17.0629056	0.33742
7	12.86602	6.1328875	0.94803
8	7.8254299	1.4123836	1.1094
9	−2.4998058	3.88543	0.4506
10	−2.5378407	8.3695419	0.2559
11	−6.6016795	9.9428645	0.77241
12	9.5246104	13.2109367	0.3707
13	−4.6222565	12.7082078	0.35477
14	10.3590898	9.4368107	0.29737
15	6.5936256	10.0758751	0.62692
16	6.0172696	15.4801344	0.54801
17	5.4927833	−2.122805	0.17332
18	−6.4449116	6.3504329	0.97968
19	11.6189996	1.7535332	1.3688
20	−0.5090778	11.9470749	0.37997
21	4.0738014	2.2208489	0.67802
22	5.3790547	6.0552358	0.53311
23	2.1775499	8.9743649	0.9368
24	10.3294307	−2.5482636	0.55764
25	12.4568628	16.0298895	0.47352
26	3.6786224	12.5873912	0.90425
27	−5.8302737	2.1316645	1.0141
28	−4.5715688	−0.8569324	0.3142
29	14.2820471	9.4893215	0.30181
30	−1.9906776	15.3468051	0.2656

. An estimation of the spatial mean in the area and standard deviation was obtained from the stratified samples. In the case of the epoch of DOY 275 of the year 2013 in Figure 8, the ROTI spatial mean was 0.59 TECu/min, and the standard deviation was 0.32 TECu/min.

Similarly, in Figure 9 and

Table 3. Coordinates of the cells in Figure 6 with ROTI values at centroids.

#Cluster	Lon.	Lat.	Roti
1	5.33433035	17.935405	0.25196
2	11.44710797	9.1741088	0.26178
3	10.34243442	5.8323767	0.31078
4	−1.85233973	13.635944	0.4807
5	−2.71429565	0.4553038	0.25837
6	−0.11120521	10.2096185	0.38294
7	1.4351383	13.4653392	0.21679
8	6.47425093	−2.8970902	0.15069
9	1.64522135	17.1000017	0.34119
10	2.02938754	0.8038183	0.50119
11	12.1004933	12.6976959	0.16265
12	10.18263684	1.8554506	0.43152
13	−0.06530547	3.510308	0.84188
14	−3.46841862	10.1579548	0.24123
15	7.92117177	10.954396	0.074431
16	11.95834614	16.4556678	0.23902
17	3.62740587	4.1560713	0.30639
18	6.34823896	0.5981264	0.2985
19	9.45749638	−2.020902	0.38298
20	−1.95639684	6.8142173	0.30587
21	6.22859991	7.6175812	0.1379
22	4.84957444	14.1012447	0.13958
23	3.26280007	−2.1870801	0.26696
24	−3.61956196	3.6408247	0.33199
25	−5.20300959	6.5937315	0.39425
26	3.74976675	10.4215318	0.13719
27	8.4389697	14.5183123	0.072586
28	7.04080151	4.1116067	0.21658
29	8.87521784	17.9407346	0.22048
30	1.67284558	7.114281	0.3569

, the results that illustrate another case during a quiet time corresponding to an epoch on DOY 301 of the year 2013 are given. In this example, the spatial mean found was 0.29 TECu/min, and the standard deviation was 0.15 TECu/min, corresponding to a low occurrence of irregular structures and low ROTI values in all the region.

4. Conclusions

This study used an unsupervised machine learning technique to provide a spatial mean of the rate of change of TEC index (ROTI) to characterize the irregularities of the ionosphere over the West African region. Due to limitation of the availability of GNSS ground-based stations in the region, we performed ordinary Kriging interpolation to monitor and map the evolution of the ionospheric plasma irregularities. Two cases representing geomagnetically disturbed and quiet conditions during a period of high occurrence of ionospheric plasma irregularities in the equinoctial month of October 2013 were analyzed. A sampling technique based on the K-means algorithm was applied to generate an optimal design of clusters of the ROTI maps to provide a representative regional spatial mean of the ionospheric irregularities’ distributions. This method is recognized to perform better in identifying clusters with a spherical shape as in the case of the ROTI regional maps over the Western African region.

The conclusions derived from the cases under study are as follows: (1) the results show that the ordinary Kriging algorithm produced low values of variance of approximately 0.5 TECu/min² in about 99% of the study area, implying trust in the results obtained; (2) it has been shown that the instantaneous ROTI maps produced by the ordinary Kriging algorithm could be a valid representation of the spatio-temporal plasma electron density distribution in regions with a low density of receivers such as the study region; and (3) the partitioning of the region into compact clusters allowed us to obtain a spatial mean of the ROTI values that could be useful for monitoring the irregularities’ occurrence over a region. If applied to denser GNSS stations regions, where no interpolation is needed, the clusters could be determined directly from the ROTI values at the IPPs, allowing instantaneous values for monitoring and mitigation purposes in real time.