Skip to content

Tackle the -px... -py... walls plot option #9047

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.

Sign up for GitHub

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

Open
joa-quim wants to merge 10 commits into
base: master
Choose a base branch
from
+108 −14
Open

Tackle the -px... -py... walls plot option #9047

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
10 commits
Select commit Hold shift + click to select a range
9ff015e
Try fixing again #529
joa-quim May 28, 2026
2741c2c
One more step towards the goal
joa-quim May 28, 2026
1faab34
Finally (?) fix the different scales issues
joa-quim May 28, 2026
8316a39
Fix the non-closing cune (+b) at the base.
joa-quim May 28, 2026
8fa3ab7
Do not skip Z-annotation label equal to plane level.
joa-quim May 28, 2026
e113c3d
Fix for last the Y-scales/annotations
joa-quim May 28, 2026
b8b89a5
Remove commented lines
joa-quim May 29, 2026
52d8aba
Add a test
joa-quim May 29, 2026
4b2520f
Planes test
joa-quim May 29, 2026
a402ba3
Add execut permission
seisman May 29, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
20 changes: 15 additions & 5 deletions src/gmt_map.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6367,6 +6367,7 @@ GMT_LOCAL int gmtmap_init_three_D (struct GMT_CTRL *GMT) {
unsigned int i, i_min = 0, i_max = 0;
bool easy, positive;
double x, y, zmin = 0.0, zmax = 0.0, z_range;
double _zmin = 0.0, _zmax = 0.0;

GMT->current.proj.three_D = (GMT->current.proj.z_project.view_azimuth != 180.0 || GMT->current.proj.z_project.view_elevation != 90.0);
GMT->current.proj.scale[GMT_Z] = GMT->current.proj.z_pars[0];
Expand All @@ -6389,6 +6390,9 @@ GMT_LOCAL int gmtmap_init_three_D (struct GMT_CTRL *GMT) {
case GMT_LINEAR: /* Regular scaling */
zmin = (GMT->current.proj.xyz_pos[GMT_Z]) ? GMT->common.R.wesn[i_min] : GMT->common.R.wesn[i_max];
zmax = (GMT->current.proj.xyz_pos[GMT_Z]) ? GMT->common.R.wesn[i_max] : GMT->common.R.wesn[i_min];
_zmin = (GMT->current.proj.xyz_pos[GMT_Z]) ? GMT->common.R.wesn[ZLO] : GMT->common.R.wesn[ZHI];
_zmax = (GMT->current.proj.xyz_pos[GMT_Z]) ? GMT->common.R.wesn[ZHI] : GMT->common.R.wesn[ZLO];

GMT->current.proj.fwd_z = &gmtlib_translin;
GMT->current.proj.inv_z = &gmtlib_itranslin;
break;
Expand All @@ -6397,21 +6401,27 @@ GMT_LOCAL int gmtmap_init_three_D (struct GMT_CTRL *GMT) {
GMT_Report (GMT->parent, GMT_MSG_ERROR, "Option -Jz -JZ: limits must be positive for log10 projection\n");
return GMT_PROJECTION_ERROR;
}
zmin = (GMT->current.proj.xyz_pos[GMT_Z]) ? d_log10 (GMT, GMT->common.R.wesn[i_min]) : d_log10 (GMT, GMT->common.R.wesn[i_max]);
zmax = (GMT->current.proj.xyz_pos[GMT_Z]) ? d_log10 (GMT, GMT->common.R.wesn[i_max]) : d_log10 (GMT, GMT->common.R.wesn[i_min]);
zmin = (GMT->current.proj.xyz_pos[GMT_Z]) ? d_log10(GMT, GMT->common.R.wesn[i_min]) : d_log10(GMT, GMT->common.R.wesn[i_max]);
zmax = (GMT->current.proj.xyz_pos[GMT_Z]) ? d_log10(GMT, GMT->common.R.wesn[i_max]) : d_log10(GMT, GMT->common.R.wesn[i_min]);
_zmin = (GMT->current.proj.xyz_pos[GMT_Z]) ? d_log10(GMT, GMT->common.R.wesn[ZLO]) : d_log10(GMT, GMT->common.R.wesn[ZHI]);
_zmax = (GMT->current.proj.xyz_pos[GMT_Z]) ? d_log10(GMT, GMT->common.R.wesn[ZHI]) : d_log10(GMT, GMT->common.R.wesn[ZLO]);

GMT->current.proj.fwd_z = &gmtproj_translog10;
GMT->current.proj.inv_z = &gmtproj_itranslog10;
break;
case GMT_POW: /* x^y transformation */
GMT->current.proj.xyz_pow[GMT_Z] = GMT->current.proj.z_pars[1];
GMT->current.proj.xyz_ipow[GMT_Z] = 1.0 / GMT->current.proj.z_pars[1];
positive = !((GMT->current.proj.xyz_pos[GMT_Z] + (GMT->current.proj.xyz_pow[GMT_Z] > 0.0)) % 2);
zmin = (positive) ? pow (GMT->common.R.wesn[ZLO], GMT->current.proj.xyz_pow[GMT_Z]) : pow (GMT->common.R.wesn[i_max], GMT->current.proj.xyz_pow[GMT_Z]);
zmax = (positive) ? pow (GMT->common.R.wesn[i_max], GMT->current.proj.xyz_pow[GMT_Z]) : pow (GMT->common.R.wesn[i_min], GMT->current.proj.xyz_pow[GMT_Z]);
zmin = (positive) ? pow(GMT->common.R.wesn[i_min], GMT->current.proj.xyz_pow[GMT_Z]) : pow(GMT->common.R.wesn[i_max], GMT->current.proj.xyz_pow[GMT_Z]);
zmax = (positive) ? pow(GMT->common.R.wesn[i_max], GMT->current.proj.xyz_pow[GMT_Z]) : pow(GMT->common.R.wesn[i_min], GMT->current.proj.xyz_pow[GMT_Z]);
_zmin = (positive) ? pow(GMT->common.R.wesn[ZLO], GMT->current.proj.xyz_pow[GMT_Z]) : pow(GMT->common.R.wesn[ZHI], GMT->current.proj.xyz_pow[GMT_Z]);
_zmax = (positive) ? pow(GMT->common.R.wesn[ZHI], GMT->current.proj.xyz_pow[GMT_Z]) : pow(GMT->common.R.wesn[ZLO], GMT->current.proj.xyz_pow[GMT_Z]);

GMT->current.proj.fwd_z = &gmtproj_transpowz;
GMT->current.proj.inv_z = &gmtproj_itranspowz;
}
z_range = zmax - zmin;
z_range = _zmax - _zmin; /* This one is for the Z-axis scale. zmin & zmax are for -p[x|y] levels */
if (z_range == 0.0 && GMT->current.proj.compute_scale[GMT_Z])
GMT->current.proj.scale[GMT_Z] = 0.0; /* No range given, just flat projected map */
else if (GMT->current.proj.compute_scale[GMT_Z])
Expand Down
1 change: 1 addition & 0 deletions src/gmt_parse.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -98,6 +98,7 @@ GMT_LOCAL int gmtparse_B_arg_inspector (struct GMT_CTRL *GMT, char *in) {
k = (in[0] == 'p' || in[0] == 's') ? 1 : 0; /* Skip p|s in -Bp|s */
if (strchr ("xyz", in[k])) gmt5++; /* Definitively GMT5 */
if (k == 0 && !isdigit (in[0]) && strchr ("WESNwesn", in[1])) gmt5++; /* Definitively GMT5 */
if (in[0] == 's' && in[1] && strchr ("WESNZwenzlrbtu", in[1])) gmt5++; /* -Bs<axis-letter> is GMT5 frame (south + other axes), not -B[s] secondary prefix */
j = k;
while (j < last && (in[j] == 'x' || in[j] == 'y' || in[j] == 'z')) j++;
custom = (in[j] == 'c'); /* Got -B[p|s][xyz]c<customfile> */
Expand Down
82 changes: 76 additions & 6 deletions src/gmt_plot.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3186,7 +3186,7 @@ GMT_LOCAL void gmtplot_cube_box (struct GMT_CTRL *GMT, struct PSL_CTRL *PSL, int
xx[0] = xx[1] = nesw[(quadrant+1)%4]; xx[2] = xx[3] = nesw[(quadrant+3)%4];
yy[0] = yy[3] = GMT->current.proj.zmin; yy[1] = yy[2] = GMT->current.proj.zmax;
gmt_setpen (GMT, &GMT->current.map.frame.pen);
PSL_plotline (PSL, xx, yy, 4, PSL_MOVE|PSL_STROKE);
PSL_plotline (PSL, xx, yy, 4, PSL_MOVE|PSL_STROKE|PSL_CLOSE);
}

GMT_LOCAL void gmtplot_timestamp(struct GMT_CTRL *GMT, struct PSL_CTRL *PSL, double x, double y, unsigned int justify, char *U_label) {
Expand Down Expand Up @@ -5909,10 +5909,10 @@ void gmt_xy_axis (struct GMT_CTRL *GMT, double x0, double y0, double length, dou

for (i = 0; i < nx1; i++) {
if (gmtlib_annot_pos (GMT, val0, val1, T, &knots[i], &t_use)) continue; /* Outside range */
if (axis == GMT_Z && fabs (knots[i] - GMT->current.proj.z_level) < GMT_CONV8_LIMIT) continue; /* Skip z annotation coinciding with z-level plane */
if (axis == GMT_Z && (GMT->current.proj.z_project.view_plane % 3) == GMT_Z && fabs (knots[i] - GMT->current.proj.z_level) < GMT_CONV8_LIMIT) continue; /* Skip z annotation coinciding with z-level plane (only for -pz) */
x = (*xyz_fwd) (GMT, knots[i]); /* Convert to inches on the page */
if (skip_center_annot && doubleAlmostEqualZero (knots[i], skip_val)) continue; /* Don't want annotations exactly at intersection between graph axes */
if (gmtplot_skip_end_annotation (GMT, axis, x, length)) continue; /* Don't want annotations exactly at one or both ends of the axis */
if (skip_center_annot && doubleAlmostEqualZero(knots[i], skip_val)) continue; /* Don't want annotations exactly at intersection between graph axes */
if (gmtplot_skip_end_annotation(GMT, axis, x, length)) continue; /* Don't want annotations exactly at one or both ends of the axis */
if (!is_interval && gmtplot_skip_second_annot (k, knots[i], knots_p, np, primary)) continue; /* Secondary annotation skipped when coinciding with primary annotation */
if (label_c && label_c[i] && label_c[i][0]) {
strncpy (string, label_c[i], GMT_LEN256-1);
Expand Down Expand Up @@ -5949,7 +5949,7 @@ void gmt_xy_axis (struct GMT_CTRL *GMT, double x0, double y0, double length, dou

for (i = 0; i < nx1; i++) {
if (gmtlib_annot_pos (GMT, val0, val1, T, &knots[i], &t_use)) continue; /* Outside range */
if (axis == GMT_Z && fabs (knots[i] - GMT->current.proj.z_level) < GMT_CONV8_LIMIT) continue; /* Skip z annotation coinciding with z-level plane */
if (axis == GMT_Z && (GMT->current.proj.z_project.view_plane % 3) == GMT_Z && fabs(knots[i] - GMT->current.proj.z_level) < GMT_CONV8_LIMIT) continue; /* Skip z annotation coinciding with z-level plane (only for -pz) */
x = (*xyz_fwd) (GMT, t_use); /* Convert to inches on the page */
if (skip_center_annot && doubleAlmostEqualZero (knots[i], skip_val)) continue; /* Don't want annotations exactly at intersection between graph axes */
if (gmtplot_skip_end_annotation (GMT, axis, x, length)) continue; /* Don't want annotations exactly at one or both ends of the axis */
Expand Down Expand Up @@ -6648,6 +6648,48 @@ void gmt_map_basemap (struct GMT_CTRL *GMT) {
bool clip_on = false, too_crowded = false;
char *order[2] = {"before", "after"};
struct PSL_CTRL *PSL= GMT->PSL;
/* For -px/-py world-coord wall plan, vertical of basemap is the Z axis. Swap only Y annotation state with Z
* (labels and range), scale Y so labels still land within plan's projected y-box [0,map.height]; matrix d-stretch
* then maps that visually to Z height. Restore before gmt_vertical_axis. */
bool swap_yz = (GMT->current.proj.z_project.plane >= GMT_ZW && (GMT->current.proj.z_project.plane % 3) != GMT_Z);
bool swap_xy = (GMT->current.proj.z_project.plane >= GMT_ZW && (GMT->current.proj.z_project.plane % 3) == GMT_X);
double saved_y_wesn[2] = {0.0, 0.0}, saved_y_scale = 0.0, saved_y_origin = 0.0;
double saved_x_wesn[2] = {0.0, 0.0}, saved_x_scale = 0.0, saved_x_origin = 0.0;
double saved_rect_x[2] = {0.0, 0.0}, saved_map_width = 0.0;
struct GMT_PLOT_AXIS saved_axis_y, saved_axis_x;
if (swap_yz) {
double z_range = GMT->common.R.wesn[ZHI] - GMT->common.R.wesn[ZLO];
saved_y_wesn[0] = GMT->common.R.wesn[YLO];
saved_y_wesn[1] = GMT->common.R.wesn[YHI];
saved_y_scale = GMT->current.proj.scale[GMT_Y];
saved_y_origin = GMT->current.proj.origin[GMT_Y];
saved_axis_y = GMT->current.map.frame.axis[GMT_Y];
GMT->common.R.wesn[YLO] = GMT->common.R.wesn[ZLO];
GMT->common.R.wesn[YHI] = GMT->common.R.wesn[ZHI];
GMT->current.proj.scale[GMT_Y] = (z_range > 0.0) ? GMT->current.map.height / z_range : saved_y_scale;
GMT->current.proj.origin[GMT_Y] = -GMT->common.R.wesn[ZLO] * GMT->current.proj.scale[GMT_Y];
GMT->current.map.frame.axis[GMT_Y] = GMT->current.map.frame.axis[GMT_Z];
GMT->current.map.frame.axis[GMT_Y].id = GMT_Y; /* So gmt_xy_axis uses gmt_y_to_yy(with our swapped scale[Y]) not gmt_z_to_zz */
}
if (swap_xy) { /* -px: plan horizontal is the Y direction; swap X annotation state with original Y so plan top/bottom show Y range */
saved_x_wesn[0] = GMT->common.R.wesn[XLO];
saved_x_wesn[1] = GMT->common.R.wesn[XHI];
saved_x_scale = GMT->current.proj.scale[GMT_X];
saved_x_origin = GMT->current.proj.origin[GMT_X];
saved_axis_x = GMT->current.map.frame.axis[GMT_X];
saved_rect_x[0] = GMT->current.proj.rect[XLO];
saved_rect_x[1] = GMT->current.proj.rect[XHI];
saved_map_width = GMT->current.map.width;
GMT->common.R.wesn[XLO] = saved_y_wesn[0];
GMT->common.R.wesn[XHI] = saved_y_wesn[1];
GMT->current.proj.scale[GMT_X] = saved_y_scale;
GMT->current.proj.origin[GMT_X] = saved_y_origin;
GMT->current.map.frame.axis[GMT_X] = saved_axis_y;
GMT->current.map.frame.axis[GMT_X].id = GMT_X;
GMT->current.proj.rect[XLO] = GMT->current.proj.rect[YLO];
GMT->current.proj.rect[XHI] = GMT->current.proj.rect[YHI];
GMT->current.map.width = GMT->current.map.height;
}

/* 0. Determine if we need to be here and set a few parameters */

Expand Down Expand Up @@ -6712,6 +6754,24 @@ void gmt_map_basemap (struct GMT_CTRL *GMT) {

if (GMT->current.proj.got_azimuths) gmt_M_uint_swap (GMT->current.map.frame.side[E_SIDE], GMT->current.map.frame.side[W_SIDE]); /* Undo temporary swap */

if (swap_xy) { /* Restore X-axis state we swapped with Y for the plan horizontal */
GMT->common.R.wesn[XLO] = saved_x_wesn[0];
GMT->common.R.wesn[XHI] = saved_x_wesn[1];
GMT->current.proj.scale[GMT_X] = saved_x_scale;
GMT->current.proj.origin[GMT_X] = saved_x_origin;
GMT->current.map.frame.axis[GMT_X] = saved_axis_x;
GMT->current.proj.rect[XLO] = saved_rect_x[0];
GMT->current.proj.rect[XHI] = saved_rect_x[1];
GMT->current.map.width = saved_map_width;
}
if (swap_yz) { /* Restore Y-axis state we swapped with Z for the vertical wall plan */
GMT->common.R.wesn[YLO] = saved_y_wesn[0];
GMT->common.R.wesn[YHI] = saved_y_wesn[1];
GMT->current.proj.scale[GMT_Y] = saved_y_scale;
GMT->current.proj.origin[GMT_Y] = saved_y_origin;
GMT->current.map.frame.axis[GMT_Y] = saved_axis_y;
}

if (GMT->current.proj.three_D && GMT->current.map.frame.drawz) GMT->current.map.frame.plotted_header = false; /* Now we can plot the title [if selected via -B+t] */

if ((GMT->current.map.frame.order == GMT_BASEMAP_BEFORE && (GMT->current.plot.mode_3D & 1)) || (GMT->current.map.frame.order == GMT_BASEMAP_AFTER && (GMT->current.plot.mode_3D & 2)))
Expand Down Expand Up @@ -10582,13 +10642,21 @@ void gmt_plane_perspective (struct GMT_CTRL *GMT, int plane, double level) {
if (plane < 0) /* Reset to original matrix */
PSL_command (PSL, "PSL_GPP setmatrix\n"); /* Return to normal current transformation matrix */
else { /* New perspective plane: compute all derivatives and use full matrix */
if (plane >= GMT_ZW) level = gmt_z_to_zz (GMT, level); /* First convert world z coordinate to projected z coordinate */
if (plane >= GMT_ZW) { /* World coordinate: pick the axis perpendicular to the plane so the correct scale (X/Y/Z) is applied */
switch (plane % 3) {
case GMT_X: level = gmt_x_to_xx(GMT, level); break; /* Constant-X plane: level is in user X units */
case GMT_Y: level = gmt_y_to_yy(GMT, level); break; /* Constant-Y plane: level is in user Y units */
case GMT_Z: level = gmt_z_to_zz(GMT, level); break; /* Constant-Z plane: level is in user Z units */
}
}
switch (plane % 3) {
case GMT_X: /* Constant x, Convert y,z to x',y' */
a = GMT->current.proj.z_project.sin_az;
b = -GMT->current.proj.z_project.cos_az * GMT->current.proj.z_project.sin_el;
c = 0.0;
d = GMT->current.proj.z_project.cos_el;
if (plane >= GMT_ZW && GMT->current.map.height > 0.0) /* User -px plan: stretch in_y so plan height matches Z axis height (proj.zmax) */
d *= GMT->current.proj.zmax / GMT->current.map.height;
e = GMT->current.proj.z_project.x_off - level * GMT->current.proj.z_project.cos_az;
f = GMT->current.proj.z_project.y_off - level * GMT->current.proj.z_project.sin_az * GMT->current.proj.z_project.sin_el;
break;
Expand All @@ -10597,6 +10665,8 @@ void gmt_plane_perspective (struct GMT_CTRL *GMT, int plane, double level) {
b = -GMT->current.proj.z_project.sin_az * GMT->current.proj.z_project.sin_el;
c = 0.0;
d = GMT->current.proj.z_project.cos_el;
if (plane >= GMT_ZW && GMT->current.map.height > 0.0) /* User -py plan: stretch in_y so plan height matches Z axis height (proj.zmax) */
d *= GMT->current.proj.zmax / GMT->current.map.height;
e = GMT->current.proj.z_project.x_off + level * GMT->current.proj.z_project.sin_az;
f = GMT->current.proj.z_project.y_off - level * GMT->current.proj.z_project.cos_az * GMT->current.proj.z_project.sin_el;
break;
Expand Down
6 changes: 3 additions & 3 deletions test/baseline/psbasemap.dvc
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
outs:
- md5: 6f6177c3673dc054c5489c208a9d072a.dir
nfiles: 57
- md5: 499cd79959131cb42466afb2ba0bda90.dir
nfiles: 58
path: psbasemap
hash: md5
size: 2936773
size: 2981097
13 changes: 13 additions & 0 deletions test/psbasemap/perspective_planes.sh
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,13 @@
#!/usr/bin/env bash
#
# Test 3D perspective basemap with -px, -py, -pz view planes side by side.
# Each cube uses -BwsENZ1+b to draw all 12 cube edges. Y range 0-20, X range 0-10, Z range 0-30.
# -px: looking along X (plan = YZ); -py: looking along Y (plan = XZ); -pz: looking along Z (plan = XY).

ps=perspective_planes.ps

common="-R0/10/0/20/0/30 -JX2.5c/5c -JZ7.5c -Bxa+lx -Bya+ly -Bza+lz -BwsENZ1+b --PS_MEDIA=A4"

gmt psbasemap $common -px135/40/5 -P -K -Xf1c -Yf2c > $ps
gmt psbasemap $common -py135/40/10 -O -K -Xf8c >> $ps
gmt psbasemap $common -pz135/40/20 -O -Xf14c >> $ps
Loading