1 /*--------------------------------------------------------------------
2 *
3 * Copyright (c) 1991-2021 by the GMT Team (https://www.generic-mapping-tools.org/team.html)
4 * See LICENSE.TXT file for copying and redistribution conditions.
5 *
6 * This program is free software; you can redistribute it and/or modify
7 * it under the terms of the GNU Lesser General Public License as published by
8 * the Free Software Foundation; version 3 or any later version.
9 *
10 * This program is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 * GNU Lesser General Public License for more details.
14 *
15 * Contact info: www.generic-mapping-tools.org
16 *--------------------------------------------------------------------*/
17 /*
18 * Brief synopsis: grdinterpolate reads a 3D netcdf spatial data cube with
19 * the 3rd dimension either depth/height or time. It then interpolates
20 * the cube at arbitrary depth z (or time) values and writes either a single
21 * slice 2-D grid or another multi-level 3-D data cube. Alternatively,
22 * we can read a stack of input 2-D grids instead of the 3D cube. Finally,
23 * we may sample time-series (-S) or extract a vertical slice (-E) rather
24 * than write gridded horizontal output slice(s).
25 *
26 * Author: Paul Wessel
27 * Date: 1-AUG-2019
28 * Version: 6 API
29 *
30 */
31
32 #include "gmt_dev.h"
33
34 #define THIS_MODULE_CLASSIC_NAME "grdinterpolate"
35 #define THIS_MODULE_MODERN_NAME "grdinterpolate"
36 #define THIS_MODULE_LIB "core"
37 #define THIS_MODULE_PURPOSE "Interpolate a 3-D cube, 2-D grids or 1-D series from a 3-D data cube or stack of 2-D grids"
38 #define THIS_MODULE_KEYS "<G{+,G?},ED(,SD("
39 #define THIS_MODULE_NEEDS ""
40 #define THIS_MODULE_OPTIONS "-:>RVbdefghinoqs"
41
42 struct GRDINTERPOLATE_CTRL {
43 struct GRDINTERPOLATE_In {
44 bool active;
45 char **file;
46 unsigned int n_files;
47 } In;
48 struct GRDEDIT_D { /* -D[+x<xname>][+yyname>][+z<zname>][+v<zname>][+s<scale>][+o<offset>][+n<invalid>][+t<title>][+r<remark>] */
49 bool active;
50 char *information;
51 } D;
52 struct GRDINTERPOLATE_E { /* -E<file>|<line1>[,<line2>,...][+a<az>][+c][+g][+i<step>][+l<length>][+n<np][+o<az>][+p][+r<radius>][+x] */
53 bool active;
54 unsigned int mode;
55 char *lines;
56 double step;
57 char unit;
58 } E;
59 struct GRDINTERPOLATE_F { /* -Fl|a|c[1|2] */
60 bool active;
61 unsigned int mode;
62 unsigned int type;
63 char spline[GMT_LEN8];
64 } F;
65 struct GRDINTERPOLATE_G { /* -G<output_grdfile> */
66 bool active;
67 char *file;
68 } G;
69 struct GRDINTERPOLATE_S { /* -S<x>/<y>|<pointfile>[+h<header>] */
70 bool active;
71 double x, y;
72 char *file;
73 char *header;
74 } S;
75 struct GRDINTERPOLATE_T { /* -T<start>/<stop>/<inc> or -T<value> */
76 bool active;
77 struct GMT_ARRAY T;
78 char *string;
79 } T;
80 struct GRDINTERPOLATE_Z { /* -Z<min>/<max>/<inc>, -Z<file>, or -Z<list> */
81 bool active;
82 unsigned mode; /* 1 if need to create 0/n-1/1 array */
83 struct GMT_ARRAY T;
84 } Z;
85 };
86
New_Ctrl(struct GMT_CTRL * GMT)87 static void *New_Ctrl (struct GMT_CTRL *GMT) { /* Allocate and initialize a new control structure */
88 static char type[3] = {'l', 'a', 'c'};
89 struct GRDINTERPOLATE_CTRL *C;
90
91 C = gmt_M_memory (GMT, NULL, 1, struct GRDINTERPOLATE_CTRL);
92
93 /* Initialize values whose defaults are not 0/false/NULL */
94 sprintf (C->F.spline, "%c", type[GMT->current.setting.interpolant]); /* Set default interpolant */
95 return (C);
96 }
97
grdinterpolate_free_files(struct GMT_CTRL * GMT,char *** list,unsigned int n)98 GMT_LOCAL void grdinterpolate_free_files (struct GMT_CTRL *GMT, char ***list, unsigned int n) {
99 for (unsigned int k = 0; k < n; k++)
100 gmt_M_str_free ((*list)[k]);
101 gmt_M_free (GMT, *list);
102 }
103
Free_Ctrl(struct GMT_CTRL * GMT,struct GRDINTERPOLATE_CTRL * C)104 static void Free_Ctrl (struct GMT_CTRL *GMT, struct GRDINTERPOLATE_CTRL *C) { /* Deallocate control structure */
105 if (!C) return;
106 grdinterpolate_free_files (GMT, &(C->In.file), C->In.n_files);
107 gmt_M_str_free (C->D.information);
108 gmt_M_str_free (C->G.file);
109 gmt_M_str_free (C->S.file);
110 gmt_M_str_free (C->S.header);
111 gmt_M_str_free (C->T.string);
112 gmt_free_array (GMT, &(C->T.T));
113 gmt_free_array (GMT, &(C->Z.T));
114 gmt_M_free (GMT, C);
115 }
116
usage(struct GMTAPI_CTRL * API,int level)117 static int usage (struct GMTAPI_CTRL *API, int level) {
118 static char type[3] = {'l', 'a', 'c'};
119 const char *name = gmt_show_name_and_purpose (API, THIS_MODULE_LIB, THIS_MODULE_CLASSIC_NAME, THIS_MODULE_PURPOSE);
120 if (level == GMT_MODULE_PURPOSE) return (GMT_NOERROR);
121 GMT_Usage (API, 0, "usage: %s <cube> | <grd1> <grd2> <grd3> ... -G<outfile> [%s] "
122 "[-E<file>|<line1>[,<line2>,...][+a<az>][+g][+i<step>][+l<length>][+n<np][+o<az>][+p][+r<radius>][+x]] "
123 "[-Fl|a|c|n[+1|2]] [%s] [-S<x>/<y>|<table>[+h<header>]] [-T[<min>/<max>/]<inc>[+i|n]] [%s] "
124 "[-Z[<levels>]] [%s] [%s] [%s] [%s] [%s] [%s] [%s] [%s] [%s] [%s] [%s] [%s]\n",
125 name, GMT_GRDEDIT3D, GMT_Rgeo_OPT, GMT_V_OPT, GMT_b_OPT, GMT_d_OPT, GMT_e_OPT, GMT_f_OPT, GMT_g_OPT,
126 GMT_h_OPT, GMT_i_OPT, GMT_n_OPT, GMT_o_OPT, GMT_q_OPT, GMT_s_OPT, GMT_colon_OPT, GMT_PAR_OPT);
127
128 if (level == GMT_SYNOPSIS) return (GMT_MODULE_SYNOPSIS);
129
130 GMT_Message (API, GMT_TIME_NONE, " REQUIRED ARGUMENTS:\n");
131
132 GMT_Usage (API, 1, "\n<cube> is the name of the input 3D netCDF data cube. However, with -Z we instead expect a "
133 "series of 2-D grids.");
134 GMT_Usage (API, 1, "\n-G<outfile>");
135 GMT_Usage (API, -2, "Specify a single output file name (or a filename format template; also see -S) To write a "
136 "series of 2-D grids instead of a cube, include a floating-point C-format statement in <outfile> set via "
137 "-G for embedding the level in the file name.");
138 GMT_Message (API, GMT_TIME_NONE, "\n OPTIONAL ARGUMENTS:\n");
139 gmt_cube_info_syntax (API->GMT, 'D');
140 GMT_Usage (API, 1, "\n-E<file>|<line1>[,<line2>,...][+a<az>][+g][+i<step>][+l<length>][+n<np][+o<az>][+p][+r<radius>][+x]");
141 GMT_Usage (API, -2, "Set up a single cross-section based on <file> or on the given <line1>[,<line2>,...]. Give "
142 "start and stop coordinates for each line segment. The format of each <line> is <start>/<stop>, where "
143 "<start> or <stop> are coordinate pairs, e.g., <lon1/lat1>/<lon2>/<lat2>. Append +i<inc> to set the "
144 "sampling increment [Default is 0.5 x min of grid's (x_inc, y_inc)] Instead of <start/stop>, give <origin> "
145 "and append +a|o|l|n|r as required:");
146 GMT_Usage (API, -3, "+a Define a profiles from <origin> in <az> direction. Add +l<length>.");
147 GMT_Usage (API, -3, "+g Use gridline coordinates (degree longitude or latitude) if <line> is so aligned [great circle].");
148 GMT_Usage (API, -3, "+o Define a profile centered on <origin> in <az> direction. Add +l<length>.");
149 GMT_Usage (API, -3, "+p Sample along the parallel if <line> has constant latitude.");
150 GMT_Usage (API, -3, "+r Define a circle about <origin> with given <radius>. Add +i<inc> or +n<np>.");
151 GMT_Usage (API, -3, "+n Set the number of output points as <np> and computes <inc> from <length>.");
152 GMT_Usage (API, -3, "+x Follow a loxodrome (rhumbline) [great circle].");
153 GMT_Usage (API, -2, "Note: A unit is optional. Only ONE unit type from %s can be used throughout this option, so "
154 "mixing of units is not allowed [Default unit is km, if grid is geographic].");
155 GMT_Usage (API, 1, "\n-Fl|a|c|n][+1|2]");
156 GMT_Usage (API, -2, "Set the grid interpolation mode. Choose from:");
157 GMT_Usage (API, -3, "l: Linear interpolation.");
158 GMT_Usage (API, -3, "a: Akima spline interpolation.");
159 GMT_Usage (API, -3, "c: Cubic spline interpolation.");
160 GMT_Usage (API, -3, "n: No interpolation (nearest point).");
161 GMT_Usage (API, -2, "Optionally, append +1 for 1st derivative or +2 for 2nd derivative. [Default is -F%c].",
162 type[API->GMT->current.setting.interpolant]);
163 GMT_Option (API, "R");
164 GMT_Usage (API, 1, "\n-S<x>/<y>|<table>[+h<header>]");
165 GMT_Usage (API, -2, "Give a fixed point for across-stack sampling [Default] or interpolation [with -T]. For "
166 "multiple points, give a <table> of points instead (one point per record). Output is a multi-segment table "
167 "written to standard output unless -G is used to set a file name. To write each series to separate files, let "
168 "-G<outfile> contain a C-format integer specifier (e.g, %%d) for embedding the running point number. "
169 "Append a fixed header via +h<header> [trailing text per record in <table>].");
170 GMT_Usage (API, 1, "\n-T[<file>|<list>|<min>/<max>/<inc>[+b|i|l|n]]");
171 GMT_Usage (API, -2, "Interpolate the 3-D grid at given levels across the 3rd dimension. Make evenly spaced output "
172 "level steps from <min> to <max> by <inc>. Control setup via modifiers:");
173 GMT_Usage (API, 3, "+b Select log2 spacing in <inc>");
174 GMT_Usage (API, 3, "+i Indicate <inc> is the reciprocal of desired <inc> (e.g., 3 for 0.3333.....).");
175 GMT_Usage (API, 3, "+l Select log10 spacing via <inc> = 1,2,3.");
176 GMT_Usage (API, 3, "+n Let <inc> mean the number of points instead. of increment");
177 GMT_Usage (API, -2, "Alternatively, give a <file> with output times in the first column, or a comma-separated <list>.");
178 GMT_Option (API, "V");
179 GMT_Usage (API, 1, "\n-Z[<levels>]");
180 GMT_Usage (API, -2, "Read or write 2-D grids that make up a virtual 3-D data cube. To read a series of 2-D grids, "
181 "give -Z<levels>, where <levels> for each grid is set via <min>/<max>/<inc>, <zfile>, or a comma-separated "
182 "list. No argument means let levels be 0, 1, 2, ...");
183 GMT_Usage (API, -2, "Note: If -Z and no -T, -E, -S then we simply write the stack as a 3-D data cube.");
184 GMT_Option (API, "a,bi2,bo,d,e,f,g,h,i,n,o,q,s,:,.");
185
186 return (GMT_MODULE_USAGE);
187 }
188
parse(struct GMT_CTRL * GMT,struct GRDINTERPOLATE_CTRL * Ctrl,struct GMT_OPTION * options)189 static int parse (struct GMT_CTRL *GMT, struct GRDINTERPOLATE_CTRL *Ctrl, struct GMT_OPTION *options) {
190 /* This parses the options provided to grdcut and sets parameters in CTRL.
191 * Any GMT common options will override values set previously by other commands.
192 * It also replaces any file names specified as input or output with the data ID
193 * returned when registering these sources/destinations with the API.
194 */
195
196 unsigned int n_errors = 0, n_files = 0, n_alloc = 0, k = 0;
197 char *c = NULL;
198 struct GMT_OPTION *opt = NULL;
199 struct GMTAPI_CTRL *API = GMT->parent;
200
201 for (opt = options; opt; opt = opt->next) {
202 switch (opt->option) {
203
204 case '<': /* Input files */
205 if (n_alloc <= Ctrl->In.n_files) {
206 n_alloc += GMT_SMALL_CHUNK;
207 Ctrl->In.file = gmt_M_memory (GMT, Ctrl->In.file, n_alloc, char *);
208 }
209 Ctrl->In.file[Ctrl->In.n_files] = strdup (opt->arg);
210 if (GMT_Get_FilePath (API, GMT_IS_GRID, GMT_IN, GMT_FILE_REMOTE, &(Ctrl->In.file[Ctrl->In.n_files]))) n_errors++;
211 Ctrl->In.n_files++;
212 break;
213
214 /* Processes program-specific parameters */
215
216 case 'D': /* Give grid information */
217 n_errors += gmt_M_repeated_module_option (API, Ctrl->D.active);
218 Ctrl->D.active = true;
219 Ctrl->D.information = strdup (opt->arg);
220 break;
221 case 'E': /* Create or read an equidistant profile for slicing */
222 n_errors += gmt_M_repeated_module_option (API, Ctrl->E.active);
223 Ctrl->E.active = true;
224 Ctrl->E.lines = strdup (opt->arg);
225 break;
226 case 'F': /* Set the spline type */
227 n_errors += gmt_M_repeated_module_option (API, Ctrl->F.active);
228 Ctrl->F.active = true;
229 switch (opt->arg[0]) {
230 case 'l': Ctrl->F.mode = GMT_SPLINE_LINEAR; break;
231 case 'a': Ctrl->F.mode = GMT_SPLINE_AKIMA; break;
232 case 'c': Ctrl->F.mode = GMT_SPLINE_CUBIC; break;
233 case 'n': Ctrl->F.mode = GMT_SPLINE_NN; break;
234 default:
235 GMT_Report (API, GMT_MSG_ERROR, "Option -F: Bad spline selector %c\n", opt->arg[0]);
236 n_errors++;
237 break;
238 }
239 if (opt->arg[1] == '+') Ctrl->F.type = (opt->arg[2] - '0'); /* Want first or second derivative */
240 strncpy (Ctrl->F.spline, opt->arg, GMT_LEN8-1); /* Keep track of what was given since it may need to be passed verbatim to other modules */
241 break;
242 case 'G': /* Output file or name template */
243 n_errors += gmt_M_repeated_module_option (API, Ctrl->G.active);
244 if (n_files++ > 0) { n_errors++; continue; }
245 Ctrl->G.file = strdup (opt->arg);
246 if (strchr (Ctrl->G.file, '%') == NULL) { /* Gave a fixed output file, can check */
247 if (GMT_Get_FilePath (API, GMT_IS_GRID, GMT_OUT, GMT_FILE_LOCAL, &(Ctrl->G.file))) n_errors++;
248 }
249 break;
250 case 'T': /* Set level sampling spacing */
251 n_errors += gmt_M_repeated_module_option (API, Ctrl->T.active);
252 Ctrl->T.active = true;
253 Ctrl->T.string = strdup (opt->arg);
254 n_errors += gmt_parse_array (GMT, 'T', opt->arg, &(Ctrl->T.T), GMT_ARRAY_TIME | GMT_ARRAY_SCALAR | GMT_ARRAY_RANGE | GMT_ARRAY_UNIQUE, GMT_Z);
255 break;
256
257 case 'S': /* Sample vertically across the grid stack at one or more points */
258 n_errors += gmt_M_repeated_module_option (API, Ctrl->S.active);
259 Ctrl->S.active = true;
260 if ((c = strstr (opt->arg, "+h"))) { /* Got a fixed header string for output segment headers */
261 if (c[2])
262 Ctrl->S.header = strdup (&c[2]);
263 else {
264 GMT_Report (API, GMT_MSG_ERROR, "Option -S: Modifier +h not given a header string\n");
265 n_errors++;
266 }
267 c[0] = '\0'; /* Chop off all modifiers */
268 }
269 if (access (opt->arg, F_OK) && strchr (opt->arg, '/')) { /* Got a single point and not a valid path */
270 char txt_a[GMT_LEN256] = {""}, txt_b[GMT_LEN256] = {""};
271 if (sscanf (opt->arg, "%[^/]/%s", txt_a, txt_b) != 2) {
272 GMT_Report (API, GMT_MSG_ERROR, "Option -S: Cannot parse point coordinates from %s\n", opt->arg);
273 n_errors++;
274 }
275 else { /* OK to convert to numbers */
276 n_errors += gmt_verify_expectations (GMT, gmt_M_type (GMT, GMT_IN, GMT_X),
277 gmt_scanf_arg (GMT, txt_a, gmt_M_type (GMT, GMT_IN, GMT_X), false, &Ctrl->S.x), txt_a);
278 n_errors += gmt_verify_expectations (GMT, gmt_M_type (GMT, GMT_IN, GMT_Y),
279 gmt_scanf_arg (GMT, txt_b, gmt_M_type (GMT, GMT_IN, GMT_Y), false, &Ctrl->S.y), txt_b);
280 }
281 }
282 else {
283 Ctrl->S.file = strdup (opt->arg);
284 if (GMT_Get_FilePath (API, GMT_IS_DATASET, GMT_IN, GMT_FILE_REMOTE, &(Ctrl->S.file))) n_errors++;
285 }
286 if (c) c[0] = '+'; /* Restore modifiers */
287 break;
288
289 case 'Z': /* Control input/output grid management */
290 n_errors += gmt_M_repeated_module_option (API, Ctrl->Z.active);
291 Ctrl->Z.active = true;
292 if (opt->arg[0] == 'i') k = 1; /* Skip the now deprecated 'i' in -Zi */
293 if (opt->arg[k])
294 n_errors += gmt_parse_array (GMT, 'Z', &opt->arg[k], &(Ctrl->Z.T), GMT_ARRAY_TIME | GMT_ARRAY_SCALAR | GMT_ARRAY_RANGE | GMT_ARRAY_UNIQUE, GMT_Z);
295 else
296 Ctrl->Z.mode = 1;
297 break;
298
299 default: /* Report bad options */
300 n_errors += gmt_default_error (GMT, opt->option);
301 break;
302 }
303 }
304 if (Ctrl->In.n_files) Ctrl->In.file = gmt_M_memory (GMT, Ctrl->In.file, Ctrl->In.n_files + 1, char *); /* One extra so we have a NULL-terminated array */
305
306 n_errors += gmt_M_check_condition (GMT, Ctrl->In.n_files < 1, "Error: No input grid(s) specified.\n");
307 n_errors += gmt_M_check_condition (GMT, !Ctrl->Z.active && Ctrl->In.n_files != 1, "Must specify a single input 3D grid cube file unless -Z is set\n");
308 n_errors += gmt_M_check_condition (GMT, Ctrl->F.type > 2, "Option -F: Only 1st or 2nd derivatives may be requested\n");
309 if (!(Ctrl->S.active || Ctrl->E.active)) { /* Under -S and -E, the -T and -G are optional */
310 n_errors += gmt_M_check_condition (GMT, !Ctrl->G.file, "Option -G: Must specify output grid file\n");
311 n_errors += gmt_M_check_condition (GMT, n_files != 1, "Must specify only one output file name\n");
312 }
313 n_errors += gmt_M_check_condition (GMT, Ctrl->E.active && Ctrl->S.active, "Options -E and -S cannot be used together\n");
314 if (Ctrl->E.active) {
315 n_errors += gmt_M_check_condition (GMT, strstr (Ctrl->E.lines, "+d"), "Option -E: Unrecognized modifier +d\n");
316 }
317 n_errors += gmt_M_check_condition (GMT, Ctrl->Z.active && !(Ctrl->T.active || Ctrl->E.active || Ctrl->S.active) && strchr (Ctrl->G.file, '%'), "Options -Z: If -T not given the we must write a single data cube\n");
318 n_errors += gmt_M_check_condition (GMT, Ctrl->G.active && API->external && strchr (Ctrl->G.file, '%'), "Option -G: Cannot contain format-specifiers when not used on the command line\n");
319
320 return (n_errors ? GMT_PARSE_ERROR : GMT_NOERROR);
321 }
322
grdinterpolate_equidistant_levels(struct GMT_CTRL * GMT,double * z,unsigned int nz)323 GMT_LOCAL bool grdinterpolate_equidistant_levels (struct GMT_CTRL *GMT, double *z, unsigned int nz) {
324 /* Return true if spacing between layers is constant */
325 unsigned int k;
326 double dz;
327
328 if (gmt_M_is_verbose (GMT, GMT_MSG_INFORMATION)) { /* Report levels */
329 for (k = 0; k < nz; k++)
330 GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Level %d: %g\n", k, z[k]);
331 }
332 if (nz < 3) return true; /* Special case of a single layer */
333 dz = z[1] - z[0]; /* Get first increment, then compare to the rest */
334 for (k = 2; k < nz; k++)
335 if (!doubleAlmostEqual (dz, z[k]-z[k-1])) return false;
336 return true;
337 }
338
339 #define bailout(code) {gmt_M_free_options (mode); return (code);}
340 #define Return(code) {Free_Ctrl (GMT, Ctrl); gmt_end_module (GMT, GMT_cpy); bailout (code);}
341
GMT_grdinterpolate(void * V_API,int mode,void * args)342 EXTERN_MSC int GMT_grdinterpolate (void *V_API, int mode, void *args) {
343 char file[PATH_MAX] = {""}, cube_layer[GMT_LEN64] = {""}, *nc_z_named = NULL;
344 bool equi_levels, convert_to_cube = false, z_is_abstime = false, got_cube = false;
345 int error = 0;
346 unsigned int int_mode, row, col, level_type, dtype = 0;
347 uint64_t n_layers = 0, k, node, start_k, stop_k, n_layers_used, *this_dim = NULL, dims[3] = {0, 0, 0};
348 double wesn[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, inc[3] = {0.0, 0.0, 0.0}, w_range[2] = {0.0, 0.0};
349 double *level = NULL, *i_value = NULL, *o_value = NULL;
350 struct GMT_GRID *Grid = NULL;
351 struct GMT_CUBE *C[2] = {NULL, NULL}; /* Structures to hold input/output cubes */
352
353 struct GMT_DATASET *In = NULL, *Out = NULL;
354 struct GRDINTERPOLATE_CTRL *Ctrl = NULL;
355 struct GMT_CTRL *GMT = NULL, *GMT_cpy = NULL;
356 struct GMT_OPTION *options = NULL;
357 struct GMTAPI_CTRL *API = gmt_get_api_ptr (V_API); /* Cast from void to GMTAPI_CTRL pointer */
358
359 /*----------------------- Standard module initialization and parsing ----------------------*/
360
361 if (API == NULL) return (GMT_NOT_A_SESSION);
362 if (mode == GMT_MODULE_PURPOSE) return (usage (API, GMT_MODULE_PURPOSE)); /* Return the purpose of program */
363 options = GMT_Create_Options (API, mode, args); if (API->error) return (API->error); /* Set or get option list */
364
365 if ((error = gmt_report_usage (API, options, 0, usage)) != GMT_NOERROR) bailout (error); /* Give usage if requested */
366
367 /* Parse the command-line arguments */
368
369 if ((GMT = gmt_init_module (API, THIS_MODULE_LIB, THIS_MODULE_CLASSIC_NAME, THIS_MODULE_KEYS, THIS_MODULE_NEEDS, NULL, &options, &GMT_cpy)) == NULL) bailout (API->error); /* Save current state */
370 if (GMT_Parse_Common (API, THIS_MODULE_OPTIONS, options)) Return (API->error);
371 Ctrl = New_Ctrl (GMT); /* Allocate and initialize a new control structure */
372 if ((error = parse (GMT, Ctrl, options)) != 0) Return (error);
373
374 /*---------------------------- This is the grdinterpolate main code ----------------------------*/
375
376 gmt_grd_set_datapadding (GMT, true); /* Turn on gridpadding when reading a subset */
377
378 if (Ctrl->Z.active) { /* Create the input level array */
379 if (Ctrl->Z.mode == 1) /* No levels given, auto-make a list from 0, 1, ... */
380 n_layers = gmt_make_equidistant_array (GMT, 0.0, Ctrl->In.n_files-1.0, 1.0, &level);
381 else { /* Gave information to build a level array */
382 if (gmt_create_array (GMT, 'Z', &(Ctrl->Z.T), NULL, NULL)) {
383 GMT_Report (API, GMT_MSG_ERROR, "Option -Z: Unable to set up input level array\n");
384 Return (API->error);
385 }
386 if (Ctrl->In.n_files != Ctrl->Z.T.n) {
387 GMT_Report (API, GMT_MSG_ERROR, "Option -Z: Number of input 2-D grids does not match number of levels given via -Z\n");
388 Return (API->error);
389 }
390 n_layers = Ctrl->Z.T.n; /* Set number of layers anticipated */
391 level = Ctrl->Z.T.array; /* Pointer to allocated array with the level values */
392 z_is_abstime = Ctrl->Z.T.temporal; /* In case we parsed abs time for levels */
393 }
394 if (!(Ctrl->T.active || Ctrl->E.active || Ctrl->S.active)) convert_to_cube = true; /* Just want to build cube from input stack */
395 }
396 else if (gmt_M_file_is_memory (Ctrl->In.file[0])) { /* Got a memory reference */
397 if (Ctrl->In.n_files == 1) { /* Got one memory reference so it must be a cube */
398 if (Ctrl->In.file[0][GMTAPI_OBJECT_FAMILY_START] != 'U') {
399 GMT_Report (API, GMT_MSG_ERROR, "Input memory reference is not a cube!\n");
400 Return (API->error);
401 }
402 C[GMT_IN] = GMT_Read_VirtualFile (API, Ctrl->In.file[0]);
403 n_layers = C[GMT_IN]->header->n_bands;
404 level = C[GMT_IN]->z;
405 got_cube = true;
406 }
407 else { /* If many references given then clearly we are missing -Z */
408 GMT_Report (API, GMT_MSG_ERROR, "Multiple grid memory references requires -Z!\n");
409 Return (API->error);
410 }
411 }
412 else { /* See if we got a 3D netCDF data cube; if so return number of layers and and the levels array */
413 nc_z_named = strchr (Ctrl->In.file[0], '?'); /* Maybe given a specific variable? */
414 if (nc_z_named) { /* Gave a specific variable. Keep variable name and remove from filename */
415 strcpy (cube_layer, &nc_z_named[1]);
416 nc_z_named[0] = '\0'; /* Chop off layer name for now */
417 }
418 if ((error = gmt_nc_read_cube_info (GMT, Ctrl->In.file[0], w_range, &n_layers, &level, NULL))) {
419 Return (error);
420 }
421 got_cube = true;
422 }
423
424 if (n_layers == 1) {
425 GMT_Report (API, GMT_MSG_ERROR, "Only one layer given - need at least two to interpolate across levels\n");
426 Return (GMT_RUNTIME_ERROR);
427 }
428
429 if (got_cube && !(Ctrl->E.active || Ctrl->S.active || Ctrl->T.active)) { /* If given a cube we requires other options */
430 GMT_Report (API, GMT_MSG_ERROR, "Data cube read but none of -E, -S -T were given\n");
431 Return (GMT_RUNTIME_ERROR);
432 }
433 /* Create output level array, if selected */
434 if (Ctrl->T.active && gmt_create_array (GMT, 'T', &(Ctrl->T.T), NULL, NULL)) {
435 GMT_Report (API, GMT_MSG_ERROR, "Option -T: Unable to set up output level array\n");
436 Return (API->error);
437 }
438 equi_levels = grdinterpolate_equidistant_levels (GMT, level, n_layers); /* Are levels equidistant? */
439 level_type = gmt_M_type (GMT, GMT_IN, GMT_Z); /* Either time or floating point values like depth */
440
441 /* Determine the range of input layers needed for interpolation */
442
443 start_k = 0; stop_k = n_layers - 1; /* We first assume all layers are needed */
444 if (Ctrl->T.active) {
445 if (Ctrl->T.T.array[0] > level[stop_k] || Ctrl->T.T.array[Ctrl->T.T.n-1] < level[start_k]) {
446 GMT_Report (API, GMT_MSG_ERROR, "Option -T: Specified range outside that of the data cube\n");
447 Return (GMT_RUNTIME_ERROR);
448 }
449 while (start_k < n_layers && Ctrl->T.T.array[0] > level[start_k]) /* Find the first layer that is inside the output time range */
450 start_k++;
451 if (start_k && Ctrl->T.T.array[0] < level[start_k]) start_k--; /* Go back one if start time is less than first layer */
452 if (start_k && (Ctrl->F.mode == GMT_SPLINE_AKIMA || Ctrl->F.mode == GMT_SPLINE_CUBIC))
453 start_k--; /* One more to define the spline coefficients */
454 while (stop_k && Ctrl->T.T.array[Ctrl->T.T.n-1] < level[stop_k]) /* Find the last layer that is inside the output time range */
455 stop_k--;
456 if (stop_k < n_layers && Ctrl->T.T.array[Ctrl->T.T.n-1] > level[stop_k]) stop_k++; /* Go forward one if stop time is larger than last layer */
457 if (stop_k < (n_layers-1) && (Ctrl->F.mode == GMT_SPLINE_AKIMA || Ctrl->F.mode == GMT_SPLINE_CUBIC))
458 stop_k++; /* One more to define the spline coefficients */
459 }
460 n_layers_used = stop_k - start_k + 1; /* Total number of input layers needed */
461 if (n_layers_used == 1) { /* Might have landed exactly on one of the grid levels, but gmt_intpol needs at least 2 inputs */
462 if (start_k) start_k--;
463 else stop_k++; /* We know there are at least 2 input grids at this point in the code */
464 n_layers_used = 2;
465 }
466
467 if (Ctrl->E.active) { /* Create or read profile */
468 if (!gmt_access (GMT, Ctrl->E.lines, F_OK)) { /* Gave a file with a pre-computed crossection */
469 struct GMT_DATASEGMENT *Si = NULL;
470 if ((In = GMT_Read_Data (API, GMT_IS_DATASET, GMT_IS_FILE, GMT_IS_LINE, GMT_READ_NORMAL, NULL, Ctrl->E.lines, NULL)) == NULL) {
471 GMT_Report (API, GMT_MSG_ERROR, "Unable to open or read file %s.\n", Ctrl->E.lines);
472 Return (API->error);
473 }
474 if (In->n_segments > 1) {
475 GMT_Report (API, GMT_MSG_ERROR, "File %s contains more than one segment\n", Ctrl->E.lines);
476 Return (API->error);
477 }
478 if (In->n_columns < 2) {
479 GMT_Report (API, GMT_MSG_ERROR, "File %s does not contain coordinate pairs\n", Ctrl->E.lines);
480 Return (API->error);
481 }
482 Si = In->table[0]->segment[0]; /* Short-hand to only segment */
483 if (In->n_columns < 3) { /* Need to create distance column array */
484 gmt_adjust_dataset (GMT, In, 3); /* Add one more output column */
485 gmt_M_free (GMT, Si->data[GMT_Z]); /* Free it so we can add it next */
486 Si->data[GMT_Z] = gmt_dist_array (GMT, Si->data[GMT_X], Si->data[GMT_Y], Si->n_rows, true);
487 }
488 if (!grdinterpolate_equidistant_levels (GMT, Si->data[GMT_Z], Si->n_rows)) {
489 GMT_Report (API, GMT_MSG_ERROR, "File %s does not contain equidistant coordinates\n", Ctrl->E.lines);
490 Return (API->error);
491 }
492 }
493 else { /* Create profile */
494 char prof_args[GMT_LEN128] = {""};
495 if (!(equi_levels || Ctrl->T.active)) {
496 GMT_Report (API, GMT_MSG_ERROR, "Option -E requires either equidistant levels or resampling via -T\n");
497 Return (GMT_RUNTIME_ERROR);
498 }
499 if (gmt_get_dist_units (GMT, Ctrl->E.lines, &Ctrl->E.unit, &Ctrl->E.mode)) { /* Bad mixing of units in -E specification */
500 GMT_Report (API, GMT_MSG_ERROR, "Cannot mix units in -E (%s)\n", Ctrl->E.lines);
501 Return (GMT_RUNTIME_ERROR);
502 }
503 if (gmt_init_distaz (GMT, Ctrl->E.unit, Ctrl->E.mode, GMT_MAP_DIST) == GMT_NOT_A_VALID_TYPE) /* Initialize the distance unit and scaling */
504 Return (GMT_NOT_A_VALID_TYPE);
505
506 /* Need to get dx,dy from one grid */
507 if (Ctrl->Z.active) /* Get the first file */
508 sprintf (file, "%s", Ctrl->In.file[0]);
509 else /* Get the first layer from 3D cube possibly via a selected variable */
510 sprintf (file, "%s?%s[0]", Ctrl->In.file[0], cube_layer);
511 if ((Grid = GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_ONLY, NULL, file, NULL)) == NULL) {
512 GMT_Report (API, GMT_MSG_ERROR, "Unable to read header from file %s.\n", file);
513 Return (API->error);
514 }
515 /* Set default spacing to half the min grid spacing: */
516 Ctrl->E.step = 0.5 * MIN (Grid->header->inc[GMT_X], Grid->header->inc[GMT_Y]);
517 if (GMT_Destroy_Data (API, &Grid)) {
518 GMT_Report (API, GMT_MSG_ERROR, "Unable to destroy temporary grid?\n");
519 Return (API->error);
520 }
521 if (gmt_M_is_geographic (GMT, GMT_IN)) { /* Convert to km if geographic or degrees if arc-units */
522 if (!GMT->current.map.dist[GMT_MAP_DIST].arc) Ctrl->E.step *= GMT->current.proj.DIST_M_PR_DEG; /* Convert from degrees to meters or from min/secs to degrees */
523 Ctrl->E.step *= GMT->current.map.dist[GMT_MAP_DIST].scale; /* Scale to chosen unit */
524 GMT_Report (API, GMT_MSG_INFORMATION, "Default sampling interval in -E is %g %c (may be overridden by -E modifiers).\n", Ctrl->E.step, Ctrl->E.unit);
525 }
526 if (strstr (Ctrl->E.lines, "+c"))
527 strncpy (prof_args, Ctrl->E.lines, GMT_LEN128-1);
528 else /* Make sure we request continuous if possible */
529 sprintf (prof_args, "%s+c", Ctrl->E.lines);
530 In = gmt_make_profiles (GMT, 'E', prof_args, true, false, true, Ctrl->E.step, 0, NULL, &dtype);
531 if (In->table[0] == NULL)
532 Return (GMT_RUNTIME_ERROR);
533 }
534 }
535 else if (Ctrl->S.active) { /* Create time/depth-series and not grid output */
536 /* Since we let grdtrack read the grids we do a separate branch here and the return from the module */
537 uint64_t seg, row;
538 uint64_t dim[4] = {1, 1, 1, 2}; /* Dataset dimension for one point */
539 char header[GMT_LEN256] = {""};
540 struct GMT_DATASEGMENT *Si = NULL;
541
542 if (Ctrl->S.file) { /* Read a list of points, the list may have trailing text which may be activated by -S...+h */
543 if ((In = GMT_Read_Data (API, GMT_IS_DATASET, GMT_IS_FILE, GMT_IS_POINT, GMT_READ_NORMAL, NULL, Ctrl->S.file, NULL)) == NULL) {
544 GMT_Report (API, GMT_MSG_ERROR, "Unable to open or read file %s.\n", Ctrl->S.file);
545 Return (API->error);
546 }
547 if (Ctrl->S.header) { /* Want to use this fixed text to add to the trailing text of all the points */
548 for (seg = 0; seg < In->n_segments; seg++) {
549 Si = In->table[0]->segment[seg]; /* Short hand to this segment */
550 if (Si->text == NULL) /* Input file did not have any trailing text so add array now */
551 Si->text = gmt_M_memory (GMT, NULL, Si->n_rows, char *);
552 for (row = 0; row < Si->n_rows; row++) {
553 if (Si->text[row]) { /* Already has trailing text, combine with user argument */
554 sprintf (header, "%s %s", Si->text[row], Ctrl->S.header);
555 gmt_M_str_free (Si->text[row]);
556 Si->text[row] = strdup (header);
557 }
558 else /* Just use user argument */
559 Si->text[row] = strdup (Ctrl->S.header);
560 }
561 }
562 }
563 }
564 else { /* Single point, simplify logic below by making a 1-point dataset */
565 if ((In = GMT_Create_Data (API, GMT_IS_DATASET, GMT_IS_POINT, (Ctrl->S.header) ? GMT_WITH_STRINGS : 0, dim, NULL, NULL, 0, 0, NULL)) == NULL) {
566 GMT_Report (API, GMT_MSG_ERROR, "Unable to open or read file %s.\n", Ctrl->S.file);
567 Return (API->error);
568 }
569 Si = In->table[0]->segment[0]; /* Short hand to the first and only segment */
570 Si->data[GMT_X][0] = Ctrl->S.x;
571 Si->data[GMT_Y][0] = Ctrl->S.y;
572 if (Ctrl->S.header) /* Set the fixed header here via training text */
573 Si->text[0] = strdup (Ctrl->S.header);
574 Si->n_rows = 1;
575 gmt_set_dataset_minmax (GMT, In);
576 }
577 }
578
579 if (Ctrl->E.active || Ctrl->S.active) { /* Vertical profiles or slice */
580 unsigned int io_mode = GMT_WRITE_NORMAL;
581 uint64_t seg, row, rec, col;
582 uint64_t dim[4] = {1, 1, 1, 2}; /* Dataset dimension for one point */
583 char i_file[GMT_VF_LEN] = {""}, o_file[GMT_VF_LEN] = {""}, grid[GMT_LEN128] = {""}, header[GMT_LEN256] = {""}, cmd[GMT_LEN128] = {""};
584 struct GMT_DATASET *D = NULL;
585 struct GMT_DATASEGMENT *Si = NULL, *So = NULL;
586
587 dim[GMT_SEG] = In->n_records; /* One output time-series per input data location */
588 dim[GMT_ROW] = n_layers_used; /* Length of each sampled time-series per input data location */
589 dim[GMT_COL] = In->n_columns + 2; /* x, y [,col3, col4...], time|z, value */
590 if ((Out = GMT_Create_Data (API, GMT_IS_DATASET, GMT_IS_LINE, GMT_WITH_STRINGS, dim, NULL, NULL, 0, 0, NULL)) == NULL) {
591 GMT_Report (API, GMT_MSG_ERROR, "Unable to create output dataset for time-series\n");
592 Return (API->error);
593 }
594 if (GMT_Open_VirtualFile (API, GMT_IS_DATASET, GMT_IS_POINT, GMT_IN|GMT_IS_REFERENCE, In, i_file) != GMT_NOERROR) {
595 GMT_Report (API, GMT_MSG_ERROR, "Unable to create virtual dataset for time-series\n");
596 Return (API->error);
597 }
598
599 gmt_disable_bghio_opts (GMT); /* Do not want any -b -g -h -i -o to affect the workings of grdtrack calls */
600
601 for (k = start_k; k <= stop_k; k++) { /* For all selected input levels k */
602 GMT_Init_VirtualFile (API, 0, i_file); /* Reset so it can be read again */
603 if (Ctrl->Z.active) /* Get the k'th file */
604 sprintf (grid, "%s", Ctrl->In.file[k]);
605 else /* Get the k'th layer from 3D cube */
606 sprintf (grid, "%s?%s[%" PRIu64 "]", Ctrl->In.file[0], cube_layer, k);
607 if (GMT_Open_VirtualFile (API, GMT_IS_DATASET, GMT_IS_POINT, GMT_OUT|GMT_IS_REFERENCE, NULL, o_file) == GMT_NOTSET) {
608 GMT_Report (API, GMT_MSG_ERROR, "Unable to create virtual dataset for time-series\n");
609 Return (API->error);
610 }
611
612 sprintf (cmd, "%s -G%s ->%s", i_file, grid, o_file);
613 if (GMT->common.R.active[RSET]) { /* Gave a subregion, so pass -R along */
614 strcat (cmd, " -R");
615 strcat (cmd, GMT->common.R.string);
616 }
617 GMT_Report (API, GMT_MSG_DEBUG, "Sampling the grid: %s\n", cmd);
618 if (GMT_Call_Module (API, "grdtrack", GMT_MODULE_CMD, cmd) != GMT_NOERROR) { /* Sample this layer at the points */
619 Return (API->error);
620 }
621 if ((D = GMT_Read_VirtualFile (API, o_file)) == NULL) { /* Get the results from grdtrack */
622 GMT_Report (API, GMT_MSG_ERROR, "Unable to read virtual dataset for time-series created by grdtrack\n");
623 Return (API->error);
624 }
625 if (D->n_tables == 0 || D->table[0]->n_segments == 0) {
626 GMT_Report (API, GMT_MSG_ERROR, "No data time-series created by grdtrack\n");
627 Return (API->error);
628 }
629
630 for (seg = rec = 0; seg < D->table[0]->n_segments; seg++) { /* For each point we sampled at */
631 Si = D->table[0]->segment[seg]; /* Short hand to this segment */
632
633 for (row = 0; row < Si->n_rows; row++, rec++) { /* For each selected point which matches each output segment */
634 So = Out->table[0]->segment[rec]; /* Short hand to this output segment */
635 if (k == start_k) { /* Set the segment header just once */
636 if (Si->text && Si->text[row])
637 sprintf (header, "Location %g,%g %s", Si->data[GMT_X][row], Si->data[GMT_Y][row], Si->text[row]);
638 else
639 sprintf (header, "Location %g,%g", Si->data[GMT_X][row], Si->data[GMT_Y][row]);
640 So->header = strdup (header);
641 }
642 if (Ctrl->S.active) { /* Want x,y,z[,.....],value output */
643 for (col = 0; col < GMT_Z; col++) /* Copy over x,y */
644 So->data[col][k] = Si->data[col][row];
645 So->data[GMT_Z][k] = level[k]; /* Add level as the z column */
646 while (col < Si->n_columns) { /* Copy over the rest */
647 So->data[col+1][k] = Si->data[col][row];
648 col++;
649 }
650 }
651 else { /* Format for -E is x,y[,....],value */
652 for (col = 0; col < Si->n_columns; col++) /* Copy over the various columns */
653 So->data[col][k] = Si->data[col][row];
654 So->data[col][k] = level[k]; /* Add level as the last data column */
655 }
656 }
657 }
658 for (col = GMT_Z; col < Si->n_columns; col++)
659 gmt_set_column_type (GMT, GMT_OUT, col, GMT_IS_FLOAT); /* These are data columns */
660 if (GMT_Close_VirtualFile (API, i_file) != GMT_NOERROR) {
661 GMT_Report (API, GMT_MSG_ERROR, "Unable to close input virtual dataset for time-series\n");
662 Return (API->error);
663 }
664 if (GMT_Close_VirtualFile (API, o_file) != GMT_NOERROR) {
665 GMT_Report (API, GMT_MSG_ERROR, "Unable to close output virtual dataset for time-series\n");
666 Return (API->error);
667 }
668 if (GMT_Destroy_Data (API, &D) != GMT_OK) {
669 GMT_Report (API, GMT_MSG_ERROR, "Unable to delete dataset returned by grdtrack\n");
670 Return (API->error);
671 }
672 }
673 gmt_set_column_type (GMT, GMT_OUT, col, level_type); /* This is the grid-level data type which on output is in this column */
674
675 if (Ctrl->T.active) { /* Want to interpolate through the sampled points using the specified spline */
676 if (GMT_Open_VirtualFile (API, GMT_IS_DATASET, GMT_IS_LINE, GMT_IN|GMT_IS_REFERENCE, Out, i_file) != GMT_NOERROR) {
677 GMT_Report (API, GMT_MSG_ERROR, "Unable to create virtual dataset for sampled time-series\n");
678 Return (API->error);
679 }
680 if (GMT_Open_VirtualFile (API, GMT_IS_DATASET, GMT_IS_LINE, GMT_OUT|GMT_IS_REFERENCE, NULL, o_file) == GMT_NOTSET) {
681 GMT_Report (API, GMT_MSG_ERROR, "Unable to create virtual dataset for sampled time-series\n");
682 Return (API->error);
683 }
684 sprintf (cmd, "%s -F%s -N%d -T%s ->%s", i_file, Ctrl->F.spline, (int)(Out->n_columns - 1), Ctrl->T.string, o_file);
685 if (GMT_Call_Module (API, "sample1d", GMT_MODULE_CMD, cmd) != GMT_NOERROR) { /* Interpolate each profile per -T */
686 Return (API->error);
687 }
688 if (GMT_Destroy_Data (API, &Out) != GMT_OK) {
689 GMT_Report (API, GMT_MSG_ERROR, "Unable to delete virtual dataset for time-series\n");
690 Return (API->error);
691 }
692 if ((Out = GMT_Read_VirtualFile (API, o_file)) == NULL) { /* Get the updated results from sample1d */
693 GMT_Report (API, GMT_MSG_ERROR, "Unable to read virtual dataset for time-series created by sample1d\n");
694 Return (API->error);
695 }
696 dim[GMT_ROW] = Out->table[0]->segment[0]->n_rows; /* Update new row dim */
697 }
698
699 gmt_reenable_bghio_opts (GMT); /* Recover settings provided by user (if -b -g -h -i were used at all) */
700
701 if (Ctrl->S.active) { /* Write the table(s) */
702 if (Ctrl->G.file && strchr (Ctrl->G.file, '%')) { /* Want separate files per series, so change mode and build file names per segment */
703 struct GMT_DATASEGMENT_HIDDEN *SH = NULL;
704 io_mode = GMT_WRITE_SEGMENT; /* Flag that we want individual segment files */
705 for (seg = 0; seg < Out->table[0]->n_segments; seg++) { /* Set the output file names */
706 SH = gmt_get_DS_hidden (Out->table[0]->segment[seg]);
707 sprintf (file, Ctrl->G.file, seg);
708 SH->file[GMT_OUT] = strdup (file);
709 }
710 }
711 /* Time to write out the results to stdout, a file, or individual segment files */
712 gmt_set_segmentheader (GMT, GMT_OUT, true); /* Here we always want to write the headers we built */
713 if (GMT_Write_Data (API, GMT_IS_DATASET, GMT_IS_FILE, GMT_IS_LINE, io_mode, NULL, Ctrl->G.file, Out) != GMT_NOERROR) {
714 GMT_Report (API, GMT_MSG_ERROR, "Unable to write dataset for time-series\n");
715 Return (API->error);
716 }
717 }
718 else { /* Here we return the slice as a grid */
719 uint64_t ij;
720 char unit[GMT_LEN64] = {""};
721 double wesn[4], inc[2];
722
723 wesn[XLO] = Out->table[0]->segment[0]->data[2][0]; wesn[XHI] = Out->table[0]->segment[dim[GMT_SEG]-1]->data[2][0];
724 if (dtype == GMT_IS_LON && wesn[XHI] < wesn[XLO]) wesn[XHI] += 360.0; /* Watch out for 360 wraps */
725 wesn[YLO] = Out->table[0]->segment[0]->data[4][0]; wesn[YHI] = Out->table[0]->segment[0]->data[4][dim[GMT_ROW]-1];
726 inc[GMT_X] = gmt_M_get_inc (GMT, wesn[XLO], wesn[XHI], Out->n_segments, GMT_GRID_NODE_REG);
727 inc[GMT_Y] = gmt_M_get_inc (GMT, wesn[YLO], wesn[YHI], Out->table[0]->segment[0]->n_rows, GMT_GRID_NODE_REG);
728 gmt_set_column_type (GMT, GMT_OUT, GMT_X, dtype); /* This is data type of the distances */
729 gmt_set_column_type (GMT, GMT_OUT, GMT_Y, GMT_IS_FLOAT); /* This is data type of depth or time */
730 gmt_set_column_type (GMT, GMT_OUT, GMT_Z, GMT_IS_FLOAT); /* This is data type of the cube values */
731 if ((Grid = GMT_Create_Data (API, GMT_IS_GRID, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA, NULL, wesn, inc, GMT_GRID_NODE_REG, GMT_NOTSET, NULL)) == NULL) {
732 GMT_Report (API, GMT_MSG_ERROR, "Unable to create a grid for slice\n");
733 Return (API->error);
734 }
735 if (GMT_Set_Comment (API, GMT_IS_GRID, GMT_COMMENT_IS_REMARK, "Grid slide through 3-D data cube", Grid)) Return (API->error);
736 if (GMT_Set_Comment (API, GMT_IS_GRID, GMT_COMMENT_IS_OPTION | GMT_COMMENT_IS_COMMAND, options, Grid)) Return (API->error);
737 if (dtype == GMT_IS_LON)
738 strcpy (Grid->header->x_units, "longitude [degrees_east]");
739 else if (dtype == GMT_IS_LAT)
740 strcpy (Grid->header->x_units, "latitude [degrees_north]");
741 else {
742 sprintf (unit, "Distance (%c)", Ctrl->E.unit);
743 strcpy (Grid->header->x_units, unit);
744 }
745
746 for (seg = 0; seg < Out->n_segments; seg++) { /* Each segment represents one x-coordinate */
747 So = Out->table[0]->segment[seg]; /* Short hand to this output segment */
748 for (row = 0; row < So->n_rows; row++) {
749 ij = gmt_M_ijp (Grid->header, So->n_rows-row-1, seg); /* Must flip order since rows in grid goes down */
750 Grid->data[ij] = (float)So->data[3][row];
751 }
752 }
753 if (GMT_Destroy_Data (API, &Out) != GMT_OK) {
754 GMT_Report (API, GMT_MSG_ERROR, "Unable to delete final dataset for time-series\n");
755 Return (API->error);
756 }
757 GMT_Report (API, GMT_MSG_INFORMATION, "Write sliced grid to output file %s\n", Ctrl->G.file);
758 if (GMT_Write_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA, NULL, Ctrl->G.file, Grid) != GMT_NOERROR) {
759 Return (API->error);
760 }
761 }
762 if (!Ctrl->Z.active)
763 gmt_M_free (GMT, level);
764 Return (GMT_NOERROR);
765 }
766
767 /* Get here if neither -E nor -S were selected: We want to interpolate for one or more horizontal slices in the cube and need to read/write cubes */
768
769 int_mode = gmt_set_interpolate_mode (GMT, Ctrl->F.mode, Ctrl->F.type); /* What mode we pass to the interpolation */
770
771 if (GMT->common.R.active[RSET]) /* Use current -R setting for subsets, if given */
772 gmt_M_memcpy (wesn, GMT->common.R.wesn, 4, double);
773 wesn[ZLO] = level[start_k]; wesn[ZHI] = level[stop_k]; /* Then add the zmin/zmax limitation */
774
775 GMT_Report (API, GMT_MSG_INFORMATION, "Will read %" PRIu64 " layers (%" PRIu64 " - %" PRIu64 ") for levels %g to %g.\n", n_layers_used, start_k, stop_k, level[start_k], level[stop_k]);
776
777 /* Read the selected subset of the cube into C[GMT_IN] */
778
779 if (Ctrl->Z.active) { /* Need to read in individual grids and convert to cube first */
780 size_t here = 0;
781 unsigned int N = n_layers;
782 struct GMT_GRID **G = NULL;
783
784 if ((G = GMT_Read_Group (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA, wesn, Ctrl->In.file, &N, NULL)) == NULL)
785 Return (EXIT_FAILURE);
786 if (!GMT->common.R.active[RSET]) /* Use current -R setting for subsets, if given */
787 gmt_M_memcpy (wesn, G[0]->header->wesn, 4, double);
788 dims[GMT_Z] = n_layers; /* Number of input levels */
789 gmt_M_memcpy (inc, G[0]->header->inc, 2U, double);
790 if ((C[GMT_IN] = GMT_Create_Data (API, GMT_IS_CUBE, GMT_IS_VOLUME, GMT_CONTAINER_AND_DATA, dims, wesn, inc, G[0]->header->registration, GMT_NOTSET, NULL)) == NULL)
791 Return (GMT_MEMORY_ERROR);
792 for (k = 0; k < n_layers; k++, here += G[0]->header->size)
793 gmt_M_memcpy (&(C[GMT_IN]->data[here]), G[k]->data, G[0]->header->size, gmt_grdfloat);
794 if (GMT_Destroy_Group (API, &G, n_layers))
795 Return (API->error);
796 if (C[GMT_IN]->z == NULL && GMT_Put_Levels (API, C[GMT_IN], level, n_layers))
797 Return (API->error);
798 C[GMT_IN]->mode = GMT_CUBE_IS_STACK; /* Flag that the source was a stack of grids and not a cube */
799 if (gmt_M_is_geographic (GMT, GMT_IN))
800 gmt_set_geographic (GMT, GMT_OUT);
801 else
802 gmt_set_cartesian (GMT, GMT_OUT);
803 if (z_is_abstime) gmt_set_column_type (GMT, GMT_OUT, GMT_Z, GMT_IS_ABSTIME);
804 }
805 else if (C[GMT_IN] == NULL) { /* Read the cube */
806 gmt_M_free (GMT, level); /* Free this one now since it will be re-read by GMT_Read_Data */
807 if ((C[GMT_IN] = GMT_Read_Data (API, GMT_IS_CUBE, GMT_IS_FILE, GMT_IS_VOLUME, GMT_CONTAINER_AND_DATA, wesn, Ctrl->In.file[0], NULL)) == NULL)
808 Return (GMT_DATA_READ_ERROR);
809 }
810
811 if (convert_to_cube) { /* Just want to build cube from input stack */
812 if (Ctrl->D.active && gmt_decode_cube_h_info (GMT, Ctrl->D.information, C[GMT_IN])) {
813 GMT_Destroy_Data (API, &C[GMT_IN]);
814 Return (GMT_PARSE_ERROR);
815 }
816 GMT_Report (API, GMT_MSG_INFORMATION, "Convert %" PRIu64 " grid layers to a single data cube %s.\n", n_layers, Ctrl->G.file);
817 if (GMT_Set_Comment (API, GMT_IS_CUBE, GMT_COMMENT_IS_OPTION | GMT_COMMENT_IS_COMMAND, options, C[GMT_IN]))
818 Return (EXIT_FAILURE);
819 if (GMT_Write_Data (API, GMT_IS_CUBE, GMT_IS_FILE, GMT_IS_VOLUME, GMT_CONTAINER_AND_DATA, NULL, Ctrl->G.file, C[GMT_IN]))
820 Return (EXIT_FAILURE);
821 Return (GMT_NOERROR);
822 }
823
824 GMT_Report (API, GMT_MSG_INFORMATION, "Interpolate %" PRIu64 " new layers (%g to %g in steps of %g).\n", Ctrl->T.T.n, Ctrl->T.T.array[0], Ctrl->T.T.array[Ctrl->T.T.n-1]);
825
826 gmt_set_column_type (GMT, GMT_OUT, GMT_Z, GMT_IS_FLOAT); /* The 3rd dimension is not time in the grids, but we may have read time via -Z with -f2T */
827
828 /* Create cube with layers for each output level */
829
830 gmt_M_memcpy (wesn, C[GMT_IN]->header->wesn, 4, double); /* This is the output common x/y region now */
831 inc[GMT_X] = C[GMT_IN]->header->inc[GMT_X]; inc[GMT_Y] = C[GMT_IN]->header->inc[GMT_Y]; /* And common x/y increments */
832 if (Ctrl->T.T.var_inc || Ctrl->T.T.n == 1) { /* Not equidistant output levels selected via -T so must pass the number of output levels instead of increment */
833 dims[GMT_Z] = Ctrl->T.T.n; /* Number of output levels */
834 this_dim = dims; /* Pointer to the dims instead of NULL */
835 inc[GMT_Z] = 0.0;
836 }
837 else { /* Normal equidistant output levels lets us pass z-inc */
838 inc[GMT_Z] = Ctrl->T.T.inc;
839 wesn[ZLO] = Ctrl->T.T.min;
840 wesn[ZHI] = Ctrl->T.T.max;
841 }
842
843 if ((C[GMT_OUT] = GMT_Create_Data (API, GMT_IS_CUBE, GMT_IS_VOLUME, GMT_CONTAINER_AND_DATA, this_dim, wesn, inc, C[GMT_IN]->header->registration, GMT_NOTSET, NULL)) == NULL)
844 Return (GMT_MEMORY_ERROR);
845
846 if (Ctrl->D.active && gmt_decode_cube_h_info (GMT, Ctrl->D.information, C[GMT_OUT])) {
847 GMT_Destroy_Data (API, &C[GMT_OUT]);
848 Return (GMT_PARSE_ERROR);
849 }
850
851 /* If not equidistant we must add in the level array manually */
852 if (C[GMT_OUT]->z == NULL && GMT_Put_Levels (API, C[GMT_OUT], Ctrl->T.T.array, Ctrl->T.T.n))
853 Return (API->error);
854
855 /* Allocate input and output arrays for the 1-D spline */
856 if ((i_value = gmt_M_memory (GMT, NULL, C[GMT_IN]->header->n_bands, double)) == NULL) Return (GMT_MEMORY_ERROR);
857 if ((o_value = gmt_M_memory (GMT, NULL, Ctrl->T.T.n, double)) == NULL) Return (GMT_MEMORY_ERROR);
858
859 /* Loop over all coregistered x/y nodes then drill through the cube to resample at new layer values */
860 for (row = 0; row < C[GMT_IN]->header->n_rows; row++) {
861 node = gmt_M_ijp (C[GMT_IN]->header, row, 0); /* Relative node numbers in the input and output layers */
862 for (col = 0; col < C[GMT_IN]->header->n_columns; col++, node++) {
863 for (k = 0; k < C[GMT_IN]->header->n_bands; k++) /* For all available input levels, extract data[x,y,z(k)] */
864 i_value[k] = C[GMT_IN]->data[node+k*C[GMT_IN]->header->size];
865 gmt_intpol (GMT, C[GMT_IN]->z, i_value, NULL, C[GMT_IN]->header->n_bands, Ctrl->T.T.n, Ctrl->T.T.array, o_value, 0.0, int_mode); /* Resample at requested output levels */
866 for (k = 0; k < Ctrl->T.T.n; k++) /* For all output levels, place the interpolated values at this (x,y) across all levels */
867 C[GMT_OUT]->data[node+k*C[GMT_OUT]->header->size] = (float)o_value[k];
868 }
869 }
870 GMT_Destroy_Data (API, &C[GMT_IN]); /* Done with the input cube */
871
872 if (GMT_Set_Comment (API, GMT_IS_CUBE, GMT_COMMENT_IS_OPTION | GMT_COMMENT_IS_COMMAND, options, C[GMT_OUT]))
873 Return (EXIT_FAILURE);
874 if (GMT_Write_Data (API, GMT_IS_CUBE, GMT_IS_FILE, GMT_IS_VOLUME, GMT_CONTAINER_AND_DATA, NULL, Ctrl->G.file, C[GMT_OUT]))
875 Return (EXIT_FAILURE);
876
877 GMT_Destroy_Data (API, &C[GMT_OUT]); /* Done with the output cube */
878
879 /* Done with everything; free up remaining memory */
880
881 gmt_M_free (GMT, i_value);
882 gmt_M_free (GMT, o_value);
883
884 Return (GMT_NOERROR);
885 }
886