Efficient clustering of GNSS stations for processing using double differences

Efficient clustering of GNSS stations for processing using double differences | Research Square window.SnipcartSettings = { analytics: { enabled: false } }; (function() { var accessVector = localStorage.getItem('access_vector') || ''; window.dataLayer = window.dataLayer || []; if (accessVector) { window.dataLayer.push({ user: { profile: { profileInfo: { snid: accessVector } } } }); } })(); (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0],j=d.createElement(s),dl=l!='dataLayer'?'&l='+l:'';j.async=true;j.src='https://www.googletagmanager.com/gtm.js?id='+i+dl;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-K279D39R'); Browse Preprints In Review Journals COVID-19 Preprints AJE Video Bytes Research Tools Research Promotion AJE Professional Editing AJE Rubriq About Preprint Platform In Review Editorial Policies Our Team Advisory Board Help Center Sign In Submit a Preprint Cite Share Download PDF Research Article Efficient clustering of GNSS stations for processing using double differences Shane P. Grigsby, Demián D. Gómez This is a preprint; it has not been peer reviewed by a journal. https://doi.org/ 10.21203/rs.3.rs-7096364/v1 This work is licensed under a CC BY 4.0 License Status: Published Journal Publication published 21 Jan, 2026 Read the published version in GPS Solutions → Version 1 posted 8 You are reading this latest preprint version Abstract The rapid growth of GNSS networks poses significant challenges for efficiently processing large datasets using double-difference techniques. In this study, we introduce a novel clustering algorithm, qmeans , which is based on bisecting k-means, to partition GNSS networks into smaller, manageable subnetworks or clusters for double-difference processing. We explore the trade-offs between cluster size, computational cost, and solution quality using a comprehensive dataset of approximately 1,200 stations distributed across México, the United States, and Canada. Our results demonstrate that partitioning the network into clusters of 20 to 30 stations with 6 overlap stations between clusters can reduce processing time by ~20%, while larger clusters of 40-50 stations with 10 overlap stations slightly improve solution precision. We show that the number of shared stations between clusters impacts both the computational efficiency and the precision of the final solution, with higher counts leading to better precision but also increased processing time. The qmeans algorithm is integrated into the open-source Parallel.GAMIT software, offering a scalable, flexible solution that can be applied to large GNSS networks. Our work sets a foundation for selecting optimal subnetwork sizes based on specific needs of a GNSS processing project, enabling faster processing without significantly sacrificing solution quality. GNSS processing double-difference processing clustering GNSS network GNSS solution scatter segmentation Figures Figure 1 Figure 2 Figure 3 Figure 4 Figure 5 Figure 6 1. Introduction Global GNSS networks providing open data have grown geometrically since the early 1990s from a few hundred continuously operating stations to over 20,000 known stations today. This number is even larger when considering intermittently occupied sites for survey campaigns and private stations that do not share their data openly. The vast quantity of GNSS observations generated daily is essential for numerous scientific, engineering, and mapping applications, requiring automated processing to manage the computational burden efficiently (Gómez et al. 2024 ). GNSS double-difference processing techniques remain a standard approach for high-precision positioning and network adjustments. To our knowledge, network design for large scale differential processing is a topic that has not been thoroughly discussed in geodesy. Yet, one of the main challenges in processing large GNSS networks at Ohio State University (OSU) has been efficiently dividing the network into subnetworks for processing in GAMIT/GLOBK (Herring et al. 2018 ). Currently, our GNSS database has over 6,100 stations adding up to ~ 19M station-days, with a typical GNSS processing project at OSU containing about 2,000 simultaneous stations. Current computational power restricts double-difference processing to a maximum of approximately 80 simultaneous stations per session. Additionally, processing time increases as s 2 , where s is the number of sites in the processing session, creating a significant bottleneck in the computation of large network solutions. As GNSS networks continue to expand, this constraint poses real challenges for large-scale analyses, requiring partitioning strategies, hierarchical processing schemes (International Organization for Standardization 2020 ), or alternative methodologies to integrate observations from thousands of stations without excessive computational cost. The increasing size of GNSS networks requires scalable solutions that balance computational feasibility with the need for accurate and consistent geodetic products. Addressing these limitations is crucial for maintaining the reliability and improving GNSS-derived positioning, velocity fields, and global reference frame realizations, such as the International Terrestrial Reference Frame (ITRF, Altamimi et al. 2023 ), as well as regional realizations like the Geodetic Reference System for the Americas (SIRGAS, Alves Costa et al. 2022 ). In this paper, we present a new algorithm we term qmeans clustering, based on the sklearn (Pedregosa et al. 2011 ) bisecting-kmeans, to partition a network of GNSS stations into subnetworks of fewer stations for efficient processing. Our results show that subnetworks of 20 to 30 stations provide an optimal partitioning to minimize the processing time while slightly increasing the scatter level. In contrast, subnetworks with 40 to 50 stations increase the computation time with a slight improvement in scatter. These results establish a foundation for selecting subnetwork sizes based on the objectives of a GNSS processing project, enabling a balance between processing speed and solution precision. Throughout this work, we refer to subnetworks of GNSS stations as ‘clusters’, and stations that are shared across clusters as overlap stations , in contrast to the more typical tie station language frequently used in GNSS literature. Typically, tie stations refer to an exchange of stations between clusters, while our use of overlap refers to expansion of a given cluster to include stations from neighboring clusters. The distinction here is in the mutual reciprocity implied by tie stations –i.e., ‘tying’ clusters A and B together with tie stations would involve adding N stations from A into B, while also adding N stations from B back into A. In contrast, overlap stations may or may not have any reciprocal relationship between specific clusters: if cluster A ‘overlaps’ into cluster B by N stations, cluster B may then instead overlap into another cluster that isn’t A. We use overlap stations because we are able to define guaranteed algorithmic performance when discussing the processing behavior in the context of single clusters and their expansion, and further note that the number of overlap stations for a single cluster can be conceptually thought of as half of the number of analogous tie stations . The remaining sections of this manuscript are divided as follows: Methods describes the qmeans clustering algorithm and the methodology used to connect neighboring clusters. Results show a one-year ~ 1,200 station network GAMIT run using two cluster sizes and three levels of station overlap, and the impact of cluster size and overlap in the weighted root mean square (wrms) scatter of the solutions. Discussion develops and discusses a model to predict the execution time of a double-difference processing session of s stations– and allows us to understand the trade-off between cluster size, shared stations between neighboring clusters, and session execution time. 2. Methods Segmentation of points in space is a problem typically addressed using unsupervised learning, specifically unsupervised clustering. The problem is unsupervised because we do not have an a priori set of labels to train from, and need to accomplish the segmentation into clusters using only the information inherent in the location data of our stations at run-time. This also means that our cluster labeling and assignment is independent between runs, although the inherent spatial structure may lead to some persistent clusters that are forced to reform due to geometry constraints imposed by coasts, national borders, and other human or natural factors that influence station placement and density. For our particular use case in defining subnetworks of GNSS networks, we have an additional constraint that is not typically present in other clustering problems: we need to define our partitions with overlapping point membership between the subnetworks, so that the reference frame maintains consistency between adjacent subnetworks. Existing clustering techniques include k-means variants, spectral methods, agglomerative techniques, and density-based methods. None of these techniques address our problem completely. The k-means techniques split points into groupings of equal variance, minimizing the within-cluster sum of squares (termed ‘inertia’) across the dataset (Forgy 1965 ; Lloyd 1982 )—however, k-means requires that the number of output clusters be specificized in advance as the ‘k’ parameter. Spectral clustering similarly requires an a priori number of clusters to be provided in advance, and additionally struggles as the number of clusters to segment grows (Ng et al. 2001 ). Agglomerative techniques initiate each observation as a cluster and iteratively ‘link’ neighboring clusters together as merges; while these methods are hierarchical, they also tend to lead to unbalanced cluster membership, a behavior that is at odds with our downstream processing goals. Density based methods include popular algorithms such as DBSCAN (Ester et al. 1996 ; Schubert et al. 2017 ) and OPTICS (Ankerst et al. 1999 ), and can be conceptually thought of as convolving a density kernel over the input data points, and then passing a plane through the output heatmap topology of density where the isolated ‘peaks’ form distinct cluster objects—however, this leaves a subset of points (in the valleys) that are not assigned to any cluster and are simply labeled as noise. Other clustering algorithms such as Mean Shift and Affinity Propagation have poor computational complexity, and none of the above methods address the need for overlapping membership between clusters. As our base clustering model, we selected bisecting-kmeans clustering, a hierarchical partitioning variant of traditional kmeans (Steinbach 2000 ), which we modified and improved to adjust it to our needs. More specifically, we choose the variant from sklearn, which itself uses and modifies kmeans++ (Arthur and Vassilvitskii 2006 ). Our algorithm modifies the existing sklearn implementation with a new set of termination conditions which eliminates the k parameter requirement of specifying the number of output clusters at algorithm run-time. Since the ‘ k’ in ‘kmeans’ refers to the number of output clusters, we term our modified clustering algorithm qmeans to distinguish it from its parent algorithms. Our qmeans algorithm conceptually operates opposite to the agglomerative techniques: rather than starting with all points as independent clusters that are iteratively merged, we instead start with a single cluster of all stations that is iteratively bisected. We then supplement our new qmeans algorithm with two post processing steps: an overcluster routine to expand cluster membership to adjacent station clusters, and then a prune step to filter subnetwork clusters that are fully redundant after cluster expansion. Note that we adopt the convention of italicizing names of functions , classes , and parameters that we use, inherit from, or modify. 2.1 Qmeans cluster formation Since qmeans is a modification of bisecting-kmeans , which is itself a variant of kmeans + + and kmeans clustering, it is useful to review the parent algorithms to understand how qmeans works. As mentioned above, kmeans works by dividing N points (i.e., our stations) into K non-overlapping clusters C , and minimizing inertia which is defined as: Where µ j is the centroid of proposed cluster j , and x i is a point coordinate from the set of all coordinates proposed for assignment within cluster j . In the original algorithm, this is done by selecting k centroid centers, and then perturbing those centers (and adjacent point membership) to iteratively minimize inertia; the sklearn implementation uses kmeans + + which accelerates convergence by using random seeding, an enhancement which is of marginal impact when k is low. For both qmeans and bisecting-kmeans , there are two nested loops iterating: an inner-loop iterating to converge on the centroid and membership proposal that will optimally bisect a given cluster ‘node’ within the recursion node tree, and an outer-loop iterating over nodes to bisect. For k = N , the bisecting-kmeans algorithm sets k inner = 2 , and splits the dataset after optimizing cluster centroid location and membership to minimize inertia as measured in (1). This bisection process is repeated recursively in the outer-loop by selecting the cluster with either the largest inertia or largest point membership, and splitting that cluster (using k inner = 2 ), with the recursion terminating when the total number of clusters is equal to the overall number of output clusters requested by k . As a consequence, the bisecting-kmeans algorithm is hierarchical, reduces run-time complexity for high values of k , and produces clusters with more uniform membership. These latter two properties make bisecting-kmeans especially appealing for our GAMIT processing pipeline, as having lower variance in cluster membership number allows us to further optimize for run-time. The base kmeans algorithm must effectively recalculate the entire clustering structure from scratch whenever the k parameter is incremented. In contrast, the bisecting variant maintains the prior i- th outer-loop iteration cluster centroids and membership labels in a recursion tree while swapping the current highest inertia (or largest membership) cluster node within a given outer-loop iteration for two new subcluster nodes that bisect that parent. This outer-loop iterative structure for the bisecting algorithm allows us to remove the k parameter entirely, and instead replace it with our own boundary conditions. Our resulting qmeans algorithm replaces the k parameter with a single new qmax parameter for maximum cluster size. The qmax parameter is a hard boundary condition that provides an algorithmic guarantee that the output clustering will include only clusters of size at or below the qmax parameter setting– if a node within the qmeans recursion tree is larger than qmax it is bisected, if the node is smaller than qmax it is not, and the algorithm terminates when there are no remaining nodes left in the tree to bisect. Our qmeans algorithm produces clusters with different properties and streamlines the implementation code when compared with the parent kmeans + + and bisecting-kmeans algorithms. As mentioned above, bisecting-kmeans chooses the next node to bisect by selecting the node with either the largest inertia or the largest point membership, and will produce a different output clustering depending on which of the two metrics is selected. For qmeans , the distinction is meaningless and does not impact the output clustering at all. Deciding on which cluster to bisect matters in bisecting-kmeans because the algorithm increments the total number of clusters at each outer-loop iteration and will therefore terminate after k − 1 total outer-loop iterations–this means that even though the bisection of any particular node will produce the same output subnodes from applying (1), the choice of which order to apply the fixed k bisections will produce different results. In contrast, qmeans is not iteration bound for the outer-loop and will apply bisections, using (1) in the inner-loop to determine output nodes until there are no nodes above qmax remaining to bisect. One drawback of kmeans algorithms in general is that they perform best when clusters are convex; because the technique bisects points into groupings of equal variance, this family of algorithms handles irregular shapes poorly. Low membership bisections can occur on our data inputs because of variable density and sparse station coverage–for example, across oceans with widely separated island stations in a geometry that yields a small member cluster for the tightly grouped GNSS stations over islands such as Hawaiʻi, and single member clusters at distant atoll stations such as Guam. While both low cluster membership and high cluster membership present problems for our GAMIT runs, cluster expansion for overlap stations ensures that even single station clusters such as Guam will have multiple stations to form a subnetwork for a GAMIT session. Since clusters with high station membership can cause non-recoverable crashes in our processing runs due to memory errors, by design we provide algorithmic guarantees for a ceiling on the station count membership during both cluster formation from qmeans and expansion via overcluster . A final additional post processing by the prune function (Section 2.3 ) is implemented to remove any redundant clusters that may be formed entirely by stations that overlap with other clusters. 2.2 Cluster expansion When processing subnetworks within GAMIT, we require there to be overlapping stations between the subnetworks, and we use the overcluster function to ensure this overlap. The overcluster function takes three required parameter arguments: overlap , nmax , and rejection threshold . The rejection threshold is a distance set by default to 5,000 km, and rejects any points that fall outside of that distance—this is mainly present as a sanity check for very isolated stations in polar regions or oceans, and ensures that GAMIT can form double differences on adjacent stations. While the rejection threshold was enabled to be user modified, we have had no reason to do so and it remains fixed for all of our processing runs. The overlap parameter specifies how many additional stations we want to expand our cluster by, while nmax sets a maximum number of stations to include from any given adjacent cluster. The cluster expansion is recursive, and uses a balltree spatial index from sklearn to calculate nearest neighbors (Omohundro 1989 ). More explicitly, we take all of the GNSS stations except for the cluster being expanded and build a spatial index, using the index to calculate the k = 1 nearest neighbors from each station outside of the cluster, to each member of the cluster. This intersection of cluster and non-cluster members results in a sorted list of station distances to our cluster’s member stations, and we add the external station with the shortest distance as an overlap point. We then update the cluster and non-cluster membership to include/exclude this new station and rerun the nearest neighbor analysis. As we add points from neighbors to our cluster, we also track the original label of the added stations; when the sum of a given label’s occurrence is equal to the nmax parameter, we remove that entire neighboring cluster from future consideration. This ensures that we get broad connectivity and representation from neighboring clusters. The overcluster expansion algorithm terminates when the number of stations added to the cluster under expansion reaches that value set by the overlap parameter—unless there are no stations to add under the distance threshold parameter and the algorithm terminates early, a condition which was not triggered for any of the results we present given our other parameter choices. The overlap parameter provides us with a simple algorithmic guarantee for our processing pipeline. The maximum size for any cluster provided to GAMIT will be: Max cluster size < = qmax + overlap (2) The nmax parameter ensures that the cluster overlap doesn’t heavily favor a single neighboring cluster. If nmax = 1 , then the number of clusters that are expanded over will be equal to what the overlap parameter is set to. Using the default nmax = 2, the number of clusters that are expanded over will be bounded between overlap and overlap / 2. It is important to note that while we are guaranteed these minimum overlaps, in practice many central clusters will have more redundancy than this. Asymmetry in reciprocity means that some central clusters will have substantial additional redundancy when the subnetworks are reconciled. While redundancy improves the scatter of the GNSS solutions, some clusters are either so small or so central that their membership is entirely redundant, with the stations fully represented across other clusters. Thus, we implement a final postprocessing step to prune redundant clusters and improve run-time. 2.3 Cluster pruning The output from the overcluster function is a Boolean M by N matrix, where the M rows correspond to M clusters, and the N columns correspond to N total stations. Each row has columns marked with True values to indicate which stations belong to that cluster after the cluster has been expanded. Depending on the GNSS station distribution geometry and overcluster parameters, there will be varying redundancy in the membership of the clusters. We require that each station is represented at least once; while we expect that some stations will be represented multiple times due to overlaps, we proactively remove cases where the full cluster itself is redundant because that entire subnetwork membership is fully present in one or more other subnetwork clusters. To check for this redundancy, we take the sum of each column to verify if that GNSS station is captured one or more times, and then iterate through the matrix by temporarily removing a given row. As we remove a row / cluster, we recalculate the column totals; if removing that row causes us to lose representation of any station in our input dataset, then we pass on to the next row and cluster. However, if removing a cluster leaves all GNSS stations present in at least one other subnetwork, then we delete the cluster and update the M by N matrix encoding before continuing to recursively check the remaining rows and clusters for if they have any fully redundant coverage. Using our production parameters, a linear scan through the cluster encodings trims approximately 20–30 percent of overall GNSS stations that are ultimately fed into GAMIT, with a concurrent drop in overall processing time. We expect the impact to processing time is disproportionately impacted by the extreme tail of smaller sized clusters as explained in Section 4 , therefore, while we maintain the option to execute the prune function as a simple linear scan, our current default method will monotonically sort the rows for the Boolean overcluster matrix according to the row sum, and then preferentially eliminate redundant rows by removing the smallest redundant clusters first. The successive steps of qmeans , overcluster , and prune are highlighted for a part of a global processing run in Fig. 1 . 2.4 Computational considerations and software availability The bulk of the processing time is in the GAMIT runs, which is why we seek an efficient cluster structure—so we can more efficiently process daily GNSS solutions on the time scale of years. The algorithms presented here are computationally efficient; the processing of a full day of ~ 1,000 GNSS stations may take ~ 20 hours in GAMIT, while our clustering processing and postprocessing adds only a few seconds per day, depending on station count. The bisecting qmeans algorithm is stochastic, initiating random centroid locations as it searches for an optimum solution. However, to ensure that our results are reproducible, we typically fix the seed in the random number generator so that we have deterministic behavior when comparing runs with different parameter values. The qmeans algorithm, as well as the overcluster and prune functions, are open-source functions integrated into Parallel.GAMIT, a BSD licensed open-source library maintained by the authors to parallelize GAMIT runs. Our library is available through GitHub ( https://github.com/demiangomez/Parallel.GAMIT ) , and is additionally installable from the pypi community repository under the package name pgamit when using the python ‘pip install’ command. While GAMIT ( http://www-gpsg.mit.edu/gg/ ) is available from its authors on request, distribution is limited to non-commercial research use, which also means that it is distributed separately from the pgamit package. The clustering functions we develop are freely available, and can be run both within the Parallel.GAMIT parallelization routines, with clustering functions (within the pgamit.cluster module) following the sklearn documentation and API style so that they can be used separately as needed. 3. Results To test qmeans and evaluate the time-precision trade-off, we processed one year of data (2008) for a GNSS network in North America—composed of ~ 1,200 unique stations in México, the United States, and Canada—and compared the run-times and wrms of four different cluster configurations: qmax = 20 with 4, 6, and 10 overlap stations respectively, and qmax = 40 with 10 overlap stations which is our baseline case for comparisons. Since qmeans cannot guarantee clusters with a uniform number of stations, Fig. 2 a shows the stations per cluster distribution for the test year using qmax = 40 with 10 overlap stations, which results in a membership floor for clusters (at least 11 stations), as well as a ceiling (no more than 50 stations). Figure 2 b shows that for qmax = 20 with 10 overlap stations, the total cluster count per day rises as expected, and the clusters are distributed across a tighter range between 11 and 30 stations per cluster. For a graphical representation of the clusters using qmax = 40 , see Fig S1 in the Supporting Information, and Fig S2 for cluster size and overlap station count histograms for a single day. 3.1 Clustering performance of the qmeans algorithm In terms of the number of clusters produced for a Parallel.GAMIT session, Fig. 3 a shows that with 10 overlap stations, qmax = 20 results in nearly double the number of clusters compared to the qmax = 40 configuration. It is important to note that the number of clusters produced by qmeans is invariant with respect to the overlap parameter specified, as the partitioning in the qmeans process is independent of the overcluster function, which does not add or remove clusters. Although the prune function is impacted by structure and number of overlap stations when determining redundant clusters to remove, the output number of clusters is primarily determined by the qmax parameter of qmeans. Figure 3 b shows the observed cumulative time difference between our reference case of qmax = 40 with 10 overlap stations, and qmax = 20 with 10, 6, and 4 overlap stations respectively. For the case where only qmax differed and overlap stations were held at 10, our results showed that qmax = 40 took ~ 250 more hours to run relative to qmax = 20 . Based on the total processing time for qmax = 40 (~ 8,500 hours), 250 hours represents only a 3% reduction of run-time. A more significant run-time reduction was observed for qmax = 20 when varying overlap stations to 6 and 4, where the total run-time difference from our qmax = 40 reference reached 2,000 and ~ 2,600 hours (23% and 30% reduction), respectively. Table 1 Average total stations to process after overcluster and prune . Overlap = 4 Overlap = 6 Overlap = 10 qmax = 40 1,320 1,408 1,608 qmax = 20 1,480 1,617 2,010 Table 1 provides a summary of the mean number of stations per cluster configuration for an average network size of ~ 1,200 stations, and highlights the fundamental tradeoff between cluster size, total number of stations to process, and run-time impact. Using 10 overlap stations, the total number of input stations for GAMIT to process when qmax = 40 is 2,010 versus 1,608 for qmax = 20 ; however, the total cumulative run-time is 3% less for qmax = 20 , despite processing 402 additional stations, which is 25% more stations than the reference qmax = 40 . This is expected given the polynomial behavior of run-time as a function of cluster size (discussed below in Section 3.3 ), but optimization is complicated because of how overall station count varies as a function of both number of clusters (itself a function of qmax ) and number of overlap stations. As Table 1 shows, since smaller clusters are more numerous, increasing the number of overlap stations has an outsized effect on overall number of cumulative stations to process– going from 10 overlap stations to 6 overlap stations reduces cumulative station count by only 200 for qmax = 40 , but by over 400 for qmax = 20 since there are approximately double the clusters produced to add overlap stations to. Given that qmax = 40 with 10 overlap stations and qmax = 20 with 6 overlap stations have approximately the same cumulative number of stations to process (1,608 versus 1,617), the red line from Fig. 3 b presents an interesting comparison on run-time, a difference driven entirely by the size distribution of the clusters. 3.2 Solution quality impact Cluster size and number of overlap stations is expected to imply a trade-off both between execution time, and solution quality. To assess this trade-off in solution quality, we analyzed the daily GNSS repeatability using qmax = 40 with 10 overlap stations as a reference solution and compared it against solutions using qmax = 20 with station overlaps of 4, 6, and 10 for the same test year and the same network. Daily repeatabilities are obtained by transforming consecutive daily GNSS solutions onto one another using six-parameter Helmert transformations (Bevis and Brown 2014 ). After all solutions have been aligned, the differences between individual station coordinates are computed and a wrms value in north, east, and up provides an estimate of the scatter level of the solution for each individual site. Figure 4 a-c show the cumulative empirical distribution function for the daily repeatability wrms scatter for the north, east, and up components for all stations in the analysis. A clear pattern emerges from these plots, which is most pronounced in the east component, where we note that the base reference case ( qmax = 40 with 10 overlap stations) outperforms the other solutions using qmax = 20 across all overlap parameter values, with overlap = 4 producing the worst results, while solutions for overlap = 6 and overlap = 10 closely follow each other. Interestingly, the vertical component in Fig. 4 c shows that solutions with qmax = 40 or 20 and overlap = 10 closely follow each other. Although the smallest subnetwork clusters with the least overlap are consistently outperformed by other configurations, as shown by the 0.1 and 0.3 mm shaded regions of scatter for the horizontal and vertical components in Fig. 4 a-c, all tested cluster configurations fall within these tight regions and therefore can be considered statistically equal in terms of expected GNSS solution noise. Despite these results, the level of agreement between the wrms of these experiments cannot be always guaranteed, given that other effects such as network geometry or latitude location can affect the overall solution quality. Hence, whenever minimum scatter is desired, cluster size and number of overlap stations should be increased to achieve needed precision. 3.3 Processing time as function of station count from historical data As part of our parallel processing system for GAMIT, we save all of the relevant statistics for each run to a Postgres database, including execution time, overlap station count, and solution quality indicators. Since the deployment of Parallel.GAMIT at OSU in 2019, we have performed a total of ~ 2.5M individual GAMIT runs for various projects. Using this Postgres database, we selected a subset of GAMIT runs for execution time statistics over an average GNSS project of ~ 1,000 stations. We did not use the entire 2.5M runs in order to avoid double counting the same network distributions and duplicate runs, since many of the sessions are densifications or reruns using different models and orbits. After removing any ‘outlier’ runs (those that took more than 80 minutes to complete) that could alter our statistics, the dataset used in this work is composed of ~ 360k GAMIT executions. The removed long execution times account for less than 0.5% of the dataset, and were mostly related to external events not associated with the GNSS processing, such as computer network outages and other technological difficulties encountered during the processing. Figure 5 a shows the histogram of execution times for all the selected GNSS sessions from 1994 to 2022. For simplicity and consistency of the results, we performed all statistics and analysis using GPS-only solutions (rather than multi-constellation). The station clusters for these runs were created using an ad-hoc, non-optimized script (deployed from 2019 to 2024, prior to developing qmeans ) which creates station clusters by proximity and splits the resulting subnetworks if the number of stations in the processing session is above a threshold of ~ 60. Once each cluster is determined, the clusters contain sets of unique stations to process and overlap stations from neighboring clusters, which are used to merge the solutions with GLOBK upon finishing all the processing jobs. Additionally, the clusters are also tied together using a ‘backbone’ network of evenly distributed sites across the entire area of interest to provide a ‘rigid support’ for the chainmail-like network of smaller clusters (see Appendix A for a description of the backbone network algorithm) and guarantee a robust combination. Figure 5 b shows the cluster size histogram where we note two peaks, one around 15 and a second one near 40 stations. This distribution is produced due to the impossibility of splitting the stations into clusters with an exact number of stations. The data for execution time as a function of station count (Fig. 5 c) can be fitted using a quadratic expression to obtain an estimate of the execution time; to account for other factors affecting the execution time at lower station counts, such as overhead time to read observation and metadata files, models, and others, we fit a second-degree polynomial rather than a pure quadratic. We performed this fit and found that the execution time t in minutes can be expressed as: $$\:t=6.04\:-0.10\cdot\:s\:+0.021\cdot\:{s}^{2}$$ 3 where s is the number of stations in the processing session. Given the large number of data points used in the fit, the weighted root mean square error of our polynomial fit (95% confidence interval) was ~ 2 minutes. 4. Discussion Although determining the exact reason behind the solution quality improvement with increasing cluster size and overlap station count is outside the scope of this work, we surmise that clusters with more stations have, on average, a larger spatial aperture, and also larger number of double differences. Larger aperture and more double differences will improve the determination of the zenith tropospheric delay parameters which, in turn, improves both horizontal and vertical component scatter. Additionally, increasing the overlap station count increases the redundancy of the network, which helps mitigate the solution scatter when combining all clusters together. 4.1 Execution time model Due to memory limits and also because execution time increases quadratically, GNSS sessions with N > 60–80 stations need first to be partitioned into C clusters of s stations each. As seen in Section 3 , our GAMIT runs occur on datasets where the cluster sizes are not uniform; the compute run-time for these production sessions can be approximated by modifying Eq. 3 : $$\:T={\sum\:}_{i=1}^{C}\left(6.04\:-0.10\cdot\:({s}_{i}+e)\:+0.021\cdot\:{({s}_{i}\:+e)}^{2}\right)$$ 4 Where e is the overlap parameter value (constant for all clusters), C is the actual number of clusters formed, and s i is the cluster size of cluster at index i . Eq. ( 4 ) does not account for geometry considerations within GAMIT or other background processes that impact run-time; still, we find that the time estimate T has only slight bias in predicting overall run-time when compared to metadata within our postgres database which records the empirical run-time. Equation 4 is of little practical value by itself given that we already log the actual run-time of the GAMIT processing runs, and don’t need predictive estimation of individual subnetwork run-time for our operational processing. More useful is generalizing Eq. 4 to visualize the behavior of the subnetworks and help illustrate parameter choices that can guide us in optimizing our compute pipeline. The ratio between parameters e and s , \(\:R=\frac{e}{s}\) , can be defined as the ‘redundancy’ of a given Parallel.GAMIT session, where \(\:R=0\) indicates that there are no shared stations between clusters (this is an undesired condition) and \(\:R=1\) indicates that all stations are processed separately in two or more clusters. To idealize the compute run-time of a given GNSS network, the total number of clusters to be processed is \(\:C=\frac{N}{s}\) , assuming N total session stations are divisible by s , where stations per cluster s is now fixed at a single value for all clusters rather than varying per cluster. Assuming this idealized uniform distribution of clusters, the total execution time T , including the redundancy, can be estimated as: $$\:T={\sum\:}_{i=1}^{C}{t}_{i}=C\cdot\:\left(6.04\:-0.10\cdot\:s\cdot\:(1+R)\:+0.021\cdot\:{\left[s\cdot\:(1+R)\right]}^{2}\right)$$ 5 Equation ( 5 ) uses R instead of e , as redundancy is defined as the ratio between e and s , which simplifies the interpretation of contour plot in Fig. 6 . Figure 6 shows a contour plot of T (for N = 1,000) for Eq. ( 5 ) as a function of s and R , also including the contours of equal overlap stations, which vary with s for a constant R . For clarity, we plotted a continuous contour field regardless of the value of \(\:C=\frac{N}{s}\) , although the idealized execution time model is only valid when C is an integer. We also show, as a black dashed contour, the minimum run-time for each value of R and its corresponding cluster size. It should be noted that this line marks the limit of the cluster mediated computation efficiency: cluster sizes to the left of this dashed line will increase the run-time, decreasing efficiency rather than improving it. In other words, the execution time needed to complete a multi-GAMIT run session decreases concurrently with cluster size (i.e., as the number of clusters formed per session goes up) until the limit denoted by the intersection of the black dashed contour. Thus, if for example one uses 4 overlap stations, a session divided into uniform clusters of, say, 20 stations, takes ~ 400 minutes less to finish than if it were divided into uniform clusters of 50 stations, although the number of clusters is larger for the former case. This tendency of decreasing time with lower cluster size, however, is less significant with increasing number of overlap stations. Figure 6 shows that for increasing overlap stations, the overlap contours become increasingly parallel to the execution time contours, meaning that moving along overlap contours (and changing the cluster size) does not significantly change the execution time. We find that Fig. 6 is mostly useful to examine the execution time behavior purely as a function of the overlap parameter, and help select an ideal cluster size to target. As we mentioned before, run-time is not the only consideration that a user has to account for when choosing a given partitioning. From Fig. 6 it is clear that the choice of cluster configuration can have a considerable impact in the time needed to compute the solutions. This time difference will translate into a more or less significant ‘wall time’, i.e., the actual time needed to compute the solutions, depending on the number of compute nodes used to process the data. For instance, to process the totality of the GNSS stations in México, the United States, and Canada from 1994 to 2025 using clusters of size 20 with 6 overlap station clusters requires a total computation time of 8,223 days, or ~ 33 days at 250 cores. If the same project is partitioned into clusters of size 50 with the same 6 station overlap, then the required run-time is 9,600 days, or ~ 38 days at 250 cores, a wall time difference of 5 days. The wall time difference becomes more significant with a smaller compute cluster, and in the case of using only 100 cores, the wall-time difference between these two example runs is close to 14 days. 4.2 Conclusions This study presents a novel approach for efficiently partitioning large GNSS networks into subnetworks or clusters using a modified bisecting k-means algorithm, which we called qmeans . The primary goal of this work was to develop a method for reducing the computational burden of processing large GNSS datasets while maintaining the precision of the results. Through extensive testing, we demonstrated that qmeans provides a robust and scalable solution for dividing GNSS stations into clusters of manageable sizes, facilitating faster double-difference processing in GAMIT/GLOBK. The qmeans algorithm itself offers several advantages. First, the total clustering processing time was found to add only a few seconds per day, a minimal overhead relative to the time savings achieved by efficiently dividing the network. Second, by eliminating the need to predefine the number of clusters (as in traditional k-means clustering), qmeans provides a flexible and dynamic solution that can adapt to varying network sizes and spatial configurations. This feature is particularly beneficial when processing GNSS networks where station distributions are irregular, such as in sparsely populated regions or oceanic islands. Moreover, the hierarchical nature of qmeans , combined with the post-processing steps ( overcluster and prune ), ensures that the final clusters are both computationally efficient and spatially coherent, with sufficient overlap for robust reference frame realization. As expected, the increase in the station count caused by adding additional overlap stations leads to a decrease in computational efficiency, as some stations are processed more than once. However, we have also demonstrated that this reduction in efficiency results in a corresponding increase in precision, which is a desired outcome when processing GNSS data. Our results show that partitioning the network into clusters of approximately 20 to 40 stations provides an optimal balance between reducing computational time and preserving solution accuracy. Smaller clusters (20 stations) lead to faster processing times, with slightly increased solution scatter, while larger clusters (40 stations) offer a modest reduction in scatter at the cost of increased processing time. In our experiments, the qmax = 20 with overlap = 6 configuration was shown to cut processing time by nearly 23% compared to the qmax = 40 with overlap = 10 configuration despite processing a similar number of total stations within the session. Similarly, using qmax = 20 with overlap = 4 reduced the processing time by 30%, despite generating twice the clusters– both cases with a trade-off of a slight increase in solution scatter. These findings suggest that for large-scale GNSS projects, careful consideration of cluster size and overlap station count is critical for achieving efficient processing without compromising the quality of the result. By contrast, if one is seeking to obtain a fast solution for testing purposes, reducing the number of stations and overlap stations per cluster significantly reduces the processing time for GNSS solutions. Declarations The research leading to these results received funding from the National Geodetic Survey (NGS), under Grant Agreement AWD-115866. Conflicts of interest/Competing interests All authors certify that they have no affiliations with or involvement in any organization or entity with any financial interest or non-financial interest in the subject matter or materials discussed in this manuscript. The authors have no financial or proprietary interests in any material discussed in this article. Author Contribution S.G. and D.G. contributed equally to the manuscript. S.G. designed and developed the qmeans, overcluster, and prune functions described in the manuscript, and led integration for adding them to the ParallelGAMIT library. D.G. envisioned experimental design, interpreted the results, and developed the run time predictive model. S.G. wrote the initial draft of the Methods and Discussion sections; D.G. wrote the initial draft of the Introduction and Results sections; both S.G. and D.G. reviewed, edited, and revised the entire manuscript. D.G. added supplementary material, and produced figures 2, 3, 4, 5 and 6. S.G. produced figure 1, and authored and tested the pgamit.cluster module used to create all figures except figure 5. S.G. and D.G. maintain the open source project Parallel.GAMIT (pgamit), and D.G. is the project founder of that project. Acknowledgement We acknowledge Eric Kendrick for his role in maintaining the postgres database of GAMIT run statistics, which was queried by the authors to develop the predictive runtime model. Data Availability Data availability: Data is public and available through various websites, including the International GNSS Service data repository, through the EarthScope Facility Archive, and the NOAA Continuously Operating Reference Station (CORS) Network (NCN), managed by NOAA/National Geodetic Survey. References Altamimi Z, Rebischung P, Collilieux X, et al (2023) ITRF2020: an augmented reference frame refining the modeling of nonlinear station motions. J Geod 97:47. https://doi.org/10.1007/s00190-023-01738-w Alves Costa SM, Sánchez L, Piñón D, et al (2022) Status of the SIRGAS reference frame: recent developments and new challenges. In: IAG International Symposium on Reference Frames for Applications in Geosciences. Springer Nature Switzerland Cham, pp 153–165 Ankerst M, Breunig MM, Kriegel H-P, Sander J (1999) OPTICS: ordering points to identify the clustering structure. ACM SIGMOD Rec 28:49–60. https://doi.org/10.1145/304181.304187 Arthur D, Vassilvitskii S (2006) k-means++: The advantages of careful seeding. Stanford Bevis M, Brown A (2014) Trajectory models and reference frames for crustal motion geodesy. J Geod 88:283–311. https://doi.org/10.1007/s00190-013-0685-5 Costa SMA, Sánchez L, Piñon D, et al (2023) Status of the SIRGAS reference frame: Recent developments and new challenges. In: International Association of Geodesy, Reference Frames Symposium REFAG2022 Ester M, Kriegel H-P, Sander J, Xu X (1996) A density-based algorithm for discovering clusters in large spatial databases with noise. In: kdd. pp 226–231 Forgy EW (1965) Cluster analysis of multivariate data: efficiency versus interpretability of classifications. biometrics 21:768–769 Gómez DD, Bevis MG, Caccamise DJ, et al (2024) An empirical tool for predicting the presence or absence of coseismic displacements at GNSS stations. GPS Solut 28:214. https://doi.org/10.1007/s10291-024-01758-9 Herring TA, King RW, Floyd MA, McClusky SC (2018) Introduction to GAMIT/GLOBK International Organization for Standardization (2020) Geographic information - Geodetic references - Part 1: International terrestrial reference system (ITRS) (ISO Standard No. 19161-1:2020(E)) Lloyd S (1982) Least squares quantization in PCM. IEEE Trans Inf Theory 28:129–137 Ng A, Jordan M, Weiss Y (2001) On spectral clustering: Analysis and an algorithm. Adv Neural Inf Process Syst 14: Omohundro SM (1989) Five balltree construction algorithms Pedregosa F, Varoquaux G, Gramfort A, et al (2011) Scikit-learn: Machine learning in Python. J Mach Learn Res 12:2825–2830 Schubert E, Sander J, Ester M, et al (2017) DBSCAN Revisited, Revisited: Why and How You Should (Still) Use DBSCAN. ACM Trans Database Syst 42:1–21. https://doi.org/10.1145/3068335 Steinbach M (2000) A Comparison of Document Clustering Techniques Michael Steinbach, George Karypis, and Vipin Kumar Additional Declarations No competing interests reported. Supplementary Files SupportingInformation.docx Cite Share Download PDF Status: Published Journal Publication published 21 Jan, 2026 Read the published version in GPS Solutions → Version 1 posted Editorial decision: Revision requested 10 Nov, 2025 Reviews received at journal 03 Nov, 2025 Reviewers agreed at journal 15 Sep, 2025 Reviewers agreed at journal 10 Sep, 2025 Reviewers invited by journal 13 Jul, 2025 Editor assigned by journal 13 Jul, 2025 Submission checks completed at journal 11 Jul, 2025 First submitted to journal 10 Jul, 2025 You are reading this latest preprint version Research Square lets you share your work early, gain feedback from the community, and start making changes to your manuscript prior to peer review in a journal. As a division of Research Square Company, we’re committed to making research communication faster, fairer, and more useful. We do this by developing innovative software and high quality services for the global research community. Our growing team is made up of researchers and industry professionals working together to solve the most critical problems facing scientific publishing. Also discoverable on Platform About Our Team In Review Editorial Policies Advisory Board Help Center Resources Author Services Accessibility API Access RSS feed Manage Cookie Preferences © Research Square 2026 | ISSN 2693-5015 (online) Privacy Policy Terms of Service Do Not Sell My Personal Information {"props":{"pageProps":{"initialData":{"identity":"rs-7096364","acceptedTermsAndConditions":true,"allowDirectSubmit":false,"archivedVersions":[],"articleType":"Research Article","associatedPublications":[],"authors":[{"id":485654116,"identity":"e4a69b77-3310-4f0c-91c1-0f49b6646f4f","order_by":0,"name":"Shane P. Grigsby","email":"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAZAAAAAyAQMAAABI0h/eAAAABlBMVEX///8AAABVwtN+AAAACXBIWXMAAA7EAAAOxAGVKw4bAAABG0lEQVRIiWNgGAWjYDACCSjNBkYVDAx8IB4PiDhAlJYzYJJILWBdjG1EaJGf3fzswQeG2nw+/uXPHvPOs5NjY+99+OFNBYMc340ErFoM7hwzN5zBcNyyTeKNuTHvtmRjNp7jxpJzzjAYS+LSIpFgJs3DcMyATeIMmzTvtgOJbRJpDNK8bQyJG3BokZ+R/k36D1jL8WfSvHPAWph/A7XU49LCcCPHTJqBocaAjb/BTJq3AayFDWRLggEuh93IKZPsMTgAtIXH3HDOMZBfjrFZzjkjYTjzzANcDtsm8aOizkC+//izB29q7OT42duYb7ypsJHnO47DYRC7DgMjKIGBiQchJIFbNQTUMTDwH2Bg/EFI3SgYBaNgFIxIAADgaFge4xle0AAAAABJRU5ErkJggg==","orcid":"","institution":"University of Colorado","correspondingAuthor":true,"prefix":"","firstName":"Shane","middleName":"P.","lastName":"Grigsby","suffix":""},{"id":485654117,"identity":"74fb95f4-045d-4739-84ed-8c402f487d8b","order_by":1,"name":"Demián D. Gómez","email":"","orcid":"","institution":"Ohio State University","correspondingAuthor":false,"prefix":"","firstName":"Demián","middleName":"D.","lastName":"Gómez","suffix":""}],"badges":[],"createdAt":"2025-07-10 22:53:14","currentVersionCode":1,"declarations":"","doi":"10.21203/rs.3.rs-7096364/v1","doiUrl":"https://doi.org/10.21203/rs.3.rs-7096364/v1","draftVersion":[],"editorialEvents":[{"content":"https://doi.org/10.1007/s10291-025-02020-6","type":"published","date":"2026-01-21T15:58:30+00:00"}],"editorialNote":"","failedWorkflow":false,"files":[{"id":87012133,"identity":"204baaa7-2f27-48b6-8e2f-48ed195a11ac","added_by":"auto","created_at":"2025-07-18 09:29:31","extension":"png","order_by":1,"title":"Figure 1","display":"","copyAsset":false,"role":"figure","size":632882,"visible":true,"origin":"","legend":"\u003cp\u003eSuccessive steps of \u003cem\u003eqmeans\u003c/em\u003e, \u003cem\u003eovercluster\u003c/em\u003e, and \u003cem\u003eprune\u003c/em\u003e(a) Initial \u003cem\u003eqmeans\u003c/em\u003e clustering using \u003cem\u003eqmax=16 \u003c/em\u003efor a subset of 1,008 stations from January 5th, 2022; the degenerate \u003cem\u003ered\u003c/em\u003e cluster of two stations was split off from the \u003cem\u003eyellow\u003c/em\u003e cluster, which had 17 stations prior to the split. (b) \u003cem\u003eOvercluster\u003c/em\u003e using \u003cem\u003eoverlap=4\u003c/em\u003e and \u003cem\u003enmax=2\u003c/em\u003e, shown as the hull of original \u0026amp; added stations (cyan cluster hull excluded for readability). The \u003cem\u003eblue\u003c/em\u003ecluster adds 2 \u003cem\u003epurple\u003c/em\u003e and 2 \u003cem\u003ered\u003c/em\u003e stations; \u003cem\u003ered\u003c/em\u003e cluster adds 2 \u003cem\u003epurple\u003c/em\u003eand 2 \u003cem\u003eyellow\u003c/em\u003e; \u003cem\u003epurple \u003c/em\u003eadds 2 \u003cem\u003eyellow\u003c/em\u003e and 2 \u003cem\u003ecyan\u003c/em\u003e; \u003cem\u003eyellow \u003c/em\u003eadds 2 \u003cem\u003epurple\u003c/em\u003e and 2 \u003cem\u003egreen\u003c/em\u003e; \u003cem\u003egreen \u003c/em\u003eadds 2 \u003cem\u003eyellow\u003c/em\u003e, 1 \u003cem\u003epurple, \u003c/em\u003eand 1 \u003cem\u003ecyan\u003c/em\u003e. (c) Final post \u003cem\u003eprune\u003c/em\u003e output, with degenerate \u003cem\u003ered\u003c/em\u003e cluster removed since it had complete overlap with the \u003cem\u003eblue\u003c/em\u003ecluster. Note the symmetry of overlap between the \u003cem\u003eyellow\u003c/em\u003e–\u003cem\u003egreen\u003c/em\u003e and \u003cem\u003epurple\u003c/em\u003e–\u003cem\u003eyellow\u003c/em\u003e clusters, and the asymmetric overlap for the \u003cem\u003eblue\u003c/em\u003e–\u003cem\u003epurple\u003c/em\u003eand \u003cem\u003egreen\u003c/em\u003e–\u003cem\u003epurple\u003c/em\u003e clusters.\u003c/p\u003e","description":"","filename":"floatimage2.png","url":"https://assets-eu.researchsquare.com/files/rs-7096364/v1/096c6d1484f0e2a5885249a5.png"},{"id":87012136,"identity":"07a7b6d8-e324-4696-8f29-e96d050abe30","added_by":"auto","created_at":"2025-07-18 09:29:31","extension":"png","order_by":2,"title":"Figure 2","display":"","copyAsset":false,"role":"figure","size":192765,"visible":true,"origin":"","legend":"\u003cp\u003eCluster size distribution for test year 2008 for two configurations. (a) \u003cem\u003eqmax=40 \u003c/em\u003ewith 10 overlap stations. (b) \u003cem\u003eqmax=20\u003c/em\u003ewith 10 overlap stations. The line of clusters with 45 stations corresponds to the backbone networks which are always calculated for each day regardless of \u003cem\u003eqmax\u003c/em\u003e parameter settings.\u003c/p\u003e","description":"","filename":"floatimage3.png","url":"https://assets-eu.researchsquare.com/files/rs-7096364/v1/d59fbc3dbc940d98898eafb4.png"},{"id":87012137,"identity":"ae61214d-ac4e-4f4e-aa6a-881d1e444c58","added_by":"auto","created_at":"2025-07-18 09:29:31","extension":"png","order_by":3,"title":"Figure 3","display":"","copyAsset":false,"role":"figure","size":102384,"visible":true,"origin":"","legend":"\u003cp\u003e(a) subnet count for \u003cem\u003eqmax=20 \u003c/em\u003eand \u003cem\u003eqmax=40\u003c/em\u003e configurations, with 10 overlap stations. (b) Cumulative run-time difference between reference \u003cem\u003eqmax=40\u003c/em\u003ewith 10 overlap stations, and \u003cem\u003eqmax=20\u003c/em\u003ewith 10 (blue line), 6 (red line), and 4 (yellow line) overlap stations respectively.\u003c/p\u003e","description":"","filename":"floatimage4.png","url":"https://assets-eu.researchsquare.com/files/rs-7096364/v1/f93422e35454392db67ef190.png"},{"id":87012447,"identity":"381041cd-5ac5-48c7-b2ad-c6115df1f8c1","added_by":"auto","created_at":"2025-07-18 09:37:31","extension":"png","order_by":4,"title":"Figure 4","display":"","copyAsset":false,"role":"figure","size":155374,"visible":true,"origin":"","legend":"\u003cp\u003eEmpirical cumulative distribution function for repeatability analysis of four cluster configurations (a) Weighted root mean square scatter for the north component. Red shaded area represents the ± 0.1 mm region. (b) same as (a) for the east component. (c) same as (b) for the up component. Red shaded area represents the ± 0.3 mm region.\u003c/p\u003e","description":"","filename":"floatimage5.png","url":"https://assets-eu.researchsquare.com/files/rs-7096364/v1/300252fe6186da32dbc63c62.png"},{"id":87012147,"identity":"647e94cb-34b5-4fb2-b8b6-96c7f756cd55","added_by":"auto","created_at":"2025-07-18 09:29:31","extension":"png","order_by":5,"title":"Figure 5","display":"","copyAsset":false,"role":"figure","size":227192,"visible":true,"origin":"","legend":"\u003cp\u003eEmpirical GAMIT execution time (a) Histogram of execution times used to obtain the polynomial quadratic parameters. (b) Same as (a) but for the cluster sizes. (c) Density of data points and fitted polynomial quadratic curve (3), execution time as a function of station count.\u003c/p\u003e","description":"","filename":"floatimage6.png","url":"https://assets-eu.researchsquare.com/files/rs-7096364/v1/68b61b0ce1b4245d9a73439c.png"},{"id":87012452,"identity":"95898c4f-3200-4e13-bc72-f73e47c10cda","added_by":"auto","created_at":"2025-07-18 09:37:31","extension":"png","order_by":6,"title":"Figure 6","display":"","copyAsset":false,"role":"figure","size":194824,"visible":true,"origin":"","legend":"\u003cp\u003eColor contour plot with 500-minute intervals of total execution time for a network of 1,000 stations. Black dashed line shows the minimum run-time for each redundancy value as a function of stations per cluster. Black contours show the overlap station count per cluster, with redundancy values on the y-axis.\u003c/p\u003e","description":"","filename":"floatimage7.png","url":"https://assets-eu.researchsquare.com/files/rs-7096364/v1/baab0e8410375972dfd51a14.png"},{"id":101152490,"identity":"e200ebe4-70d2-490a-9606-20cfc0def772","added_by":"auto","created_at":"2026-01-26 16:12:03","extension":"pdf","order_by":0,"title":"","display":"","copyAsset":false,"role":"manuscript-pdf","size":2061666,"visible":true,"origin":"","legend":"","description":"","filename":"manuscript.pdf","url":"https://assets-eu.researchsquare.com/files/rs-7096364/v1/9fdd3973-80cb-43f5-8544-575057ee02bb.pdf"},{"id":87012141,"identity":"5e23280b-c5bf-4b1a-95b6-a3b7080eaf35","added_by":"auto","created_at":"2025-07-18 09:29:31","extension":"docx","order_by":0,"title":"","display":"","copyAsset":false,"role":"supplement","size":254188,"visible":true,"origin":"","legend":"","description":"","filename":"SupportingInformation.docx","url":"https://assets-eu.researchsquare.com/files/rs-7096364/v1/d577553459a45ccf37465e27.docx"}],"financialInterests":"No competing interests reported.","formattedTitle":"Efficient clustering of GNSS stations for processing using double differences","fulltext":[{"header":"1. Introduction","content":"\u003cp\u003eGlobal GNSS networks providing open data have grown geometrically since the early 1990s from a few hundred continuously operating stations to over 20,000 known stations today. This number is even larger when considering intermittently occupied sites for survey campaigns and private stations that do not share their data openly. The vast quantity of GNSS observations generated daily is essential for numerous scientific, engineering, and mapping applications, requiring automated processing to manage the computational burden efficiently (G\u0026oacute;mez et al. \u003cspan citationid=\"CR9\" class=\"CitationRef\"\u003e2024\u003c/span\u003e).\u003c/p\u003e\u003cp\u003eGNSS double-difference processing techniques remain a standard approach for high-precision positioning and network adjustments. To our knowledge, network design for large scale differential processing is a topic that has not been thoroughly discussed in geodesy. Yet, one of the main challenges in processing large GNSS networks at Ohio State University (OSU) has been efficiently dividing the network into subnetworks for processing in GAMIT/GLOBK (Herring et al. \u003cspan citationid=\"CR10\" class=\"CitationRef\"\u003e2018\u003c/span\u003e). Currently, our GNSS database has over 6,100 stations adding up to ~\u0026thinsp;19M station-days, with a typical GNSS processing project at OSU containing about 2,000 simultaneous stations. Current computational power restricts double-difference processing to a maximum of approximately 80 simultaneous stations per session. Additionally, processing time increases as \u003cem\u003es\u003c/em\u003e\u003csup\u003e2\u003c/sup\u003e, where \u003cem\u003es\u003c/em\u003e is the number of sites in the processing session, creating a significant bottleneck in the computation of large network solutions. As GNSS networks continue to expand, this constraint poses real challenges for large-scale analyses, requiring partitioning strategies, hierarchical processing schemes (International Organization for Standardization \u003cspan citationid=\"CR11\" class=\"CitationRef\"\u003e2020\u003c/span\u003e), or alternative methodologies to integrate observations from thousands of stations without excessive computational cost.\u003c/p\u003e\u003cp\u003eThe increasing size of GNSS networks requires scalable solutions that balance computational feasibility with the need for accurate and consistent geodetic products. Addressing these limitations is crucial for maintaining the reliability and improving GNSS-derived positioning, velocity fields, and global reference frame realizations, such as the International Terrestrial Reference Frame (ITRF, Altamimi et al. \u003cspan citationid=\"CR1\" class=\"CitationRef\"\u003e2023\u003c/span\u003e), as well as regional realizations like the Geodetic Reference System for the Americas (SIRGAS, Alves Costa et al. \u003cspan citationid=\"CR2\" class=\"CitationRef\"\u003e2022\u003c/span\u003e). In this paper, we present a new algorithm we term \u003cem\u003eqmeans\u003c/em\u003e clustering, based on the sklearn (Pedregosa et al. \u003cspan citationid=\"CR15\" class=\"CitationRef\"\u003e2011\u003c/span\u003e) bisecting-kmeans, to partition a network of GNSS stations into subnetworks of fewer stations for efficient processing. Our results show that subnetworks of 20 to 30 stations provide an optimal partitioning to minimize the processing time while slightly increasing the scatter level. In contrast, subnetworks with 40 to 50 stations increase the computation time with a slight improvement in scatter. These results establish a foundation for selecting subnetwork sizes based on the objectives of a GNSS processing project, enabling a balance between processing speed and solution precision.\u003c/p\u003e\u003cp\u003eThroughout this work, we refer to subnetworks of GNSS stations as \u0026lsquo;clusters\u0026rsquo;, and stations that are shared across clusters as \u003cem\u003eoverlap stations\u003c/em\u003e, in contrast to the more typical \u003cem\u003etie station\u003c/em\u003e language frequently used in GNSS literature. Typically, \u003cem\u003etie stations\u003c/em\u003e refer to an exchange of stations between clusters, while our use of \u003cem\u003eoverlap\u003c/em\u003e refers to expansion of a given cluster to include stations from neighboring clusters. The distinction here is in the mutual reciprocity implied by \u003cem\u003etie stations\u003c/em\u003e\u0026ndash;i.e., \u0026lsquo;tying\u0026rsquo; clusters A and B together with \u003cem\u003etie stations\u003c/em\u003e would involve adding \u003cem\u003eN\u003c/em\u003e stations from A into B, while also adding \u003cem\u003eN\u003c/em\u003e stations from B back into A. In contrast, \u003cem\u003eoverlap stations\u003c/em\u003e may or may not have any reciprocal relationship between specific clusters: if cluster A \u0026lsquo;overlaps\u0026rsquo; into cluster B by \u003cem\u003eN\u003c/em\u003e stations, cluster B may then instead overlap into another cluster that isn\u0026rsquo;t A. We use \u003cem\u003eoverlap stations\u003c/em\u003e because we are able to define guaranteed algorithmic performance when discussing the processing behavior in the context of single clusters and their expansion, and further note that the number of \u003cem\u003eoverlap stations\u003c/em\u003e for a single cluster can be conceptually thought of as half of the number of analogous \u003cem\u003etie stations\u003c/em\u003e.\u003c/p\u003e\u003cp\u003eThe remaining sections of this manuscript are divided as follows: Methods describes the \u003cem\u003eqmeans\u003c/em\u003e clustering algorithm and the methodology used to connect neighboring clusters. Results show a one-year\u0026thinsp;~\u0026thinsp;1,200 station network GAMIT run using two cluster sizes and three levels of station overlap, and the impact of cluster size and overlap in the weighted root mean square (wrms) scatter of the solutions. Discussion develops and discusses a model to predict the execution time of a double-difference processing session of \u003cem\u003es\u003c/em\u003e stations\u0026ndash; and allows us to understand the trade-off between cluster size, shared stations between neighboring clusters, and session execution time.\u003c/p\u003e"},{"header":"2. Methods","content":"\u003cp\u003eSegmentation of points in space is a problem typically addressed using unsupervised learning, specifically unsupervised clustering. The problem is unsupervised because we do not have an \u003cem\u003ea priori\u003c/em\u003e set of labels to train from, and need to accomplish the segmentation into clusters using only the information inherent in the location data of our stations at run-time. This also means that our cluster labeling and assignment is independent between runs, although the inherent spatial structure may lead to some persistent clusters that are forced to reform due to geometry constraints imposed by coasts, national borders, and other human or natural factors that influence station placement and density. For our particular use case in defining subnetworks of GNSS networks, we have an additional constraint that is not typically present in other clustering problems: we need to define our partitions with overlapping point membership between the subnetworks, so that the reference frame maintains consistency between adjacent subnetworks.\u003c/p\u003e\n\u003cp\u003eExisting clustering techniques include k-means variants, spectral methods, agglomerative techniques, and density-based methods. None of these techniques address our problem completely. The k-means techniques split points into groupings of equal variance, minimizing the within-cluster sum of squares (termed \u0026lsquo;inertia\u0026rsquo;) across the dataset (Forgy \u003cspan class=\"CitationRef\"\u003e1965\u003c/span\u003e; Lloyd \u003cspan class=\"CitationRef\"\u003e1982\u003c/span\u003e)\u0026mdash;however, k-means requires that the number of output clusters be specificized in advance as the \u0026lsquo;k\u0026rsquo; parameter. Spectral clustering similarly requires an \u003cem\u003ea priori\u003c/em\u003e number of clusters to be provided in advance, and additionally struggles as the number of clusters to segment grows (Ng et al. \u003cspan class=\"CitationRef\"\u003e2001\u003c/span\u003e). Agglomerative techniques initiate each observation as a cluster and iteratively \u0026lsquo;link\u0026rsquo; neighboring clusters together as merges; while these methods are hierarchical, they also tend to lead to unbalanced cluster membership, a behavior that is at odds with our downstream processing goals. Density based methods include popular algorithms such as DBSCAN (Ester et al. \u003cspan class=\"CitationRef\"\u003e1996\u003c/span\u003e; Schubert et al. \u003cspan class=\"CitationRef\"\u003e2017\u003c/span\u003e) and OPTICS (Ankerst et al. \u003cspan class=\"CitationRef\"\u003e1999\u003c/span\u003e), and can be conceptually thought of as convolving a density kernel over the input data points, and then passing a plane through the output heatmap topology of density where the isolated \u0026lsquo;peaks\u0026rsquo; form distinct cluster objects\u0026mdash;however, this leaves a subset of points (in the valleys) that are not assigned to any cluster and are simply labeled as noise. Other clustering algorithms such as Mean Shift and Affinity Propagation have poor computational complexity, and none of the above methods address the need for overlapping membership between clusters.\u003c/p\u003e\n\u003cp\u003eAs our base clustering model, we selected \u003cem\u003ebisecting-kmeans\u003c/em\u003e clustering, a hierarchical partitioning variant of traditional kmeans (Steinbach \u003cspan class=\"CitationRef\"\u003e2000\u003c/span\u003e), which we modified and improved to adjust it to our needs. More specifically, we choose the variant from sklearn, which itself uses and modifies \u003cem\u003ekmeans++\u003c/em\u003e (Arthur and Vassilvitskii \u003cspan class=\"CitationRef\"\u003e2006\u003c/span\u003e). Our algorithm modifies the existing sklearn implementation with a new set of termination conditions which eliminates the \u003cem\u003ek\u003c/em\u003e parameter requirement of specifying the number of output clusters at algorithm run-time. Since the \u0026lsquo;\u003cem\u003ek\u0026rsquo;\u003c/em\u003e in \u0026lsquo;kmeans\u0026rsquo; refers to the number of output clusters, we term our modified clustering algorithm \u003cem\u003eqmeans\u003c/em\u003e to distinguish it from its parent algorithms. Our \u003cem\u003eqmeans\u003c/em\u003e algorithm conceptually operates opposite to the agglomerative techniques: rather than starting with all points as independent clusters that are iteratively merged, we instead start with a single cluster of all stations that is iteratively bisected. We then supplement our new \u003cem\u003eqmeans\u003c/em\u003e algorithm with two post processing steps: an \u003cem\u003eovercluster\u003c/em\u003e routine to expand cluster membership to adjacent station clusters, and then a \u003cem\u003eprune\u003c/em\u003e step to filter subnetwork clusters that are fully redundant after cluster expansion. Note that we adopt the convention of italicizing names of \u003cem\u003efunctions\u003c/em\u003e, \u003cem\u003eclasses\u003c/em\u003e, and \u003cem\u003eparameters\u003c/em\u003e that we use, inherit from, or modify.\u003c/p\u003e\n\u003cdiv id=\"Sec3\" class=\"Section2\"\u003e\n \u003ch2\u003e2.1 Qmeans cluster formation\u003c/h2\u003e\n \u003cp\u003eSince \u003cem\u003eqmeans\u003c/em\u003e is a modification of \u003cem\u003ebisecting-kmeans\u003c/em\u003e, which is itself a variant of \u003cem\u003ekmeans\u0026thinsp;+\u0026thinsp;+\u003c/em\u003e\u0026thinsp;and kmeans clustering, it is useful to review the parent algorithms to understand how \u003cem\u003eqmeans\u003c/em\u003e works. As mentioned above, kmeans works by dividing \u003cem\u003eN\u003c/em\u003e points (i.e., our stations) into \u003cem\u003eK\u003c/em\u003e non-overlapping clusters \u003cem\u003eC\u003c/em\u003e, and minimizing inertia which is defined as:\u003c/p\u003e\n \u003cdiv class=\"BlockQuote\"\u003e\n \u003cp\u003e\u003cimg src=\"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAoEAAACACAYAAAB0kLM0AAAAAXNSR0IArs4c6QAAAARnQU1BAACxjwv8YQUAAAAJcEhZcwAADsMAAA7DAcdvqGQAACVUSURBVHhe7d0HvFP1+cfxn3+kVhFQa5EhUmSPFlAQxWKFapkCFmSIgyJFEBxVFFCmioCMgiJLVMDBUtBWcCBSLWW5FyhgxQooArZQarUq55/PwzkQQhKSe5PcXPJ9v17h3pycrJNc8uT5/Z7nd4wX4kREREQkp/yf/1NEREREcoiCQBEREZEcpCBQREREJAcpCBQRERHJQQoCRURERHKQgkARERGRHKQgUERERCQHKQgUERERyUEKAkVERERyUFasGLJ69Wp3//33u08++cR17NjR9ezZ061atcotWLDA7d27123evNk1bdrU9evXz5UoUcK/loiIiIjkVZFhIf7vBWLbtm0W7A0aNMjt27fPDR061H399ddu+/btbvjw4a5Tp06udOnS7pZbbnE1atRwtWvX9q8pIiIiInlV4MPBn376qdu5c6crVaqU/dyzZ49l/wYMGOBOOOEE26dkyZLuuOOOc1u3brXzIiIiIpI/BR4Ennbaaa5r166OUen169fb+f79+7uiRYv6e+wPFHft2uXOOOMMf4uIiIiI5EeBB4Fnnnmma9iwoVu3bp0Fgc2aNXMnn3yyf6mzIeI33njDsoEVK1b0t4qIiIhIfmRNdfCGDRssCLzgggvciSee6G/dnwVctmyZa9KkiatUqZJty4JaFhEREZFCLWuCQDKBqF69uv0MvP/+++6dd95xzZs3twzho48+6t59913/UhERERHJi6wIAnfv3u3efvttq/ylEjjw7bffWhaQIeNzzz3XhoaXLl3qypQp4+8hIiIiInmRFUHgli1b3Nq1a21uYNmyZf2tzhUpUsQVL17cHX/88e7LL790EydOdK1bt3annHKKv4eIiIiI5EWB9wkEGb///Oc/rl27dq5q1ar+1lCE+n//Z/MAmQO4Y8cOV65cOdeyZUsLCkVEREQk77JixRDwMI455hj/3KEIEv/3v/+5YsWKWWAoIiIiIvmTNUGgiIiIiGSO0moiIiIiOUhBoIiIiEgOUhAoIiIikoMUBIqIiIjkIAWBIiIiIjlIQaCIiIhIDkpJi5gffvjB+vxlI/oKqregiIiIyKFSEgT26tXLlnWLbPZMk2fW+02F4Lb5GXkiyONnpO+++85Vr17d3XjjjbbaiIiIiIjsl5IgsFatWm7dunX+uYPOOeccd+yxx+Y7S8j1CSYJ6ggsgxNLze3du9e2x0KAOG/ePNehQwd/i4iIiIikJAh89tln3eWXX+7+/e9/+1ucBX9TpkxxPXr0yHcQGASA//3vf92ePXvstHv3bltPeNu2bW7Lli3uH//4h/vwww/dBx98YAFiuCuuuMLdf//97qSTTvK3iEi68XcbZOtFRCT7pCQI/P77793YsWPdwIED/S37nXbaae5Pf/qTZQTTjbWF//73v7v33nvPrVixwi1cuNCCQ5xwwgnu5Zdfdg0bNrTzIpI+b775pv3d88XtX//6lytRooS79NJLM/L/gIiIJK7IsBD/9zxjyJUh4U8//dS9//77/lZnw7XvvPOO++1vf2uBWDoVKVLEnXrqqfY4LrroItekSRN33HHH2eP5+uuv7YPowgsvtAyliKQe3yeXLFliX8AuvvhiV7t2bVe5cmX30ksvuUmTJllQ2KhRI39vEREpaCnJBAYI+K688krLxgUYCurevbubMWOGvyVzqFpetGiRu+WWW9w///lPGyouX768f6mIpBKZd4qwHnjgAVe6dGl/q7NsYOPGjd3WrVvdzJkzXZs2bfxLRESkIKUkExjgP/6TTz7ZvvmHz8vbsGGDzcdr0KCBvyUzyFDWrFnT7vf555+3bRdccIHmKCXh888/t+P4ox/9yN+yH98dyPxyGRnXWBiiL1q06GHXl+z02Wef2Zen448/3t9yKF53OgGQ5T/xxBP9rfs9+eST7r777nMlS5Z0v/zlL/2tzv34xz8+8IWMbGDHjh39S/bjPYZ47yMREUm9lDfQ69y5s+vZs6cFBwE+MP74xz+6VatW+Vsyiw+kiRMnWgELBSWSmDVr1rjbb7/dffLJJ/6Wg3bu3Ok6depkc7/iadmypVVnS+FAEdWYMWP8c4f75ptvLLM+fPhwf8uh+PLH+yZSpUqV7CcBZKSlS5e6IUOGuC+++MLfIiIimZCWLsr9+/e3eXnhNm3a5Eg6Us1bEBiCOv/8891TTz3lb5F41q5d60aNGmVBPXO7IhEMUIxDxXYsZH/YR4F34cH8WarvYyETyPAuX+witWvXzk2ePNm+OEQia4xTTjnFfoYj8KSIbMCAAW7Xrl3+VhERSbe0BIE/+clP3MiRI90ZZ5zhb9nvxRdftIwcAURB4EOmQoUK/rnsxLxKjtHbb7/tb8k82u0QsFNc06xZM3/r4Y40rB60B9Hwe+GRyOsVax8CvGuvvfawaR8EjnPmzHHFixd3Xbt29bcexKhB3759bQh66NCh9uVBRETSL23rqZ111lmWSYqsxmXO0JGGENOFbENkhjKbMBzGh+BNN91kxTTxMjLpQmaPAL5YsWK2EoxIfj300EPu9ddfd3369LFOAdEwv/CGG26w9jKzZ8/2t4qISDqlLQgEQ4kENOHIAjJcTGNnORTNtoPh8vXr19tQaqa9+uqr7rnnnrMPbBVzSH698sor7q677rJM3+DBg62VUyw1atRwl112mTV2p5hMRETSK61BIENGd9xxhw0rhtu8ebN962fJNzmoSpUq7g9/+INr3ry5mzp1atT5U+lEgH7nnXdaU216KorkB1k95gfyt3733Xcn1Cu0VatWNhw8f/58f4uIiKRLWoNA0Bpm/Pjx7vTTT/e37EdF4D333OOfk0CXLl0sE3f11VcfUmGdCWQBWQOa+xbJj7feeutABpBqYrLKVA7TMiieqlWrWqPpWbNmHXFfERHJn4xEGXXr1nUjRow4rA8YbWNUrZs9aPJbpkyZuMUgIkfCF4lHHnnEAkC+1ARoJk2m+UjIQjM/dvXq1f4WERFJh4ylmjp06OB+//vf++f2Y/iRil0qYqVgMReRPo7nnXde3HlbIvEwl5UhYCqBqfZlze5ly5ZZZwCaSVNwdCQUlVWvXt2+IEZrRSMiIqmRsSCQ+UC33nqrrdgRjv6BNIplWbdsQ2EGqxnwwbZy5Ur7MKNJcoDg9fHHH3cLFiywfaJ57bXXrD0G66l+/PHH/tbDcV8EYmRRuC+GhMNXXQGVu2RTWA+Zodvly5cfKB7hJ9d77LHH3NNPPx21wXM8PDeCcobiCgIf9qw0w+Nneb99+/b5l+y/jBVfnnjiCSsYSOFKh1mNoDzW+yobffTRR+6qq66y9y7DufQNpDiMbCCNxRkerlatmr93bGXLlrXm0i+88IKCQBGRNMropDP6BjIcFNk/kJU8JkyY4J/LHrSyad26tQ2Ptm/f3ia40/SWPnoDBw60oIVl8tjGUlg9evRw27dvt+sSyFDkQRDIvEjOt23b1l1//fVRV03gg5OVTZgYz23Rb41q4XAM1zJUxiocZFZ5DAREtN8go8qazTwe7pPiErbFa+Yc7i9/+Yst6XX22Wf7WzKH+WPBsWWJMR73oEGD7PGwjTllBOMlSpSwLwxjx44tkMrpTPrrX/9qbYJ47rEQ6BO033vvvf6WgrVjxw7L4vF3wN8N7WCCE+9pWg4lWnDE3EACQFUJi4ikUSiIyLhx48Z5oQ970jkHTqHgxVu0aJG/R3YIfQh5n332mbdgwQLv2GOP9cqWLes98sgj3qhRo7y//e1v3jfffOP98MMPXijQ8iZNmmTPI/Sh7a1fv94bPHiwFwpubJ99+/Z53377rTds2DDbZ+jQoV4oiPHvZT/ua8uWLXbbxYsXt/sLBZT+pfvt3r3bC33we6EgyG6nRo0a3rJly+z2uE9uk/tiv759+9o+d999t3/t+Bo1amT3y+NNRCgQ9urVq+dNnz7d33I4js2ZZ57pjR8/3t9yuLffftu78sorvU2bNh04Jk8++aRXqlQpb+rUqV63bt28VatW2fYZM2Z4xx9/vD3OL774wrYdrW699VbvmGOO8bp06eJvOVzoC5W9xhyjVGnQoIHXp08f/9zheJ+GAjzv6quv9rccxOvH+4f3eqwT789EzJo1yytWrJg3ceJEf4uIiKRaZstPfTfffLO79NJL/XP7MRxMljCbhr8Ywqaqmazbqaeeapm5UPBny6iFgiYrdKGCl+wVGTzOP/PMM/Y8OE9BDNtolUN1ZNOmTd1Pf/pTG0Jmof5w3Fe5cuVsCK106dL+1kORCfvZz37mfve739l9Mnmex8OxZA5V0aJF7b7YLxRY2XWYkxUs0B8LrXo4/jTT5jYyhSXKKBhidRKG/4L7JpvE3DEyqaEg0p177rm2naFyMpscV57/0YzsaOjv07Vo0cLfcjgy1cy7C45PQeP14/3Oez3WKdpKI9HwXuRvQv1ERUTSp0CCQFAZ/Itf/MI/tx8ffKNHjz5sGDQbEHTwuDgR4EXicvr8MceRYO6cc87xLzmIoVoCPPokxprrxAorR2oNwwcp+7GGKwFcnTp1/EsOKlmypO3H0HMwRB0LlzMfkLlYmTRlyhQbHiTQC8e8S4Z76RcX3q6GwJahT1ag4PkdrWiNEgTukT02A7xmzMHjOETOsz0a8GWJv6nIL0siIpI6BRYE8k2flQEiGyIzoXz69OmHFAZkA7IyBF5kAaMh4AoyWWTlomU8qLpln++//z7f66PyeMj4xZrDRyDJ/RFMcYqHwJb5d2TfEs3UpAKZzGjLiFFwQ4D785//3JUvX97fur/VEMVFBNux/PnPf7ZMa36Pb0FijidFQrVq1bJ1uKMhA0xWlC8crLSRCG6XJuRfffWVvyV7sYwcfyvhhVgiIpJaBRYEgkII1sqNxPDgmjVr/HPZgaCL4axSpUr5Ww7HPmQvCM6i4XJOyG+Qy+0w9Bbv8SD8PmMhC8jjYWgxkyh6IOMTif5wZEp/85vfJBWUkjW67bbbrHgk2SCQzBrTFLp27WoZx/ycqIilQTqBbF4QrDE8TxYw+GIRiWppniND54liLe/evXtHLUzKNjxvvsjw3hQRkfQo0CCQ/+SpgIzsH8iwKXPwsg0BSbweegRbXE7GMBN4PKm8ryMFi6nG0Hjk8SQruXHjRvs92abVZM2oKqZRcazgKRYCajK0BJ/5PTFEz2uT7GPAt99+a1MKwHzAWK8vFeDcB/NME8XfGcuxVahQwd+S/Y40NUJERPKuwP+HZdiHIb5gDh3tVGbOnBl3yC+bJZJ5S6VU3BcZQIIxCjUyHQhGos0NvRDJpsYaeo+FQgIyecyRSyaDCN53zDekFyG9CvNzmjdvnq2WkUhj5EgEwMwJJPiL1VOPdkNk83i+F110kb/1yAiqL7vssoxnfPOCKQxkp8m+i4hIemTF12yyHwz78ME9adKkpLIbkn+s7kDQQRBY0AgC6Q3YuHHjqMEKmTYydpGYH8dl+QliGconsMrvieAvL1lAUB3PiS9FsYpfaArOUDNzBoOMOX9DnKLh2JBhzbZ5tvGQTWWearTpAiIikhoFHgRu3brV9evXz7377rtu5MiRtrKAZBbrBRNw0QS7oDOBrIbCXLdf/epXNkQbjsCHoV7eM+FYiYVikDFjxliT5cIU7EQiACQL1qBBg5iZxBUrVljAHp4FXLp0qTX8DsdxpOJ+7ty5NkeR9kWFZY4dBSE81sI0dC0iUtgUaBDIt31Wh2B5KFbSYAgtU/PpkpXs8GJhQgDIfDpWfIiVTUo1KoMJTqhyDRCEEgSCIChyvuDixYstqxU+X5TrTJs2zfoMEjiwbFlhKHyIZvfu3db2BbTNITMZiecY9M779a9/bT937dpl2yKDZo4tw9tMraCAiFVoCJILA94fBLqJVj6LiEjyCiwIJOPEQvPMoaJNyODBg/M0hypTGJpCrGCQ7UcKFJnkHkx0j1VgwuXB7cQKiNkeZOxi3U6wPZHHBYYWyaDRwzAT7rnnHltT9qabbvK3OMsGc0Jk6yAqf2kdwzBx8D4hY0aVLHPnqKTledJaJdpwcWHAMDiZOxAARWY0ec3pkciQOcPNQYDEPEKCw/Cm0RwH9u3Tp49V4ZNRI8Bne2FAcQzD+9H6bYqISGoUWBDIqhr0CeSDi8bR2Tj3hw9dgj+G6JhTRRaKD2CCj6AFCR/U7MOQJAEUH96sdxreC5DbYT+2B0OZ3A4fysEHfbAPQQDtQbj+2rVrD7RvAT+571WrVtl9kkl944037H6CffidfQimuA2CA7JL4ftEw/ArgUUmWvPwuIJgL2i8zYc+z/fGG2+090JwOchQku2jJ2J4oMOx5nnRloViEo4LQ6SFtZE0759gxZyFCxdahi/Ac73vvvvs9/r161vAS9EE7wPeA6zHHZ45JKBs06aNrTADMqx79uw5cD7b8Z6ll2hhebwiIoVSKPjIuClTptgasBUqVLA1b7PVnDlzvNAHrlelShUvFFh4JUqU8MqUKWPnWZuXNXpnz57t1alTx6tcubKtacupfPnyXrVq1Wyf0Ie39+ijj9p51tHldtgn9OHm1apVy/YBa6WyrWLFit5JJ51k98XxOeuss+z6GDlypFe1alXbh8s5cV+cZz1mDBgwwKtdu7Zt4364P26XxzxhwgQvFDTYfpFYp5j1euOtVRsuv2sHP/HEE94555zjDRkyxNZL7tGjh7dixQq7LBTweaGAz9bGZe3jq666ylu8ePFhj53bDwU29vvDDz9MatSbP3++nS9sQgG6d++999pz6Nixo9e8eXOvSZMm9rryunfv3t3eI6EvBbbeMseeYzZw4EDvjjvusHV5w/G+4/ggFHTbMQwF197GjRttWyz5WTs4VUIBrD0/HjPPQ0RE0uMY/rFoMENY77Rnz56W4WIoL1vWPY2GTB3ZP7IuwdAs2TQyWWRhqAQl68Y+YB/2DbJuZNaCfcjkMNQb3E6wD+dpk8N9MfzF+WAol/vh5WHOHvO9uB+OG5cH+3A77Mc+nLgNskM8DvYJHg8nLo82zyzAsHwowLBsVOT8skgMz7Zt29aaD0f2eQzw/JiPxlxP1gEOx/Mik0lGFKwOQpUyuIwsGK1QOH6swMIQcKyecRzbUEBkBROslRwKyO02eO6FBat4XHPNNe7pp5+2eXwMlTPMy3EmsxkK/u1YBK87rzOvFceM4xOvlQpZ6lBgadnCULBsfThjYfiVE1X60XCsKd5iDimtnNKBY8BygaEvWPYeExGR9MjocDANbikEYbiT+UrZHACCQIj+cXwIE4Rw4kOXD1F+J8gI9uFEMMd2etxxPnwfrhN+O8E+XCe4L+bBsT3Yh/3ZJwjICOK4nfB9+J3rcRnYxnW4bvjj4XrxAkAwf4xghKXX0o3jwrFs1KiRnYIAEFxG8QdD1BSIcFmsABDMFeQLBcEJw4f02QsKLAoLhu1fffVVe/wVK1a0bQTQtEtiGJxjEASA4HU9//zzbf3tI/XSY6idY8TxjBcAZotly5bZMYi1JKKIiKRGxoJAPphZJowsE6089A0/+9BkmaCL9ZsLE+YC0j6GZebIpE6YMMGyjIUJ8wEJwMmIEvylEnNRyQSTMcx2ZIYJAskEnn766f5WERFJh4wEgQztMRxIe4r+/fvb0F0wLJppDKdKdAxfU6VNuxEya4UBw74ETwQ4DAO/9NJLlhWtU6eOv0f2I0AjSw4qfhlqTZVg2JhjU758eX9r9qIghkw2UxNERCS90h4EEnSxLBxzAaniZDiYuU0Fgb5pM2bM8M8d3ZhjSK81qpojEXTEct5557nOnTtb5Xa062Ybho4vvfRSV7duXWs2TnU1bWfyumJHQWCu54svvmjDurGWissr3gN8+WKeXzDMnK3efPNN+3+C109NokVE0i/tQeDw4cNtlQca27I2K9/yCwJFCnfddddR3XKCZsNMpm/Xrp1lUoYNG+ZuuOEG161bNyuaAMPyFCDEQvBEIQfzCSdPnuxvzW7MHWMImOdKkQoroBQmzJHl9SFIS3UGk9Y5FJgwd7Cgvnwlgvfu2LFjbUifIhYREUm/tAaBBBGjRo1ytWvXdlOnTnWlS5f2L8k8HgcrSSSz4H5hQXUwTbdppPzss8/a0Dvz+kaPHu0mTpzoBg0aZJWcZEGpOj3SfDn6sw0ZMsStXr3aLViwwN96OALrRDBkm26siMH7rKC+ZOQHx6dcuXLu4osvzncxBNnuJUuW2HuCE9XS9F1s2LChv0d8ZImP9Lomsk8yGC0YN26cffHgi5qIiGRG2oJA2jww9Mvk7unTp9ucpILAByHtLlihokOHDgcqbbMRj5UqUdqC0Ow3OHGeuV3RkEW67rrrrOiGRstz5syxKlCqa6kQpjKYY0+WhUCcuWfNmjXzrx0bK4iQxWU92miVtsy7Y45Z5Moe4ajo5fUvyOC/MKB1C69z0Aw6r8jyUuFNNpT3BQUWNNlm24UXXujvFR+Z8nhzB3lNWdKOoDVVGAonE8iXxoKaKywikovS0ieQoa3LL7/cKjbJSDFnqyCwrixz2wh+CLBYNYE1ZrMNff3IuvGhTZVoeMBHxoVAjg92WqmE4/gSADKPiiH3yy67zL8kOgLBu+++21aYSPQ4EJQyjBhtKJHjS8uR8PYukT755BNXtmzZrA6+jxa8VjfffLMFaiwVt3TpUmsjw3sk0eCKZeV4PWO9pvx3QYDJ+zJ8Def8YMiaoJL5nSIikjkpDwJp8dC1a1f7SRaC4CXT/7mTEaG69dFHH7XgCvSQYzg06MuXLRgKY64kVa0MZ9I7kTltQe83Xh4CKJoFh1eNEjgOHDjQhtEY7k1kGI0gk+psjokyLkcnAjQKLCg2OeussywAFxERiSalQSBZLNq/EGwwH42sRCaCDQIpAj8yXFRC0hKDE4FSgGzZJZdc4p9LHD3oaNLLeq3hDYv5kKXyMj9ZLrIp/fr1c+vWrbMMHfeRKIJG5pBxHYbeExme43V54YUX7LVR1kVERCS3pSwIpFcb1ZmPP/64tRmhzQPDmAzDpgIPk6COtiXMH+L+GI5kLhXDomxjOIxF8iOfEhWXixYtSrpFBv0NGWKlUIJh7fCVGZi/xVzHZ555Js9DzCwPNn78eHtsyfRwI3ikwIX5erfffrsFkIngOLGEHUUUIiIikttSEgTSk27o0KGWYQLLk9FqJEXx5QEEP8GJCsV4/e7C0SqFJsjxlh6LhmwbQ9vXXnutPb9g2S4CW+Y8UolJppCJ/bEQmG7fvt0KLcLxHOjdRhCXbGNclgGrWbOmFWU89dRTNu9LREREJBkpqQ5+8MEHDwSAYHiWjB2tSFJ5YvF6bptsVqIBIBWyTJJPNgAElbS0lWnevPkh67bSd40h3ETWYmW+HlnRSB988IENKbNUG88r3on9CBoDPC6yokzMZ5k3ERERkWSlJBNIlo1h2WybZ0agRLEFxSn0SksW2b7nnnvOFt8Pz/YxB4+KZ6qOyRLGQ7BHAMpyYOHo58ecSeYphgeYkQj+KGZp3769rSsLGiPTC5AglCHhRDB/keFsglbNBxQREZGUVwcfLShyYT4gwSNz98KbELM8GcO4r7zyimXyYiGAi5WBpDcarTu6d+8eN0vJbdCTj5UUgiHlhx9+2ApwmBfI7SSCQJ1CkniPV0RERHKHgsAYFi9ebM2lKbogYxdkzxjm7tWr14F1Tik2IctGU+QAzZWpUqawhOIYsoWRgR6ZU5o2czvJLudFv8N69erZerkMDR8J2UgCR4alkylAERERkaNXSuYEHo3Wr19v8w8pwAgfPmX1DgpCqIBmdQX6EdKaJrBhwwYbQmZVBfq0jRgxwq1Zs8a/9CAWyKf3H0O7yeIxtWjRwh7LypUr/a3RUZhCZTNrNysAFBERyRwWz+jcubN/LjoSRvQLZnWvYJ3/aObOnWv7sX+qKAiMgkIMgjlQhBJgJQ8aUJPFIxCjApoWMSzOH2D+IM2emzZtakUpVE4znBsNFccEaAwrJ4Os4pgxY6wK+84777TilWg+/vhjWy6Px0dxi4iIiKQfgRrBH4kagrdo6G9MUFelShVbyIHP7Hi4Peoc+Dwn1kiFIsPonyKHIMibMmWK+/zzz63vIEO9O3bscDNnzrTltGjNQnBIsEgzbLJsQVNssnsEiARqBHhkFIcMGRK1aTaNpqnwveOOO+x2yQ5SvJFI4QbXY54gBSYLFixwxYoVs/ukXyLzGZnHuHDhQnfllVe6Vq1aqRhEREQkAwgACdSY8x+reJSM3/z58y2Zw++MPILrxGs7R8DYsGFD16RJE4sZCArzQ3MCo3j55ZctsKMAhGXceIEI0lq3bm2Noen7R4BFYEUfwVjrrFK4QUEJQVq8IIzijuHDh1vlLu1swlch4eUhMGRlkGhvDDKNr776qluyZIkFrVQaM+xLgMibRMGfiIhIZgQBILFCrAxgJLKBZAJBPQFxwJFw2126dHFz5sw54nBzPAoCo3jggQdc375987zUHMgcEkDywvbs2dPfGhtrvpLVo/k0b6KgLyA/yfJRDazKXhERkexFALh27VrrJxy+3n88eQkCEdzX8uXLbWW0vDiq5gSmIp6l+nf16tXW14/ijryiapdGz/TySwRZQIZuJ0+e7ObNm2fZQ06sCDJ79mwFgCIiIlls2rRptj4/mblEA8D8YDEKEkgkifIq7UEgc+po1pzo+rZ58eGHH7rbbrvNDRgwwPr30dKFNip5QRZu2bJl1n6lWrVq/tbk0cSZHoOsOywiIiJHL2KHgQMH2u8sNJEJZAwZdn799dcTHnqOlPYgcNu2bW7GjBkWqKUDlbW9e/e2yZL8JCJmTV4mY9KqJVnMxWvbtq3r1KlT1GKOWCjIIC1LNpKiEdrCMCePoVwRERE5ej300EOWlUOiw7mpwPx/kBXMi7QHgUSpmzZtsnl2qUYhBJW3rKRBtpG+fQzhssQaw7C0YNm6dau/d2JKly7t7rvvvqTmAhL4PfbYY65bt27uvffecxMnTrS1jZlXSBsZEREROXpNnz7dfrIIRCYFU8VoL5OXtjFpDwIJgggE6ZmXasy7o3KXAoxw3CdLpHF5stlAqmm5fryl3CJxnQsvvNCif94IlHpT+k12UkRERI5eBF9Bj7+zzz7bfmZK+GplxB3JSmsQSI88soDp8N1331kBB8uylSpVyt96EL388NZbb9nPdCMbSQA4adIkW6eXNi8iIiJydCMWCaQj4RUPRayB8NXLEpWWIJDhUSpa6bc3btw4GytPNRo1E2Qy5y7aihxk86i4JTpnbV8RERGRVKMmIBCv0XM6hFchU5mcrLT0CWRpNfrk0UiZql1631E0QeNk0H/vySeftAKKZJoZ0zOPOX80Vt67d681b96yZYutjhE5EZNCFAozKlasaPP18tPuRURERCQa+vUFAVgyff4Cee0TGAiPo5IN6VKeCST4IwDjoNAehUIJMnUnnniiv8f+ZU+oaGGfZE4tWrSwSZBU7RIQMveOuXvRAkm2sx/DxhRpiIiIiMhBKc8EsiAyw69MjmSyJAUT119/vRs9erS/R2rs3LnTgsIvv/zSPf7444dFzh999JFlAknNcjnFKSIiIiKpRJJKmUAfw66NGjWyfnsPPvig/WQljFQjy0fxB1m+aE+a7WQKmTPIAs0iIiIiqda0aVP/N+dWrlzp/1Y4pK06mAzd888/7+rVq+dq165tw7IEZSBTSIdr5gkmc6IB87p16yzoI7ikLyBLs1EkEon5hjRupHw6E8u3iIiISO7JdDFILPXr1/d/S1zagsA5c+bYMip9+vSxoI0ALmjXws8JEya4MWPGuLFjxyZ8Yn+KQAjwmGfIyiAEetxPJJarQ82aNZUJFBERkbQgFgls3rzZ/y0zwhtEB6uHJCMt1cGgWTNZO/oEbt++3U2ePNn16tXLhou/+uorW+2D4o5k8FApMCEDSOEHLWK4nyuuuMKNGjXK32t/FTFj7FOnTrUq5LxExyIiIiKJqFy5srWkY8UQRkGTkZ85gUuWLHGtWrWy3xcvXuxatmxpvyeMIDAdLrnkEq9Ro0be7t27vdGjR3szZ870L0mt6dOne9WrV/dCB93f4nkrVqzwqlSp4k2aNMnfIiIiIpIeU6dOJaFmp2SEAkevfv36B67bu3dvb+fOnf6lR3b77bfb9SpVquRvSU7aMoEbN260IVzmArZt29Z6+hUpUsS/NHXI+hEJz5o160DrGHTp0sW1adMmqeXfRERERJLFtDTa3zFFLZGMHEveNm7c2D93uEQzig0aNLAaC6bgde7c2d+auLQFgSIiIiK5Ytq0aTbtrXfv3jYFLt1oyUf7O6a8vfbaa/7W5CgIFBEREUmBoGcg8wPTvVLZdddd5+bOneuWL1/u6tSp429NjoJAERERkRRgWJhAkNZ0yRaIJIOq4Lp16+Z5GDigCXMiIiIiKRAEfwSDZOrSgdtu3759vgNAKAgUERERSZEgEDz55JNTHgiSAezataubOXNmvgNAaDhYREREJA2oAn7uuefciBEj/C15x22xLN0111yTspXQFASKiIiI5CANB4uIiIjkIAWBIiIiIjlIQaCIiIhIDlIQKCIiIpKDFASKiIiI5CAFgSIiIiI5SEGgiIiISA5SECgiIiKSgxQEioiIiOQgBYEiIiIiOUhBoIiIiEgOUhAoIiIikoMUBIqIiIjkIAWBIiIiIjlIQaCIiIhIDlIQKCIiIpKDFASKiIiI5CAFgSIiIiI5SEGgiIiISA5SECgiIiKSgxQEioiIiOQgBYEiIiIiOUhBoIiIiEjOce7/Aem7tj8pLP1fAAAAAElFTkSuQmCC\" width=\"641\" height=\"128\"\u003e\u003c/p\u003e\n \u003c/div\u003e\n \u003cp\u003eWhere \u003cem\u003e\u0026micro;\u003c/em\u003e\u003csub\u003e\u003cem\u003ej\u003c/em\u003e\u003c/sub\u003e is the centroid of proposed cluster \u003cem\u003ej\u003c/em\u003e, and \u003cem\u003ex\u003c/em\u003e\u003csub\u003e\u003cem\u003ei\u003c/em\u003e\u003c/sub\u003e is a point coordinate from the set of all coordinates proposed for assignment within cluster \u003cem\u003ej\u003c/em\u003e. In the original algorithm, this is done by selecting \u003cem\u003ek\u003c/em\u003e centroid centers, and then perturbing those centers (and adjacent point membership) to iteratively minimize inertia; the sklearn implementation uses \u003cem\u003ekmeans\u0026thinsp;+\u0026thinsp;+\u003c/em\u003e\u0026thinsp;which accelerates convergence by using random seeding, an enhancement which is of marginal impact when \u003cem\u003ek\u003c/em\u003e is low.\u003c/p\u003e\n \u003cp\u003eFor both \u003cem\u003eqmeans\u003c/em\u003e and \u003cem\u003ebisecting-kmeans\u003c/em\u003e, there are two nested loops iterating: an inner-loop iterating to converge on the centroid and membership proposal that will optimally bisect a given cluster \u0026lsquo;node\u0026rsquo; within the recursion node tree, and an outer-loop iterating over nodes to bisect. For \u003cem\u003ek\u0026thinsp;=\u0026thinsp;N\u003c/em\u003e, the \u003cem\u003ebisecting-kmeans\u003c/em\u003e algorithm sets \u003cem\u003ek\u003c/em\u003e\u003csub\u003e\u003cem\u003einner\u003c/em\u003e\u003c/sub\u003e\u003cem\u003e= 2\u003c/em\u003e, and splits the dataset after optimizing cluster centroid location and membership to minimize inertia as measured in (1). This bisection process is repeated recursively in the outer-loop by selecting the cluster with either the largest inertia or largest point membership, and splitting that cluster (using \u003cem\u003ek\u003c/em\u003e\u003csub\u003e\u003cem\u003einner\u003c/em\u003e\u003c/sub\u003e\u003cem\u003e= 2\u003c/em\u003e), with the recursion terminating when the total number of clusters is equal to the overall number of output clusters requested by \u003cem\u003ek\u003c/em\u003e. As a consequence, the \u003cem\u003ebisecting-kmeans\u003c/em\u003e algorithm is hierarchical, reduces run-time complexity for high values of \u003cem\u003ek\u003c/em\u003e, and produces clusters with more uniform membership. These latter two properties make \u003cem\u003ebisecting-kmeans\u003c/em\u003e especially appealing for our GAMIT processing pipeline, as having lower variance in cluster membership number allows us to further optimize for run-time.\u003c/p\u003e\n \u003cp\u003eThe base kmeans algorithm must effectively recalculate the entire clustering structure from scratch whenever the \u003cem\u003ek\u003c/em\u003e parameter is incremented. In contrast, the bisecting variant maintains the prior \u003cem\u003ei-\u003c/em\u003eth outer-loop iteration cluster centroids and membership labels in a recursion tree while swapping the current highest inertia (or largest membership) cluster node within a given outer-loop iteration for two new subcluster nodes that bisect that parent. This outer-loop iterative structure for the bisecting algorithm allows us to remove the \u003cem\u003ek\u003c/em\u003e parameter entirely, and instead replace it with our own boundary conditions. Our resulting \u003cem\u003eqmeans\u003c/em\u003e algorithm replaces the \u003cem\u003ek\u003c/em\u003e parameter with a single new \u003cem\u003eqmax\u003c/em\u003e parameter for maximum cluster size. The \u003cem\u003eqmax\u003c/em\u003e parameter is a hard boundary condition that provides an algorithmic guarantee that the output clustering will include only clusters of size at or below the \u003cem\u003eqmax\u003c/em\u003e parameter setting\u0026ndash; if a node within the \u003cem\u003eqmeans\u003c/em\u003e recursion tree is larger than \u003cem\u003eqmax\u003c/em\u003e it is bisected, if the node is smaller than \u003cem\u003eqmax\u003c/em\u003e it is not, and the algorithm terminates when there are no remaining nodes left in the tree to bisect.\u003c/p\u003e\n \u003cp\u003eOur \u003cem\u003eqmeans\u003c/em\u003e algorithm produces clusters with different properties and streamlines the implementation code when compared with the parent \u003cem\u003ekmeans\u0026thinsp;+\u0026thinsp;+\u003c/em\u003e\u0026thinsp;and \u003cem\u003ebisecting-kmeans\u003c/em\u003e algorithms. As mentioned above, \u003cem\u003ebisecting-kmeans\u003c/em\u003e chooses the next node to bisect by selecting the node with either the largest inertia or the largest point membership, and will produce a different output clustering depending on which of the two metrics is selected. For \u003cem\u003eqmeans\u003c/em\u003e, the distinction is meaningless and does not impact the output clustering at all. Deciding on which cluster to bisect matters in \u003cem\u003ebisecting-kmeans\u003c/em\u003e because the algorithm increments the total number of clusters at each outer-loop iteration and will therefore terminate after \u003cem\u003ek \u0026minus;\u0026thinsp;1\u003c/em\u003e total outer-loop iterations\u0026ndash;this means that even though the bisection of any particular node will produce the same output subnodes from applying (1), the choice of which order to apply the fixed \u003cem\u003ek\u003c/em\u003e bisections will produce different results. In contrast, \u003cem\u003eqmeans\u003c/em\u003e is not iteration bound for the outer-loop and will apply bisections, using (1) in the inner-loop to determine output nodes until there are no nodes above \u003cem\u003eqmax\u003c/em\u003e remaining to bisect.\u003c/p\u003e\n \u003cp\u003eOne drawback of kmeans algorithms in general is that they perform best when clusters are convex; because the technique bisects points into groupings of equal variance, this family of algorithms handles irregular shapes poorly. Low membership bisections can occur on our data inputs because of variable density and sparse station coverage\u0026ndash;for example, across oceans with widely separated island stations in a geometry that yields a small member cluster for the tightly grouped GNSS stations over islands such as Hawaiʻi, and single member clusters at distant atoll stations such as Guam. While both low cluster membership and high cluster membership present problems for our GAMIT runs, cluster expansion for overlap stations ensures that even single station clusters such as Guam will have multiple stations to form a subnetwork for a GAMIT session. Since clusters with high station membership can cause non-recoverable crashes in our processing runs due to memory errors, by design we provide algorithmic guarantees for a ceiling on the station count membership during both cluster formation from \u003cem\u003eqmeans\u003c/em\u003e and expansion via \u003cem\u003eovercluster\u003c/em\u003e. A final additional post processing by the \u003cem\u003eprune\u003c/em\u003e function (Section \u003cspan class=\"InternalRef\"\u003e2.3\u003c/span\u003e) is implemented to remove any redundant clusters that may be formed entirely by stations that overlap with other clusters.\u003c/p\u003e\n\u003c/div\u003e\n\u003cdiv id=\"Sec4\" class=\"Section2\"\u003e\n \u003ch2\u003e2.2 Cluster expansion\u003c/h2\u003e\n \u003cp\u003eWhen processing subnetworks within GAMIT, we require there to be overlapping stations between the subnetworks, and we use the \u003cem\u003eovercluster\u003c/em\u003e function to ensure this overlap. The \u003cem\u003eovercluster\u003c/em\u003e function takes three required parameter arguments: \u003cem\u003eoverlap\u003c/em\u003e, \u003cem\u003enmax\u003c/em\u003e, and \u003cem\u003erejection threshold\u003c/em\u003e. The \u003cem\u003erejection threshold\u003c/em\u003e is a distance set by default to 5,000 km, and rejects any points that fall outside of that distance\u0026mdash;this is mainly present as a sanity check for very isolated stations in polar regions or oceans, and ensures that GAMIT can form double differences on adjacent stations. While the \u003cem\u003erejection threshold\u003c/em\u003e was enabled to be user modified, we have had no reason to do so and it remains fixed for all of our processing runs. The \u003cem\u003eoverlap\u003c/em\u003e parameter specifies how many additional stations we want to expand our cluster by, while \u003cem\u003enmax\u003c/em\u003e sets a maximum number of stations to include from any given adjacent cluster. The cluster expansion is recursive, and uses a \u003cem\u003eballtree\u003c/em\u003e spatial index from sklearn to calculate nearest neighbors (Omohundro \u003cspan class=\"CitationRef\"\u003e1989\u003c/span\u003e). More explicitly, we take all of the GNSS stations except for the cluster being expanded and build a spatial index, using the index to calculate the \u003cem\u003ek\u0026thinsp;=\u0026thinsp;1\u003c/em\u003e nearest neighbors from each station outside of the cluster, to each member of the cluster. This intersection of cluster and non-cluster members results in a sorted list of station distances to our cluster\u0026rsquo;s member stations, and we add the external station with the shortest distance as an overlap point. We then update the cluster and non-cluster membership to include/exclude this new station and rerun the nearest neighbor analysis. As we add points from neighbors to our cluster, we also track the original label of the added stations; when the sum of a given label\u0026rsquo;s occurrence is equal to the \u003cem\u003enmax\u003c/em\u003e parameter, we remove that entire neighboring cluster from future consideration. This ensures that we get broad connectivity and representation from neighboring clusters. The \u003cem\u003eovercluster\u003c/em\u003e expansion algorithm terminates when the number of stations added to the cluster under expansion reaches that value set by the \u003cem\u003eoverlap\u003c/em\u003e parameter\u0026mdash;unless there are no stations to add under the \u003cem\u003edistance threshold\u003c/em\u003e parameter and the algorithm terminates early, a condition which was not triggered for any of the results we present given our other parameter choices.\u003c/p\u003e\n \u003cp\u003eThe \u003cem\u003eoverlap\u003c/em\u003e parameter provides us with a simple algorithmic guarantee for our processing pipeline. The maximum size for any cluster provided to GAMIT will be:\u003c/p\u003e\n \u003cp\u003e\u003cem\u003eMax cluster size\u0026thinsp;\u0026lt;\u0026thinsp;=\u0026thinsp;qmax\u0026thinsp;+\u0026thinsp;overlap\u003c/em\u003e (2)\u003c/p\u003e\n \u003cp\u003eThe \u003cem\u003enmax\u003c/em\u003e parameter ensures that the cluster overlap doesn\u0026rsquo;t heavily favor a single neighboring cluster. If \u003cem\u003enmax\u0026thinsp;=\u0026thinsp;1\u003c/em\u003e, then the number of clusters that are expanded over will be equal to what the \u003cem\u003eoverlap\u003c/em\u003e parameter is set to. Using the default \u003cem\u003enmax\u003c/em\u003e\u0026thinsp;=\u0026thinsp;2, the number of clusters that are expanded over will be bounded between \u003cem\u003eoverlap\u003c/em\u003e and \u003cem\u003eoverlap\u003c/em\u003e / 2. It is important to note that while we are guaranteed these minimum overlaps, in practice many central clusters will have more redundancy than this. Asymmetry in reciprocity means that some central clusters will have substantial additional redundancy when the subnetworks are reconciled. While redundancy improves the scatter of the GNSS solutions, some clusters are either so small or so central that their membership is entirely redundant, with the stations fully represented across other clusters. Thus, we implement a final postprocessing step to \u003cem\u003eprune\u003c/em\u003e redundant clusters and improve run-time.\u003c/p\u003e\n\u003c/div\u003e\n\u003cdiv id=\"Sec5\" class=\"Section2\"\u003e\n \u003ch2\u003e2.3 Cluster pruning\u003c/h2\u003e\n \u003cp\u003eThe output from the \u003cem\u003eovercluster\u003c/em\u003e function is a Boolean \u003cem\u003eM\u003c/em\u003e by \u003cem\u003eN\u003c/em\u003e matrix, where the \u003cem\u003eM\u003c/em\u003e rows correspond to \u003cem\u003eM\u003c/em\u003e clusters, and the \u003cem\u003eN\u003c/em\u003e columns correspond to \u003cem\u003eN\u003c/em\u003e total stations. Each row has columns marked with True values to indicate which stations belong to that cluster after the cluster has been expanded. Depending on the GNSS station distribution geometry and \u003cem\u003eovercluster\u003c/em\u003e parameters, there will be varying redundancy in the membership of the clusters. We require that each station is represented at least once; while we expect that some stations will be represented multiple times due to overlaps, we proactively remove cases where the full cluster itself is redundant because that entire subnetwork membership is fully present in one or more other subnetwork clusters. To check for this redundancy, we take the sum of each column to verify if that GNSS station is captured one or more times, and then iterate through the matrix by temporarily removing a given row. As we remove a row / cluster, we recalculate the column totals; if removing that row causes us to lose representation of any station in our input dataset, then we pass on to the next row and cluster. However, if removing a cluster leaves all GNSS stations present in at least one other subnetwork, then we delete the cluster and update the \u003cem\u003eM\u003c/em\u003e by \u003cem\u003eN\u003c/em\u003e matrix encoding before continuing to recursively check the remaining rows and clusters for if they have any fully redundant coverage. Using our production parameters, a linear scan through the cluster encodings trims approximately 20\u0026ndash;30 percent of overall GNSS stations that are ultimately fed into GAMIT, with a concurrent drop in overall processing time. We expect the impact to processing time is disproportionately impacted by the extreme tail of smaller sized clusters as explained in Section \u003cspan class=\"InternalRef\"\u003e4\u003c/span\u003e, therefore, while we maintain the option to execute the \u003cem\u003eprune\u003c/em\u003e function as a simple linear scan, our current default method will monotonically sort the rows for the Boolean \u003cem\u003eovercluster\u003c/em\u003e matrix according to the row sum, and then preferentially eliminate redundant rows by removing the smallest redundant clusters first. The successive steps of \u003cem\u003eqmeans\u003c/em\u003e, \u003cem\u003eovercluster\u003c/em\u003e, and \u003cem\u003eprune\u003c/em\u003e are highlighted for a part of a global processing run in Fig. \u003cspan class=\"InternalRef\"\u003e1\u003c/span\u003e.\u003c/p\u003e\n\u003c/div\u003e\n\u003cdiv id=\"Sec6\" class=\"Section2\"\u003e\n \u003ch2\u003e2.4 Computational considerations and software availability\u003c/h2\u003e\n \u003cp\u003eThe bulk of the processing time is in the GAMIT runs, which is why we seek an efficient cluster structure\u0026mdash;so we can more efficiently process daily GNSS solutions on the time scale of years. The algorithms presented here are computationally efficient; the processing of a full day of ~\u0026thinsp;1,000 GNSS stations may take\u0026thinsp;~\u0026thinsp;20 hours in GAMIT, while our clustering processing and postprocessing adds only a few seconds per day, depending on station count. The bisecting \u003cem\u003eqmeans\u003c/em\u003e algorithm is stochastic, initiating random centroid locations as it searches for an optimum solution. However, to ensure that our results are reproducible, we typically fix the seed in the random number generator so that we have deterministic behavior when comparing runs with different parameter values. The \u003cem\u003eqmeans\u003c/em\u003e algorithm, as well as the \u003cem\u003eovercluster\u003c/em\u003e and \u003cem\u003eprune\u003c/em\u003e functions, are open-source functions integrated into Parallel.GAMIT, a BSD licensed open-source library maintained by the authors to parallelize GAMIT runs. Our library is available through GitHub (\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://github.com/demiangomez/Parallel.GAMIT\u003c/span\u003e\u003c/span\u003e\u003cspan type=\"Underline\" class=\"Underline\" name=\"Emphasis\"\u003e)\u003c/span\u003e, and is additionally installable from the pypi community repository under the package name \u003cem\u003epgamit\u003c/em\u003e when using the python \u0026lsquo;pip install\u0026rsquo; command. While GAMIT (\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttp://www-gpsg.mit.edu/gg/\u003c/span\u003e\u003c/span\u003e\u003cspan type=\"Underline\" class=\"Underline\" name=\"Emphasis\"\u003e)\u003c/span\u003e is available from its authors on request, distribution is limited to non-commercial research use, which also means that it is distributed separately from the \u003cem\u003epgamit\u003c/em\u003e package. The clustering functions we develop are freely available, and can be run both within the Parallel.GAMIT parallelization routines, with clustering functions (within the pgamit.cluster module) following the sklearn documentation and API style so that they can be used separately as needed.\u003c/p\u003e\n\u003c/div\u003e"},{"header":"3. Results","content":"\u003cp\u003eTo test \u003cem\u003eqmeans\u003c/em\u003e and evaluate the time-precision trade-off, we processed one year of data (2008) for a GNSS network in North America\u0026mdash;composed of ~\u0026thinsp;1,200 unique stations in M\u0026eacute;xico, the United States, and Canada\u0026mdash;and compared the run-times and wrms of four different cluster configurations: \u003cem\u003eqmax\u0026thinsp;=\u0026thinsp;20\u003c/em\u003e with 4, 6, and 10 overlap stations respectively, and \u003cem\u003eqmax\u0026thinsp;=\u0026thinsp;40\u003c/em\u003e with 10 overlap stations which is our baseline case for comparisons. Since \u003cem\u003eqmeans\u003c/em\u003e cannot guarantee clusters with a uniform number of stations, Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003ea shows the stations per cluster distribution for the test year using \u003cem\u003eqmax\u0026thinsp;=\u0026thinsp;40\u003c/em\u003e with 10 overlap stations, which results in a membership floor for clusters (at least 11 stations), as well as a ceiling (no more than 50 stations). Figure\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003eb shows that for \u003cem\u003eqmax\u0026thinsp;=\u0026thinsp;20\u003c/em\u003e with 10 overlap stations, the total cluster count per day rises as expected, and the clusters are distributed across a tighter range between 11 and 30 stations per cluster. For a graphical representation of the clusters using \u003cem\u003eqmax\u0026thinsp;=\u0026thinsp;40\u003c/em\u003e, see Fig \u003cspan refid=\"MOESM1\" class=\"InternalRef\"\u003eS1\u003c/span\u003e in the Supporting Information, and Fig S2 for cluster size and overlap station count histograms for a single day.\u003c/p\u003e\u003cp\u003e\u003c/p\u003e\u003cdiv id=\"Sec8\" class=\"Section2\"\u003e\u003ch2\u003e3.1 Clustering performance of the \u003cem\u003eqmeans\u003c/em\u003e algorithm\u003c/h2\u003e\u003cp\u003eIn terms of the number of clusters produced for a Parallel.GAMIT session, Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003ea shows that with 10 overlap stations, \u003cem\u003eqmax\u0026thinsp;=\u0026thinsp;20\u003c/em\u003e results in nearly double the number of clusters compared to the \u003cem\u003eqmax\u0026thinsp;=\u0026thinsp;40\u003c/em\u003e configuration. It is important to note that the number of clusters produced by \u003cem\u003eqmeans\u003c/em\u003e is invariant with respect to the \u003cem\u003eoverlap\u003c/em\u003e parameter specified, as the partitioning in the \u003cem\u003eqmeans\u003c/em\u003e process is independent of the \u003cem\u003eovercluster\u003c/em\u003e function, which does not add or remove clusters. Although the \u003cem\u003eprune\u003c/em\u003e function is impacted by structure and number of overlap stations when determining redundant clusters to remove, the output number of clusters is primarily determined by the \u003cem\u003eqmax\u003c/em\u003e parameter of \u003cem\u003eqmeans.\u003c/em\u003e Figure\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003eb shows the observed cumulative time difference between our reference case of \u003cem\u003eqmax\u0026thinsp;=\u0026thinsp;40\u003c/em\u003e with 10 overlap stations, and \u003cem\u003eqmax\u0026thinsp;=\u0026thinsp;20\u003c/em\u003e with 10, 6, and 4 overlap stations respectively. For the case where only \u003cem\u003eqmax\u003c/em\u003e differed and overlap stations were held at 10, our results showed that \u003cem\u003eqmax\u0026thinsp;=\u0026thinsp;40\u003c/em\u003e took\u0026thinsp;~\u0026thinsp;250 more hours to run relative to \u003cem\u003eqmax\u0026thinsp;=\u0026thinsp;20\u003c/em\u003e. Based on the total processing time for \u003cem\u003eqmax\u0026thinsp;=\u0026thinsp;40\u003c/em\u003e (~\u0026thinsp;8,500 hours), 250 hours represents only a 3% reduction of run-time. A more significant run-time reduction was observed for \u003cem\u003eqmax\u0026thinsp;=\u0026thinsp;20\u003c/em\u003e when varying overlap stations to 6 and 4, where the total run-time difference from our \u003cem\u003eqmax\u0026thinsp;=\u0026thinsp;40\u003c/em\u003e reference reached 2,000 and ~\u0026thinsp;2,600 hours (23% and 30% reduction), respectively.\u003c/p\u003e\u003cp\u003e\u003c/p\u003e\u003cp\u003e\u003cdiv class=\"gridtable\"\u003e\u003ctable float=\"Yes\" id=\"Tab1\" border=\"1\"\u003e\u003ccaption language=\"En\"\u003e\u003cdiv class=\"CaptionNumber\"\u003eTable 1\u003c/div\u003e\u003cdiv class=\"CaptionContent\"\u003e\u003cp\u003eAverage total stations to process after \u003cem\u003eovercluster\u003c/em\u003e and \u003cem\u003eprune\u003c/em\u003e.\u003c/p\u003e\u003c/div\u003e\u003c/caption\u003e\u003ccolgroup cols=\"4\"\u003e\u003cdiv align=\"left\" class=\"colspec\" colname=\"c1\" colnum=\"1\"\u003e\u003c/div\u003e\u003cdiv align=\"char\" char=\".\" class=\"colspec\" colname=\"c2\" colnum=\"2\"\u003e\u003c/div\u003e\u003cdiv align=\"char\" char=\".\" class=\"colspec\" colname=\"c3\" colnum=\"3\"\u003e\u003c/div\u003e\u003cdiv align=\"char\" char=\".\" class=\"colspec\" colname=\"c4\" colnum=\"4\"\u003e\u003c/div\u003e\u003cthead\u003e\u003ctr\u003e\u003cth align=\"left\" colname=\"c1\"\u003e\u0026nbsp;\u003c/th\u003e\u003cth align=\"left\" colname=\"c2\"\u003e\u003cp\u003eOverlap\u0026thinsp;=\u0026thinsp;4\u003c/p\u003e\u003c/th\u003e\u003cth align=\"left\" colname=\"c3\"\u003e\u003cp\u003eOverlap\u0026thinsp;=\u0026thinsp;6\u003c/p\u003e\u003c/th\u003e\u003cth align=\"left\" colname=\"c4\"\u003e\u003cp\u003eOverlap\u0026thinsp;=\u0026thinsp;10\u003c/p\u003e\u003c/th\u003e\u003c/tr\u003e\u003c/thead\u003e\u003ctbody\u003e\u003ctr\u003e\u003ctd align=\"left\" colname=\"c1\"\u003e\u003cp\u003e\u003cb\u003eqmax\u0026thinsp;=\u0026thinsp;40\u003c/b\u003e\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"char\" char=\".\" colname=\"c2\"\u003e\u003cp\u003e1,320\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"char\" char=\".\" colname=\"c3\"\u003e\u003cp\u003e1,408\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"char\" char=\".\" colname=\"c4\"\u003e\u003cp\u003e1,608\u003c/p\u003e\u003c/td\u003e\u003c/tr\u003e\u003ctr\u003e\u003ctd align=\"left\" colname=\"c1\"\u003e\u003cp\u003e\u003cb\u003eqmax\u0026thinsp;=\u0026thinsp;20\u003c/b\u003e\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"char\" char=\".\" colname=\"c2\"\u003e\u003cp\u003e1,480\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"char\" char=\".\" colname=\"c3\"\u003e\u003cp\u003e1,617\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"char\" char=\".\" colname=\"c4\"\u003e\u003cp\u003e2,010\u003c/p\u003e\u003c/td\u003e\u003c/tr\u003e\u003c/tbody\u003e\u003c/colgroup\u003e\u003c/table\u003e\u003c/div\u003e\u003c/p\u003e\u003cp\u003eTable\u0026nbsp;\u003cspan refid=\"Tab1\" class=\"InternalRef\"\u003e1\u003c/span\u003e provides a summary of the mean number of stations per cluster configuration for an average network size of ~\u0026thinsp;1,200 stations, and highlights the fundamental tradeoff between cluster size, total number of stations to process, and run-time impact. Using 10 overlap stations, the total number of input stations for GAMIT to process when \u003cem\u003eqmax\u0026thinsp;=\u0026thinsp;40\u003c/em\u003e is 2,010 versus 1,608 for \u003cem\u003eqmax\u0026thinsp;=\u0026thinsp;20\u003c/em\u003e; however, the total cumulative run-time is 3% less for \u003cem\u003eqmax\u0026thinsp;=\u0026thinsp;20\u003c/em\u003e, despite processing 402 additional stations, which is 25% more stations than the reference \u003cem\u003eqmax\u0026thinsp;=\u0026thinsp;40\u003c/em\u003e. This is expected given the polynomial behavior of run-time as a function of cluster size (discussed below in Section \u003cspan refid=\"Sec10\" class=\"InternalRef\"\u003e3.3\u003c/span\u003e), but optimization is complicated because of how overall station count varies as a function of both number of clusters (itself a function of \u003cem\u003eqmax\u003c/em\u003e) and number of overlap stations. As Table\u0026nbsp;\u003cspan refid=\"Tab1\" class=\"InternalRef\"\u003e1\u003c/span\u003e shows, since smaller clusters are more numerous, increasing the number of overlap stations has an outsized effect on overall number of cumulative stations to process\u0026ndash; going from 10 overlap stations to 6 overlap stations reduces cumulative station count by only 200 for \u003cem\u003eqmax\u0026thinsp;=\u0026thinsp;40\u003c/em\u003e, but by over 400 for \u003cem\u003eqmax\u0026thinsp;=\u0026thinsp;20\u003c/em\u003e since there are approximately double the clusters produced to add overlap stations to. Given that \u003cem\u003eqmax\u0026thinsp;=\u0026thinsp;40\u003c/em\u003e with 10 overlap stations and \u003cem\u003eqmax\u0026thinsp;=\u0026thinsp;20\u003c/em\u003e with 6 overlap stations have approximately the same cumulative number of stations to process (1,608 versus 1,617), the red line from Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003eb presents an interesting comparison on run-time, a difference driven entirely by the size distribution of the clusters.\u003c/p\u003e\u003c/div\u003e\u003cdiv id=\"Sec9\" class=\"Section2\"\u003e\u003ch2\u003e3.2 Solution quality impact\u003c/h2\u003e\u003cp\u003eCluster size and number of overlap stations is expected to imply a trade-off both between execution time, and solution quality. To assess this trade-off in solution quality, we analyzed the daily GNSS repeatability using \u003cem\u003eqmax\u0026thinsp;=\u0026thinsp;40\u003c/em\u003e with 10 overlap stations as a reference solution and compared it against solutions using \u003cem\u003eqmax\u0026thinsp;=\u0026thinsp;20\u003c/em\u003e with station overlaps of 4, 6, and 10 for the same test year and the same network. Daily repeatabilities are obtained by transforming consecutive daily GNSS solutions onto one another using six-parameter Helmert transformations (Bevis and Brown \u003cspan citationid=\"CR5\" class=\"CitationRef\"\u003e2014\u003c/span\u003e). After all solutions have been aligned, the differences between individual station coordinates are computed and a wrms value in north, east, and up provides an estimate of the scatter level of the solution for each individual site.\u003c/p\u003e\u003cp\u003e\u003c/p\u003e\u003cp\u003eFigure\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003ea-c show the cumulative empirical distribution function for the daily repeatability wrms scatter for the north, east, and up components for all stations in the analysis. A clear pattern emerges from these plots, which is most pronounced in the east component, where we note that the base reference case (\u003cem\u003eqmax\u0026thinsp;=\u0026thinsp;40\u003c/em\u003e with 10 overlap stations) outperforms the other solutions using \u003cem\u003eqmax\u0026thinsp;=\u0026thinsp;20\u003c/em\u003e across all overlap parameter values, with \u003cem\u003eoverlap\u0026thinsp;=\u0026thinsp;4\u003c/em\u003e producing the worst results, while solutions for \u003cem\u003eoverlap\u0026thinsp;=\u0026thinsp;6\u003c/em\u003e and \u003cem\u003eoverlap\u0026thinsp;=\u0026thinsp;10\u003c/em\u003e closely follow each other. Interestingly, the vertical component in Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003ec shows that solutions with \u003cem\u003eqmax\u0026thinsp;=\u0026thinsp;40\u003c/em\u003e or \u003cem\u003e20\u003c/em\u003e and \u003cem\u003eoverlap\u0026thinsp;=\u0026thinsp;10\u003c/em\u003e closely follow each other. Although the smallest subnetwork clusters with the least overlap are consistently outperformed by other configurations, as shown by the 0.1 and 0.3 mm shaded regions of scatter for the horizontal and vertical components in Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003ea-c, all tested cluster configurations fall within these tight regions and therefore can be considered statistically equal in terms of expected GNSS solution noise. Despite these results, the level of agreement between the wrms of these experiments cannot be always guaranteed, given that other effects such as network geometry or latitude location can affect the overall solution quality. Hence, whenever minimum scatter is desired, cluster size and number of overlap stations should be increased to achieve needed precision.\u003c/p\u003e\u003c/div\u003e\u003cdiv id=\"Sec10\" class=\"Section2\"\u003e\u003ch2\u003e3.3 Processing time as function of station count from historical data\u003c/h2\u003e\u003cp\u003eAs part of our parallel processing system for GAMIT, we save all of the relevant statistics for each run to a Postgres database, including execution time, overlap station count, and solution quality indicators. Since the deployment of Parallel.GAMIT at OSU in 2019, we have performed a total of ~\u0026thinsp;2.5M individual GAMIT runs for various projects. Using this Postgres database, we selected a subset of GAMIT runs for execution time statistics over an average GNSS project of ~\u0026thinsp;1,000 stations. We did not use the entire 2.5M runs in order to avoid double counting the same network distributions and duplicate runs, since many of the sessions are densifications or reruns using different models and orbits. After removing any \u0026lsquo;outlier\u0026rsquo; runs (those that took more than 80 minutes to complete) that could alter our statistics, the dataset used in this work is composed of ~\u0026thinsp;360k GAMIT executions. The removed long execution times account for less than 0.5% of the dataset, and were mostly related to external events not associated with the GNSS processing, such as computer network outages and other technological difficulties encountered during the processing.\u003c/p\u003e\u003cp\u003eFigure\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003ea shows the histogram of execution times for all the selected GNSS sessions from 1994 to 2022. For simplicity and consistency of the results, we performed all statistics and analysis using GPS-only solutions (rather than multi-constellation). The station clusters for these runs were created using an ad-hoc, non-optimized script (deployed from 2019 to 2024, prior to developing \u003cem\u003eqmeans\u003c/em\u003e) which creates station clusters by proximity and splits the resulting subnetworks if the number of stations in the processing session is above a threshold of ~\u0026thinsp;60. Once each cluster is determined, the clusters contain sets of unique stations to process and overlap stations from neighboring clusters, which are used to merge the solutions with GLOBK upon finishing all the processing jobs. Additionally, the clusters are also tied together using a \u0026lsquo;backbone\u0026rsquo; network of evenly distributed sites across the entire area of interest to provide a \u0026lsquo;rigid support\u0026rsquo; for the chainmail-like network of smaller clusters (see Appendix A for a description of the backbone network algorithm) and guarantee a robust combination.\u003c/p\u003e\u003cp\u003e\u003c/p\u003e\u003cp\u003eFigure\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003eb shows the cluster size histogram where we note two peaks, one around 15 and a second one near 40 stations. This distribution is produced due to the impossibility of splitting the stations into clusters with an exact number of stations. The data for execution time as a function of station count (Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003ec) can be fitted using a quadratic expression to obtain an estimate of the execution time; to account for other factors affecting the execution time at lower station counts, such as overhead time to read observation and metadata files, models, and others, we fit a second-degree polynomial rather than a pure quadratic. We performed this fit and found that the execution time \u003cem\u003et\u003c/em\u003e in minutes can be expressed as:\u003cdiv id=\"Equ1\" class=\"Equation\"\u003e\u003cdiv format=\"TEX\" class=\"mathdisplay\" id=\"FileID_Equ1\" name=\"EquationSource\"\u003e\n$$\\:t=6.04\\:-0.10\\cdot\\:s\\:+0.021\\cdot\\:{s}^{2}$$\u003c/div\u003e\u003cdiv class=\"EquationNumber\"\u003e3\u003c/div\u003e\u003c/div\u003e\u003c/p\u003e\u003cp\u003ewhere \u003cem\u003es\u003c/em\u003e is the number of stations in the processing session. Given the large number of data points used in the fit, the weighted root mean square error of our polynomial fit (95% confidence interval) was ~\u0026thinsp;2 minutes.\u003c/p\u003e\u003c/div\u003e"},{"header":"4. Discussion","content":"\u003cp\u003eAlthough determining the exact reason behind the solution quality improvement with increasing cluster size and overlap station count is outside the scope of this work, we surmise that clusters with more stations have, on average, a larger spatial aperture, and also larger number of double differences. Larger aperture and more double differences will improve the determination of the zenith tropospheric delay parameters which, in turn, improves both horizontal and vertical component scatter. Additionally, increasing the overlap station count increases the redundancy of the network, which helps mitigate the solution scatter when combining all clusters together.\u003c/p\u003e\u003cdiv id=\"Sec12\" class=\"Section2\"\u003e\u003ch2\u003e4.1 Execution time model\u003c/h2\u003e\u003cp\u003eDue to memory limits and also because execution time increases quadratically, GNSS sessions with \u003cem\u003eN\u003c/em\u003e\u0026thinsp;\u0026gt;\u0026thinsp;60\u0026ndash;80 stations need first to be partitioned into \u003cem\u003eC\u003c/em\u003e clusters of \u003cem\u003es\u003c/em\u003e stations each. As seen in Section \u003cspan refid=\"Sec7\" class=\"InternalRef\"\u003e3\u003c/span\u003e, our GAMIT runs occur on datasets where the cluster sizes are not uniform; the compute run-time for these production sessions can be approximated by modifying Eq.\u0026nbsp;\u003cspan refid=\"Equ1\" class=\"InternalRef\"\u003e3\u003c/span\u003e:\u003cdiv id=\"Equ2\" class=\"Equation\"\u003e\u003cdiv format=\"TEX\" class=\"mathdisplay\" id=\"FileID_Equ2\" name=\"EquationSource\"\u003e\n$$\\:T={\\sum\\:}_{i=1}^{C}\\left(6.04\\:-0.10\\cdot\\:({s}_{i}+e)\\:+0.021\\cdot\\:{({s}_{i}\\:+e)}^{2}\\right)$$\u003c/div\u003e\u003cdiv class=\"EquationNumber\"\u003e4\u003c/div\u003e\u003c/div\u003e\u003c/p\u003e\u003cp\u003eWhere \u003cem\u003ee\u003c/em\u003e is the \u003cem\u003eoverlap\u003c/em\u003e parameter value (constant for all clusters), \u003cem\u003eC\u003c/em\u003e is the actual number of clusters formed, and \u003cem\u003es\u003c/em\u003e\u003csub\u003e\u003cem\u003ei\u003c/em\u003e\u003c/sub\u003e is the cluster size of cluster at index \u003cem\u003ei\u003c/em\u003e. Eq.\u0026nbsp;(\u003cspan refid=\"Equ2\" class=\"InternalRef\"\u003e4\u003c/span\u003e) does not account for geometry considerations within GAMIT or other background processes that impact run-time; still, we find that the time estimate \u003cem\u003eT\u003c/em\u003e has only slight bias in predicting overall run-time when compared to metadata within our postgres database which records the empirical run-time.\u003c/p\u003e\u003cp\u003eEquation \u003cspan refid=\"Equ2\" class=\"InternalRef\"\u003e4\u003c/span\u003e is of little practical value by itself given that we already log the actual run-time of the GAMIT processing runs, and don\u0026rsquo;t need predictive estimation of individual subnetwork run-time for our operational processing. More useful is generalizing Eq.\u0026nbsp;\u003cspan refid=\"Equ2\" class=\"InternalRef\"\u003e4\u003c/span\u003e to visualize the behavior of the subnetworks and help illustrate parameter choices that can guide us in optimizing our compute pipeline. The ratio between parameters \u003cem\u003ee\u003c/em\u003e and \u003cem\u003es\u003c/em\u003e, \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:R=\\frac{e}{s}\\)\u003c/span\u003e\u003c/span\u003e, can be defined as the \u0026lsquo;redundancy\u0026rsquo; of a given Parallel.GAMIT session, where \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:R=0\\)\u003c/span\u003e\u003c/span\u003e indicates that there are no shared stations between clusters (this is an undesired condition) and \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:R=1\\)\u003c/span\u003e\u003c/span\u003e indicates that all stations are processed separately in two or more clusters. To idealize the compute run-time of a given GNSS network, the total number of clusters to be processed is \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:C=\\frac{N}{s}\\)\u003c/span\u003e\u003c/span\u003e, assuming \u003cem\u003eN\u003c/em\u003e total session stations are divisible by \u003cem\u003es\u003c/em\u003e, where stations per cluster \u003cem\u003es\u003c/em\u003e is now fixed at a single value for all clusters rather than varying per cluster. Assuming this idealized uniform distribution of clusters, the total execution time \u003cem\u003eT\u003c/em\u003e, including the redundancy, can be estimated as:\u003cdiv id=\"Equ3\" class=\"Equation\"\u003e\u003cdiv format=\"TEX\" class=\"mathdisplay\" id=\"FileID_Equ3\" name=\"EquationSource\"\u003e\n$$\\:T={\\sum\\:}_{i=1}^{C}{t}_{i}=C\\cdot\\:\\left(6.04\\:-0.10\\cdot\\:s\\cdot\\:(1+R)\\:+0.021\\cdot\\:{\\left[s\\cdot\\:(1+R)\\right]}^{2}\\right)$$\u003c/div\u003e\u003cdiv class=\"EquationNumber\"\u003e5\u003c/div\u003e\u003c/div\u003e\u003c/p\u003e\u003cp\u003eEquation (\u003cspan refid=\"Equ3\" class=\"InternalRef\"\u003e5\u003c/span\u003e) uses \u003cem\u003eR\u003c/em\u003e instead of \u003cem\u003ee\u003c/em\u003e, as redundancy is defined as the ratio between \u003cem\u003ee\u003c/em\u003e and \u003cem\u003es\u003c/em\u003e, which simplifies the interpretation of contour plot in Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003e.\u003c/p\u003e\u003cp\u003e\u003c/p\u003e\u003cp\u003eFigure\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003e shows a contour plot of \u003cem\u003eT\u003c/em\u003e (for \u003cem\u003eN\u003c/em\u003e\u0026thinsp;=\u0026thinsp;1,000) for Eq.\u0026nbsp;(\u003cspan refid=\"Equ3\" class=\"InternalRef\"\u003e5\u003c/span\u003e) as a function of \u003cem\u003es\u003c/em\u003e and \u003cem\u003eR\u003c/em\u003e, also including the contours of equal overlap stations, which vary with \u003cem\u003es\u003c/em\u003e for a constant \u003cem\u003eR\u003c/em\u003e. For clarity, we plotted a continuous contour field regardless of the value of \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:C=\\frac{N}{s}\\)\u003c/span\u003e\u003c/span\u003e, although the idealized execution time model is only valid when \u003cem\u003eC\u003c/em\u003e is an integer. We also show, as a black dashed contour, the minimum run-time for each value of \u003cem\u003eR\u003c/em\u003e and its corresponding cluster size. It should be noted that this line marks the limit of the cluster mediated computation efficiency: cluster sizes to the left of this dashed line will increase the run-time, decreasing efficiency rather than improving it. In other words, the execution time needed to complete a multi-GAMIT run session decreases concurrently with cluster size (i.e., as the number of clusters formed per session goes up) until the limit denoted by the intersection of the black dashed contour. Thus, if for example one uses 4 overlap stations, a session divided into uniform clusters of, say, 20 stations, takes ~\u0026thinsp;400 minutes less to finish than if it were divided into uniform clusters of 50 stations, although the number of clusters is larger for the former case. This tendency of decreasing time with lower cluster size, however, is less significant with increasing number of overlap stations. Figure\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003e shows that for increasing overlap stations, the overlap contours become increasingly parallel to the execution time contours, meaning that moving along overlap contours (and changing the cluster size) does not significantly change the execution time. We find that Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003e is mostly useful to examine the execution time behavior purely as a function of the overlap parameter, and help select an ideal cluster size to target.\u003c/p\u003e\u003cp\u003eAs we mentioned before, run-time is not the only consideration that a user has to account for when choosing a given partitioning. From Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003e it is clear that the choice of cluster configuration can have a considerable impact in the time needed to compute the solutions. This time difference will translate into a more or less significant \u0026lsquo;wall time\u0026rsquo;, i.e., the actual time needed to compute the solutions, depending on the number of compute nodes used to process the data. For instance, to process the totality of the GNSS stations in M\u0026eacute;xico, the United States, and Canada from 1994 to 2025 using clusters of size 20 with 6 overlap station clusters requires a total computation time of 8,223 days, or ~\u0026thinsp;33 days at 250 cores. If the same project is partitioned into clusters of size 50 with the same 6 station overlap, then the required run-time is 9,600 days, or ~\u0026thinsp;38 days at 250 cores, a wall time difference of 5 days. The wall time difference becomes more significant with a smaller compute cluster, and in the case of using only 100 cores, the wall-time difference between these two example runs is close to 14 days.\u003c/p\u003e\u003c/div\u003e\u003cdiv id=\"Sec13\" class=\"Section2\"\u003e\u003ch2\u003e4.2 Conclusions\u003c/h2\u003e\u003cp\u003eThis study presents a novel approach for efficiently partitioning large GNSS networks into subnetworks or clusters using a modified bisecting k-means algorithm, which we called \u003cem\u003eqmeans\u003c/em\u003e. The primary goal of this work was to develop a method for reducing the computational burden of processing large GNSS datasets while maintaining the precision of the results. Through extensive testing, we demonstrated that \u003cem\u003eqmeans\u003c/em\u003e provides a robust and scalable solution for dividing GNSS stations into clusters of manageable sizes, facilitating faster double-difference processing in GAMIT/GLOBK.\u003c/p\u003e\u003cp\u003eThe \u003cem\u003eqmeans\u003c/em\u003e algorithm itself offers several advantages. First, the total clustering processing time was found to add only a few seconds per day, a minimal overhead relative to the time savings achieved by efficiently dividing the network. Second, by eliminating the need to predefine the number of clusters (as in traditional k-means clustering), \u003cem\u003eqmeans\u003c/em\u003e provides a flexible and dynamic solution that can adapt to varying network sizes and spatial configurations. This feature is particularly beneficial when processing GNSS networks where station distributions are irregular, such as in sparsely populated regions or oceanic islands. Moreover, the hierarchical nature of \u003cem\u003eqmeans\u003c/em\u003e, combined with the post-processing steps (\u003cem\u003eovercluster\u003c/em\u003e and \u003cem\u003eprune\u003c/em\u003e), ensures that the final clusters are both computationally efficient and spatially coherent, with sufficient overlap for robust reference frame realization. As expected, the increase in the station count caused by adding additional overlap stations leads to a decrease in computational efficiency, as some stations are processed more than once. However, we have also demonstrated that this reduction in efficiency results in a corresponding increase in precision, which is a desired outcome when processing GNSS data.\u003c/p\u003e\u003cp\u003eOur results show that partitioning the network into clusters of approximately 20 to 40 stations provides an optimal balance between reducing computational time and preserving solution accuracy. Smaller clusters (20 stations) lead to faster processing times, with slightly increased solution scatter, while larger clusters (40 stations) offer a modest reduction in scatter at the cost of increased processing time. In our experiments, the \u003cem\u003eqmax\u0026thinsp;=\u0026thinsp;20\u003c/em\u003e with \u003cem\u003eoverlap\u0026thinsp;=\u0026thinsp;6\u003c/em\u003e configuration was shown to cut processing time by nearly 23% compared to the \u003cem\u003eqmax\u0026thinsp;=\u0026thinsp;40\u003c/em\u003e with \u003cem\u003eoverlap\u0026thinsp;=\u0026thinsp;10\u003c/em\u003e configuration despite processing a similar number of total stations within the session. Similarly, using \u003cem\u003eqmax\u0026thinsp;=\u0026thinsp;20\u003c/em\u003e with \u003cem\u003eoverlap\u0026thinsp;=\u0026thinsp;4\u003c/em\u003e reduced the processing time by 30%, despite generating twice the clusters\u0026ndash; both cases with a trade-off of a slight increase in solution scatter. These findings suggest that for large-scale GNSS projects, careful consideration of cluster size and overlap station count is critical for achieving efficient processing without compromising the quality of the result. By contrast, if one is seeking to obtain a fast solution for testing purposes, reducing the number of stations and overlap stations per cluster significantly reduces the processing time for GNSS solutions.\u003c/p\u003e\u003c/div\u003e"},{"header":"Declarations","content":"\u003cp\u003eThe research leading to these results received funding from the National Geodetic Survey (NGS), under Grant Agreement AWD-115866.\u003c/p\u003e\n\u003ch2\u003eConflicts of interest/Competing interests\u003c/h2\u003e\n\u003cp\u003eAll authors certify that they have no affiliations with or involvement in any organization or entity with any financial interest or non-financial interest in the subject matter or materials discussed in this manuscript. The authors have no financial or proprietary interests in any material discussed in this article.\u003c/p\u003e\n\u003ch2\u003eAuthor Contribution\u003c/h2\u003e\n\u003cp\u003eS.G. and D.G. contributed equally to the manuscript. S.G. designed and developed the qmeans, overcluster, and prune functions described in the manuscript, and led integration for adding them to the ParallelGAMIT library. D.G. envisioned experimental design, interpreted the results, and developed the run time predictive model. S.G. wrote the initial draft of the Methods and Discussion sections; D.G. wrote the initial draft of the Introduction and Results sections; both S.G. and D.G. reviewed, edited, and revised the entire manuscript. D.G. added supplementary material, and produced figures 2, 3, 4, 5 and 6. S.G. produced figure 1, and authored and tested the pgamit.cluster module used to create all figures except figure 5. S.G. and D.G. maintain the open source project Parallel.GAMIT (pgamit), and D.G. is the project founder of that project.\u003c/p\u003e\n\u003ch2\u003eAcknowledgement\u003c/h2\u003e\n\u003cp\u003eWe acknowledge Eric Kendrick for his role in maintaining the postgres database of GAMIT run statistics, which was queried by the authors to develop the predictive runtime model.\u003c/p\u003e\n\u003ch2\u003eData Availability\u003c/h2\u003e\n\u003cp\u003eData availability: Data is public and available through various websites, including the International GNSS Service data repository, through the EarthScope Facility Archive, and the NOAA Continuously Operating Reference Station (CORS) Network (NCN), managed by NOAA/National Geodetic Survey.\u003c/p\u003e"},{"header":"References","content":"\u003col\u003e\n \u003cli\u003eAltamimi Z, Rebischung P, Collilieux X, et al (2023) ITRF2020: an augmented reference frame refining the modeling of nonlinear station motions. J Geod 97:47. https://doi.org/10.1007/s00190-023-01738-w\u003c/li\u003e\n \u003cli\u003eAlves Costa SM, S\u0026aacute;nchez L, Pi\u0026ntilde;\u0026oacute;n D, et al (2022) Status of the SIRGAS reference frame: recent developments and new challenges. In: IAG International Symposium on Reference Frames for Applications in Geosciences. Springer Nature Switzerland Cham, pp 153\u0026ndash;165\u003c/li\u003e\n \u003cli\u003eAnkerst M, Breunig MM, Kriegel H-P, Sander J (1999) OPTICS: ordering points to identify the clustering structure. ACM SIGMOD Rec 28:49\u0026ndash;60. https://doi.org/10.1145/304181.304187\u003c/li\u003e\n \u003cli\u003eArthur D, Vassilvitskii S (2006) k-means++: The advantages of careful seeding. Stanford\u003c/li\u003e\n \u003cli\u003eBevis M, Brown A (2014) Trajectory models and reference frames for crustal motion geodesy. J Geod 88:283\u0026ndash;311. https://doi.org/10.1007/s00190-013-0685-5\u003c/li\u003e\n \u003cli\u003eCosta SMA, S\u0026aacute;nchez L, Pi\u0026ntilde;on D, et al (2023) Status of the SIRGAS reference frame: Recent developments and new challenges. In: International Association of Geodesy, Reference Frames Symposium REFAG2022\u003c/li\u003e\n \u003cli\u003eEster M, Kriegel H-P, Sander J, Xu X (1996) A density-based algorithm for discovering clusters in large spatial databases with noise. In: kdd. pp 226\u0026ndash;231\u003c/li\u003e\n \u003cli\u003eForgy EW (1965) Cluster analysis of multivariate data: efficiency versus interpretability of classifications. biometrics 21:768\u0026ndash;769\u003c/li\u003e\n \u003cli\u003eG\u0026oacute;mez DD, Bevis MG, Caccamise DJ, et al (2024) An empirical tool for predicting the presence or absence of coseismic displacements at GNSS stations. GPS Solut 28:214. https://doi.org/10.1007/s10291-024-01758-9\u003c/li\u003e\n \u003cli\u003eHerring TA, King RW, Floyd MA, McClusky SC (2018) Introduction to GAMIT/GLOBK\u003c/li\u003e\n \u003cli\u003eInternational Organization for Standardization (2020) Geographic information - Geodetic references - Part 1: International terrestrial reference system (ITRS) (ISO Standard No. 19161-1:2020(E))\u003c/li\u003e\n \u003cli\u003eLloyd S (1982) Least squares quantization in PCM. IEEE Trans Inf Theory 28:129\u0026ndash;137\u003c/li\u003e\n \u003cli\u003eNg A, Jordan M, Weiss Y (2001) On spectral clustering: Analysis and an algorithm. Adv Neural Inf Process Syst 14:\u003c/li\u003e\n \u003cli\u003eOmohundro SM (1989) Five balltree construction algorithms\u003c/li\u003e\n \u003cli\u003ePedregosa F, Varoquaux G, Gramfort A, et al (2011) Scikit-learn: Machine learning in Python. J Mach Learn Res 12:2825\u0026ndash;2830\u003c/li\u003e\n \u003cli\u003eSchubert E, Sander J, Ester M, et al (2017) DBSCAN Revisited, Revisited: Why and How You Should (Still) Use DBSCAN. ACM Trans Database Syst 42:1\u0026ndash;21. https://doi.org/10.1145/3068335\u003c/li\u003e\n \u003cli\u003eSteinbach M (2000) A Comparison of Document Clustering Techniques Michael Steinbach, George Karypis, and Vipin Kumar\u003c/li\u003e\n\u003c/ol\u003e"}],"fulltextSource":"","fullText":"","funders":[],"hasAdminPriorityOnWorkflow":false,"hasManuscriptDocX":true,"hasOptedInToPreprint":true,"hasPassedJournalQc":"","hasAnyPriority":false,"hideJournal":false,"highlight":"","institution":"","isAcceptedByJournal":true,"isAuthorSuppliedPdf":false,"isDeskRejected":"","isHiddenFromSearch":false,"isInQc":false,"isInWorkflow":false,"isPdf":false,"isPdfUpToDate":true,"isWithdrawnOrRetracted":false,"journal":{"display":true,"email":"[email protected]","identity":"gps-solutions","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":false,"externalIdentity":"gpss","sideBox":"Learn more about [GPS Solutions](http://link.springer.com/journal/10291)","snPcode":"10291","submissionUrl":"https://submission.nature.com/new-submission/10291/3","title":"GPS Solutions","twitterHandle":"","acdcEnabled":true,"dfaEnabled":true,"editorialSystem":"em","reportingPortfolio":"Springer Hybrid","inReviewEnabled":true,"inReviewRevisionsEnabled":false},"keywords":"GNSS processing, double-difference processing, clustering, GNSS network, GNSS solution scatter, segmentation","lastPublishedDoi":"10.21203/rs.3.rs-7096364/v1","lastPublishedDoiUrl":"https://doi.org/10.21203/rs.3.rs-7096364/v1","license":{"name":"CC BY 4.0","url":"https://creativecommons.org/licenses/by/4.0/"},"manuscriptAbstract":"\u003cp\u003eThe rapid growth of GNSS networks poses significant challenges for efficiently processing large datasets using double-difference techniques. In this study, we introduce a novel clustering algorithm, \u003cem\u003eqmeans\u003c/em\u003e, which is based on bisecting k-means, to partition GNSS networks into smaller, manageable subnetworks or clusters for double-difference processing. We explore the trade-offs between cluster size, computational cost, and solution quality using a comprehensive dataset of approximately 1,200 stations distributed across México, the United States, and Canada. Our results demonstrate that partitioning the network into clusters of 20 to 30 stations with 6 overlap stations between clusters can reduce processing time by ~20%, while larger clusters of 40-50 stations with 10 overlap stations slightly improve solution precision. We show that the number of shared stations between clusters impacts both the computational efficiency and the precision of the final solution, with higher counts leading to better precision but also increased processing time. The \u003cem\u003eqmeans\u003c/em\u003e algorithm is integrated into the open-source Parallel.GAMIT software, offering a scalable, flexible solution that can be applied to large GNSS networks. Our work sets a foundation for selecting optimal subnetwork sizes based on specific needs of a GNSS processing project, enabling faster processing without significantly sacrificing solution quality.\u003c/p\u003e","manuscriptTitle":"Efficient clustering of GNSS stations for processing using double differences","msid":"","msnumber":"","nonDraftVersions":[{"code":1,"date":"2025-07-18 09:29:26","doi":"10.21203/rs.3.rs-7096364/v1","editorialEvents":[{"type":"communityComments","content":0},{"type":"decision","content":"Revision requested","date":"2025-11-11T02:34:36+00:00","index":"","fulltext":""},{"type":"editorInvitedReview","content":"","date":"2025-11-03T23:05:42+00:00","index":"hide","fulltext":""},{"type":"reviewerAgreed","content":"26624931078527562094406038102128810505","date":"2025-09-15T16:44:58+00:00","index":"hide","fulltext":""},{"type":"reviewerAgreed","content":"179510822520525232727396998312969688924","date":"2025-09-10T22:12:00+00:00","index":"hide","fulltext":""},{"type":"reviewersInvited","content":"","date":"2025-07-13T06:46:42+00:00","index":"","fulltext":""},{"type":"editorAssigned","content":"","date":"2025-07-13T06:41:22+00:00","index":"","fulltext":""},{"type":"checksComplete","content":"","date":"2025-07-11T05:20:06+00:00","index":"","fulltext":""},{"type":"submitted","content":"GPS Solutions","date":"2025-07-10T22:47:15+00:00","index":"","fulltext":""}],"status":"published","journal":{"display":true,"email":"[email protected]","identity":"gps-solutions","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":false,"externalIdentity":"gpss","sideBox":"Learn more about [GPS Solutions](http://link.springer.com/journal/10291)","snPcode":"10291","submissionUrl":"https://submission.nature.com/new-submission/10291/3","title":"GPS Solutions","twitterHandle":"","acdcEnabled":true,"dfaEnabled":true,"editorialSystem":"em","reportingPortfolio":"Springer Hybrid","inReviewEnabled":true,"inReviewRevisionsEnabled":false}}],"origin":"","ownerIdentity":"d63cbff9-cba7-4b5c-8148-bba1e07ce558","owner":[],"postedDate":"July 18th, 2025","published":true,"recentEditorialEvents":[],"rejectedJournal":[],"revision":"","amendment":"","status":"published-in-journal","subjectAreas":[],"tags":[],"updatedAt":"2026-01-26T16:08:49+00:00","versionOfRecord":{"articleIdentity":"rs-7096364","link":"https://doi.org/10.1007/s10291-025-02020-6","journal":{"identity":"gps-solutions","isVorOnly":false,"title":"GPS Solutions"},"publishedOn":"2026-01-21 15:58:30","publishedOnDateReadable":"January 21st, 2026"},"versionCreatedAt":"2025-07-18 09:29:26","video":"","vorDoi":"10.1007/s10291-025-02020-6","vorDoiUrl":"https://doi.org/10.1007/s10291-025-02020-6","workflowStages":[]},"version":"v1","identity":"rs-7096364","journalConfig":"researchsquare"},"__N_SSP":true},"page":"/article/[identity]/[[...version]]","query":{"redirect":"/article/rs-7096364","identity":"rs-7096364","version":["v1"]},"buildId":"8U1c8b4HqxoKbykW_rLl7","isFallback":false,"isExperimentalCompile":false,"dynamicIds":[84888],"gssp":true,"scriptLoader":[]}