From 369ac9402c8414a09c89637506cd8f2af3bc0738 Mon Sep 17 00:00:00 2001 From: Joaquim Date: Wed, 15 Apr 2026 22:04:34 +0100 Subject: [PATCH 1/2] pscoast: fix offset after a -p90 MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit 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/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. --- src/gmt_map.c | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/src/gmt_map.c b/src/gmt_map.c index 2c59f2f93e7..6856a87f632 100644 --- a/src/gmt_map.c +++ b/src/gmt_map.c @@ -6524,6 +6524,22 @@ GMT_LOCAL int gmtmap_init_three_D (struct GMT_CTRL *GMT) { GMT->current.proj.z_project.x_off = GMT->current.proj.z_project.view_x - x; GMT->current.proj.z_project.y_off = GMT->current.proj.z_project.view_y - y; } + else if (doubleAlmostEqualZero(GMT->current.proj.z_project.view_elevation, 90.0)) { + /* Pure z-rotation (elevation 90 = looking straight down). Place the map centre + * at (map.width/2, map.height/2) on page, so that the rotated disc/map centre + * stays at the same location a non-rotated plot would produce. We cannot just + * centre the bbox because perimeter sampling inflates it asymmetrically. + * See issue #8794. */ + 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); /* 2-D projected coords of map centre */ + /* Apply the rotation part of gmt_xyz_to_xy (x_off=y_off=0) to get post-rotation position of 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; + } else { GMT->current.proj.z_project.x_off = -GMT->current.proj.z_project.xmin; GMT->current.proj.z_project.y_off = -GMT->current.proj.z_project.ymin; From 67ffe26013c32093228540ff8aab9ae3de8dd54e Mon Sep 17 00:00:00 2001 From: Joaquim Date: Thu, 16 Apr 2026 00:27:48 +0100 Subject: [PATCH 2/2] Apply this pathch only if -p has been used. This at least prevents breaking other tests. --- src/gmt_map.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/gmt_map.c b/src/gmt_map.c index 6856a87f632..63c29802cd6 100644 --- a/src/gmt_map.c +++ b/src/gmt_map.c @@ -6524,7 +6524,7 @@ GMT_LOCAL int gmtmap_init_three_D (struct GMT_CTRL *GMT) { GMT->current.proj.z_project.x_off = GMT->current.proj.z_project.view_x - x; GMT->current.proj.z_project.y_off = GMT->current.proj.z_project.view_y - y; } - else if (doubleAlmostEqualZero(GMT->current.proj.z_project.view_elevation, 90.0)) { + else if (GMT->common.p.active && doubleAlmostEqualZero(GMT->current.proj.z_project.view_elevation, 90.0)) { /* Pure z-rotation (elevation 90 = looking straight down). Place the map centre * at (map.width/2, map.height/2) on page, so that the rotated disc/map centre * stays at the same location a non-rotated plot would produce. We cannot just