Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

plot/plot3d: Simplify the documentation for the -A option #8645

Open
seisman opened this issue Dec 25, 2024 · 9 comments · May be fixed by #8674
Open

plot/plot3d: Simplify the documentation for the -A option #8645

seisman opened this issue Dec 25, 2024 · 9 comments · May be fixed by #8674
Labels
documentation Improve documentation
Milestone

Comments

@seisman
Copy link
Member

seisman commented Dec 25, 2024

https://docs.generic-mapping-tools.org/latest/plot.html#a

The -A option has the following syntax:

-A[m|p|x|y|r|t]

in which m|p works for geographic map, x|y works for Cartensian plot, and r|t works for polar projections.

But in the source code, m/y/r works the same way, and p/x/t works the same way.

gmt/src/psxy.c

Lines 838 to 848 in b34b125

case 'A': /* Turn off draw_arc mode */
n_errors += gmt_M_repeated_module_option (API, Ctrl->A.active);
switch (opt->arg[0]) {
case 'm': case 'y': case 'r': Ctrl->A.mode = GMT_STAIRS_Y; break;
case 'p': case 'x': case 't': Ctrl->A.mode = GMT_STAIRS_X; break;
#ifdef DEBUG
default: Ctrl->A.step = atof (opt->arg); break; /* Undocumented test feature */
#endif
}
break;

So, instead of showing the -A[m|p|x|y|r|t], we can simplify it to:

-A[x|y]

then the documentation can be something like:

x: Draw the line along x first then along y.
y: Draw the line along y first then along x.

Here, x and y have different meanings in different coordinate systems.

  • For Cartesian data, x and y means X- and Y-axis, respectively
  • For polar projection, x means theta and y means r
  • For geographic projections, x means parallel and y means meridian

I'd like to work on it if approved.

@seisman seisman added the documentation Improve documentation label Dec 25, 2024
@joa-quim
Copy link
Member

Yes, it will look better that way.

@seisman
Copy link
Member Author

seisman commented Dec 31, 2024

I'm updating the documentation of this option on the PyGMT side first, then @yvonnefroehlich found an issue at GenericMappingTools/pygmt#3720 (comment), which I think is a GMT bug. Here is the CLI version just in case you want to try it:

gmt begin map

gmt basemap -JH10c -R0/360/-90/90 -Bafg -B+t"-Am"
echo 90 -50 | gmt plot -Sc0.2c -Gred 
gmt plot -W1p -Am << EOF
90  -50
150 50
EOF

gmt basemap -JH10c -R0/360/-90/90 -Bafg -B+t"-Ap" -Xw+1c
echo 90 -50 | gmt plot -Sc0.2c -Gred 
gmt plot -W1p -Ap << EOF
90  -50
150 50
EOF
gmt end show

Here are the related codes:

gmt/src/gmt_vector.c

Lines 1607 to 1653 in b877d72

if (mode == GMT_STAIRS_Y) { /* First follow meridian, then parallel */
gmt_M_set_delta_lon (lon[i-1], lon[i], dlon); /* Beware of jumps due to sign differences */
lon_i = lon[i-1] + dlon; /* Use lon_i instead of lon[i] in the marching since this avoids any jumping */
theta = fabs (dlon) * cosd (lat[i-1]);
n_step = lrint (theta / step);
for (j = 1; j <= n_step; j++) {
c = j / (double)n_step;
gmt_prep_tmp_arrays (GMT, GMT_NOTSET, n_new, 2); /* Init or reallocate tmp vectors */
GMT->hidden.mem_coord[GMT_X][n_new] = lon[i-1] * (1 - c) + lon_i * c;
GMT->hidden.mem_coord[GMT_Y][n_new] = lat[i-1];
n_new++;
}
theta = fabs (lat[i]-lat[i-1]);
n_step = lrint (theta / step);
for (j = k; j < n_step; j++) { /* Start at 0 to make sure corner point is saved */
c = j / (double)n_step;
gmt_prep_tmp_arrays (GMT, GMT_NOTSET, n_new, 2); /* Init or reallocate tmp vectors */
GMT->hidden.mem_coord[GMT_X][n_new] = lon[i];
GMT->hidden.mem_coord[GMT_Y][n_new] = lat[i-1] * (1 - c) + lat[i] * c;
n_new++;
}
k = 0;
}
else if (mode == GMT_STAIRS_X) { /* First follow parallel, then meridian */
theta = fabs (lat[i]-lat[i-1]);
n_step = lrint (theta / step);
for (j = 1; j <= n_step; j++) {
c = j / (double)n_step;
gmt_prep_tmp_arrays (GMT, GMT_NOTSET, n_new, 2); /* Init or reallocate tmp read vectors */
GMT->hidden.mem_coord[GMT_X][n_new] = lon[i-1];
GMT->hidden.mem_coord[GMT_Y][n_new] = lat[i-1] * (1 - c) + lat[i] * c;
n_new++;
}
gmt_M_set_delta_lon (lon[i-1], lon[i], dlon); /* Beware of jumps due to sign differences */
lon_i = lon[i-1] + dlon; /* Use lon_i instead of lon[i] in the marching since this avoids any jumping */
theta = fabs (dlon) * cosd(lat[i]);
n_step = lrint (theta / step);
for (j = k; j < n_step; j++) { /* Start at 0 to make sure corner point is saved */
c = j / (double)n_step;
gmt_prep_tmp_arrays (GMT, GMT_NOTSET, n_new, 2); /* Init or reallocate tmp read vectors */
GMT->hidden.mem_coord[GMT_X][n_new] = lon[i-1] * (1 - c) + lon_i * c;
GMT->hidden.mem_coord[GMT_Y][n_new] = lat[i];
n_new++;
}
k = 0;
}

For the GMT_STAIR_Y case, it should follow Y (meridian) first, then X (parallel), so we should create the temporary array for the line with constant longitude first, then another line with consistent latitude. However, the source code is doing it in the wrong way. The fix is straightforward.

@joa-quim @Esteban82 Could you please verify if it's indeed a bug, since the codes have been here for more than ten years and nobody complains about it.

@Esteban82
Copy link
Member

Yes @seisman, it is a bug.

map

@joa-quim
Copy link
Member

joa-quim commented Jan 1, 2025

Indeed, it's a bug.

@seisman
Copy link
Member Author

seisman commented Jan 4, 2025

Test sample1d/sample.sh started to fail after #8648.

Here is the comparison between the baseline image and the new image.

Baseline Generated
sample sample

The green circles are the input data, and the first point is at (-5, 74). When -Am (blue lines) is used, the data should be resampled following meridian (Y) first, then parallel (X). Apparently, in the baseline image, the lines for -Am/-Ap are incorrect.

In the image generated by the current master branch, the geographic one (upper panel) looks correct, but the Cartesian one (lower panel) is still wrong. Am I understanding it correctly?

@joa-quim
Copy link
Member

joa-quim commented Jan 4, 2025

Hmm, aren't you swapping what's wrong/right? With master (confirmed with a local run). "Generated" geog is right. "BLUE LINE: -Am mode" is right on master but not on Baseline. The inverse ... hold on ...

Lower panels show different behavior. The wither, at the top, are wrong but the narrower at the the bottom are correct.

@seisman
Copy link
Member Author

seisman commented Jan 4, 2025

The "Baseline" one on the left is the current sample1d/sample.ps file, in which both panels are wrong.

After #8648, we get the "Generated" one on the right, in which the "geographic" one is correct but the Cartesian one is wrong.

With PR #8667, we will get the following image, which is correct for both "geographic" and "Cartesian"
image

@joa-quim
Copy link
Member

joa-quim commented Jan 4, 2025

I find the Cartesian example too complicate for me to understand what is right or wrong.

@remkos
Copy link
Contributor

remkos commented Jan 4, 2025

I agree that the cartesian example is a bit too complicated. But in any case with @seisman's #8667, both the cartesian and geographic versions are correct. Thus I approved #8667

@seisman seisman linked a pull request Jan 7, 2025 that will close this issue
@seisman seisman added this to the 6.6.0 milestone Jan 7, 2025
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
documentation Improve documentation
Projects
None yet
Development

Successfully merging a pull request may close this issue.

4 participants