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
Proj="-JM14c"
Limit="-R95/107/-5/5"
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
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;
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?
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.
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?
I have applied with -Q and it done the task. I think it solved the issue. Files are attached;
Input_data: 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.
I think there is some problem in Output map, problem is;
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.
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).
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;
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.
