Skip to content

pscoast: fix offset after a -p90#8981

Merged
joa-quim merged 2 commits intomasterfrom
fix-8794
Apr 16, 2026
Merged

pscoast: fix offset after a -p90#8981
joa-quim merged 2 commits intomasterfrom
fix-8794

Conversation

@joa-quim
Copy link
Copy Markdown
Member

@joa-quim joa-quim commented Apr 15, 2026

Assisted-by: Claude Opus 4.6

Fix #8794

Claude explanations that may be useful in future related issues:

Problem

When overlaying a rotated plot on top of a non-rotated one using -p/90 (with no +v/+w modifiers), the two plots end up misaligned by a
small diagonal offset, even though elevation = 90 is mathematically equivalent to a pure z-rotation that should preserve the map centre.

Reproducer from the issue:
gmt pscoast -Rg -JF180/0/60/3 -Baf -BWSen -Ggray -Sblue -Da -P -K > lixo.ps
gmt pscoast -Rg -JF180/0/60/3 -B -X3 -p90/90 -Ggray -Sblue -Da -O >> lixo.ps

With -X3 equal to exactly one disc diameter, the two circles were expected to be tangent at the equator, but instead showed a few-millimetre
offset in both X and Y.

Root cause

In gmtmap_init_three_D (src/gmt_map.c), when z_project.fixed is false (i.e. no explicit +v/+w pivot), x_off/y_off are computed as -xmin/-ymin,
pinning the rotated bounding box to (0,0) of the local coordinate system.

The bounding box is built by sampling points along the lon/lat perimeter and transforming each through the 3D projection, so it is inflated
asymmetrically relative to the disc/map itself. Pinning the inflated bbox corner to (0,0) therefore shifts the map centre away from where a
non-rotated plot would place it — the source of the diagonal misalignment on overlays.

Fix

Add a new branch between fixed and the existing else that handles the pure-rotation case (view_elevation == 90, no +v/+w). Instead of centring
the (asymmetric) bbox, place the map's geographic centre point (central meridian, mid-latitude) at (map.width/2, map.height/2) on the page.
This guarantees the disc/map centre ends up at exactly the same page location as the non-rotated plot, so -X overlays align as expected.

Implementation:

  else if (doubleAlmostEqualZero (GMT->current.proj.z_project.view_elevation, 90.0)) {
      double wx, wy, xc, yc;
      double world_x = (gmt_M_is_geographic (GMT, GMT_IN)) ? GMT->current.proj.central_meridian : 0.5 * (GMT->common.R.wesn[XLO] +
  GMT->common.R.wesn[XHI]);
      double world_y = 0.5 * (GMT->common.R.wesn[YLO] + GMT->common.R.wesn[YHI]);
      gmt_geo_to_xy (GMT, world_x, world_y, &wx, &wy);
      /* Apply the rotational part of gmt_xyz_to_xy (with x_off=y_off=0) to obtain the post-rotation position of the centre */
      xc = -wx * GMT->current.proj.z_project.cos_az + wy * GMT->current.proj.z_project.sin_az;
      yc = -(wx * GMT->current.proj.z_project.sin_az + wy * GMT->current.proj.z_project.cos_az) * GMT->current.proj.z_project.sin_el;
      GMT->current.proj.z_project.x_off = 0.5 * GMT->current.map.width  - xc;
      GMT->current.proj.z_project.y_off = 0.5 * GMT->current.map.height - yc;
  }

The formula reuses the rotational part of gmt_xyz_to_xy (gmt_map.c:6357-6358); cos_az/sin_az/sin_el are already populated at line 6341, well
before this point.

Scope

  • Only triggered when elevation == 90 and no +v/+w — i.e. the user asked for what is semantically a pure rotation.
  • Genuine 3-D perspective plots (el < 90), explicit pivots (+v/+w), and -p with no elevation (handled via do_z_rotation) all follow the original code paths unchanged.

Assisted-by: Claude Opus 4.6

Claude explanations that may be useful in future related issues:

  Problem

  When overlaying a rotated plot on top of a non-rotated one using -p<az>/90 (with no +v/+w modifiers), the two plots end up misaligned by a
  small diagonal offset, even though elevation = 90 is mathematically equivalent to a pure z-rotation that should preserve the map centre.

  Reproducer from the issue:
  gmt pscoast -Rg -JF180/0/60/3 -Baf -BWSen -Ggray -Sblue -Da -P -K > lixo.ps
  gmt pscoast -Rg -JF180/0/60/3 -B -X3 -p90/90 -Ggray -Sblue -Da -O >> lixo.ps

  With -X3 equal to exactly one disc diameter, the two circles were expected to be tangent at the equator, but instead showed a few-millimetre
  offset in both X and Y.

  Root cause

  In gmtmap_init_three_D (src/gmt_map.c), when z_project.fixed is false (i.e. no explicit +v/+w pivot), x_off/y_off are computed as -xmin/-ymin,
  pinning the rotated bounding box to (0,0) of the local coordinate system.

  The bounding box is built by sampling points along the lon/lat perimeter and transforming each through the 3D projection, so it is inflated
  asymmetrically relative to the disc/map itself. Pinning the inflated bbox corner to (0,0) therefore shifts the map centre away from where a
  non-rotated plot would place it — the source of the diagonal misalignment on overlays.

  Fix

  Add a new branch between fixed and the existing else that handles the pure-rotation case (view_elevation == 90, no +v/+w). Instead of centring
  the (asymmetric) bbox, place the map's geographic centre point (central meridian, mid-latitude) at (map.width/2, map.height/2) on the page.
  This guarantees the disc/map centre ends up at exactly the same page location as the non-rotated plot, so -X overlays align as expected.

  Implementation:

  else if (doubleAlmostEqualZero (GMT->current.proj.z_project.view_elevation, 90.0)) {
      double wx, wy, xc, yc;
      double world_x = (gmt_M_is_geographic (GMT, GMT_IN)) ? GMT->current.proj.central_meridian : 0.5 * (GMT->common.R.wesn[XLO] +
  GMT->common.R.wesn[XHI]);
      double world_y = 0.5 * (GMT->common.R.wesn[YLO] + GMT->common.R.wesn[YHI]);
      gmt_geo_to_xy (GMT, world_x, world_y, &wx, &wy);
      /* Apply the rotational part of gmt_xyz_to_xy (with x_off=y_off=0) to obtain the post-rotation position of the centre */
      xc = -wx * GMT->current.proj.z_project.cos_az + wy * GMT->current.proj.z_project.sin_az;
      yc = -(wx * GMT->current.proj.z_project.sin_az + wy * GMT->current.proj.z_project.cos_az) * GMT->current.proj.z_project.sin_el;
      GMT->current.proj.z_project.x_off = 0.5 * GMT->current.map.width  - xc;
      GMT->current.proj.z_project.y_off = 0.5 * GMT->current.map.height - yc;
  }

  The formula reuses the rotational part of gmt_xyz_to_xy (gmt_map.c:6357-6358); cos_az/sin_az/sin_el are already populated at line 6341, well
  before this point.

  Scope

  - Only triggered when elevation == 90 and no +v/+w — i.e. the user asked for what is semantically a pure rotation.
  - Genuine 3-D perspective plots (el < 90), explicit pivots (+v/+w), and -p<az> with no elevation (handled via do_z_rotation) all follow the
  original code paths unchanged.
@joa-quim joa-quim requested a review from a team April 15, 2026 21:06
Copy link
Copy Markdown
Member

@Esteban82 Esteban82 left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It works.

@Esteban82 Esteban82 added the add-changelog Add PR to the changelog label Apr 15, 2026
@joa-quim
Copy link
Copy Markdown
Member Author

But breaks several other tests.

@joa-quim joa-quim merged commit a5f8a31 into master Apr 16, 2026
10 of 13 checks passed
@joa-quim joa-quim deleted the fix-8794 branch April 16, 2026 01:35
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

add-changelog Add PR to the changelog

Projects

None yet

Development

Successfully merging this pull request may close these issues.

Two images of same size offset with -X do not align. -p bug?

2 participants