diff --git a/doc/rst/source/gmt.conf.rst b/doc/rst/source/gmt.conf.rst index 8a306e238f6..36f0111065e 100644 --- a/doc/rst/source/gmt.conf.rst +++ b/doc/rst/source/gmt.conf.rst @@ -180,7 +180,11 @@ FORMAT Parameters without leading zeros (default uses fixed width formats). As examples, try hh:mm, hh.mm.ss, hh:mm:ss.xxxx, hha.m., etc. [hh:mm:ss]. If the format is simply - then no clock is output and - the ISO T divider between date and clock is omitted. + the ISO T divider between date and clock is omitted. **Note**: When + high-precision time-series are written to ASCII output the default + format may not be adequate. Many modules automatically handle + this by extending the format, but you should be alert of unusual + situations where data may appear truncated to nearest second. **FORMAT_DATE_IN** Formatting template that indicates how an input date string is diff --git a/src/filter1d.c b/src/filter1d.c index c24a632bcf1..0136a16e934 100644 --- a/src/filter1d.c +++ b/src/filter1d.c @@ -1028,6 +1028,8 @@ 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 */ + 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"); for (tbl = tseg = 0; tbl < D->n_tables; ++tbl) { /* For each input table */ diff --git a/src/gmt_io.c b/src/gmt_io.c index dcad576ed64..58ee548c1af 100644 --- a/src/gmt_io.c +++ b/src/gmt_io.c @@ -4572,6 +4572,82 @@ int gmtlib_nc_get_att_text (struct GMT_CTRL *GMT, int ncid, int varid, char *nam return status; } +GMT_LOCAL int gmtio_get_precision_width (struct GMT_CTRL *GMT, double x) { + /* Here, x < 1 seconds and we want to determine an approximate precision with max truncation error of GMT_CONV4_LIMIT */ + int k = -1; + double trunc_err, inv_x = 1.0 / x; + static double power_U[9] = {1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9}; + static double power_D[9] = {1e-1, 1e-2, 1e-3, 1e-4, 1e-5, 1e-6, 1e-7, 1e-8, 1e-9}; + + /* Here we compute |D-T|/D, where D is the incoming dt and T is the value recovered if printed using k decimals. + * We want this relative truncation error to be < 10^-4 */ + do { + k++; /* First time k becomes 0 */ + trunc_err = fabs (x - rint (x * power_U[k]) * power_D[k]) * inv_x; + // fprintf (stderr, "k = %d x = %g trunc_err = %g\n", k, x, trunc_err); + } while (k < 8 && trunc_err > GMT_CONV4_LIMIT); /* Stop at nanoseconds if things go off the rails */ + return (k + 1); +} + +GMT_LOCAL void gmtio_check_abstime_format (struct GMT_CTRL *GMT, struct GMT_DATASET *D) { + bool abstime_found = false; + unsigned int col, row; + int w_max = 0, this_w; + double sub; + struct GMT_DATASEGMENT *S = NULL; + + if (GMT->common.b.active[GMT_OUT]) return; /* Nothing to do if using binary i/o */ + if (D == NULL || D->table == NULL || D->table[0]->segment == NULL || D->table[0]->segment[0] == NULL) return; /* Nothing to work with */ + if (GMT->current.setting.time_system.unit != 's') return; /* Current unit is not second */ + if (strcmp (GMT->current.setting.format_clock_out, "hh:mm:ss")) return; /* User has changed from default format - do nothing */ + for (col = 0; !abstime_found && col < D->n_columns; col++) + if (GMT->current.io.col_type[GMT_OUT][col] == GMT_IS_ABSTIME) abstime_found = true; + if (!abstime_found) return; /* No absolute time columns to worry about */ + + /* If we get here we have ASCII absolute time output and the default FORMAT_CLOCK_OUT remains in effect. + * To assist the user we will scan the first segment's rows up to MIN (20, n_rows) and make our best guess + * from that subset. We will give information messages about our decision but not a warning. */ + + S = D->table[0]->segment[0]; /* The first segment */ + for (col = 0; col < D->n_columns; col++) { + if (GMT->current.io.col_type[GMT_OUT][col] != GMT_IS_ABSTIME) continue; /* Not an abstime column */ + for (row = 0; row < MIN (S->n_rows, 20); row++) { /* Maximum 20 rows are examined */ + sub = S->data[col][row] - floor (S->data[col][row]); /* Fractional second */ + if (gmt_M_is_zero (sub)) continue; /* Not checking zeros for width */ + if ((this_w = gmtio_get_precision_width (GMT, sub)) > w_max) + w_max = this_w; + } + } + if (w_max) { /* Need to append w_max "x" to the default format */ + strcat (GMT->current.setting.format_clock_out, "."); + while (w_max) { + strcat (GMT->current.setting.format_clock_out, "x"); + w_max--; + } + GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "To prevent loss of time-series precision we have changed FORMAT_CLOCK_OUT to %s\n", GMT->current.setting.format_clock_out); + gmtlib_clock_C_format (GMT, GMT->current.setting.format_clock_out, &GMT->current.io.clock_output, 1); + } +} + +void gmt_increase_abstime_format_precision (struct GMT_CTRL *GMT, unsigned int col, double dt) { + /* Adjust output format for absolute time if we need more precision in the seconds */ + double sub; + int w; + if (GMT->current.io.col_type[GMT_OUT][col] != GMT_IS_ABSTIME) return; /* Not an abstime column */ + if (GMT->current.setting.time_system.unit != 's') return; /* Current unit is not second */ + if (strcmp (GMT->current.setting.format_clock_out, "hh:mm:ss")) return; /* User has changed from default format - do nothing */ + sub = dt - floor (dt); /* Fractional second */ + if (gmt_M_is_zero (sub)) return; /* Not checking zeros for width */ + w = gmtio_get_precision_width (GMT, sub); /* Get desired precision */ + strcat (GMT->current.setting.format_clock_out, "."); + while (w) { + strcat (GMT->current.setting.format_clock_out, "x"); + w--; + } + GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "To prevent loss of time-series precision we have changed FORMAT_CLOCK_OUT to %s\n", GMT->current.setting.format_clock_out); + gmtlib_clock_C_format (GMT, GMT->current.setting.format_clock_out, &GMT->current.io.clock_output, 1); +} + /*! . */ int gmtlib_write_dataset (struct GMT_CTRL *GMT, void *dest, unsigned int dest_type, struct GMT_DATASET *D, bool use_GMT_io, int table) { /* Writes an entire data set to file or stream */ @@ -4591,6 +4667,8 @@ int gmtlib_write_dataset (struct GMT_CTRL *GMT, void *dest, unsigned int dest_ty strcpy (open_mode, (append) ? "a" : "w"); GMT->current.io.record_type[GMT_OUT] = D->type; + gmtio_check_abstime_format (GMT, D); /* If some columns are absolute time destined for ASCII formatting, we may need to change FORMAT_CLOCK_OUT template */ + /* Convert any destination type to stream */ switch (dest_type) { diff --git a/src/gmt_prototypes.h b/src/gmt_prototypes.h index b9d1cfde69b..b7983e2dde6 100644 --- a/src/gmt_prototypes.h +++ b/src/gmt_prototypes.h @@ -273,6 +273,7 @@ EXTERN_MSC int gmt_set_psfilename (struct GMT_CTRL *GMT); /* gmt_io.c: */ +EXTERN_MSC void gmt_increase_abstime_format_precision (struct GMT_CTRL *GMT, unsigned int col, double dt); EXTERN_MSC char *gmt_get_filename (struct GMTAPI_CTRL *API, const char* filename, const char *mods); EXTERN_MSC void gmt_quit_bad_record (struct GMTAPI_CTRL *API, struct GMT_RECORD *In); EXTERN_MSC int gmt_get_ogr_id (struct GMT_OGR *G, char *name); diff --git a/src/greenspline.c b/src/greenspline.c index d8a05a16b29..129dc230029 100644 --- a/src/greenspline.c +++ b/src/greenspline.c @@ -2036,6 +2036,9 @@ EXTERN_MSC int GMT_greenspline (void *V_API, int mode, void *args) { } #endif + + if (dimension == 1) gmt_increase_abstime_format_precision (GMT, GMT_X, Ctrl->I.inc[GMT_X]); /* In case we need more sub-second precision output */ + if (Ctrl->E.active) { /* Need to duplicate the data since SVD destroys it */ orig_obs = gmt_M_memory (GMT, NULL, nm, double); gmt_M_memcpy (orig_obs, obs, nm, double);