Currently I am searching for a convenient approach to generate topographic swath profiles solely be applying GMT tools. In the past I computed swath profiles via home brewed codes and plotted the results with GMT. This was fine as long as the profile just followed a straight line.
However, now I need a profile consisting of several segements and before I start extending my buggy code, I’m just wondering if GMT offers a solution out of the box. I found a similar request but not a working solution at the forum. I’m just wondering if the deluxe solution of grdtrack already exists: https://flint.soest.hawaii.edu/t/grdtrack-vertical-profile-along-track/1547
For my approach I would require: max, min, mean, 1 std (median, 25% and 75% Percentiles) of the cross profiles next to the distance, lat, lon of the original profile line. By reading the grdtrack manual I understand that I could first produce the cross profiles and compute the statistics in a second step. However, it would make me quite happy if I could avoid generating thousands of temporary files.
I was thinking to compute min, max, mean, … values along a profile line within a circular sliding window of a given radius. By reading the grdtrack manual I found that -E in combination with +r might be a solution but I’m not sure how to use it and how to pass a long list of coordinates to -E.
Currently, I just can think of a brute force method but computing spatial statistics (min, max,…) for the entire grid via grdfilter and using grdtrack to extract the desired values along the profile line. Anyhow, I’m quite sure that there are more elegant solutions for such a task.
Hi! I am also working on the similar thing with you recently, about Swath Profile. However, it is so sorry that I have not found the specialized commands to compute and plot swath profiles in GMT yet. I do not know whether you have found some solution. Anyway, now I have a indirect way to it. If it is still helpful for you, you can see the following. (I am so sorry my strategy is a little complex and long … )
Okay, I choose to create swath profile based on the function stacking (-S option) of multiple cross-profiles (-C option) in GMT command grdtrack.
Next let me introduce my method to you. I recommend you glance through the output figure in the last part when you read my bash script. And I have written the necessary interpretation in this bash script, for clarity, as following. Also, the following blod font represents the variable using in the bash script or the important element in the method (shown in the figure).
Input files:
relief_TibetanPlateau_forplot.nc, a example topography file.
mychem.txt, a data file including data which are expected to project to the same swath profile, with the format of lon, lat, value:
95 32 5000
##Define swath profile #Here you should select your desired survey line, with start point and end point or start point and azimuth clockwise from north and length of survey lines, as the central thick black line shown in the figure. However, (1). This line cannot accurately be the real survey line for swath profile, and the nearby thick red line will be the real survey line used to project topography and data. That may be because of the existence of grid spacing of DEM data and spherical earth. For example, I define a line with GMT, then I compute its orthogonal line. However, if I compute cross profile which is orthogonal to its orthogonal line using grdtrack, they can not be consistent with the original line. As you see, the thick red line can not consistent with the black line. (2) Whether you define survey line with start point and end point or start point, azimuth and length, both of selection need to compute the azimuth and length (using mapproject) because we will need those information to decide the position of midpoint using project.
#If you define the desired survey line with azimuth and length rather than end point.
#endlon=$( gmt project -C${orilon}/${orilat} -A${azimuth} -G${length} -L0/${length} -Q | tail -1 | awk ‘{print $1}’ )
#endlat=$( gmt project -C${orilon}/${orilat} -A${azimuth} -G${length} -L0/${length} -Q | tail -1 | awk ‘{print $2}’ )
gmt plot -W2p,black << EOF
${orilon} ${orilat}
${endlon} ${endlat}
EOF
#We must need to find the position (cenlon,cenlat) of midpoint of the desired survey line, as the central blue circle in the figure.
cendist=$( echo “${length}/2.0” | bc -l )
cenlon=$( gmt project -C${orilon}/${orilat} -A${azimuth} -G${cendist} -L0/${cendist} -Q | tail -1 | awk ‘{print $1}’ )
cenlat=$( gmt project -C${orilon}/${orilat} -A${azimuth} -G${cendist} -L0/${cendist} -Q | tail -1 | awk ‘{print $2}’ )
echo ${cenlon} ${cenlat} | gmt plot -Sc0.2 -Gblue
#We find one endpoint (cen1lon,cen1lat) of the orthogonal line, as the blue circle on the border in the figure.
cen1azimuth=$( echo “${azimuth}-90” | bc -l )
cen1lon=$( gmt project -C${cenlon}/${cenlat} -A${cen1azimuth} -G${swathw} -L0/${swathw} -Q | tail -1 | awk ‘{print $1}’ )
cen1lat=$( gmt project -C${cenlon}/${cenlat} -A${cen1azimuth} -G${swathw} -L0/${swathw} -Q | tail -1 | awk ‘{print $2}’ )
echo ${cen1lon} ${cen1lat} | gmt plot -Sc0.2 -Gblue
#We find another endpoint (cen2lon,cen2lat) of the orthogonal line, as the blue circle on the other border in the figure.
cen2azimuth=$( echo “${azimuth}+90” | bc -l )
cen2lon=$( gmt project -C${cenlon}/${cenlat} -A${cen2azimuth} -G${swathw} -L0/${swathw} -Q | tail -1 | awk ‘{print $1}’ )
cen2lat=$( gmt project -C${cenlon}/${cenlat} -A${cen2azimuth} -G${swathw} -L0/${swathw} -Q | tail -1 | awk ‘{print $2}’ )
echo ${cen2lon} ${cen2lat} | gmt plot -Sc0.2 -Gblue
##track and stack along the cross-profiles of CENTRAL NORMAL LINE #Here we have find a new line which is formed by connecting the two blue endpoints, orthogonal to our desired survey line. Next we will determine many cross profiles (i.e. the thin lines in the figure) of the new line using -C option and -S option of grdtrack. Unfortunately, there is a so minor angle shift between cross profiles (the thin black line) and our desired survey line (the thick black line). Then the analysis for stacking cross profiles will produce average value, median value, lower value, upper value and others with -S option of grdtrack.
#Now we have known that the desired survey line (the thick black line) can not be used, and the cross profiles (the thin black line) will be plotted. Here I choose the median cross profile (the central thin black line) from the output file of grdtrack, named stack_profile.txt, which includes the sampling point data for each cross profile. Then the selected cross profile will be the true coordinate system used to plot. So we can find the information about the true line (the central thin black line), such as longitude (newlon) and latitude (newlat) of its midpoint (i.e. the central red circle in the figure) and azimuth (newazimuth). Also, the central cross profile is shown as the central red thick line in the figure, whose data is stored in the file stack_oneProfile.txt. Note: the newazimuth is counterclockwise from South!
#Also, because we will project some data in the input file mychem.txt to the (swath) profile, here we calculate its two endpoints (i.e. the red circles on the endpoints of the red thick lines in the figure). The two endpoints new1 (new1lon,new1lat) and new2 (new2lon,new2lat) will be used in the following project.
##The Real Track #Here we project our data in mychem.txt to the true swath profile (the central thin black line or the red thick line), based on the two endpoint new1 (new1lon,new1lat) and new2 (new2lon,new2lat).
Output files in addition to the figure (.pdf) like in the last part:
swath_chemistry.txt: the projected data along the true survey line.
swath_env.txt: the envolve line which is generated from swath_stack.txt
swath_oneProfile.txt: the cross profile as the final projecting profile.
swath_profile.txt: all cross profiles as represented by the black thin lines in the figure.
swath_stackLOW.txt: distance, lower value of stacking cross profiles (swath profile), and others.
swath_stack.txt: distance, average value of stacking cross profiles (swath profile), and others.
swath_stackUP.txt: distance, upper value of stacking cross profiles (swath profile), and others.
swath_topo.int: the gradient file generate by grdgradient.
This script is a little long and complex. I still want to caution:
There is a angle shift between the really used survey line (the central thin black line or the red thick line) and the desired survey lines (the black thick line).
The bash script will echo some important information: (1) the longitude and latitude of end point (endlon,endlat), azimuth (azimuth) and length (length) of the desired survey line, if defining the desired survey line with start point and end point. (2) the longitude and latitude of midpoint (newlon,newlat), and azimuth (newazimuth) of true survey line counterclockwise from South. (3) the longitude and latitude (new1lon,new1lat) of one endpoint of the true survey line. (4) the longitude and latitude (new2lon,new2lat) of another endpoint of the true survey line.
If you want to plot multiple connecting survey lines, you should use the endpoints of true survey lines (new1lon,new1lat and new2lon,new2lat) rather than the endpoints of your desired survey lines.
)