Fix thin polar sectors that end up parallel to an axis (#7059)

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

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<scale> 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;

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<scale> 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