diff --git a/doc/rst/source/filter1d.rst b/doc/rst/source/filter1d.rst index 192a9f40ce5..5959e3f573a 100644 --- a/doc/rst/source/filter1d.rst +++ b/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 -------- diff --git a/src/filter1d.c b/src/filter1d.c index 2d5c7053f5e..d73d36673c2 100644 --- a/src/filter1d.c +++ b/src/filter1d.c @@ -71,13 +71,16 @@ struct FILTER1D_CTRL { bool active; double value; } L; - struct FILTER1D_N { /* -N or -Nc|[+a] */ + struct FILTER1D_N { /* -N[+e] or -Nc|[+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 */ 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 [] -F[+h] [-D] [-E] " - "[-L] [-N] [-Q] [-S] [-T[//]||[+a][+e|i|n]] " + "[-L] [-N[+e]] [-Q] [-S] [-T[//]||[+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"); GMT_Usage (API, -2, "Check for lack of data condition. If input data has a gap exceeding " " then no output will be given at that point [Default does not check Lack]."); - GMT_Usage (API, 1, "\n-N"); + GMT_Usage (API, 1, "\n-N[+e]"); GMT_Usage (API, -2, "Set the column that contains the independent variable (time) [0]. " "The left-most column is 0, the right-most is ( - 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 . " + "Requires output at the input abscissae (i.e., do not combine with -T)."); GMT_Usage (API, 1, "\n-Q"); GMT_Usage (API, -2, "Assess quality of output value by checking mean weight in convolution. " "Append 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 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"); diff --git a/src/longopt/filter1d_inc.h b/src/longopt/filter1d_inc.h index 177088c0157..36cb1dc8f4d 100644 --- a/src/longopt/filter1d_inc.h +++ b/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", diff --git a/test/filter1d/except_columns.sh b/test/filter1d/except_columns.sh new file mode 100644 index 00000000000..c5a39c09b08 --- /dev/null +++ b/test/filter1d/except_columns.sh @@ -0,0 +1,65 @@ +#!/usr/bin/env bash +# Testing gmt filter1d -N+e: 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 diff --git a/test/filter1d/filter1d-l2s.sh b/test/filter1d/filter1d-l2s.sh index e012e729561..c71cb5bd987 100755 --- a/test/filter1d/filter1d-l2s.sh +++ b/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