WIP: grdcontour gdal

src/gmt_prototypes.h

@@ -719,6 +719,8 @@ EXTERN_MSC int gmt_polygon_centroid (struct GMT_CTRL *GMT, double *x, double *y,
 EXTERN_MSC int gmt_get_distance (struct GMT_CTRL *GMT, char *line, double *dist, char *unit);
 EXTERN_MSC unsigned int *gmt_contour_edge_init (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *header, unsigned int *n_edges);
 EXTERN_MSC int64_t gmt_contours (struct GMT_CTRL *GMT, struct GMT_GRID *Grid, unsigned int smooth_factor, unsigned int int_scheme, int orient, unsigned int *edge, bool *first, double **x, double **y);
+EXTERN_MSC int64_t gmt_smooth_contour(struct GMT_CTRL *GMT, double **x_in, double **y_in, uint64_t n, int sfactor, int stype);
+EXTERN_MSC void gmt_orient_contour(struct GMT_GRID *G, double *x, double *y, uint64_t n, int orient);
 EXTERN_MSC int gmt_get_format (struct GMT_CTRL *GMT, double interval, char *unit, char *prefix, char *format);
 EXTERN_MSC int gmt_get_index (struct GMT_CTRL *GMT, struct GMT_PALETTE *P, double *value);
 EXTERN_MSC int gmt_get_rgb_from_z (struct GMT_CTRL *GMT, struct GMT_PALETTE *P, double value, double *rgb);

src/gmt_support.c

@@ -11890,6 +11890,19 @@ unsigned int *gmt_contour_edge_init (struct GMT_CTRL *GMT, struct GMT_GRID_HEADE
 	return (edge);
 }
 
+/*! . */
+int64_t gmt_smooth_contour(struct GMT_CTRL *GMT, double **x_in, double **y_in, uint64_t n, int sfactor, int stype) {
+	/* Public front to gmtsupport_smooth_contour so that modules that get their contours from
+	 * elsewhere (e.g., grdcontour, which traces them with GDAL) can still honor -S */
+	return (gmtsupport_smooth_contour(GMT, x_in, y_in, n, sfactor, stype));
+}
+
+/*! . */
+void gmt_orient_contour(struct GMT_GRID *G, double *x, double *y, uint64_t n, int orient) {
+	/* Public front to gmtsupport_orient_contour, for the same reason as gmt_smooth_contour above */
+	gmtsupport_orient_contour(G, x, y, n, orient);
+}
+
 /*! . */
 int64_t gmt_contours (struct GMT_CTRL *GMT, struct GMT_GRID *G, unsigned int smooth_factor, unsigned int int_scheme, int orient, unsigned int *edge, bool *first, double **x, double **y) {
 	/* The routine finds the zero-contour in the grd dataset.  it assumes that

test/grdcontour/gdal_circles.sh

@@ -0,0 +1,95 @@
+#!/usr/bin/env bash
+#
+# Contour z = x^2 + y^2, whose level L is exactly the circle of radius sqrt(L),
+# and check the traced lines against that analytic answer.  This pins down the
+# things a marching-squares tracer can get subtly wrong on a curved contour:
+# the radius of every vertex, whether an interior contour comes back as a single
+# closed ring, and whether -S densifies the line without moving it.
+#
+# The same grid is then punched with a NaN disk centered at (3.5,0) with radius
+# 0.8, which the circles of levels 10 and 15 (radius 3.162 and 3.873) run
+# straight through while those of levels 5 and 20 (2.236 and 4.472) miss.  The
+# two affected levels must come back as a single open arc each, the other two
+# must still be single closed rings, and no vertex may fall inside the hole.
+#
+# Interpolating a quadratic along a grid edge is off by at most h^2/4 = 6.3e-4
+# in z, i.e. 1.4e-4 in radius on the tightest circle here, so 1e-3 is a tight
+# tolerance that still catches an error of half a cell (0.025).
+
+rm -f fail
+
+gmt grdmath -R-5/5/-5/5 -I0.05 X SQR Y SQR ADD = r2.nc
+gmt grdmath r2.nc X 3.5 SUB SQR Y SQR ADD SQRT 0.8 GE 0 NAN MUL = r2nan.nc
+
+# radii <dump> <tag>: emit one line per level -
+#   <tag> <level> <n segments> <closure pattern> <worst radial error> <vertices in the hole>
+radii () {
+	$AWK -v tag="$2" -v holed=$3 '
+	function report () {
+		if (n == 0) return
+		nseg[lev]++
+		gap = sqrt ((x0-xp)*(x0-xp) + (y0-yp)*(y0-yp))
+		closed[lev] = closed[lev] (gap < 0.005 ? "C" : "O")
+	}
+	/^>/ { report(); n = 0; next }
+	NF >= 3 {
+		lev = $3
+		r = sqrt ($1*$1 + $2*$2)
+		d = r - sqrt (lev); if (d < 0) d = -d
+		if (d > worst[lev]) worst[lev] = d
+		if (holed && sqrt (($1-3.5)*($1-3.5) + $2*$2) < 0.79) hole[lev]++
+		if (n == 0) { x0 = $1; y0 = $2 }
+		xp = $1; yp = $2
+		n++
+	}
+	END {	report()
+		for (L in nseg)
+			printf "%s %g %d %s %.6f %d\n", tag, L, nseg[L], closed[L], worst[L]+0, hole[L]+0
+	}' "$1"
+}
+
+# expect <tag> <level> <n segments> <closure pattern> <radial tolerance>
+expect () {
+	$AWK -v tag="$1" -v L="$2" -v ns="$3" -v cl="$4" -v tol="$5" '
+	$1 == tag && $2 == L {
+		seen = 1
+		if ($3 != ns)  printf "%s: level %s has %s segments, expected %s\n", tag, L, $3, ns
+		if ($4 != cl)  printf "%s: level %s closure is %s, expected %s\n", tag, L, $4, cl
+		if ($5 > tol)  printf "%s: level %s is off the analytic radius by %s (tol %s)\n", tag, L, $5, tol
+		if ($6 != 0)   printf "%s: level %s has %s vertices inside the NaN hole\n", tag, L, $6
+	}
+	END { if (!seen) printf "%s: level %s is missing from the dump\n", tag, L }' summary.txt
+}
+
+# 1. Clean grid: four interior circles, one closed ring each
+gmt grdcontour r2.nc -C5 -L5/20 -JX10c -Dclean.txt
+radii clean.txt clean 0 > summary.txt
+
+# 2. Same, smoothed 4x: must have far more points but stay on the same circles
+gmt grdcontour r2.nc -C5 -L5/20 -S4 -JX10c -Dsmooth.txt
+radii smooth.txt smooth 0 >> summary.txt
+
+# 3. NaN-holed grid
+gmt grdcontour r2nan.nc -C5 -L5/20 -JX10c -Dnan.txt
+radii nan.txt nan 1 >> summary.txt
+
+for t in clean smooth; do
+	expect $t 5  1 C 0.001 >> fail
+	expect $t 10 1 C 0.001 >> fail
+	expect $t 15 1 C 0.001 >> fail
+	expect $t 20 1 C 0.001 >> fail
+done
+
+# Cutting a closed ring once leaves ONE open arc.  Its two ends stop on the last
+# cell that still had four valid corners, so they can sit up to half a cell off
+# the analytic radius - hence the looser tolerance on the two cut levels.
+expect nan 5  1 C 0.001 >> fail
+expect nan 10 1 O 0.03  >> fail
+expect nan 15 1 O 0.03  >> fail
+expect nan 20 1 C 0.001 >> fail
+
+# 4. -S4 must densify the lines about fourfold
+$AWK 'END { print NR }' clean.txt  > n_clean.txt
+$AWK 'END { print NR }' smooth.txt > n_smooth.txt
+$AWK 'NR == FNR { a = $1; next } { if ($1 < 3 * a) printf "-S4 produced %d lines from %d, expected at least 3x\n", $1, a }' \
+	n_clean.txt n_smooth.txt >> fail

test/grdcontour/gdal_level_accuracy.sh

@@ -0,0 +1,57 @@
+#!/usr/bin/env bash
+#
+# GMT_KNOWN_FAILURE
+#
+# A contour must be drawn at the level it is labelled with, whatever the local
+# gradient of the grid is.  This test fails on the GDAL tracer path and passes
+# with gmt_contours().
+#
+# grdcontour computes a noise threshold
+#
+#	small = (C interval > 0) ? MIN (C interval, z_range) * 1e-6 : z_range * 1e-6
+#
+# and gmt_contours() used it as
+#
+#	G = G_orig - cval;  if (G == 0) G += small;
+#
+# i.e. it nudged only the nodes whose value was EXACTLY the contour level, so
+# that no contour ran through a node, and left the contour's position alone.
+# The GDAL path instead hands the tracer cval - small as the level, which moves
+# the whole contour by small / |grad z|.  That is invisible on a steep grid and
+# unbounded on a flat one.
+#
+# Here z = exp(-x) over 0/20, so z_range is 1 and, because the level comes from
+# a contour FILE, the C interval is zero and small is the full z_range * 1e-6 =
+# 1e-6.  The wanted level is 3e-6, whose analytic position is
+#
+#	x = -ln (3e-6) = 12.7168983
+#
+# while the level actually traced, 3e-6 - 1e-6 = 2e-6, sits at
+#
+#	x = -ln (2e-6) = 13.1223634
+#
+# 0.405 away, or eight grid cells.  Linear interpolation across one 0.05 cell of
+# an exponential is good to 3e-4 in x, so half a cell is a generous tolerance.
+
+rm -f fail
+
+echo "0.000003 C" > cont.txt
+gmt grdmath -R0/20/0/1 -I0.05 X NEG EXP = decay.nc
+gmt grdcontour decay.nc -Ccont.txt -JX10c/2c -Ddecay.txt
+
+$AWK '
+BEGIN { want = 12.7168983; tol = 0.025 }
+/^>/ { nseg++; next }
+NF >= 3 {
+	d = $1 - want; if (d < 0) d = -d
+	if (d > worst) { worst = d; xworst = $1 }
+	if (n == 0 || $2 < ymin) ymin = $2
+	if (n == 0 || $2 > ymax) ymax = $2
+	n++
+}
+END {
+	if (n == 0) { print "no contour was traced at all"; exit }
+	if (nseg != 1) printf "%d segments, expected 1\n", nseg
+	if (worst > tol) printf "contour traced at x = %.4f, wanted %.4f (off by %.4f)\n", xworst, want, worst
+	if (ymin > 0.01 || ymax < 0.99) printf "contour spans only y = %g to %g, expected 0 to 1\n", ymin, ymax
+}' decay.txt >> fail

test/grdcontour/gdal_phase.sh

@@ -0,0 +1,60 @@
+#!/usr/bin/env bash
+#
+# Check that -Z+p (phase data) still honors the 360-degree wrap.
+#
+# grdcontour traces everything with GDAL except -Z+p, which stays on
+# gmt_contours().  The reason is gmtsupport_setcontjump(): on every cell it
+# unwraps that cell's corner values onto the branch nearest its first corner, so
+# a 179/-179 pair of nodes reads as 179/181 and no contour is drawn along the
+# wrap.  That decision is made per cell, i.e. the field is multivalued, while a
+# GDAL band carries one value per node, so the rule cannot be handed to GDAL.
+# If -Z+p ever gets routed through the GDAL path this test is what fires.
+#
+# The grid is a pure sawtooth, z = wrap(x + 90) in -180/180, so it has a branch
+# cut at x = 90 where the value jumps from 179 to -179.  Its contours are single
+# meridians, analytically at
+#
+#	level  -120  -60    0   60  120
+#	x       150 -150  -90  -30   30
+#
+# and there must be NOTHING near the cut: without the wrap handling the cell at
+# x = 90 spans -180 to 180 and so crosses every single level, painting a
+# spurious contour along that meridian for all of them.
+#
+# -L-150/150 keeps the -180 level, which genuinely does sit on the cut, out of
+# the comparison.
+
+rm -f fail
+
+gmt grdmath -R-180/180/-60/60 -I1 X 90 ADD 180 ADD 360 FMOD 180 SUB = saw.nc
+gmt grdcontour saw.nc -C60 -L-150/150 -Z+p -JX10c/4c -Dsaw.txt
+
+$AWK '
+function report () {
+	if (n == 0) return
+	nseg[lev]++
+	if (xmax - xmin > 0.01) printf "level %g is not a meridian: x from %g to %g\n", lev, xmin, xmax
+	if (ymin > -59.99 || ymax < 59.99) printf "level %g spans only y = %g to %g\n", lev, ymin, ymax
+	d = xmin - want[lev]; if (d < 0) d = -d
+	if (d > 0.51) printf "level %g is at x = %g, expected %g\n", lev, xmin, want[lev]
+	if (xmin > 88 && xmin < 92) printf "level %g put a contour on the branch cut at x = %g\n", lev, xmin
+}
+BEGIN {	want[-120] = 150; want[-60] = -150; want[0] = -90; want[60] = -30; want[120] = 30 }
+/^>/ { report(); n = 0; next }
+NF >= 3 {
+	lev = $3
+	if (n == 0) { xmin = xmax = $1; ymin = ymax = $2 }
+	if ($1 < xmin) xmin = $1;  if ($1 > xmax) xmax = $1
+	if ($2 < ymin) ymin = $2;  if ($2 > ymax) ymax = $2
+	n++
+}
+END {	report()
+	nl = 0
+	for (L in nseg) {
+		nl++
+		if (!(L in want)) printf "unexpected level %s in the dump\n", L
+		if (nseg[L] != 1) printf "level %s has %d segments, expected 1\n", L, nseg[L]
+	}
+	if (nl != 5) printf "%d levels traced, expected 5\n", nl
+	for (L in want) if (!(L in nseg)) printf "level %s is missing from the dump\n", L
+}' saw.txt >> fail

test/grdcontour/gdal_registration.sh

@@ -0,0 +1,68 @@
+#!/usr/bin/env bash
+#
+# Check that grdcontour puts the contours of a linear ramp exactly where the
+# analytic answer says they are, for both gridline- and pixel-registered grids.
+#
+# The tracer reads the grid through an in-memory GDAL raster, and GDAL samples a
+# raster at its cell centers while GMT keeps its values at the grid nodes.  The
+# wrapper's geotransform therefore has to be shifted by half a cell for gridline
+# registration and not at all for pixel registration.  A wrong shift, a wrong
+# sign on the y increment (rows run north to south), or a line offset that
+# ignores the grid pad all move every contour by half a cell or flip the grid,
+# and that is what this test looks for.
+#
+#   z = x  ->  contours at x = 1,2,...,10, each spanning the full height
+#   z = y  ->  contours at y = 1,2,...,10, each spanning the full width
+#
+# The gridline-registered grid runs from 0 to 10 in z so the 10 contour lands on
+# its far edge and the 0 contour does not exist (every node is >= 0), which is 10
+# lines.  The pixel-registered one only reaches 0.25 to 9.75, so it has 9.
+#
+# It also checks that no dumped point lies outside the area covered by the grid
+# nodes: GDAL pads a raster with a virtual half-cell border repeating the
+# outermost row/column, so a contour reaching the edge of the grid comes back
+# with a half-cell stub sticking out of the domain unless it is trimmed off.
+
+rm -f fail
+
+# check <dump> <level column> <span column> <span min> <span max> <n segments>
+check () {
+	$AWK -v c=$2 -v o=$3 -v omin=$4 -v omax=$5 -v want=$6 -v tag="$1" '
+	function report () {
+		if (n == 0) return
+		nseg++
+		if (wmin > omin + 0.01 || wmax < omax - 0.01)
+			printf "%s: level %g spans only %g/%g, expected %g/%g\n", tag, lev, wmin, wmax, omin, omax
+		if (wmin < omin - 0.01 || wmax > omax + 0.01)
+			printf "%s: level %g runs outside the node area: %g/%g\n", tag, lev, wmin, wmax
+	}
+	/^>/ { report(); n = 0; next }
+	NF >= 3 {
+		lev = $3
+		d = $c - lev; if (d < 0) d = -d
+		if (d > 0.005) printf "%s: level %g has a point %g off the level\n", tag, lev, d
+		if (n == 0 || $o < wmin) wmin = $o
+		if (n == 0 || $o > wmax) wmax = $o
+		n++
+	}
+	END {	report()
+		if (nseg != want) printf "%s: %d segments, expected %d\n", tag, nseg, want
+	}' "$1"
+}
+
+# 1. Gridline registration: the nodes cover the whole -R, so the contours must too
+gmt grdmath -R0/10/0/10 -I0.5 X = ramp_x.nc
+gmt grdmath -R0/10/0/10 -I0.5 Y = ramp_y.nc
+gmt grdcontour ramp_x.nc -C1 -JX10c -Dgrid_x.txt
+gmt grdcontour ramp_y.nc -C1 -JX10c -Dgrid_y.txt
+check grid_x.txt 1 2 0 10 10 >> fail
+check grid_y.txt 2 1 0 10 10 >> fail
+
+# 2. Pixel registration: the nodes sit half a cell inside -R, so the contours stop
+#    half a cell short of it, which is exactly where gmt_contours also stopped
+gmt grdmath -R0/10/0/10 -I0.5 -r X = pix_x.nc
+gmt grdmath -R0/10/0/10 -I0.5 -r Y = pix_y.nc
+gmt grdcontour pix_x.nc -C1 -JX10c -Dpix_x.txt
+gmt grdcontour pix_y.nc -C1 -JX10c -Dpix_y.txt
+check pix_x.txt 1 2 0.25 9.75 9 >> fail
+check pix_y.txt 2 1 0.25 9.75 9 >> fail

test/grdcontour/gdal_selectors.sh

@@ -0,0 +1,65 @@
+#!/usr/bin/env bash
+#
+# Check the option that decide WHICH levels get traced and in which direction.
+#
+# On the GDAL path grdcontour pre-filters the level list before handing it to
+# the tracer, so that levels -L or -Q+z would throw away are not traced at all.
+# That filter has to agree exactly with the tests at the top of the contour loop
+# or levels go missing (or reappear) without any other symptom, so here we pin
+# the level set each selector must produce.
+#
+# z = x on -R-9.5/9.5 with -C2 gives the nine levels -8,-6,...,8, none of them
+# sitting on the domain boundary, and each one is a single vertical line.
+#
+# -F is checked on a closed ring, where flipping the orientation has to flip the
+# sign of the signed area while leaving its magnitude alone.
+
+rm -f fail
+
+gmt grdmath -R-9.5/9.5/-5/5 -I0.5 X = ramp.nc
+
+# levels <dump>: the sorted, unique level set of a dump, space separated.  The
+# segment header is "> <level> contour -Z<level>", so the level is the -Z field.
+levels () { $AWK '/^>/ { s = $NF; sub (/^-Z/, "", s); print s }' "$1" | sort -n -u | tr '\n' ' '; }
+
+# want <tag> <dump> <expected level set>
+want () {
+	got=$(levels "$2")
+	test "$got" = "$3 " || echo "$1: levels [$got] expected [$3]" >> fail
+	n=$($AWK '/^>/ { n++ } END { print n+0 }' "$2")
+	m=$($AWK -v s="$3" 'BEGIN { print split (s, a, " ") }')
+	test "$n" -eq "$m" || echo "$1: $n segments for $m levels, expected one each" >> fail
+}
+
+gmt grdcontour ramp.nc -C2       -JX10c -Dall.txt
+gmt grdcontour ramp.nc -C2 -Ln   -JX10c -Dln.txt
+gmt grdcontour ramp.nc -C2 -LN   -JX10c -DlN.txt
+gmt grdcontour ramp.nc -C2 -LP   -JX10c -DlP.txt
+gmt grdcontour ramp.nc -C2 -Lp   -JX10c -Dlp.txt
+gmt grdcontour ramp.nc -C2 -L-5/5 -JX10c -Drange.txt
+gmt grdcontour ramp.nc -C2 -Q+z  -JX10c -Dqz.txt
+
+want all   all.txt   "-8 -6 -4 -2 0 2 4 6 8"
+want Ln    ln.txt    "-8 -6 -4 -2"
+want LN    lN.txt    "-8 -6 -4 -2 0"
+want LP    lP.txt    "0 2 4 6 8"
+want Lp    lp.txt    "2 4 6 8"
+want L-5/5 range.txt "-4 -2 0 2 4"
+want Q+z   qz.txt    "-8 -6 -4 -2 2 4 6 8"
+
+# -F: same ring, opposite winding
+gmt grdmath -R-5/5/-5/5 -I0.05 X SQR Y SQR ADD = r2.nc
+gmt grdcontour r2.nc -C5 -L5/5 -Fl -JX10c -Dleft.txt
+gmt grdcontour r2.nc -C5 -L5/5 -Fr -JX10c -Dright.txt
+
+# area <dump>: twice the signed area of the first segment, by the shoelace formula
+area () { $AWK '/^>/ { if (n) exit; next } NF >= 3 { if (n) a += xp*$2 - $1*yp; xp = $1; yp = $2; n++ } END { printf "%.6f\n", a }' "$1"; }
+
+al=$(area left.txt)
+ar=$(area right.txt)
+$AWK -v l="$al" -v r="$ar" 'BEGIN {
+	if (l == 0 || r == 0) { print "-F: a signed area came out zero"; exit }
+	if (l * r > 0) printf "-F: -Fl and -Fr wind the same way (%g and %g)\n", l, r
+	d = (l < 0 ? -l : l) - (r < 0 ? -r : r); if (d < 0) d = -d
+	if (d > 0.001 * (l < 0 ? -l : l)) printf "-F: -Fl and -Fr enclose different areas (%g and %g)\n", l, r
+}' >> fail