GMT Project projecting all data points to cross-section

I have plotted a map with data points. Purpose is to extract data points using a cross-section from map. For this I used GMT Project to project all data points within specified width to cross-section line on map. The code and map are;

gmt begin 
gmt figure Map png
gmt basemap $Proj $Limit -Bxa1 -Bya1 -BWESN
awk '{print ($1,$2)}' data_points | gmt psxy $Proj $Limit -Sc0.24 -Gblue
# draw line on map with values
echo 101 -4 >    Input_line
echo 101  4 >>  Input_line
gmt psxy Input_line $Proj $Limit -W2,brown
# project events to cross-section line
gmt project data_points -C101/-4 -E101/4 -W-110/110  > Projected_events
gmt end


data_points: data_points.dat (337 Bytes)

I want to project only those data points which lie within width(-110/110 km) on either side of cross-section line but “GMT project” projected all data points to cross-section. I have follwoing queries;

  1. Why “GMT project” projected all data points to cross-section line while width (-110/110) is already specified? If “Lw” is used in “GMT project”, it again project nearly all data points. How to project only specified width data points to cross-section?

  2. I want to draw a boundary/box around cross-section line which shows the projected points area on map as shown in below figure.

How user can draw such a box around cross-section line on map such as in this case to draw box with width (-100/100km) around cross-section? So it clearly shows the area of cross-section extraction for data points.

Kindly assist me in these issue. Thanks

Hard to know without seeing the input file and the result. However, I suspect that you should use -Q.


I used -Q and it has done the job to restrict projecting those events which are defined by width.

For understanding, -Q used to specify units (Q detail in GMT project ; Specify x, y, r, s are in degrees while p, q, dist, lmin, lmax, wmin, wmax are in km. If -Q is not set, then all these are assumed to be in the same units).

Do you know the logic behind this how (Q) is restricting to project events to cross-section based on defined width while it is used to specify output data units?

Sorry @Meyal. I didn’t see the attach file.

The logic is that -Q says that, in your case, the values for -C and -E correspond to geographic coordinates (longitude and latitude).

I have applied with -Q and it done the task. I think it solved the issue. Files are attached;

data_points.dat (337 Bytes)
Project data without using Q;
Command: gmt project data_points.dat -C101/-4 -E101/4 -W-100/100 -Lw > Projected_data_withoutQ.dat
Projected_data_withoutQ.dat (1.3 KB)
Project data using Q;
Command: gmt project data_points.dat -C101/-4 -E101/4 -W-100/100 -Lw -Q
Projected_data_withQ.dat (292 Bytes)

By seeing output files (with & without Q), you can see in file (without Q) it projected nearly all data but (withQ) it projected only those data (5 points) which lie within (W-100/100) on either side of of cross-section of below map.


Hi @Meyal,

Were you able to reproduce that red box indicating the area comprising the selected data points?


Moving to next part, I want to plot my projected data w.r.t cross-section on map like shown in following figure 1;

To get plot like above figure 1 for my data, I used the follwoing way;

gmt begin
gmt figure Map png
gmt basemap $Proj $Limit -Bxa1 -Bya1 -BWESN -Y12
awk '{print ($1,$2)}' data_points | gmt psxy $Proj $Limit -Sc0.24 -Gblue
# plot line 
echo 101 -4 >  7.3.Input_line
echo 101  4 >> 7.3.Input_line
gmt psxy 7.3.Input_line $Proj $Limit -W2,brown
# project data within (W-600/600) of cross-section
gmt project data_points -C101/-4 -E101/4 -W-600/600 -Lw -Q	> Projected_data
# plot p vs z (distance vs depth)
gmt basemap -JX12/-8 -R-1000/1000/0/100 -Bxafg+l"Distance (km)" -Byafg+l"Depth (km)" -Y-10
awk '{print ($4,$3)}' Projected_data   |   gmt psxy -JX -Sc0.16 -W0.5 -Gred
gmt end

data_points.dat (337 Bytes)
Projected_data.dat (1.3 KB)
Output map;

I think there is some problem in Output map, problem is;

  1. Most of such maps are plotted with x-axis scale (vary from -10/10 or -50/50 or some value like this). But in my case data (p vs z) is plotted only when I set scale (-1000/1000) see “Cross-section” plot in attached figure. Also all data is plotted on (0 to 1000) side of scale and no data on (0 to -1000) side on map.

I have not see any map which uses such big scale value (-1000/1000) in literature. Kindly have a look is there any problem in my way which I used for projecting data & map plotting?
and please suggest
How I can plot a cross-section like as I attached in above reference figure 1? Thanks.

@Meyal I shared this figure and script that could help you.



BTW, another way to do this is just using -i. Notice that GMT starts counting in 0.

gmt psxy -JX -Sc0.16 -W0.5 -Gred -i3,2

Notice that p (the 4rd column on Projected_data_withQ.dat) is the distance along the profile. It is 0 at the long and lat of -C. I.e. In your case, increased from South (at 4°S) to North (4°N).

If you want to plot distances to the profile, you should use the value q (the 5th column on your projected data).

A map like available in your shared github source can be plotted and attached below;

Thanks for sharing information about using -i option, it make work easy.

I understand your tip, if user want to plot distances along profile ( p ) value to be used while for distances to profile (q) do the job.


About this topic, 2 unsolved queries are still remaining, which are;

  1. How user can draw a box around cross-section profile which will show area of selection for data on map like shown in below figure;


  1. How to plot cross-section map for projected data of this topic (figure 1) like scale (x-axis) as shown in reference figure 2.

figure 1;

reference figure 2;

If you know about these or someone have the knowledge can share the way to solve these queries. Will be thankful for help.


Hi, not find a solution yet to plot bounding box/boundary around the cross-section profile.

I also would like to know how to do it. But I don’t know a direct way to do it.
I check gmt plot -L but I don’t realized how to do it (if possible).

After search, I have plotted a map for my data like cross-section shown in reference figure (focus on x axis scale). I used s vs z (column 7 vs column 3) in Projected_data.dat and output map are;

Projected_data.dat (1.3 KB)

I think using s values vs z for plotting a map with scale (x-axis of reference figure) is correct.

Have a look on Projected data and plotted map to see, whether using s values in such case is correct? Thanks.

I think that there is no direct way to do it. Maybe we should make a feature request. I will ping @pwessel.

Since in general this band would be oblique I think the only good solution is to add a new option to project that saves the polygon that reflects your search area. Feature requests are done on the GitHub issue page.


What could be an alternative way to achieve this? I have been experimenting with different approaches, and here’s what I have tried so far:

  • First, I calculate the azimuth of the two perpendicular segments based on the azimuth of the cross section just by adding 90°;
  • Next, I use project to obtain the coordinates of the two perpendicular segments.
    However, when I plot the rectangular polygon on a map using a Mercator projection, the resulting rectangle does not appear correct. What am I doing wrong?
    I have attached a figure that shows the resulting rectangle (dashed rectangle) obtained using this approach and the ideal rectangle that I would like to achieve.poly_tr_ale