Help with PyGMT X2SYS_CROSS

Q&A PyGMT Q&A
1

Hey there!
I am new to pygmt. I am trying to calculate the intersections between two ship tracks using x2sys_cross.

The tracks have 4 points of intersection but x2sys_cross outputs only one intersection.

track_plot
track_plot913×737 111 KB

I wrote the code to calculate the intersection following the code posted on another question in the forum.

# running x2sys_init
pygmt.x2sys_init(
    tag="mag",  
    fmtfile="xyzt",  
    suffix="xyzt",
    discontinuity="d", # geographical coordinates discontinuity at dataline
    spacing="10m",    
    units=["dk", "sn"],  # distance in kilometers, speed in knots
    gap="d1k",  # distance gap up to 1 km allowed
    force=True, # replace old files
    verbose="w"
)

# reading ship tracks as dataframes
track1: pd.DataFrame = pd.read_csv('dsdp74gc.xyzt',sep='\t',skiprows=1)  # 1st file with x, y, z, t columns
track2: pd.DataFrame = pd.read_csv('rc2711.xyzt',sep='\t',skiprows=1)  # 2nd file with x, y, z, t columns

# running x2sys_cross
df_crossover: pd.DataFrame = pygmt.x2sys_cross(
    tracks=[track1, track2],
    tag="mag",
    interpolation="l",  # linear interpolation
    coe="e",  # external crossovers
    speed=["l1", "u15", "h1"], 
    #trackvalues=False,  
    trackvalues=True  # Get track 1 height (h_1) and track 2 height (h_2)
)

# Results of Cross over
print (df_crossover)

        lon        lat                 t_1                 t_2       dist_1  \
0  6.041174 -28.864184 1980-06-25 12:18:54 1986-12-07 14:39:15  1472.200472   

       dist_2     head_1      head_2    vel_1      vel_2        z_1        z_2  
0  137.521438  20.389965  293.781108  8.29052  10.143803  79.679097  96.833434

Maybe I’m missing something. Can anyone please help?

This zip folder has the script, input data and the format definition file.
data.zip (130.4 KB)

Appreciate any leads!

2

Hi @Thor524! I haven’t had time to check out your data yet, but just to reply you real quick, could you try increasing the distance of the gap (-W in GMT) parameter from 1km to something more? You might also need to try increasing the spacing (-I in GMT) value from 10m to something higher. See https://www.pygmt.org/v0.4.0/api/generated/pygmt.x2sys_init.html and https://docs.generic-mapping-tools.org/latest/supplements/x2sys/x2sys_init.html#i for more details.

If that still doesn’t work, leave a reply and we’ll have someone take a closer look at your data to see why the crossovers are not being picked up :slight_smile:

3

Hi @weiji14! That did the trick!! x2sys is identifying all the crossovers.

I have a few related questions -

  1. The doc says -Wgap is for spline interpolation. Since I’m using linear interpolation, is this option still relevant?

  2. What is the use of the spacing parameter? I don’t completely understand how it affects the crossovers?

Thanks in advance :slightly_smiling_face:

4

Great, glad to hear it worked :grinning_face_with_smiling_eyes: To answer your questions:

Reading x2sys_cross — GMT 6.2.0 documentation, my intuition is that the documentation for the gap (-W) flag is just misworded, and that it should read as “Give the maximum number of data points on either side of the crossover to use in the interpolation [3].” (i.e. without the 'spline). But best to check with @pwessel or the GMT core team who know more.

The spacing (-I) parameter in x2sys_init shouldn’t affect your crossover calculation too much (other than in some very specific cases maybe?). According to x2sys_init — GMT 6.2.0 documentation, it says “These spacings refer to the binning used in the track bin-index data base.”. So say you choose a 50m spacing, the point data will then be grouped into 50x50m squares. This might affect the speed at which the algorithm finds crossovers? But I’ll defer again to @pwessel. See also Sec 2.2 of the X2SYS paper (https://www.sciencedirect.com/science/article/pii/S0098300409002945#sec2).

5

The x2sys_init settings for -I does not affect x2sys_cross in any way. It is for organizing line coverage via x2sys_binlist only. The -W option in x2sys_init determines when the distance (or time) between two consecutive points on a track are just too large to consider the line continuous. When it exceeds the gap we basically break the track into sections and no crossings found in those gaps are reported. So if gap is too small you will find no crossings…

6

Hi @weiji14, @pwessel ! Thanks for the clarification.

On a separate note, I used x2sys_solve to get the correction for the crossover bwetween these two tracks.

When I used constant offset correction (-Ec), I get these corrections, which seem way too high

dsdp74gc z 632.936
rc2711 z -632.936

But when I used linear drift (-Ed), I get these, which look right

dsdp74gc z 3.3087 -2.16173e-06*((dist))
rc2711 z -3.3087 0.0040923*((dist))

Is there any issue with the constant offset correction?

Also @weiji14, I was going through some of your code on Github. Most of them are calculating the crossovers. Do you have any code to solve the errors as well? :slightly_smiling_face:

Appreciate the help!!

7

If you can post the input file and options used with x2sys_solve then I can have a look.

8

Hi @pwessel! Here is the data set Data_x2sys_solve.zip (129.5 KB). It has two input files, fmt file and the shell script I used.

Thanks.

9

Looks like your x2sys_list command is asking for too much output and then feeding that to x2sys_solve means it reads the wrong information. See the example in the x2sys_solve documentation for a simple -Ec model - I think you just need

gmt x2sys_list $TAG.CROSS -T$TAG -Cz -Fnc -Qe -Vq > “${TAG}.LIST”

10

In my case, the crossover error/anomaly was actually the signal I’m interested in, so no, I don’t have code any code for solving the errors unfortunately. But I’ve put in a feature request for x2sys_solve at Wrapper for x2sys_solve · Issue #1406 · GenericMappingTools/pygmt · GitHub. Hopefully someone will give it a go!

11

Hi, I used the same code and data as Thor524 , but I have this error： ‘ No column to parse from file ’. I’ve used other examples in the forums with this problem as well. Is there a good way to deal with it? :rofl:

数据出错
数据出错1920×864 304 KB

12

Hi @jcf, the EmptyDataError means that no crossover points are found (because the gap is too small). You’ll need to increase the value of the gap (-W) parameter first as I suggested above.

13

I modified the gap value as you suggested and got the correct result. Thank you very much. One other little problem. Does PyGMT have a correlation function that outputs points’ values on the original track(s) near the crossover points? :stuck_out_tongue:

14

Not sure if I’m understanding this correctly, but if you want to get the original values on each track (say t_1 and t_2), set the parameter trackvalues=True in pygmt.x2sys_cross, which will report the height_1 and height_2 values. The default trackvalues=False would only give height_X (the crossover value interpolated from height_1 and height_2) and height_M (mean of height_1 and height_2).

From there, I assume you’ll need to run a correlation function on the height_1 and height_2 columns. That will need to be done in another package, maybe pandas.DataFrame.corr or something else.

15

Thank you for your reply. What I want to find is the observation point data on the original tracks closest to the crosscrovers. I found those points by finding the minimum value difference of longitude and latitude. :thinking:

16

Ah I see, then you definitely want to set trackvalues=True to get those original observation point values. Unfortunately the distance from the crossover to the closest original point isn’t reported, so yes, you’ll have to work it out yourself from the longitude/latitude values :sweat_smile: