From 939d9600f86d5f515c08052c3264fd184c33e610 Mon Sep 17 00:00:00 2001 From: Joaquim Date: Sun, 30 Aug 2026 17:09:15 +0100 Subject: [PATCH] Fix thin polar sectors that end up parallel to an axis (#7059) A single polar axis is a sector so thin that its bounding box has a width but no height. Rotating it by 90 degrees swaps the two and broke both the scale and the annotations: * gmtmap_setinfo matched -JP to the map width, which for a sector standing on end is nothing but round-off, so -JP12c+t90 produced a plot 12 x 3437746 cm. It now checks the dimension it was asked to match and falls back to the other one when the first is under 0.01 % of it, with a warning (and a pointer to -Jp for polar plots). A fully degenerate map no longer divides by zero. * n_lon_nodes/n_lat_nodes are derived from the map dimensions and rounded down to zero for a map with no width, so the boundary-walking loops in gmtlib_map_latcross/gmtlib_map_loncross never ran and no ticks or annotations were placed. The three places that derive them now floor them at 2. Adds test/psbasemap/polar_thin_sector.sh and documents the meaning of the -JP size for sectors, plus the -Jp recipe for radar plots. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01TZPAjvoG9oCgbMMbenuV8c --- .../reference/coordinate-transformations.rst | 18 +++++++ src/gmt_map.c | 38 +++++++++----- test/psbasemap/polar_thin_sector.sh | 51 +++++++++++++++++++ 3 files changed, 93 insertions(+), 14 deletions(-) create mode 100644 test/psbasemap/polar_thin_sector.sh diff --git a/doc/rst/source/reference/coordinate-transformations.rst b/doc/rst/source/reference/coordinate-transformations.rst index f54e12988cd..0beb703ca47 100644 --- a/doc/rst/source/reference/coordinate-transformations.rst +++ b/doc/rst/source/reference/coordinate-transformations.rst @@ -312,6 +312,24 @@ in the plane it is a **two-dimensional** projection. The transformation comes in #. The radius *r* can either be radius or inverted to mean depth from the surface, planetary radii, or even elevations in degrees. +**NOTE**: With **-JP** the given dimension is the dimension of the *plot*, not the radius: a full circle is +2\ *r* across, so **-JP**\ 12c draws a circle 12 cm wide, while the quarter sector **-R**\ 0/90/0/1 **-JP**\ 12c +is also 12 cm wide but has a 12 cm radius. Since the bounding box of a sector changes when the sector is rotated +with **+t**, the same sector given different **+t** values will come out at different radial scales, and a thin +sector that ends up parallel to the *y*-axis has no width for the given dimension to set (GMT then sizes the height +instead and warns). Furthermore, each plot is placed with the lower left corner of its bounding box at the current +plot origin, so overlays that differ in **+t** line up by their bounding boxes and not by their *r* = 0 point. Use +**-Jp**\ *scale* when several polar plots must share the same radial scale, e.g., when building a radar (spider) +plot from several rotated axes: + + :: + + gmt begin radar + gmt basemap -R-0.0001/0.0001/0/1 -Jp12c -Bya + gmt basemap -Jp12c+t-45 -Bya + gmt basemap -Jp12c+t-90 -Bya + gmt end show + **Example** As an example of this projection we will create a gridded data set in polar coordinates diff --git a/src/gmt_map.c b/src/gmt_map.c index 6145f68f8d5..09e75f125d0 100644 --- a/src/gmt_map.c +++ b/src/gmt_map.c @@ -2423,14 +2423,21 @@ GMT_LOCAL void gmtmap_setinfo (struct GMT_CTRL *GMT, double xmin, double xmax, d w = (xmax - xmin) * GMT->current.proj.scale[GMT_X]; h = (ymax - ymin) * GMT->current.proj.scale[GMT_Y]; - if (GMT->current.proj.gave_map_width == 1) /* Must rescale to given width */ - factor = scl / w; - else if (GMT->current.proj.gave_map_width == 2) /* Must rescale to given height */ - factor = scl / h; - else if (GMT->current.proj.gave_map_width == 3) /* Must rescale to max dimension */ - factor = scl / MAX (w, h); - else if (GMT->current.proj.gave_map_width == 4) /* Must rescale to min dimension */ - factor = scl / MIN (w, h); + if (GMT->current.proj.gave_map_width) { /* Gave a map dimension (1 = width, 2 = height, 3 = max, 4 = min) that the plot must match */ + static char *kind[5] = {"", "width", "height", "maximum dimension", "minimum dimension"}; + unsigned int type = GMT->current.proj.gave_map_width; + double big = MAX(w, h), dim = (type == 1) ? w : ((type == 2) ? h : ((type == 3) ? big : MIN(w, h))); + if (big > 0.0 && dim < GMT_CONV4_LIMIT * big) { + /* The dimension we were asked to match is nothing but round-off, e.g., a thin polar sector rotated to + * lie along the other axis. Matching it would blow the scale up to absurdity, so use the other one */ + GMT_Report(GMT->parent, GMT_MSG_WARNING, "Your map has no %s to speak of (%g %% of its other dimension) so the given size cannot set it; sizing the other dimension instead.\n", + kind[type], 100.0 * dim / big); + if (GMT->current.proj.projection == GMT_POLAR) + GMT_Report(GMT->parent, GMT_MSG_WARNING, "For polar plots, -Jp sets the radial scale directly and is not affected by the +t rotation.\n"); + dim = big; + } + if (dim > 0.0) factor = scl / dim; + } GMT->current.proj.scale[GMT_X] *= factor; GMT->current.proj.scale[GMT_Y] *= factor; GMT->current.proj.w_r *= factor; @@ -5383,8 +5390,9 @@ void gmt_wesn_search (struct GMT_CTRL *GMT, double xmin, double xmax, double ymi /* Search for extreme lon/lat coordinates by matching along the rectangular boundary */ - if (!GMT->current.map.n_lon_nodes) GMT->current.map.n_lon_nodes = urint (GMT->current.map.width / GMT->current.setting.map_line_step); - if (!GMT->current.map.n_lat_nodes) GMT->current.map.n_lat_nodes = urint (GMT->current.map.height / GMT->current.setting.map_line_step); + /* Need at least a couple of nodes along each side, even if the map has no width or height to speak of */ + if (!GMT->current.map.n_lon_nodes) GMT->current.map.n_lon_nodes = MAX(2, urint(GMT->current.map.width / GMT->current.setting.map_line_step)); + if (!GMT->current.map.n_lat_nodes) GMT->current.map.n_lat_nodes = MAX(2, urint(GMT->current.map.height / GMT->current.setting.map_line_step)); if (GMT->current.map.width > 400.0 && gmt_M_is_grdmapproject (GMT)) { /* ***project calling with true scale, probably. Reset to sane values */ GMT->current.map.n_lon_nodes = MIN (GMT->current.map.n_lon_nodes, 360); @@ -10126,8 +10134,9 @@ int gmt_proj_setup (struct GMT_CTRL *GMT, double wesn[]) { if (GMT->current.proj.central_meridian > GMT->common.R.wesn[XHI] && (GMT->current.proj.central_meridian - 360.0) >= GMT->common.R.wesn[XLO]) GMT->current.proj.central_meridian -= 360.0; } - if (!GMT->current.map.n_lon_nodes) GMT->current.map.n_lon_nodes = urint (GMT->current.map.width / GMT->current.setting.map_line_step); - if (!GMT->current.map.n_lat_nodes) GMT->current.map.n_lat_nodes = urint (GMT->current.map.height / GMT->current.setting.map_line_step); + /* Need at least a couple of nodes along each side, even if the map has no width or height to speak of */ + if (!GMT->current.map.n_lon_nodes) GMT->current.map.n_lon_nodes = MAX(2, urint(GMT->current.map.width / GMT->current.setting.map_line_step)); + if (!GMT->current.map.n_lat_nodes) GMT->current.map.n_lat_nodes = MAX(2, urint(GMT->current.map.height / GMT->current.setting.map_line_step)); error = gmtmap_init_three_D (GMT); @@ -10204,8 +10213,9 @@ int gmt_map_setup (struct GMT_CTRL *GMT, double wesn[]) { MAX (GMT->current.map.frame.axis[GMT_X].item[i].interval, GMT->current.map.frame.axis[GMT_Y].item[i].interval); } - if (!GMT->current.map.n_lon_nodes) GMT->current.map.n_lon_nodes = urint (GMT->current.map.width / GMT->current.setting.map_line_step); - if (!GMT->current.map.n_lat_nodes) GMT->current.map.n_lat_nodes = urint (GMT->current.map.height / GMT->current.setting.map_line_step); + /* Need at least a couple of nodes along each side, even if the map has no width or height to speak of */ + if (!GMT->current.map.n_lon_nodes) GMT->current.map.n_lon_nodes = MAX(2, urint(GMT->current.map.width / GMT->current.setting.map_line_step)); + if (!GMT->current.map.n_lat_nodes) GMT->current.map.n_lat_nodes = MAX(2, urint(GMT->current.map.height / GMT->current.setting.map_line_step)); GMT->current.map.dlon = (GMT->common.R.wesn[XHI] - GMT->common.R.wesn[XLO]) / GMT->current.map.n_lon_nodes; GMT->current.map.dlat = (GMT->common.R.wesn[YHI] - GMT->common.R.wesn[YLO]) / GMT->current.map.n_lat_nodes; diff --git a/test/psbasemap/polar_thin_sector.sh b/test/psbasemap/polar_thin_sector.sh new file mode 100644 index 00000000000..b75ca24a89e --- /dev/null +++ b/test/psbasemap/polar_thin_sector.sh @@ -0,0 +1,51 @@ +#!/usr/bin/env bash +# +# Test thin polar (-Jp|P) sectors, such as the single axes one would combine into a radar plot. +# A sector that ends up parallel to the y-axis used to blow the -JP scale up to absurdity and +# to lose all of its annotations. See https://github.com/GenericMappingTools/gmt/issues/7059 + +R=-R-0.0001/0.0001/0/1 + +# 1. The plot length, in cm, of the radial axis r = 0 to 1. +# Note: with -JP the size is the width of the plot, so rotating a thin sector by 45 degrees +# leaves it 12 cm wide but 12/cos(45) = 16.971 cm long. -Jp sets the radial scale and +# is therefore not affected by the rotation. +axis_length () { + printf "0 0\n0 1\n" | gmt mapproject $R "$1" --PROJ_LENGTH_UNIT=cm 2>/dev/null | \ + awk 'NR==1 {x=$1; y=$2} NR==2 {printf "%.3f\n", sqrt (($1-x)^2 + ($2-y)^2)}' +} +cat << EOF > answer.txt +-JP12c 12.000 +-JP12c+t45 16.971 +-JP12c+t90 12.000 +-JP12c+t-90 12.000 +-JP12c+du 12.000 +-Jp12c 12.000 +-Jp12c+t45 12.000 +-Jp12c+t90 12.000 +EOF +rm -f result.txt +for J in -JP12c -JP12c+t45 -JP12c+t90 -JP12c+t-90 -JP12c+du -Jp12c -Jp12c+t45 -Jp12c+t90; do + echo "$J $(axis_length $J)" >> result.txt +done +diff result.txt answer.txt > fail + +# 2. A full circle with -JP12c is 12 cm across, i.e., the radius is 6 cm +printf "0 0\n0 1\n" | gmt mapproject -R0/360/0/1 -JP12c --PROJ_LENGTH_UNIT=cm 2>/dev/null | \ + awk 'NR==1 {x=$1; y=$2} NR==2 {printf "%.3f\n", sqrt (($1-x)^2 + ($2-y)^2)}' > radius.txt +echo 6.000 > radius_answer.txt +diff radius.txt radius_answer.txt >> fail + +# 3. The axis keeps its 6 annotations on both radial edges no matter how it is rotated +cat << EOF > annot_answer.txt ++t0 12 ++t-45 12 ++t-90 12 ++t90 12 ++t135 12 +EOF +rm -f annot.txt +for t in 0 -45 -90 90 135; do + echo "+t$t $(gmt psbasemap $R -Jp12c+t$t -Bya -P 2>/dev/null | grep -oE '\) [a-z][a-z] Z' | wc -l | tr -d ' ')" >> annot.txt +done +diff annot.txt annot_answer.txt >> fail