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 https://docs.generic-mapping-tools.org/6.2/supplements/x2sys/x2sys_init.html and https://docs.generic-mapping-tools.org/6.2/supplements/x2sys/x2sys_cross.html
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 https://github.com/weiji14/deepicedrain/blob/v0.3.0/atlxi_lake.ipynb (scroll down to the ’ Crossover Track Analysis’ section).
#
# 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
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
pygmt.x2sys_init(
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
)
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],
tag="ICESAT2",
# 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:
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 
P.S. Remember to cite the x2sys paper at https://doi.org/10.1016/j.cageo.2009.05.009 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. 
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(
tracks=tracks,
tag=x2sys_tag,
region=[15, 55, 280, 320],
interpolation="l",
coe="e",
trackvalues=True,
)
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.
你好，我刚开始求交叉点，有接个问题想请教你，方便加微信吗？18435300280
是这样，我不知道你的数据是怎么样的，求教差点需要把上升下降轨道数据分开吗？
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.
简单来说是不需要分开的。如果你不设定coe，那么x2sys会自动找所有的交叉点（包括内在的和外在的）。但是，如果你把上升和下降的数据分开了，那你可以使用coe="e"的设定，这样可以求外在的交叉点（就是上升数据和下降数据的交叉点）。
比如说你的数据分为上升（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求交叉点原理是什么在哪里能找到啊
杨老师您好，我之前看过您的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
你好，这个图就是我的fmt里的内容；
这个就是数据文件；没有区分上升下降轨道，
这是我的sh文件，这样子算出来的结果vel1，vel2显示为NAN，
这个是cygwin提示的东西，我试过把fmt文件中，tim改成time，那样的话，ssh也就是海面高及之后的数据列全部为NAN；我把tim改成rtime，都可以显示出来，目前是这样子的，我不知道这是为什么？请教各位老师一下。
我预期的成果是，获取经纬度，时间的显示为s就可以，不用转换为日期格式，目前的话是可以获取的，但不知道为什么会这样，尤其是fmt中时间为什么要写成tim才行
你这里的意思是 rtime 计算的交叉点的的速度和时间都不是 NaN, 用 time 的话是 NaN 吗？
如果是这样的话，我想应该是这样的：
rtime 是相对时间，是相对 TIME_EPOCH 设置的 1970 年的。
如果是测高数据的话，提取的时间都是相对时间，所以用 rtime 能计算出结果。
但这还是不对的，测高的时间系统是相对 1985 年的，而不是 1970 年，也就是要加一句
gmt set TIME_SYSTEM S1985
所以实际上你只要用了上面的设置以及 fmt 中用 rtime 应该就可以了
好的，谢谢老师
嗯嗯好的，谢谢老师
各位老师好，我还有问题想请教。我用x2sys_cross ers2.dat -Ters2 -Qi -R 求出的应该是自交叉点，ers2.dat是我提取出的一个周期所有经纬度高程时间数据。那如果我要求的是同一颗卫星，不同周期间的交叉点呢，是不是写成x2sys_cross ers201.dat ers202.dat -Ters2 -Qi -R 这样子呢