Plotting Median Stacked Profile and then plotting stations coming in between

lalit.seismo · May 20, 2023, 1:41pm

Hi I am using GMT 6.0 version in my windows system. I tried to plot Median Stacked profile in bathymetry area as can be seen in the figure but I am trying to plot stations depth as well on stacked line but no success. Can anyone suggest how can I approach this? I tried to use gmt project as well but somewhere I am making mistake. Can anyone please help me? I am ready to share data files.

This below is the script I am using but no success so far! The same example is given on gmt website as example 33.

set int=C:\programs\gmt6\GMT_Course\gridfiles\GEBCO_2019.nc
::set grd=C:\programs\gmt6\GMT_Course\gridfiles\etopo2.nc

::set range=127.3/131.3/27/30.4
set range=126/132/26/31

gmt begin SFM_Profiles jpg A+m0.5c -E600
::gmt makecpt -Cwiki-2.0.cpt -T-6000/4000 -Z
gmt grdcut @earth_relief_15s -G%int%  -R%range% -V
gmt grdimage %int% -R%range% -JM18c  -I+a300 -Cglobe -V
 gmt coast -Df -Bxa1f0.5 -Bya1f0.5 -BWSEN+t"Bathymetry Depth"  -W0.8p --FONT_TITLE=35,31,black -V
::gmt basemap  -Lx1c/19c+w50k+l+f+jBL --FONT_TITLE=30,25,navy -V
     
  gmt colorbar -C%cpt%  -Bxa2000f1000+l"elevation" -By+l"m" -Dx8c/-1.8c/9c/0.5ch -I0.8 -V


:: Select two points along the ridge
echo 128 27 > ridge.txt
echo 130.5 30 >> ridge.txt
gmt plot -Sc0.1i -Gblue ridge.txt
gmt plot ridge.txt -W1p,blue
gmt grdtrack ridge.txt -G%int% -C360k/2k/10k+v -Sm+sstack.txt > table.txt
gmt plot -W0.5p table.txt
gmt plot OT_Stations.txt -i1,0 -Sa17p -W0.5p -Gyellow 
gawk "BEGIN {FS=\",\"} {print $2, $1, $3}" OT_loc.csv | gmt text -F+f17p,Times-Roman,red,thick+jBL  -Dj0.4p -N 

gmt convert stack.txt -o0,5 > env.txt
gmt convert stack.txt -o0,6 -I -T >> env.txt
::gmt plot env.txt -R-400/400/-3250/250 -JX6i/2i  -W1p,red -Bxa100f10+ukm -Bya1500f250+lElevation+um -BWnSe+tProfiles -X8i -Y3i -V

gmt plot -R-230/230/-8500/1000 -Bxafg1000+l"Distance from middle section (km)" -Byaf+l"Elevation (m)" -BWSne+t"MEDIAN STACKED PROFILE" -JX9i/4i -Glightgray env.txt -X8.5i -Y2i 
::gmt plot env.txt -R-200/200/-3500/-2000 -JX6i/3i  -W1p,gold  -V

gmt plot -W3p stack.txt
::gmt project OT_Stations_depth.txt -Cenv.txt -E$(awk '{print $1, $2}' ridge.txt) -Fpz -V | gmt plot -Sa0.15i -Gblack -W0.5p
::gmt project OT_Stations_depth.txt -Cenv.txt -Q -Fpz -V | gmt plot -Sa0.15i -Gblack -W0.5p
::gmt plot OT_Stations_depth.txt -i1,0 -Sa0.15i -Gblack -W0.5p
::gmt project stack.txt -C128/27 -Fpz | gmt plot -W3p,black

::gmt project OT_Stations_depth.txt -C128/27 -Fpz | gmt plot -Sa0.15i -Gblack -W0.5p


::gmt plot OT_Stations_depth.txt -i1,0,3 -Sa0.15i -Gblack -W0.5p

::gmt legend -DjBL+w2i/2.3i+l2+o0.25c/20c -C0.15i/0.1i -F+p+g255/228/225+p+r
gmt end  show

 del *.info *.conf *.int *.history *.grd
pause

Esteban82 · May 20, 2023, 1:28pm

Hi @lalit.seismo Welcome to the GMT forum!

I think that you should use project as you did here:

gmt project stack.txt -C128/27 -Fpz | gmt plot -W3p,black

Probably you have to select the columns (with -i) that you want to plot.

I suggest to take a look at the output of gmt project. You could share here if you want.

lalit.seismo · May 20, 2023, 1:35pm

Hi @Esteban82
Thank you for your reply.
Actually stack.txt file has multiple columns. I am sharing here please take a look. I am confused which column will be suitable because mean is calculated using these. stack.txt (17.9 KB)

Esteban82 · May 20, 2023, 1:39pm

Sorry, I got the wrong line.I mean the command with OT_Stantions_depth.txt

::gmt project OT_Stations_depth.txt -C128/27 -Fpz | gmt plot -Sa0.15i -Gblack -W0.5p

pwessel · May 20, 2023, 1:41pm

Think you use the ridge coordinates to set -C and -E in project, then project your data pints and use -Fqz where q is the cross distance.

And don’t be so hard on your self that you are using 6.0. Surely you can use 6.4 and avoid a bunch of bugs and enjoy new features?

lalit.seismo · May 20, 2023, 2:42pm

Hi I tried using -i1,0 as file has first column lat second long and then depth in 3rd:

gmt project OT_Stations_depth.txt -C128/27 -Fqz | gmt plot OT_Stations.txt -Sa0.15i -i1,0 -Gblack -W0.5p

but it has plotted star randomly. Please see the figure.

And also this file looks like thisOT_Stations_depth.txt (541 Bytes)

Please suggest something. I am also trying to plot 0 to -600m depth stacked (black line) make pink or any other color so that it can be seen clearly.

lalit.seismo · May 20, 2023, 3:07pm

Hi, Thanks for your suggestion. Actually -C option I am using as -C128/27 but unable to understand how to use -E.
Please suggest where I am making mistake:

gmt project OT_Stations_depth.txt -C128/27 -Fqz | gmt plot OT_Stations.txt -Sa0.15i -i1,0 -Gblack -W0.5p

txt file also have attached previously which has stations lat,long and depth.

pwessel · May 20, 2023, 4:27pm

Not so easy. This gives the q distances

gmt project -E128/27 -C130.5/30 -Fqz -Q OT_Stations_depth.txt

but those distances are not zero on the ridge line (like your stack). So maybe run the same project command on the C or E point and then get its q value and subtract from the star q values and that will work?

Esteban82 · May 20, 2023, 5:17pm

@pwessel is this output wrong? I was expecting the depths value first.

 echo 28 127 1127 | gmt project -C128/27 -E130.5/30 -fg -: -Fzq
1.30227717218   1127

Esteban82 · May 20, 2023, 5:23pm

I add -; because the input is Lat Long. But, it seems that -; also flip the output.

echo 28.6239 127.7732 1127 | gmt project -C128/27 -E130.5/30 -Q -Fqz
-685.31053591   1127

Esteban82 · May 24, 2023, 7:26am

@lalit.seismo does this map make sense to you?

I use -: in project becasue your data is Lat Long (and GMT expects Long Lat). In plot, I multiply the depth by -(+s-1) because in your data it was positive.

Here is the script:

gmt begin test png
gmt basemap -R-230/230/-8500/1000 -Bxafg1000+l"Distance from middle section (km)" -Byaf+l"Elevation (m)" -BWSne+t"MEDIAN STACKED PROFILE" -JX9i/4i
gmt project OT_Stations_depth.txt -C128/27 -E130.5/30 -Q -Fqz | gmt plot temp -Sa0.15i -Gblack -W0.5p -i1,0+s-1
gmt end

pwessel · May 20, 2023, 5:33pm

z is always last since may contain trailing text. Probably not well documented…

Esteban82 · May 20, 2023, 5:37pm

But here I got z first (it is the same example but I manually put the long first)

echo 127.7732 28.6239 1127 | gmt project -C128/27 -E130.5/30 -Q -Fzq
1127    123.183378965

Esteban82 · May 20, 2023, 5:39pm

I think it is very clear:

If output format is ASCII then z also includes any trailing text (which is placed at the end of the record regardless of the order of z in flags ).

pwessel · May 20, 2023, 5:39pm

Yes, I spoke before looking. Just committed a minor PR so that the usage message also mentions this.

lalit.seismo · May 21, 2023, 3:04am

Hi @Esteban82, Thank you for explaining it. Yes, plot must look like this.

gmt project OT_Stations_depth.txt -C128/27 -E130.5/30 -Q -: -Fqz | gmt plot -Sa0.15i -Gblack -W0.5p -i1,0+s-1

Only thing i changed for right depth reversed the -C and -E long/lat and removed temp

Though there is one outlier station. I did not understand how come it’s going outside the range? Can we write stations name or number as well here?

One more question i’ve, can we colored a particular section of stacked line? like from the depth o to -600 m by any how? For this i also searched a lot but no success so far!

pwessel · May 21, 2023, 3:31pm

If you want to label the stars you have to supply what you want the label to be. Probably pipe your gmt project stations via text using the +r modifier to -F to select running record number as the text.

As for color. If you mean color the solid stacked line it depends if you just want it to have a constant color between two z-values. If so, maybe run the stacked curve through sample1d to get intermediate points, then through gmt select to only output the section you want and then plot that with your color.