filter1d: Add -N+e<cols> modifier to exclude columns from filtering

doc/rst/source/filter1d.rst

@@ -17,7 +17,7 @@ Synopsis
 [ |-D|\ *increment* ]
 [ |-E| ]
 [ |-L|\ *lack\_width* ]
-[ |-N|\ *t\_col* ]
+[ |-N|\ *t\_col*\ [**+e**\ *cols*] ]
 [ |-Q|\ *q\_factor* ]
 [ |-S|\ *symmetry\_factor* ]
 [ |-T|\ [*min/max*\ /]\ *inc*\ [**+a**][**+e**\|\ **i**\|\ **n**] \|\ |-T|\ *file*\|\ *list* ]
@@ -133,9 +133,14 @@ Optional Arguments
 
 .. _-N:
 
-**-N**\ *t_col*
+**-N**\ *t_col*\ [**+e**\ *cols*]
     Indicates which column contains the independent variable (time). The
     left-most column is 0, while the right-most is (*n_cols* - 1) [Default is 0].
+    Append **+e** and a comma-separated list of extra column(s) or column
+    ranges (e.g., 1,3-5) to pass through unfiltered, in addition to *t_col*.
+    By default, every column except *t_col* is filtered; use this modifier
+    to protect specific columns (e.g., coordinates) from being altered.
+    Requires output at the input abscissae, i.e., cannot be combined with |-T|.
 
 .. _-Q:
 
@@ -232,6 +237,16 @@ and not shorten the track, and add the distances every 2km to the file, use
 
   gmt filter1d track.txt -T2k+a -E -Fg200 > smooth_track.txt
 
+Given a profile with columns X, Y, Distance, and Elevation, filtering
+normally applies to every column except the independent one, which would
+also smooth X and Y. To smooth only the Elevation (column 3) as a function
+of Distance (column 2) while leaving X and Y untouched, exclude them with
+**+e**:
+
+::
+
+  gmt filter1d profile.txt -Fb100 -N2+e0,1 -E > profile_filt.txt
+
 See Also
 --------
 

src/filter1d.c

@@ -71,13 +71,16 @@ struct FILTER1D_CTRL {
 		bool active;
 		double value;
 	} L;
-	struct FILTER1D_N {	/* -N<t_col> or -Nc|<unit>[+a] */
+	struct FILTER1D_N {	/* -N<t_col>[+e<cols>] or -Nc|<unit>[+a] */
 		bool active;
 		bool add_col;
 		char unit;
 		unsigned int mode;
 		unsigned int spatial;
 		int col;
+		bool except[GMT_MAX_COLUMNS];	/* true for extra columns (besides t_col) to exclude from filtering */
+		uint64_t n_except;		/* How many columns were flagged in except */
+		uint64_t max_except;		/* Highest column index flagged, for range-checking against the data */
 	} N;
 	struct FILTER1D_Q {	/* -Q<factor> */
 		bool active;
@@ -121,6 +124,7 @@ struct FILTER1D_INFO {	/* Control structure for all aspects of the filter setup
 	uint64_t *n_this_col;		/* Pointer to array of counters [one per column]  */
 	uint64_t *n_left;		/* Pointer to array of counters [one per column]  */
 	uint64_t *n_right;		/* Pointer to array of counters [one per column]  */
+	bool *skip_col;			/* true for t_col and any user-selected columns to pass through unfiltered  */
 	uint64_t n_cols;		/* Number of columns of input  */
 	uint64_t t_col;			/* Column of time abscissae (independent variable)  */
 	uint64_t n_f_wts;		/* Number of filter weights  */
@@ -187,7 +191,7 @@ static int usage (struct GMTAPI_CTRL *API, int level) {
 	const char *name = gmt_show_name_and_purpose (API, THIS_MODULE_LIB, THIS_MODULE_CLASSIC_NAME, THIS_MODULE_PURPOSE);
 	if (level == GMT_MODULE_PURPOSE) return (GMT_NOERROR);
 	GMT_Usage (API, 0, "usage: %s [<table>] -F<type><width>[+h] [-D<increment>] [-E] "
-		"[-L<lack_width>] [-N<t_col>] [-Q<q_factor>] [-S<symmetry>] [-T[<min>/<max>/]<inc>|<file>|<list>[+a][+e|i|n]] "
+		"[-L<lack_width>] [-N<t_col>[+e<cols>]] [-Q<q_factor>] [-S<symmetry>] [-T[<min>/<max>/]<inc>|<file>|<list>[+a][+e|i|n]] "
 		"[%s] [%s] [%s] [%s] [%s] [%s] [%s] [%s] [%s] [%s] [%s] [%s] [%s] [%s]\n",
 		name, GMT_V_OPT, GMT_a_OPT, GMT_b_OPT, GMT_d_OPT, GMT_e_OPT, GMT_f_OPT, GMT_g_OPT, GMT_h_OPT, GMT_i_OPT,
 		GMT_j_OPT, GMT_o_OPT, GMT_q_OPT, GMT_colon_OPT, GMT_PAR_OPT);
@@ -228,9 +232,12 @@ static int usage (struct GMTAPI_CTRL *API, int level) {
 	GMT_Usage (API, 1, "\n-L<lack_width>");
 	GMT_Usage (API, -2, "Check for lack of data condition.  If input data has a gap exceeding "
 		"<width> then no output will be given at that point [Default does not check Lack].");
-	GMT_Usage (API, 1, "\n-N<t_col>");
+	GMT_Usage (API, 1, "\n-N<t_col>[+e<cols>]");
 	GMT_Usage (API, -2, "Set the column that contains the independent variable (time) [0]. "
 		"The left-most column is 0, the right-most is (<n_cols> - 1).");
+	GMT_Usage (API, 3, "e: Append a comma-separated list of extra column(s) or column ranges "
+		"(e.g., 1,3-5) to pass through unfiltered, in addition to <t_col>. "
+		"Requires output at the input abscissae (i.e., do not combine with -T).");
 	GMT_Usage (API, 1, "\n-Q<q_factor>");
 	GMT_Usage (API, -2, "Assess quality of output value by checking mean weight in convolution. "
 		"Append <q_factor> between 0 and 1.  If mean weight < q_factor, output is "
@@ -254,6 +261,58 @@ static int usage (struct GMTAPI_CTRL *API, int level) {
 	return (GMT_MODULE_USAGE);
 }
 
+GMT_LOCAL unsigned int filter1d_parse_except_columns (struct GMT_CTRL *GMT, char option, char *list, struct FILTER1D_N *N) {
+	/* Parse a comma-separated list of columns and/or ranges (e.g., "1,3-5") that are to be
+	 * passed through unfiltered, and flag each one in N->except.  We validate strictly here
+	 * because a silently mis-parsed list would quietly filter the wrong columns. */
+	unsigned int n_errors = 0, pos = 0;
+	int64_t start, stop, col;
+	char *token = NULL, *sep = NULL;
+
+	if (list == NULL || list[0] == '\0') {
+		GMT_Report (GMT->parent, GMT_MSG_ERROR, "Option -%c: Modifier +e requires one or more columns\n", option);
+		return (1);
+	}
+	token = gmt_M_memory (GMT, NULL, strlen (list) + 1, char);	/* gmt_strtok needs room for the whole list */
+	while (gmt_strtok (list, ",", &pos, token)) {
+		if ((sep = strpbrk (token, "-:/"))) {	/* Got a range of columns, e.g., 3-5 */
+			char code = sep[0];
+			sep[0] = '\0';	/* Split the range into its two halves */
+			if (!gmt_is_integer (token) || !gmt_is_integer (&sep[1])) {
+				GMT_Report (GMT->parent, GMT_MSG_ERROR, "Option -%c: Modifier +e: Bad column range %s%c%s\n", option, token, code, &sep[1]);
+				++n_errors;
+				continue;
+			}
+			start = atol (token);	stop = atol (&sep[1]);
+		}
+		else {	/* Just a single column */
+			if (!gmt_is_integer (token)) {
+				GMT_Report (GMT->parent, GMT_MSG_ERROR, "Option -%c: Modifier +e: Bad column %s\n", option, token);
+				++n_errors;
+				continue;
+			}
+			start = stop = atol (token);
+		}
+		if (stop < start) {
+			GMT_Report (GMT->parent, GMT_MSG_ERROR, "Option -%c: Modifier +e: Column range %" PRId64 "-%" PRId64 " is not increasing\n", option, start, stop);
+			++n_errors;
+			continue;
+		}
+		if (stop >= GMT_MAX_COLUMNS) {
+			GMT_Report (GMT->parent, GMT_MSG_ERROR, "Option -%c: Modifier +e: Column %" PRId64 " exceeds the maximum of %d columns\n", option, stop, GMT_MAX_COLUMNS - 1);
+			++n_errors;
+			continue;
+		}
+		for (col = start; col <= stop; col++) {
+			if (!N->except[col]) N->n_except++;	/* Ignore any repeated columns */
+			N->except[col] = true;
+			if ((uint64_t)col > N->max_except) N->max_except = (uint64_t)col;
+		}
+	}
+	gmt_M_free (GMT, token);
+	return (n_errors);
+}
+
 GMT_LOCAL char filter1d_set_unit_and_mode (const char *arg, unsigned int *mode) {
 	unsigned int k = 0;
 	*mode = GMT_GREATCIRCLE;	/* Default is great circle distances */
@@ -273,7 +332,7 @@ static int parse (struct GMT_CTRL *GMT, struct FILTER1D_CTRL *Ctrl, struct GMT_O
 
 	unsigned int n_errors = 0;
 	int sval = 0;
-	char *c = NULL, p, txt[GMT_LEN64] = {""}, *t_arg = NULL;
+	char *c = NULL, *e = NULL, p, txt[GMT_LEN64] = {""}, *t_arg = NULL;
 	struct GMT_OPTION *opt = NULL;
 	struct GMTAPI_CTRL *API = GMT->parent;
 
@@ -349,8 +408,12 @@ static int parse (struct GMT_CTRL *GMT, struct FILTER1D_CTRL *Ctrl, struct GMT_O
 				n_errors += gmt_M_repeated_module_option (API, Ctrl->L.active);
 				n_errors += gmt_get_required_double (GMT, opt->arg, opt->option, 0, &Ctrl->L.value);
 				break;
-			case 'N':	/* Select column with independent coordinate [0] */
+			case 'N':	/* Select column with independent coordinate [0], optionally +e<cols> to exclude more columns from filtering */
 				n_errors += gmt_M_repeated_module_option (API, Ctrl->N.active);
+				if ((e = strstr (opt->arg, "+e")) != NULL) {	/* Gave list of extra columns to exclude from filtering */
+					n_errors += filter1d_parse_except_columns (GMT, opt->option, &e[2], &Ctrl->N);
+					e[0] = '\0';	/* Chop off the modifier so the time-column parsing below is unaffected */
+				}
 				if (gmt_M_compat_check (GMT, 4)) {	/* GMT4 LEVEL */
 					if (strchr (opt->arg, '/')) { /* Gave obsolete format */
 						int sval0;
@@ -389,6 +452,7 @@ static int parse (struct GMT_CTRL *GMT, struct FILTER1D_CTRL *Ctrl, struct GMT_O
 					n_errors += gmt_M_check_condition (GMT, sval < 0, "Option -N: Time column cannot be negative.\n");
 					Ctrl->N.col = sval;
 				}
+				if (e) e[0] = '+';	/* Restore the modifier we chopped off above */
 				break;
 			case 'Q':	/* Assess quality of output */
 				n_errors += gmt_M_repeated_module_option (API, Ctrl->Q.active);
@@ -426,6 +490,8 @@ static int parse (struct GMT_CTRL *GMT, struct FILTER1D_CTRL *Ctrl, struct GMT_O
 	n_errors += gmt_M_check_condition (GMT, Ctrl->L.active && (Ctrl->L.value < 0.0 || Ctrl->L.value > Ctrl->F.width) , "Option -L: Unreasonable lack-of-data interval\n");
 	n_errors += gmt_M_check_condition (GMT, Ctrl->S.active && (Ctrl->S.value < 0.0 || Ctrl->S.value > 1.0) , "Option -S: Enter a factor between 0 and 1\n");
 	n_errors += gmt_M_check_condition (GMT, Ctrl->Q.active && (Ctrl->Q.value < 0.0 || Ctrl->Q.value > 1.0), "Option -Q: Enter a factor between 0 and 1\n");
+	n_errors += gmt_M_check_condition (GMT, Ctrl->N.n_except > 0 && Ctrl->T.active,
+		"Option -N+e: Cannot combine with -T; excluded columns pass through at the input abscissae only\n");
 
 	return (n_errors ? GMT_PARSE_ERROR : GMT_NOERROR);
 }
@@ -658,7 +724,7 @@ GMT_LOCAL int filter1d_do_the_filter (struct GMTAPI_CTRL *C, struct FILTER1D_INF
 		for (i_col = 0; i_col < F->n_cols; ++i_col) {
 			F->n_this_col[i_col] = 0;
 			wt_sum[i_col] = data_sum[i_col] = 0.0;
-			if (i_col == F->t_col)
+			if (F->skip_col[i_col])
 				good_one[i_col] = false;
 			else if (F->check_lack)
 				good_one[i_col] = !(filter1d_lack_check (F, i_col, left, right));
@@ -719,6 +785,8 @@ GMT_LOCAL int filter1d_do_the_filter (struct GMTAPI_CTRL *C, struct FILTER1D_INF
 			for (i_col = 0; i_col < F->n_cols; ++i_col) {
 				if (i_col == F->t_col)
 					data_sum[i_col] = t_time;
+				else if (F->skip_col[i_col])
+					data_sum[i_col] = F->data[i_col][k];	/* Pass through original value unfiltered */
 				else if (good_one[i_col]) {
 					data_sum[i_col] = (F->highpass) ? F->data[i_col][k] - F->this_loc[i_col] : F->this_loc[i_col];
 					++n_good_ones;
@@ -784,6 +852,8 @@ GMT_LOCAL int filter1d_do_the_filter (struct GMTAPI_CTRL *C, struct FILTER1D_INF
 				for (i_col = 0; i_col < F->n_cols; ++i_col) {
 					if (i_col == F->t_col)
 						outval[i_col] = t_time;
+					else if (F->skip_col[i_col])
+						outval[i_col] = F->data[i_col][k];	/* Pass through original value unfiltered */
 					else if (good_one[i_col]) {
 						outval[i_col] = (F->f_operator) ? data_sum[i_col] : data_sum[i_col] / wt_sum[i_col];
 						if (F->highpass) outval[i_col] = F->data[i_col][k] - outval[i_col];
@@ -809,6 +879,7 @@ GMT_LOCAL int filter1d_do_the_filter (struct GMTAPI_CTRL *C, struct FILTER1D_INF
 
 GMT_LOCAL int filter1d_allocate_space (struct GMT_CTRL *GMT, struct FILTER1D_INFO *F) {
 	F->n_this_col = gmt_M_memory (GMT, NULL, F->n_cols, uint64_t);
+	F->skip_col = gmt_M_memory (GMT, NULL, F->n_cols, bool);
 	F->data = gmt_M_memory_aligned (GMT, NULL, F->n_cols, double *);
 
 	if (F->check_asym) F->n_left = gmt_M_memory (GMT, NULL, F->n_cols, uint64_t);
@@ -841,6 +912,7 @@ GMT_LOCAL void filter1d_free_space (struct GMT_CTRL *GMT, struct FILTER1D_INFO *
 	for (i = 0; i < F->n_cols; ++i)	gmt_M_free (GMT, F->data[i]);
 	gmt_M_free (GMT, F->data);
 	gmt_M_free (GMT, F->n_this_col);
+	gmt_M_free (GMT, F->skip_col);
 	gmt_M_free (GMT, F->n_left);
 	gmt_M_free (GMT, F->n_right);
 	gmt_M_free (GMT, F->min_loc);
@@ -1047,6 +1119,24 @@ EXTERN_MSC int GMT_filter1d (void *V_API, int mode, void *args) {
 
 	filter1d_allocate_space (GMT, &F);	/* Gets column-specific flags and uint64_t space */
 
+	F.skip_col[F.t_col] = true;
+	if (Ctrl->N.n_except) {	/* Flag the extra columns the user wants passed through unfiltered */
+		uint64_t n_to_filter = 0;
+		if (Ctrl->N.max_except >= F.n_cols) {
+			uint64_t bad_col = Ctrl->N.max_except, have_cols = F.n_cols;	/* Return() frees Ctrl before formatting the message below */
+			filter1d_free_space (GMT, &F);
+			Return (GMT_PARSE_ERROR, "Option -N+e: Column %" PRIu64 " is beyond the %" PRIu64 " columns found in the data\n",
+				bad_col, have_cols);
+		}
+		for (col = 0; col <= Ctrl->N.max_except; col++)
+			if (Ctrl->N.except[col]) F.skip_col[col] = true;
+		for (col = 0; col < F.n_cols; col++) if (!F.skip_col[col]) n_to_filter++;
+		if (n_to_filter == 0) {
+			filter1d_free_space (GMT, &F);
+			Return (GMT_PARSE_ERROR, "Option -N+e: No columns left to filter\n");
+		}
+	}
+
 	gmt_increase_abstime_format_precision (GMT, Ctrl->N.col, F.t_int);	/* In case we need more sub-second precision output */
 
 	GMT_Report (API, GMT_MSG_INFORMATION, "Filter the data columns\n");

src/longopt/filter1d_inc.h

@@ -33,7 +33,7 @@ static struct GMT_KEYWORD_DICTIONARY module_kw[] = {
 	          "h,l,u",             "highpass,lower,upper",
 		  GMT_TP_STANDARD },
 	{ 0, 'L', "gap_width",         "", "", "", "", GMT_TP_STANDARD },
-	{ 0, 'N', "time_column|time_col",          "", "", "", "", GMT_TP_STANDARD },
+	{ 0, 'N', "time_column|time_col",          "", "", "e", "exclude", GMT_TP_STANDARD },
 	{ 0, 'Q', "quality",           "", "", "", "", GMT_TP_STANDARD },
 	{ 0, 'S', "symmetry",          "", "", "", "", GMT_TP_STANDARD },
 	{ 0, 'T', "range|series",

test/filter1d/except_columns.sh

@@ -0,0 +1,65 @@
+#!/usr/bin/env bash
+# Testing gmt filter1d -N<t_col>+e<cols>: exclude extra columns from filtering
+# (e.g., keep X,Y untouched while smoothing Elevation vs Distance).
+
+AWK=${AWK:-awk}
+data=except_columns_data.txt
+out=except_columns_out.txt
+
+cat << EOF > $data
+1000	2000	0	100
+1015	1980	25	120
+1030	1960	50	90
+1045	1940	75	150
+1060	1920	100	80
+1075	1900	125	130
+1090	1880	150	95
+1105	1860	175	110
+1120	1840	200	140
+1135	1820	225	90
+1150	1800	250	105
+1165	1780	275	120
+1180	1760	300	100
+EOF
+
+# Column 2 (Distance) is independent; columns 0,1 (X,Y) are excluded from
+# filtering along with it, so only column 3 (Elevation) gets filtered.
+gmt filter1d $data -Fl250 -N2+e0,1 -E > $out
+
+# Row count must be preserved.
+n_in=$(wc -l < $data)
+n_out=$(wc -l < $out)
+if [ "$n_in" != "$n_out" ]; then
+	echo "row count mismatch: input has $n_in rows but output has $n_out" > fail
+fi
+
+# X,Y (columns 1-2) must be byte-identical to the input.
+$AWK '{print $1, $2}' $data > xy_in.txt
+$AWK '{print $1, $2}' $out > xy_out.txt
+diff xy_in.txt xy_out.txt >> fail
+
+# Without +e, gmt filter1d would also filter X,Y (this is the behavior +e
+# fixes) - confirm that still holds so we know the test data is meaningful.
+gmt filter1d $data -Fl250 -N2 -E > no_except_out.txt
+if diff -q xy_in.txt <($AWK '{print $1, $2}' no_except_out.txt) > /dev/null; then
+	echo "expected X,Y to change without +e, but they did not" >> fail
+fi
+
+# A range must mean the same thing as the equivalent explicit list.
+gmt filter1d $data -Fl250 -N2+e0-1 -E > range_out.txt
+diff $out range_out.txt >> fail
+
+# Bad column lists must be rejected, not silently mis-parsed: a non-numeric
+# token, an open-ended or reversed range, an empty list, a column beyond the
+# data, and excluding every column so nothing is left to filter.
+for bad in "+eabc" "+e1-" "+e3-1" "+e" "+e0,1,99" "+e0,1,3"; do
+	if gmt filter1d $data -Fl250 -N2$bad -E > /dev/null 2>&1; then
+		echo "-N2$bad should have failed but did not" >> fail
+	fi
+done
+
+# +e cannot be combined with -T since excluded columns only pass through at
+# the input abscissae.
+if gmt filter1d $data -Fl250 -N2+e0,1 -T0/300/25 > /dev/null 2>&1; then
+	echo "-N2+e0,1 with -T should have failed but did not" >> fail
+fi

test/filter1d/filter1d-l2s.sh

@@ -20,6 +20,7 @@ cat << EOF > $a
 --l2stranstest -Fl100 -FL100
 --l2stranstest -Fu100 -FU100
 --l2stranstest -L15 -N3 -N5
+--l2stranstest -N3+e0,1
 --l2stranstest -Q0.4 -S0.99
 --l2stranstest -T10/100/5+a -Tfile
 --l2stranstest -T5+e -T5+i -T5+n
@@ -37,6 +38,7 @@ gmt $m $l2s --filter=mlprob:75+lower --filter=rmlprob:25/5+upper >> $b
 gmt $m $l2s --filter=minall:100 --filter=minpos:100 >> $b
 gmt $m $l2s --filter=maxall:100 --filter=maxneg:100 >> $b
 gmt $m $l2s --gap_width=15 --time_col=3 --time_column=5 >> $b
+gmt $m $l2s --time_column=3+exclude:0,1 >> $b
 gmt $m $l2s --quality=0.4 --symmetry=0.99 >> $b
 gmt $m $l2s --range=10/100/5+array --range=file >> $b
 gmt $m $l2s --series=5+exact  --series=5+inverse --series=5+number >> $b