PyReX: A Recursion Based Crossover Detection Algorithm in Python for Along‐Track Geophysical Measurements

Abstract

A crossover point is the location of intersection of any two ground tracks charted by multiple platforms (ships, satellite radar and laser altimeters etc.). Detection of crossovers is of prime importance to estimate the discrepancies in the geophysical measurements at the crossover points. Usual approach of crossover detection considers consecutive data points in tracks as segments and checks for intersections between all combinations of these segments. We present a Recursion based crossover detection algorithm in Python (PyReX) for rapid detection of crossovers by avoiding redundant intersection checks. We test the performance of this algorithm using along‐track sea surface height measurements from satellite altimeters. We observe that the time taken for flagging a crossover between pair of tracks with N segments each varies as $$\log \,N$$ vis‐a‐vis the $${N}^{2}$$ dependency associated with the traditional methods. We further demonstrate that PyReX significantly improves the computation speed for high frequency along‐track measurements from satellite altimeters and ship‐borne gravity data compared to existing algorithms. PyReX is a flexible, open‐source code that could be easily customized for variety of applications involving large‐scale track‐line data sets.

Introduction

Seamless combination of measurements of geophysical parameters along intersecting tracks from multiple platforms including ocean research vessels, air and space borne instruments or terrestrial observations is often challenging due to discrepancies at the intersection points. Large errors at crossovers often resemble as artificial features and question the geophysical interpretations made from the data sets. The source of such crossover errors (COE) could be due to various reasons (Wessel, 2010) such as (a) uncertainties in the coordinates leading to erroneous locations at point of intersections, (b) a non-homogenous phenomena resulting in a non-zero COE, (c) dynamic process measured from different epochs, (d) improper calibration at instrument level, etc. Rapid detection of COE at intersecting tracks are extremely important for global and regional data analysis involving repeated measurements in space or time (Bansal et al., 2005; Sandwell et al., 2014; Sandwell & Smith, 1997). Crossover analysis has evolved as a standard technique for evaluating the performance of satellite altimeters (Dorandeu et al., 2004; Krishna et al., 2023; Prandi et al., 2015; Tai & Fu, 1986) and generation of global products such as mean sea surface, marine gravity and sea level maps (Ballarotta et al., 2023; Pujol et al., 2016; Schaeffer et al., 2012; Yuan et al., 2020). Crossover adjustments are also crucial for analysis of ship-borne, terrestrial and planetary geophysical measurements (Bansal et al., 2005; Barker et al., 2016; Bertone et al., 2021; Desprats et al., 2024; Neumann et al., 2001; Sreejith & Krishna, 2013).

Existing algorithms extract crossovers by considering consecutive data points from a pair of tracks as segments and check for possible intersections at pair-wise combinations of these segments over all tracks (Greene et al., 2017; Hsu, 1995; Li et al., 2022; Wessel, 1989, 2010). For a pair of tracks with N segments, the number of computations and time required to check for an intersection increases as $${N}^{2}$$ (Hsu, 1995). Thus, it becomes a time consuming task to check for intersections in applications dealing with large volumes of along-track data sets.

We present a Recursion based crossover detection algorithm implemented in Python (PyReX), where the time taken for flagging a crossover varies as $$\log \,N$$. We test the performance of this crossover detection method for both satellite altimeter and ship track data sets. PyReX uses Numba Just-In-Time (JIT) (Lam et al., 2015, 2024) and multiprocessing for rapid detection of crossovers in high frequency along-track data sets. We validate and demonstrate the improvement in computation speed of PyReX by comparing with other existing algorithms.

Crossover Detection

In PyReX, the crossover detection is made efficient by excluding parts of tracks which do not intersect by gradually splitting of tracks into smaller sub-tracks, first based on trend changes in latitude or longitude (Figure S1 in Supporting Information S1) and later recursively (Figure 1). Brief descriptions of major steps involved in the PyReX workflow (Vikram et al., 2024) are provided below.

[IMAGE OMITTED. SEE PDF]

Rectangle Test

For a pair of tracks of any arbitrary lengths, the two rectangular regions bounded by the endpoints of each track are checked for overlap (Figures 2a and 2b). This Rectangle Test (RT) flags for a possibility of intersection while a failure in the same rejects the chance of any intersection between the two tracks (Li et al., 2022; Wessel, 1989). A simple RT check is done initially to eliminate tracks that could not possibly intersect. The remaining tracks are split into sub-tracks at the location of trend change in latitude or longitude (Figure S1 in Supporting Information S1) and where longitude crosses 180°. These sub-tracks will be split further multiple times and RT will be implemented on them. It is efficient to have the end points of any sub-track as maxima and minima of latitude or longitude at any instant to easily find the rectangular region.

[IMAGE OMITTED. SEE PDF]

Recursion Method

After splitting the tracks based on trend change in latitude or longitude (Figure S1 in Supporting Information S1) and where longitude crosses 180°, RT is done on the sub-track combination pairs. The pairs flagged for a possible intersection after qualifying RT are further split into two sub-tracks (Figure 2c). RT is executed on all combinations of these sub-tracks and the combinations with successful RT are split further into smaller sub-tracks and tested again (Figures 2d and 2e). This procedure of splitting and testing is repeated recursively until the sub-tracks are shrunk down to individual segments between two points (Figure 2f). RT on these individual segments flags the possible segment pairs that may intersect. The need for testing intersections on all combinations of segments (Hsu, 1995; Schwarz, 2024) is thus eliminated.

For the recursive algorithm, the number of segments (N) in a track can be related with the depth of recursion (d) and number of splits or recursion type (r) to reach segments as 1 $$\frac{N}{{\mathrm{r}}^{\mathrm{d}}}\simeq 1\Rightarrow d=\frac{\log \,N}{\log \,r}$$

The time taken by the recursive algorithm to reach segment level is proportional to the depth of recursion. It is also proportional to the square of number of sub-tracks, as all the combinations of sub-tracks need to be tested. From Equation 1 it can be deduced as 2 $$t\mathit{\propto }\frac{{r}^{2}}{\log \,r}\,\log \,N$$

Thus, making the time taken proportional to $$\log \,N$$. However, if the segment combinations of non-overlapping sub-tracks are not filtered out, the number of checks (and time) required varies as $${N}^{2}$$ for track with N segments (Hsu, 1995).

Intersection Test

The segment pairs flagged by the recursion method are tested for intersections and the points of intersection or the crossover locations are calculated. The segments can be assumed as straight lines in case of smaller distances. Intersection is calculated using line equations and examined if the point of intersection lies internally that is, within both the segments (Figure 2f). However, the segments can also be defined as great circle arcs and the point of intersection is calculated by the procedure explained in Appendix A1.

Once the crossovers are detected, Haversine formula (Robusto, 1957) is used for finding the distances of crossover points along the track, and a piece-wise linear function is used to interpolate the geophysical parameters.

Implementation of PyReX

PyReX is a crossover detection tool primarily focused on the extraction of crossovers from along-track geophysical measurements (elevation, sea surface topography, sea level anomaly, gravity anomalies etc.) in an efficient manner. The configuration script (config_cross.py) is the starting point for the user working with PyReX. It is divided into multiple sections each denoted by a calc_id, and a set of initial parameters within each section to perform a specific calculation. This config_cross.py should be updated by adding new sections (with new calc_id) or editing the already existing sections based on the user's requirements.

The essence of PyReX is in the detection script (detection.py) which takes the list of track files containing along-track measurements. The user may specify the list of input tracks, a preferred intersection method (line or arc intersection), type of crossover (intersections within a track are “internal” while intersections between two tracks are “external”) as initial parameters within the calc_id section corresponding to detection in config_cross.py. Execution of detection.py with its corresponding calc_id reports the list of track pairs with intersection locations along with interpolated values as an output detection_file. The user may use this detection_file for further crossover analysis like inspecting COE statistics and crossover adjustments.

In addition, PyReX offers a crossover correction tool (Appendix A2) and a stacking tool for averaging repeated along-track measurements (Appendix A3) as auxiliary codes.

Performance Characteristics of PyReX Using Satellite Altimetric Crossovers

We use the 40 Hz along-track sea surface height (SSH) data measured by SARAL/AltiKa (AVISO+, 2020) to evaluate the performance of PyReX.

Investigation on a Pair of Satellite Tracks

We select a pair of tracks (Track #5 and #106 from cycle 15; Figure S2 in Supporting Information S1) having one crossover point and consider same number of data points on either sides of the crossover point to evaluate the dependency of length of the track on time taken to detect the crossover. Similarly, to analyze the role of along-track sampling rate on detection time, we consider every nth point from the original 40 Hz data to effectively reduce it to a (40/n) Hz data. For example, taking every second (2nd) and fourth (4th) point from 40 Hz data reduces the tracks to 20 and 10 Hz data respectively. Time taken for extraction of this crossover point by detection function (inside detection.py) is computed by averaging 1,000 runs using repeat function ($$repeat=1,000$$) from Python timeit module. These variations of time taken with track sizes (length of track and along-track sampling rates) can be seen as timing curves (Figure 3).

[IMAGE OMITTED. SEE PDF]

Tracks can be recursively split into any number of sub-tracks rather than just two (Figures 2c–2e). We change the recursion type (i.e., value of “r” in Equations 1 and 2) to assess its impact on the detection time. The time should vary logarithmically with the number of segments (N) in a track (from Equation 2). In order to capture any linear and offset components due to other steps, we fit the timing curves (Figure 3) with the function, t (N), given as 3 $$t(N)=a\,\log \,N+bN+c$$ where a, b, and c are the coefficients for logarithmic, linear and offset components respectively.

Timing curves for the detection of a single crossover point show a linear behavior for different recursion types (Figures 3a and 3b). This can be attributed to the time taken for determining the locations (indices) of trend change in latitude or longitude, and 180° longitude crossing. Moreover, when this step is not included in the timing calculation, we observe that the timing curves behave logarithmically with a significant improvement in the time taken, as confirmed by Figures 3c and 3d.

Splitting the tracks into two sub-tracks that is, recursion of type 2 is the fastest and consistent with Equation 2 (Figures 3c and 3d). This is because recursion adds overhead (contribution from $$\frac{{r}^{2}}{\log \,r}$$ factor in Equation 2) when tracks are split into larger number of sub-tracks.

Performance on a Global-Scale Analysis

We test the total time taken by PyReX detection.py script (single run without multiprocessing) from all the tracks (1,002 tracks) in cycle 15 using the bash utility/usr/bin/time (Figure 4). We change the sampling rate but the number of crossover points remain almost the same. The detection.py script calculates the locations (indices) of trend changes before testing for detection, thus avoiding the need to do this redundantly for every combination of tracks. Similarly, maxima and minima of latitude and longitudes for the preliminary RT, quality flags and the distances of each point along the track are also pre-calculated. These distances are later used to calculate the distance of the detected crossover point along the track, which is a useful information in addition to its location. Timing curves for the crossover detection from a single cycle show both linear and logarithmic behavior (Figure 4). The linear component can be reduced by pre-calculating information, whenever possible. In addition, using the Numba JIT compilation along with these precalculations, the overall time reduces even further. It should be noted that the efficiency of crossover detection is not limited by the recursion method proposed in the current study.

[IMAGE OMITTED. SEE PDF]

Comparison With the Existing Methods

The Schwarz method (Schwarz, 2024) of fast and robust curve intersection in MATLAB is the primary tool to extract crossovers in Antarctic Mapping Tools (AMT) (Greene et al., 2017) and PyXover (Bertone et al., 2020). Moreover, the Improved Rapid rejection and Straddle Test (IRST) algorithm (Li et al., 2022) is developed by modifying “tile-off” in AMT (Greene et al., 2017). PyXover (Bertone et al., 2020) subsamples the tracks at a fixed ratio (Bertone et al., 2021; Desprats et al., 2024) to get a rough/potential crossover using the Schwarz method (Schwarz, 2024) and then iteratively refines the search locally to get the exact crossover. Therefore, we utilize these crossover detection methods viz., the Schwarz method (Schwarz, 2024) and IRST (Li et al., 2022) along with the x2sys package (Wessel, 2010) from GMT (Wessel et al., 2019, 2024) to investigate the performance of PyReX (without numba and without multiprocessing). This is demonstrated by a preliminary comparision of the time taken (50 runs using the bash utility/usr/bin/time) to detect crossovers for satellite altimetric tracks subsampled at different along-track sampling rates over a local study region (40–20°S and 60–80°E).

The GMT x2sys package, written in “C”, takes a set of tracks as input. In contrast, the detection functions from Schwarz, 2024 and IRST (modified to Mercator projection and searching in ±3° of rough position), both written in MATLAB, can only accept individual tracks and had to be stitched together with matlab engine in a Python script to loop over the tracks in a data set. We observe that PyReX performs better with increasing data sizes (Figure 5) even though different implementations have different overheads and offsets. The number of crossovers (Table T1 in Supporting Information S1) along with their COE statistics (Figure S3 in Supporting Information S1) extracted from PyReX for different along-track sampling rates remains exactly the same when compared with other methods. In addition, we observe that the locations of crossovers are identical (Tables T2 and T3 in Supporting Information S1) since there is no ambiguity in the solution for intersection of two line segments.

[IMAGE OMITTED. SEE PDF]

Validation of PyReX Using Ship-Track Crossovers

Ship-borne free-air gravity measurements are known to have COE up to several tens of mGal (Talwani, 1971; Wessel & Watts, 1988). We select marine gravity measurements from 55 ship-borne tracks located in the Indian Ocean region (Table T4 and Figure S4 in Supporting Information S1) available from NCEI's Marine Trackline Geophysical Data () to compare the crossover statistics of PyReX with that obtained from GMT x2sys package (Wessel, 2010; Wessel et al., 2019, 2024). The outliers in ship track data could be due to instrumental errors, navigational uncertainties, and errors caused by incomplete correction of gravity anomalies (Wessel & Watts, 1988). These outliers are removed using statistical reporting of minimum, maximum, mean and standard deviation. Further, ship tracks are categorized into two groups to demonstrate the crossover detection capabilities of PyReX for low (Group A) and high frequency (Group B) along-track measurements. Group A consists of 53 tracks while Group B contains two extra tracks (index 13 and index 17) with a higher sampling rates of >50 samples/km in addition to the 53 tracks from Group A (Table T4 and Figure S4 in Supporting Information S1).

The time taken to detect crossovers using PyReX is measured without any cache to compare with GMT. We observe that for Group A, the time taken for detection of crossovers using PyReX and GMT were 39.6 and 27.5 s, respectively (Table T5 in Supporting Information S1). However, PyReX takes only 17.8 s for the same data using cache from previous runs. The number of crossovers is the same and the locations are identical with similar crossover statistics for PyReX and GMT (Figure 6 and Table T5 in Supporting Information S1).

[IMAGE OMITTED. SEE PDF]

In the case of Group B, PyReX extracts the crossovers within 42.6 s without cache while it takes 282.3 s for GMT (Table T5 in Supporting Information S1). In addition, the time taken to detect crossovers is just 20.6 s for PyReX using cache from previous runs for the same data. We observe that the time taken to detect crossovers using PyReX is lower by an order of magnitude compared to GMT in this case while dealing with higher sampling rate measurements. PyReX detects an extra crossover (Figure S4 in Supporting Information S1) which could be due to an edge case as described in Appendix A4. After ignoring this extra crossover, we observe that the crossover locations are identical with similar crossover statistics for PyReX and GMT even in the case of high sampling rate along-track measurements (Figure 6 and Table T5 in Supporting Information S1). This demonstrates the capability of PyReX, and its efficiency has potential implications while dealing with large data sets.

Conclusions

We present a Recursion based crossover detection algorithm in Python (PyReX) for analysis of satellite altimetric and ship track data over globe. The algorithm is designed to efficiently detect crossovers using recursion based flagging to avoid redundant intersection checks. We test this algorithm by measuring time taken to detect SSH crossovers for different satellite altimetric track sizes. We observe that it varies as $$\log \,N$$ and the recursion based flagging is not the bottle-neck. As expected splitting the tracks into two sub-tracks that is, recursion of type 2 is the fastest. We use ship track data to validate the crossover locations and observe that PyReX is efficient with increasing data sizes which has implications while dealing with large data.

Appendix

Appendix - A

A1

Intersection Test for Arc Segments

A vector joining a point with latitude (λ) and longitude (ϕ) from origin is expressed as $$P=\left[\begin{array}{@{}c@{}}x\\ y\\ z\end{array}\right]=\left[\begin{array}{@{}c@{}}\cos \,\lambda \,\cos \,\phi \\ \cos \,\lambda \,\sin \,\phi \\ \sin \,\lambda \end{array}\right]$$

For each segment along great circle arc, the cross-product between the vectors of the end points is $${N}_{1}={P}_{1a}\times {P}_{1b}$$ $${N}_{2}={P}_{2a}\times {P}_{2b}$$ where 1 and 2 are first and second arc-segments, while a and b stands for the start and endpoints of segments respectively. The direction of the intersection point is given by $$L={N}_{1}\times {N}_{2}$$ $${I}_{1}=\frac{L}{\Vert L\Vert }$$ $${I}_{2}=-{I}_{1}$$ where $${I}_{1}$$ and $${I}_{2}$$ are the two antipodal solutions of intersection points. Only one of the two antipodal solutions lies internally (Figure 2f), that is, on both the segments when they intersect. This solution is taken as the point of intersection that is, crossover point.

A2

Corrections

The difference in the measurements (Z) at crossovers, denoted by COE, are defined as $$COE={Z}_{self}-{Z}_{others}$$ where “self” is a track that is being corrected and “others” are the tracks used for correction, respectively. The main goal of this crossover correction is to reduce the COE. For each track corrections at each crossover are calculated from COE and are propagated for the entire track by curve fitting the corrections at crossovers along the distance on track using scipy.optimize.curve fit function (Virtanen et al., 2020) and subtracting the fit values. Correction at each crossover is calculated by $${X}_{correction}=\left({Z}_{self}-{Z}_{other}\right)\ast weight$$

Different values for weight can be set. The value of 1 indicates sticking “self” tracks to “others” tracks. In case of same set of tracks for both “self” and “others,” a weight of 0.5 can be used to correct for internal consistency.

The function used for curve fitting is (Yuan et al., 2020) $\begin{align*}F(d)={a}_{-1}& \text{For}\,\text{fit}\,\text{order}-1\\ ={a}_{-1}+{a}_{0}d& \text{For}\,\text{fit}\,\text{order}\,0\\ ={a}_{-1}+{a}_{0}d+\sum\limits _{i=1}^{n}\left({c}_{i}\,\cos \left(i\frac{2\pi d}{{F}_{d}}\right)+{s}_{i}\,\sin \left(i\frac{2\pi d}{{F}_{d}}\right)\right)& \text{For}\,\text{any}\,\text{fit}\,\text{order}\,\mathrm{n}\end{align*}$ where “d” is distance of a crossover along the track and “$${F}_{d}$$” is the fundamental distance. Fundamental distance can be any length greater than the length of the longest track in the data set to capture all the components. In the case of fitting a smaller track, the distance and order are scaled according to the length of the track.

The detection_file generated by detection.py is taken as input for crossover adjustments (correction.py) by adding/editing the correction calc_id section in config_cross.py. The fundamental distance should be set greater than or equal to the length of the longest track. When distance between two consecutive crossovers along a track is large, a threshold can be set to use simple linear interpolation instead of sinusoidal fit. The user should select tracks to be considered for self and other based on the track names in detection_file. After correction, a corrected detection_file will be generated based on the corrected tracks.

A3

Stacking

Along-track data contains information regarding location in geographic (longitude, latitude) or Cartesian (x, y) coordinates and a geophysical observable (z). For satellites in Exact Repeat Mode (ERM) orbit, altimetric measurement at every location along the track is repeated after a fixed time duration (10 days for Jason series and 35 days for ERS/Envisat/SARAL). Therefore, for stacking ERM measurements, a base track is generated by considering all the data corresponding to the same track number available from all cycles.

Observations lying within a given radial distance from each location on the base track are averaged to get a stacked profile. During averaging, outliers beyond 3σ are iteratively excluded until no new data is removed. The base track for a stacked profile can either be the track with highest number of observations or a user specified track of choice. The users must specify number of tracks in a cycle (e.g., SARAL has 1,002 tracks per cycle), the radial distance and the base track to be considered for stacking under the calc_id section of the config_cross.py.

A4

Edge Cases

PyReX and GMT may have some differences in the number of crossovers detected depending on how the edge cases are dealt. These edge cases can arise when segments are perfectly horizontal or vertical, segments coinciding with each other and segments intersecting exactly at the end points. The strategy of dealing with such edge cases varies between PyReX and GMT on case-to-case basis.

In the case of Group B, PyReX extracts an extra crossover between index13 and dme10 (Figure S3 in Supporting Information S1) compared to GMT. The high sampling rate measurements from index13 with insufficient precision in the locations of along-track data has an intersection with dme10 exactly on the endpoint of segment. PyReX considers this as an intersection consistent with other neighboring scenarios while GMT skips it.

Acknowledgments

This work has been carried out as part of GAP program of SAC. The satellite altimetric data products are processed on behalf of CNES SALP project and distributed by AVISO+ (). The ship track data can be accessed from NCEI's Marine Trackline Geophysical Data (). We thank the editor, Cathleen Jones, and two anonymous reviewers for their constructive suggestions. We thank Nilesh Desai and Rashmi Sharma from SAC for their support.

Data Availability Statement

PyReX code and the test data sets are available for downloading at: .

References

AVISO+. (2020). Along‐track level‐2+ (L2P) SLA; SALP‐MU‐P‐EA‐23150‐CLS; Issue 2.0. CNES.

Ballarotta, M., Ubelmann, C., Veillard, P., Prandi, P., Etienne, H., Mulet, S., et al. (2023). Improved global sea surface height and current maps from remote sensing and in situ observations. Earth System Science Data, 15(1), 295–315. [DOI: https://dx.doi.org/10.5194/essd-15-295-2023]