How to use the x2sys_init and x2sys_cross for intersections calculation?

Sorry I’m a beginner of python and GMT, even though I read the documentation from x2sys, and I also found some code from Github issues page, I still could not figure out how to implement this in my python code.

So could anyone please give some example code on how to use x2sys_init and x2sys_cross in Pygmt to find the intersections in a track dataset?

Thank you very much!

Hi @skgiter, welcome to the forum! I will have to say that x2sys is not an easy module to use, but will try to help out as much as I can. It will likely involve looking up the GMT documentation at and

To get you started, this is what I do for my work on ICESat-2 with x, y, h (height/elevation) and t (time) columns. The code is adapted from (scroll down to the ’ Crossover Track Analysis’ section).

  1. You will need to have a Format Definition File (*.fmt) in place first. Mine is at and looks like the below (!!important: put this file in X2SYS/TAGNAME where TAGNAME is the tag name you use, e.g. ICESAT2 in my case):
# Define file for X2SYS processing of ASCII tsv files
# This file applies to plain 4-column ASCII files.
#ASCII		# The input file is ASCII
#SKIP 1		# The number of header records to skip
#name	intype	NaN-proxy?	NaN-proxy	scale	offset	oformat
x	a	N		0		1	0	%g
y	a	N		0		1	0	%g
h	a	N		0		1	0	%g
t	a	N		0		1	0	%g
  1. Next, you will need to initialize an X2SYS database using x2sys_init
import os
import pygmt

x2syspath: str = "X2SYS"
os.environ["X2SYS_HOME"] = os.path.abspath(x2syspath)

# Initialize X2SYS database in the X2SYS/ICESAT2 folder
    tag="ICESAT2",  # satellite name
    fmtfile=f"{x2syspath}/ICESAT2/xyht",  # format definition file
    suffix="tsv",  # data are in tab separated format
    units=["de", "se"],  # distance in metres, speed in metres per second
    gap="d250e",  # distance gap up to 250 metres allowed
  1. Now comes the hard part, doing the actual crossover calculation using pygmt.x2sys_cross. I store my track data in pandas.DataFrame which keeps things simpler.
track1: pd.DataFrame = pd.read_csv(...)  # 1st file with x, y, h, t columns
track2: pd.DataFrame = pd.read_csv(...)  # 2nd file with x, y, h, t columns

df_crossover: pd.DataFrame = pygmt.x2sys_cross(
    tracks=[track1, track2],
    # region=[-460000, -400000, -560000, -500000],  # set a bounding box region here
    interpolation="l",  # linear interpolation
    coe="e",  # external crossovers
    trackvalues=False,  # Get crossover error (h_X) and mean height value (h_M)
    # trackvalues=True,  # Get track 1 height (h_1) and track 2 height (h_2)
    # outfile="xover_236_562.tsv"  # optionally, output to a file instead of a dataframe

and below this is what the crossover output would look like. The key things are:

  • the crossover locations (x and y)
  • the time when the two tracks crossed the crossover point (t_1 and t_2)
  • the crossover anomaly (h_X) which is the difference between h_1 and h_2
  • the mean (average) height value which is (h_1 + h_2) / 2
               x             y                 t_1                 t_2  \
0 -216453.946304 -611770.40777 2019-09-07 11:53:02 2019-09-11 21:07:28   
1 -216453.946304 -611770.40777 2019-09-07 11:53:02 2019-12-11 16:47:17   
2 -216453.946304 -611770.40777 2019-09-07 11:53:02 2020-06-10 08:06:50   
3 -216453.946304 -611770.40777 2019-09-07 11:53:02 2020-09-09 03:46:38   
4 -216453.946304 -611770.40777 2020-03-07 03:12:36 2019-09-11 21:07:28   

         dist_1        dist_2      head_1      head_2        vel_1  \
0   3296.598252    703.294816  359.881315  219.076615  6878.997551   
1   3296.598252   7107.346604  359.881315  219.076615  6878.997551   
2   3296.598252  13511.398391  359.881315  219.076615  6878.997551   
3   3296.598252  19915.450178  359.881315  219.076615  6878.997551   
4  10981.312349    703.294816  359.881315  219.076615  6906.821748   

         vel_2       h_X        h_M  
0  6879.319824  0.004869  62.719589  
1  8267.201695  0.006300  62.718874  
2  6879.319824 -0.056660  62.750354  
3  6910.273551 -0.754562  63.099305  
4  6879.319824 -0.063100  62.685604  

If you set the parameter trackvalues=True, then you will get the raw height values (h_1 and h_2).

Hopefully this will get you started. I’ve simplified things a lot, but there are some quirks with using X2SYS (e.g. data have to be in specific folder, the environment variable is confusing, etc) so let us know if you need more help (please post your error messages).

If you provide an example or sample of the data you would like to do crossover analysis on, that would also be of help, and we can give you better advice on what to do :slight_smile:

P.S. Remember to cite the x2sys paper at if you use it in your work!

Thanks for your detailed answer!I’ll read through the links and code provided and post the results here if I make any progress in days. :grin:

1 Like

Hi @weiji14, I have managed to get the cross points based on your suggestion. However, there are still some problems that I would like to ask for your help.

The tracks with the calculated cross points (black dots) are shown below:

I used the same configuration as yours:

df_crossover: pd.DataFrame = pygmt.x2sys_cross(
                                      region=[15, 55, 280, 320],  

The format of the track data I prepared is like this:

and the calculated results:

Although the positions of the cross points are successfully calculated, the calculated t_1 and t_2 are not correct

So do you have any suggestion on how to calculate the t1/t2? Thank you very much:)

Hmm, do you mean that you’re missing some precision in the t_1 and t_2 times? I remember that GMT only supports time up to microseconds according to Support various datetime types as input by seisman · Pull Request #464 · GenericMappingTools/pygmt · GitHub. If you can print out the exact value of t_1 and t_2 using df_crossover.t_1, it might be easier to inspect the actual floating point datetime value.

Also, I did notice that your h_1 and h_2 values are nan, so you might need to adjust the gap parameter in x2sys_init to be longer than 250 metres. Otherwise you won’t be getting any valid crossover heights.


你好 @lizhennan007, 你可以直接在这里提出你的问题,这样方便更多人看到然后可以帮你解答。或则你可以试看加入https://gmt-china.org的QQ群?


You can find something in :Need help for GMT x2sys adjustment

The question is: When getting crossovers, do we need to separate the ascending and descending track data.

Answer: In short, no, not necessarily. If you don’t set the coe parameter, the default with x2sys is to use find all crossovers (internal and external). However, if you do separate the ascending and descending tracks into separate files/tables, then you can set coe="e" for external crossovers only (i.e. find crossovers between ascending and descending tracks).

To be clear, if you have ascending track (A) and descending track (D), then x2sys’s default is to calculate all crossovers A-A, A-D, D-D. Using coe="i" would calculate A-A and D-D only, and using coe="e" would calculate A-D only.


比如说你的数据分为上升(A)和下降(D),那么x2sys默认会计算所有的交叉点A-A, A-D, D-D。使用coe="i"只会计算A-A和D-D,使用coe="e"只会计算A-D的交叉点。

我没有分上升下降轨道,若设置内部交叉即Qi时候就是上升(A)-A,下降(D)-D;外部交叉就是Qe即A-D这样子对吧。。默认会求所有交叉点。我用一个周期的数据算交叉点,给定R96/130/-5/30范围,fmt文件是lon,lat,ssh,time,但是结果显示时间栏是NAN,下边是图,你看一下。不知道弄,数据格式以及fmt哪里有问题,可以告诉我一下吗?谢谢,还有一个问题就是gmt cross求交叉点原理是什么在哪里能找到啊

Uploading: image.png… image
数据依次为lon lat ssh time

杨老师您好,我之前看过您的HY2B例子,我老师给我的例子就是您之前写的,但是在time设置部分有一些问题,我看介绍,要求是轨道开始时间/结束时间/轨道信息 这几个缺失,我不知道准备数据哪里出了问题

麻烦把xys.fmt文件里的文字复制出来。还有,可以的话也把您的代码和 ers_2_1.dat 里一部分的数据复制到这里,这样比较方便看出问题。x2sys算出来的结果也可以给我们看一下吗?

目前来看,你的时间(time)是GPS time吗?或许要改一下oformat,可是从图里看不出你用什么。如果你的ers2_1.dat文件里有multisegment tracks(取决于> groud track 是单独的header还是multisegment header),那你可能得在xys.fmt文件里加一行# MULTISEG。详情请参考x2sys_init — GMT 6.3.0 documentation

Wessel, P. (2010). Tools for analyzing intersecting tracks: The x2sys package. Computers & Geosciences, 36(3), 348–354. Redirecting

All x2sys modules have been translated into Chinese.
x2sys_init — GMT 中文手册 (

1 Like


你这里的意思是 rtime 计算的交叉点的的速度和时间都不是 NaN, 用 time 的话是 NaN 吗?

rtime 是相对时间,是相对 TIME_EPOCH 设置的 1970 年的。
如果是测高数据的话,提取的时间都是相对时间,所以用 rtime 能计算出结果。

但这还是不对的,测高的时间系统是相对 1985 年的,而不是 1970 年,也就是要加一句

gmt set TIME_SYSTEM S1985

所以实际上你只要用了上面的设置以及 fmt 中用 rtime 应该就可以了

参见 x2sys_init — GMT 中文手册 (



各位老师好,我还有问题想请教。我用x2sys_cross ers2.dat -Ters2 -Qi -R 求出的应该是自交叉点,ers2.dat是我提取出的一个周期所有经纬度高程时间数据。那如果我要求的是同一颗卫星,不同周期间的交叉点呢,是不是写成x2sys_cross ers201.dat ers202.dat -Ters2 -Qi -R 这样子呢