Sample long and lat from lines at specified distance in km

Q&A
1

Hi everyone!
I’d like to automize the following step in my workflow but I couldn’t find a solution in the documentation, so I am not even sure that’s something GMT can do.

I am currently using a Qgis plugin (QChainage) to sample my cross-section trace (shapefile) at a specified distance including the end point; i.e. I get a point object every 5km along the line and I force the algorithm to include the ending coordinates of the trace even if it’s not at a multiple of 5 km from the starting point. Then I export long/lat of the points to a .csv and feed it to my GMT script to sample several .grd files and plot the cross-section.

Can you point me to the right module(s) for this purpose or tell me if it’s just something currently unfeasible with GMT?

For readers interested in how I plot cross sections with GMT (click here)

Sample the z from several xyz.grd files at the coordinates in the .csv using grdtrack, then use mapproject to get incremental distances between points and finally plot the results on a cartesian graph to obtain the cross section.

In case you need to substitute commas in the .csv with white spaces:

cat csv-file | sed -e ‘s/,/ /g’ > coordinates-file

Sample z values from grid:

gmt grdtrack coordinates-file -Ggridfile -nl -Af > xyz-file

Get incremental distances in km and print “distances - z values” file (x and y to plot):

gmt mapproject xyz-file -Rregion -Gk-to | awk ‘{print $4, $3}’ > dist-z_file

2

Do you need something like example 33?

3

Thanks for that example but at a first glance I’d say no.
I actually need to sample at equidistant points a .gmt file that contains a single line (which may have several vertices) and to include the very last vertex in order to preserve the total lenght of my section.
Am I misunderstanding example 33? I didn’t get my solution in there…still looking the documentation to understand it better anyway.

edit:
what I call cross sections is not 90° from section trace, I may have used the wrong term.
I have 7 .grd files of the same area at different depths underground, which I need to sample at the same coordinates and I already have a script to plot my sections which is working just fine!
Example 33 plots cross profiles along a line, which is awesome, but it’s not what I need. I’m trying to obtain long and lat from a .shp file which I converted to .gmt in order to sample my 7 .grd files at the same location. I upload an example of section, where the black line is the topography and each of the colored lines refers to z values on my .grd files.

4
41928×1508 78.7 KB

@Esteban82 Is that sampling feasible with grdtrack using the .gmt file and specifing the distance of samples?

4

I thnk it is possible with grdtrack but I don’t know how.

You could use sample1d:

From the cookbook:

To obtain a rhumb line (loxodrome) sampled every 5 km instead, use

gmt sample1d track.txt -T5k -AR+l > new_track.txt
5

Thank you very much @Esteban82 !!!
I could have done it also with grdtrack parsing the output to get only long and lat (column 1 and 2) where column 3 was equal to zero but that would have slowed down operations a lot…

Solved with:

gmt sample1d line.gmt -T5k -AR

edit: Note that -AR adjust the spacing to fit exactly the line lenght, to instead use the exact spacing indicated by -T use -Ar

6

sample1d used to have an option to deal with the end points (-S in gmt5) but I don’t find that facility anymore. Anyway since options even after deprecation remain active you should still be able to use that -S in GMT6

7

Tried the code above with:
-S
-Slast_long
-Slast_lat

but I get an error from libgmt.so.6:

ERROR: Caught signal number 11 (Segmentation fault) at /lib/x86_64-linux-gnu/libgmt.so.6(GMT_sample1d+0x1822)[0x7f4feb8190c2] [0x0]
Stack backtrace:
/lib/x86_64-linux-gnu/libgmt.so.6(sig_handler+0x2bb)[0x7f4feb5390db]
/lib/x86_64-linux-gnu/libc.so.6(+0x41950)[0x7f4feb364950]
/lib/x86_64-linux-gnu/libgmt.so.6(GMT_sample1d+0x1822)[0x7f4feb8190c2]
/lib/x86_64-linux-gnu/libgmt.so.6(GMT_Call_Module+0x36d)[0x7f4feb55375d]
gmt(main+0xb80)[0x55a4df3b3ee0]
/lib/x86_64-linux-gnu/libc.so.6(__libc_start_main+0xf2)[0x7f4feb34bcb2]
gmt(_start+0x2e)[0x55a4df3b45ee]

using GMT 6.1 on Ubuntu 20.10. Is the syntax right if I want to keep the last vertex of line.gmt?

gmt sample1d line.gmt -T5k -Ar -Smax_long

8

The syntax should be -S0/total_dist. This assuming that your first column has the accumulated distance along the track and that first value is zero and last = total_dist

9

Ok thanks! No I compute incremental distances after sampling so I guess I can’t apply the -S flag in my case…

10

The old -S is deprecated but same effect is obtained via -T:

-Sstart[/stop] option is deprecated; use -Ttmin/tmax/tinc[+a|n] instead

If I udnerstand you correctly, you want to sample your input track every 5 km, plus output any final coordinate that did not fit into the 5-km stepping. So -Ar will do the first part but not the 2nd part, and -AR violates your 5 km spacing. Seems we would need a 3rd flavor of the -Ar|R schemes.

11

If you had computed the accumulated distances and use them as time there would be no need to the QGis visit.

12

Yes you understood it correctly and it would be totally Awesome to have that option and useful as well! :sunglasses: Just like -Ar but with the final coordinate pair appended to the list and maybe and [INFORMATION] message with the distance from the previous sample.

I can’t compute incremental distances without the equidistant sampling, can I?
Well, of course I could calculate the array measuring the total lenght of the trace and then dividing it by my desired spacing to get the number of rows needed…but then I still need to sample my coordinates adding the final one…I can’t really figure out what you meant there…

I can adjust my 5km rule anyway, this alone brings me out of GUI softwares! gmt sample1d saved me likely a thousand clicks today and who knows how many in the following years! In the end it is just sampling grids to get points to plot the frame for a stratigraphic section, so if the algorithm lowers the sampling rate to fit my lenght that’s not causing any problem at all to my research, anyway it is sampling an interpolated surface, real data is well more spaced than 5 km. It would be different if the spacing was really a strict requirement to the job, e.g. if someone was planning the placement of geophones or GCPs on a trace or systematic sampling on a route at sea. Then the new option proposed above could definitely come in very handy!

13

BTW I leave a link to a post of my blog (in Spanish) that you could fine useful.

1 Like
14

I struggled with what I think is a similar problem just last week (but did not find sample1d… :man_facepalming:). In my case, I need to divide an upper and lower fault trace in to an equal number of segments of approximately 5 km (e.g.), retaining both end points but not caring about intermediate vertices. So if one trace is 26 km long and the other is 24 km, the first will be divided into 5 segments of 5.2 km and the second will be divided into 5 segments of 4.8 km, 6 points on both.

I figured out how to do it with gmt plot -SqD...+t (see figure below) and can now get very similar results using gmt sample1d -AR -Fl -T${step}k -fg (red Xs in the figure). I haven’t had time to investigate what might be causing the differences yet, and, at least for now, they are close enough for my needs. Both methods seem to stray from the specified step length around the vertex in the middle, and the total length differs in the end. I definitely do not yet understand the nuances of -j options and spherical vs. ellipsoidal calculations.

If you want to see a version of the script I'm using, click here. 
targetStep=5 #km
# get trace lengths - maybe increase number of digits in output?
len1=$(gmt spatial -Qk trace1.gmt | cut -f 3)
len2=$(gmt spatial -Qk trace2.gmt | cut -f 3)
# calculated required number of steps and respective step lengths
nSteps=$(echo "scale=0; ($len1 + $len2) / $targetStep / 2" | bc)
step1=$(echo "scale=8; $len1 / $nSteps"  | bc)
step2=$(echo "scale=8; $len2 / $nSteps" | bc)
# sample lines, optionally pipe through mapproject for incremental and cumulative distance
gmt sample1d trace1.gmt -AR -Fl -T${step1}k -fg | gmt mapproject -G+uk+a+i > trace1-pts.gmt
gmt sample1d trace2.gmt -AR -Fl -T${step2}k -fg | gmt mapproject -G+uk+a+i > trace2-pts.gmt
# less conveniently...
gmt plot trace1.gmt $(gmt info -I.1 trace{1,2}.gmt) -Jm1:1 -N -SqD${step1}k:+LD+ttmpPts.gmt > /dev/null
awk 'NR==1{print $1,$2,0,0}' trace1.gmt > trace1-pts-alt.gmt # add angle and incremental distance
cat tmpPts.gmt >> trace1-pts-alt.gmt
awk -v len=$len1 'END{print $1,$2,0,len}' trace1.gmt >> trace1-pts-alt.gmt
rm tmpPts.gmt
# create segment boundaries (for illustration)
paste trace1-pts.gmt trace2-pts.gmt | awk '{printf(">\n%s %s\n%s %s\n",$1,$2,$5,$6)}' > segments.gmt

westFlank
westFlank1131×1235 28.9 KB

15

@jaltekruse I’m glad this post has already been useful to someone! That’s an interesting application but I don’t get the difference between the two methods from the figure alone…maybe the linear interpolation you apply with -Fl is causing small differences only noticeable comparing numbers, I’d try using -Fn or just without the -F option if you haven’t yet!

Also make sure that the faults lenghts are exactly multiples of their relative spacing if you need it to be be maintained by -AR because:

@Esteban82 your blog is awesome! Thank you for the link, I’ll definitely learn a lot from it!
I made a similar script but I just set the y axis and calculate the x axis dimension to get my desired vertical exageration instead of setting the dimension of the figure to then calculate the v.e., that’s useful if you have to compare different sections from the same study (which is actually my case).

2 Likes