1 /*--------------------------------------------------------------------
2 *
3 * Copyright (c) 1999-2021 by the GMT Team (https://www.generic-mapping-tools.org/team.html)
4 *
5 * This program is free software; you can redistribute it and/or modify
6 * it under the terms of the GNU Lesser General Public License as published by
7 * the Free Software Foundation; version 3 or any later version.
8 *
9 * This program is distributed in the hope that it will be useful,
10 * but WITHOUT ANY WARRANTY; without even the implied warranty of
11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 * GNU Lesser General Public License for more details.
13 *
14 * Contact info: www.generic-mapping-tools.org
15 *--------------------------------------------------------------------*/
16 /*
17 * GRDSPOTTER will (1) read a grid with topo/grav/whatever of seamounts,
18 * (2) read an ASCII file with stage or finite rotations, and (3)
19 * convolve the flowline of each node with its prism volume, and
20 * (4) build a cumulative volcano amplitude (CVA) grid. Optionally, it
21 * can also calculate the Data Importance (DI) associated with each node's
22 * flowline and the predicted age (PA) for each node. Finally, errors in
23 * the location of the CVA maximum can be assessed by running a bootstrap
24 * estimation. The grids are all written out in GMT format and can be
25 * processed and plotted with GMT.
26 * GRDSPOTTER is part of the SPOTTER supplemental GMT package and should
27 * be installed accordingly, see README.spotter.
28 *
29 * Author: Paul Wessel, SOEST, Univ. of Hawaii, Honolulu, HI, USA
30 * Date: 20-OCT-2007
31 * Version: 2.0
32 *
33 *-------------------------------------------------------------------------
34 * The Euler file must have following format:
35 *
36 * 1. Any number of comment lines starting with # in first column
37 * 2. Any number of blank lines (just carriage return, no spaces)
38 * 2. Any number of stage pole records which each have the format:
39 * lon(deg) lat(deg) tstart(Ma) tstop(Ma) ccw-angle(deg)
40 * 3. stage records must go from oldest to youngest rotation
41 * 4. Note tstart is larger (older) that tstop for each record
42 * 5. No gaps allowed: tstart must equal the previous records tstop
43 *
44 * Example: Duncan & Clague [1985] Pacific-Hotspot rotations:
45 *
46 * # Time in Ma, angles in degrees
47 * # lon lat tstart tend ccw-angle
48 * 165 85 150 100 24.0
49 * 284 36 100 74 15.0
50 * 265 22 74 65 7.5
51 * 253 17 65 42 14.0
52 * 285 68 42 0 34.0
53 *
54 * If total reconstruction poles are preferred,
55 * give a file like DC85 poles in finite rotation format:
56 *
57 * #longitude latitude time(My) angle(deg)
58 * 285.00000 68.00000 42.0000 34.0000
59 * 275.66205 53.05082 65.0000 43.5361
60 * 276.02501 48.34232 74.0000 50.0405
61 * 279.86436 46.30610 100.0000 64.7066
62 * 265.37800 55.69932 150.0000 82.9957
63 *
64 * Note that finite rotation files have one less column since each rotation
65 * implicitly goes to time 0.
66 *
67 * Seamount data grid is a standard GMT grid file, presumably
68 * representing topography or gravity of residual seamounts
69 * (background residual should be removed first).
70 *
71 *------------------------------------------------------------------------------
72 * REFERENCES:
73 *
74 * -> The hotspotting technique:
75 *
76 * Aslanian, D., L. Geli, and J.-L. Olivet, 1998, Hotspotting called into
77 * question, Nature, 396, 127.
78 * Wessel, P., and L. W. Kroenke, 1997, A geometric technique for
79 * relocating hotspots and refining absolute plate motions, Nature,
80 * 387, 365-369.
81 * Wessel, P., and L. W. Kroenke, 1998a, Factors influencing the locations
82 * of hot spots determined by the hot-spotting technique, Geophys. Res.
83 * Lett., 25, 55-58.
84 * Wessel, P., and L. W. Kroenke, 1998b, The geometric relationship between
85 * hot spots and seamounts: implications for Pacific hot spots, Earth
86 * Planet. Sci. Lett., 158, 1-18.
87 * Wessel, P., and L. W. Kroenke, 1998c, Hotspotting called into question
88 * - Reply, Nature, 396, 127-128.
89 *
90 * -> Suitable Seamount data set:
91 *
92 * Sandwell, D, and W. H. F. Smith, 1995, JGR, ...
93 * Smith, W. H. F., and D. Sandwell, 1997, Science, ...
94 *
95 * -> Plate motion models (stage poles):
96 *
97 * Duncan, R.A., and D. Clague, 1985, Pacific plate motion recorded by linear
98 * volcanic chains, in: A.E.M. Nairn, F. G. Stehli, S. Uyeda (eds.), The
99 * Ocean Basins and Margins, Vol. 7A, Plenum, New York, pp. 89-121.
100 * Wessel and Kroenke, 1997, (see above)
101 *
102 */
103
104 #include "gmt_dev.h"
105 #include "spotter.h"
106
107 #define THIS_MODULE_CLASSIC_NAME "grdspotter"
108 #define THIS_MODULE_MODERN_NAME "grdspotter"
109 #define THIS_MODULE_LIB "spotter"
110 #define THIS_MODULE_PURPOSE "Create CVA grid from a gravity or topography grid"
111 #define THIS_MODULE_KEYS "<G{,AG(,DG),LG),GG}"
112 #define THIS_MODULE_NEEDS "g"
113 #define THIS_MODULE_OPTIONS "-:>RVhr" GMT_OPT("F")
114
115 #define B_TO_MB (1.0 / 1048576.0)
116
117 #define TRUNC 0 /* Indices for -T */
118 #define UPPER 1
119
120 struct GRDSPOTTER_CTRL { /* All control options for this program (except common args) */
121 /* active is true if the option has been activated */
122 struct GRDSPOTTER_In {
123 bool active;
124 char *file;
125 } In;
126 struct GRDSPOTTER_A { /* -A<file> */
127 bool active;
128 char *file;
129 } A;
130 struct GRDSPOTTER_D { /* -Di<file> */
131 bool active;
132 char *file;
133 } D;
134 struct GRDSPOTTER_E { /* -E[+rotfile */
135 bool active;
136 bool mode;
137 char *file;
138 } E;
139 struct GRDSPOTTER_G { /* -Ggrdfile */
140 bool active; /* Pixel registration */
141 char *file;
142 } G;
143 struct GRDSPOTTER_I { /* -I (for checking only) */
144 bool active;
145 } I;
146 struct GRDSPOTTER_L { /* -Lfile */
147 bool active;
148 char *file;
149 } L;
150 struct GRDSPOTTER_M { /* -M */
151 bool active;
152 } M;
153 struct GRDSPOTTER_N { /* -N */
154 bool active;
155 double t_upper;
156 } N;
157 struct GRDSPOTTER_PA { /* -Dp<file> */
158 bool active;
159 char *file;
160 } PA;
161 struct GRDSPOTTER_Q { /* -Q */
162 bool active;
163 unsigned int mode;
164 unsigned int id;
165 char *file;
166 } Q;
167 struct GRDSPOTTER_S2 { /* -S2 */
168 bool active;
169 double dist;
170 } S2;
171 struct GRDSPOTTER_S { /* -S */
172 bool active;
173 } S;
174 struct GRDSPOTTER_T { /* -T */
175 bool active[2];
176 double t_fix; /* Set fixed age*/
177 } T;
178 struct GRDSPOTTER_W { /* -W */
179 bool active;
180 unsigned int n_try;
181 } W;
182 struct GRDSPOTTER_Z { /* -Z */
183 bool active;
184 bool mode;
185 double min, max, inc;
186 } Z;
187 };
188
New_Ctrl(struct GMT_CTRL * GMT)189 static void *New_Ctrl (struct GMT_CTRL *GMT) { /* Allocate and initialize a new control structure */
190 struct GRDSPOTTER_CTRL *C;
191
192 C = gmt_M_memory (GMT, NULL, 1, struct GRDSPOTTER_CTRL);
193
194 /* Initialize values whose defaults are not 0/false/NULL */
195
196 C->N.t_upper = 180.0; /* Upper age assigned to nodes on undated seafloor */
197 C->Z.max = DBL_MAX; /* Only deal with z-values inside this range [0/Inf] */
198 return (C);
199 }
200
Free_Ctrl(struct GMT_CTRL * GMT,struct GRDSPOTTER_CTRL * C)201 static void Free_Ctrl (struct GMT_CTRL *GMT, struct GRDSPOTTER_CTRL *C) { /* Deallocate control structure */
202 if (!C) return;
203 gmt_M_str_free (C->In.file);
204 gmt_M_str_free (C->A.file);
205 gmt_M_str_free (C->D.file);
206 gmt_M_str_free (C->E.file);
207 gmt_M_str_free (C->G.file);
208 gmt_M_str_free (C->L.file);
209 gmt_M_str_free (C->PA.file);
210 gmt_M_free (GMT, C);
211 }
212
usage(struct GMTAPI_CTRL * API,int level)213 static int usage (struct GMTAPI_CTRL *API, int level) {
214 const char *name = gmt_show_name_and_purpose (API, THIS_MODULE_LIB, THIS_MODULE_CLASSIC_NAME, THIS_MODULE_PURPOSE);
215 if (level == GMT_MODULE_PURPOSE) return (GMT_NOERROR);
216 GMT_Usage (API, 0, "usage: %s %s %s -G%s %s %s [-A<agegrid>] [-D[i|p]<grdfile>] "
217 "[-L<IDgrid>] [-M] [-N<upper_age>] [-Q<IDinfo>] [-S] [-Tt|u<age>] [%s] [-W<n_try>] [-Z<z_min>[/<z_max>[/<z_inc>]]] "
218 "[%s] [%s] [%s]\n", name, GMT_INGRID, SPOTTER_E_OPT, GMT_OUTGRID, GMT_I_OPT, GMT_Rgeo_OPT, GMT_V_OPT, GMT_ho_OPT, GMT_r_OPT, GMT_PAR_OPT);
219
220 if (level == GMT_SYNOPSIS) return (GMT_MODULE_SYNOPSIS);
221
222 GMT_Message (API, GMT_TIME_NONE, " REQUIRED ARGUMENTS:\n");
223 gmt_ingrid_syntax (API, 0, "Name of input grid with topo or gravity");
224 spotter_rot_usage (API);
225 gmt_outgrid_syntax (API, 'G', "Specify file name for output CVA convolution grid");
226 GMT_Option (API, "I,Rg");
227 GMT_Message (API, GMT_TIME_NONE, "\n OPTIONAL ARGUMENTS:\n");
228 GMT_Usage (API, 1, "\n-A<agegrid>");
229 GMT_Usage (API, -2, "Co-registered grid with upper ages to use [Default is flowlines for all ages].");
230 GMT_Usage (API, 1, "\n-D[i|p]<grdfile>");
231 GMT_Usage (API, -2, "Set optional output grids:");
232 GMT_Usage (API, 3, "i: Use flowlines to estimate and write data importance DI grid.");
233 GMT_Usage (API, 3, "p: Use flowlines to estimate and write predicted ages at node locations.");
234 GMT_Usage (API, 1, "\n-L<IDgrid>");
235 GMT_Usage (API, -2, "Co-registered grid with chain ID for each node [Default ignores IDs].");
236 GMT_Usage (API, 1, "\n-M Do flowline calculations as needed rather than storing in memory. "
237 "You may have to use this option if -R is too large. Cannot be used with -W or -Z-slicing.");
238 GMT_Usage (API, 1, "\n-N<upper_age>");
239 GMT_Usage (API, -2, "Set upper age in m.y. for nodes whose plate age is NaN [180].");
240 GMT_Usage (API, 1, "\n-Q<IDinfo>");
241 GMT_Usage (API, -2, "Give either single ID to use or file with list of IDs [Default uses all IDs]. "
242 "Each line would be TAG ID [w e s n] with optional zoom box.");
243 GMT_Usage (API, 1, "\n-S Normalize CVA grid to percentages of the CVA maximum.");
244 GMT_Usage (API, 1, "\n-Tt|u<age>");
245 GMT_Usage (API, -2, "Set upper ages. Repeatable, choose from:");
246 GMT_Usage (API, 3, "t: Truncate all ages to max <age> in stage pole model [Default extrapolates].");
247 GMT_Usage (API, 3, "u: After a node passes the -Z test, use this fixed <age >instead in CVA calculations.");
248 GMT_Option (API, "V");
249 GMT_Usage (API, 1, "\n-W<n_try>");
250 GMT_Usage (API, -2, "Get <n_try> bootstrap estimates of maximum CVA location [Default is no bootstrapping].");
251 GMT_Usage (API, 1, "\n-Z<z_min>[/<z_max>[/<z_inc>]]");
252 GMT_Usage (API, -2, "Ignore nodes with z-value lower than z_min [0] and optionally larger than z_max [Inf]. "
253 "Give z_min/z_max/z_inc to make CVA grids for each z-slice {Default makes 1 CVA grid].");
254 GMT_Option (API, "h,r,.");
255
256 return (GMT_MODULE_USAGE);
257 }
258
parse(struct GMT_CTRL * GMT,struct GRDSPOTTER_CTRL * Ctrl,struct GMT_OPTION * options)259 static int parse (struct GMT_CTRL *GMT, struct GRDSPOTTER_CTRL *Ctrl, struct GMT_OPTION *options) {
260 /* This parses the options provided to grdspotter and sets parameters in CTRL.
261 * Any GMT common options will override values set previously by other commands.
262 * It also replaces any file names specified as input or output with the data ID
263 * returned when registering these sources/destinations with the API.
264 */
265
266 unsigned int n_errors = 0, k, n_files = 0;
267 int m, sval;
268 struct GMT_OPTION *opt = NULL;
269 struct GMTAPI_CTRL *API = GMT->parent;
270
271 for (opt = options; opt; opt = opt->next) {
272 switch (opt->option) {
273
274 case '<': /* Input files */
275 if (n_files++ > 0) break;
276 Ctrl->In.active = true;
277 if (opt->arg[0]) Ctrl->In.file = strdup (opt->arg);
278 if (GMT_Get_FilePath (API, GMT_IS_GRID, GMT_IN, GMT_FILE_REMOTE, &(Ctrl->In.file))) n_errors++;
279 break;
280
281 /* Supplemental parameters */
282
283 case 'A':
284 n_errors += gmt_M_repeated_module_option (API, Ctrl->A.active);
285 Ctrl->A.active = true;
286 if (opt->arg[0]) Ctrl->A.file = strdup (opt->arg);
287 if (GMT_Get_FilePath (API, GMT_IS_GRID, GMT_IN, GMT_FILE_REMOTE, &(Ctrl->A.file))) n_errors++;
288 break;
289 case 'C': /* Now done automatically in spotter_init */
290 if (gmt_M_compat_check (GMT, 4))
291 GMT_Report (API, GMT_MSG_COMPAT, "-C is no longer needed as total reconstruction vs stage rotation is detected automatically.\n");
292 else
293 n_errors += gmt_default_error (GMT, opt->option);
294 break;
295 case 'E':
296 n_errors += gmt_M_repeated_module_option (API, Ctrl->E.active);
297 Ctrl->E.active = true; k = 0;
298 if (opt->arg[0] == '+') { Ctrl->E.mode = true; k = 1;}
299 if (opt->arg[k]) Ctrl->E.file = strdup (&opt->arg[k]);
300 if (GMT_Get_FilePath (API, GMT_IS_DATASET, GMT_IN, GMT_FILE_REMOTE, &(Ctrl->E.file))) n_errors++;
301 break;
302 case 'G':
303 n_errors += gmt_M_repeated_module_option (API, Ctrl->G.active);
304 Ctrl->G.active = true;
305 if (opt->arg[0]) Ctrl->G.file = strdup (opt->arg);
306 if (GMT_Get_FilePath (API, GMT_IS_GRID, GMT_OUT, GMT_FILE_LOCAL, &(Ctrl->G.file))) n_errors++;
307 break;
308 case 'D':
309 n_errors += gmt_M_repeated_module_option (API, Ctrl->D.active);
310 switch (opt->arg[0]) {
311 case 'i':
312 Ctrl->D.active = true;
313 if (opt->arg[1]) Ctrl->D.file = strdup (&opt->arg[1]);
314 if (GMT_Get_FilePath (API, GMT_IS_GRID, GMT_OUT, GMT_FILE_LOCAL, &(Ctrl->D.file))) n_errors++;
315 break;
316 case 'p':
317 Ctrl->PA.active = true;
318 if (opt->arg[1]) Ctrl->PA.file = strdup (&opt->arg[1]);
319 if (GMT_Get_FilePath (API, GMT_IS_GRID, GMT_OUT, GMT_FILE_LOCAL, &(Ctrl->PA.file))) n_errors++;
320 break;
321 default:
322 n_errors++;
323 break;
324 }
325 break;
326 case 'I':
327 n_errors += gmt_M_repeated_module_option (API, Ctrl->I.active);
328 Ctrl->I.active = true;
329 n_errors += gmt_parse_inc_option (GMT, 'I', opt->arg);
330 break;
331 case 'L':
332 n_errors += gmt_M_repeated_module_option (API, Ctrl->L.active);
333 Ctrl->L.active = true;
334 if (opt->arg[0]) Ctrl->L.file = strdup (opt->arg);
335 if (GMT_Get_FilePath (API, GMT_IS_GRID, GMT_IN, GMT_FILE_REMOTE, &(Ctrl->L.file))) n_errors++;
336 break;
337 case 's': /* Sampling interval */
338 Ctrl->S2.active = true;
339 Ctrl->S2.dist = atof (opt->arg);
340 break;
341 case 'M':
342 n_errors += gmt_M_repeated_module_option (API, Ctrl->M.active);
343 Ctrl->M.active = true;
344 break;
345 case 'N':
346 n_errors += gmt_M_repeated_module_option (API, Ctrl->N.active);
347 Ctrl->N.active = true;
348 Ctrl->N.t_upper = atof (opt->arg);
349 break;
350 case 'Q':
351 n_errors += gmt_M_repeated_module_option (API, Ctrl->Q.active);
352 Ctrl->Q.active = true;
353 if (!access (opt->arg, R_OK)) { /* The file exists */
354 Ctrl->Q.mode = 2;
355 Ctrl->Q.file = strdup (opt->arg);
356 }
357 else if ((sval = atoi (opt->arg)) > 0) { /* Got OK id value */
358 Ctrl->Q.mode = 1;
359 Ctrl->Q.id = sval;
360 }
361 else {
362 GMT_Report (API, GMT_MSG_ERROR, "Option -Q: Must give valid file or ID value\n");
363 n_errors++;
364 }
365 break;
366 case 'S':
367 n_errors += gmt_M_repeated_module_option (API, Ctrl->S.active);
368 Ctrl->S.active = true;
369 break;
370 case 'T':
371 if (opt->arg[0] == 't') {
372 n_errors += gmt_M_repeated_module_option (API, Ctrl->T.active[TRUNC]);
373 Ctrl->T.active[TRUNC] = true;
374 }
375 else if (opt->arg[0] == 'u') {
376 Ctrl->T.t_fix = atof (&opt->arg[1]);
377 n_errors += gmt_M_repeated_module_option (API, Ctrl->T.active[UPPER]);
378 Ctrl->T.active[UPPER] = true;
379 }
380 else {
381 GMT_Report (API, GMT_MSG_ERROR, "Option -T: Either use -Tt or -Tu<age>\n");
382 n_errors++;
383 }
384 break;
385 case 'W':
386 n_errors += gmt_M_repeated_module_option (API, Ctrl->W.active);
387 Ctrl->W.n_try = atoi (opt->arg);
388 Ctrl->W.active = true;
389 break;
390 case 'Z':
391 n_errors += gmt_M_repeated_module_option (API, Ctrl->Z.active);
392 Ctrl->Z.active = true;
393 m = sscanf (opt->arg, "%lf/%lf/%lf", &Ctrl->Z.min, &Ctrl->Z.max, &Ctrl->Z.inc);
394 if (m == 1) Ctrl->Z.max = 1.0e300; /* Max not specified */
395 if (m == 3) Ctrl->Z.mode = true; /* Want several slices */
396 break;
397 default: /* Report bad options */
398 n_errors += gmt_default_error (GMT, opt->option);
399 break;
400 }
401 }
402
403 n_errors += gmt_M_check_condition (GMT, !GMT->common.R.active[RSET], "Must specify -R option\n");
404 n_errors += gmt_M_check_condition (GMT, GMT->common.R.inc[GMT_X] <= 0.0 || GMT->common.R.inc[GMT_Y] <= 0.0, "Option -I: Must specify positive increment(s)\n");
405 n_errors += gmt_M_check_condition (GMT, !(Ctrl->G.active || Ctrl->G.file), "Option -G: Must specify output file\n");
406 n_errors += gmt_M_check_condition (GMT, !(Ctrl->In.active || Ctrl->In.file), "Option -Z: Must give name of topo gridfile\n");
407 n_errors += gmt_M_check_condition (GMT, Ctrl->L.file && !Ctrl->Q.mode, "Must specify both -L and -Q if one is present\n");
408 n_errors += gmt_M_check_condition (GMT, Ctrl->M.active && (Ctrl->W.active || Ctrl->Z.mode), "Cannot use -M with -W or -Z (slicing)\n");
409
410 return (n_errors ? GMT_PARSE_ERROR : GMT_NOERROR);
411 }
412
grdspotter_get_flowline(struct GMT_CTRL * GMT,double xx,double yy,double tt,struct EULER * p,unsigned int n_stages,double d_km,unsigned int step,unsigned int flag,double wesn[],double ** flow)413 GMT_LOCAL int64_t grdspotter_get_flowline (struct GMT_CTRL *GMT, double xx, double yy, double tt, struct EULER *p, unsigned int n_stages, double d_km, unsigned int step, unsigned int flag, double wesn[], double **flow) {
414 int64_t n_track, m, kx, ky, np, first, last;
415 double *c = NULL, *f = NULL;
416
417 /* Get the flowline from this point back to time tt, restricted to the given wesn box */
418 if (spotter_forthtrack (GMT, &xx, &yy, &tt, 1, p, n_stages, d_km, 0.0, flag, wesn, &c) <= 0) {
419 GMT_Report (GMT->parent, GMT_MSG_ERROR, "Nothing returned from spotter_forthtrack - skipping\n");
420 return 0LL;
421 }
422
423 n_track = lrint (c[0]); /* Number of point pairs making up this flowline */
424
425 /* Find the first point on the flowline inside the desired CVA region */
426
427 for (m = 0, ky = 2, first = GMT_NOTSET; m < n_track && first == GMT_NOTSET; m++, ky += step) { /* For each point along flowline */
428 if (c[ky] < wesn[YLO] || c[ky] > wesn[YHI]) continue; /* Latitude outside region */
429 kx = ky - 1; /* Index for the x-coordinate */
430 while (c[kx] > wesn[XHI]) c[kx] -= TWO_PI; /* Elaborate W/E test because of 360 periodicity */
431 while (c[kx] < wesn[XLO]) c[kx] += TWO_PI;
432 if (c[kx] > wesn[XHI]) continue; /* Longitude outside region */
433 first = kx; /* We are inside, this terminates the for loop */
434 }
435
436 if (first == -1) { /* Was never inside the grid, skip the entire flowline and move on */
437 gmt_M_free (GMT, c); /* Free the flowline vector */
438 return 0LL;
439 }
440
441 /* Here we know searching from the end will land inside the grid eventually so last can never exit as -1 */
442
443 for (m = n_track - 1, ky = step * m + 2, last = GMT_NOTSET; m >= 0 && last == GMT_NOTSET; m--, ky -= step) { /* For each point along flowline */
444 if (c[ky] < wesn[YLO] || c[ky] > wesn[YHI]) continue; /* Latitude outside region */
445 kx = ky - 1; /* Index for the x-coordinate */
446 while (c[kx] > wesn[XHI]) c[kx] -= TWO_PI; /* Elaborate W/E test because of 360 periodicity */
447 while (c[kx] < wesn[XLO]) c[kx] += TWO_PI;
448 if (c[kx] > wesn[XHI]) continue; /* Longitude outside region */
449 last = kx; /* We are inside, this terminates the for loop */
450 }
451
452 np = (last - first) / step + 1; /* Number of (x,y[,t]) points on this flowline inside the region */
453 if (np < n_track) { /* Just copy out the subset of points we want */
454 size_t n_alloc;
455 n_alloc = np * step; /* Number of (x,y[,t]) to copy */
456 f = gmt_M_memory (GMT, NULL, n_alloc+1, double);
457 f[0] = (double)np; /* Number of points found */
458 gmt_M_memcpy (&f[1], &c[first], n_alloc, double);
459 gmt_M_free (GMT, c); /* Free the old flowline vector */
460 *flow = f; /* Return pointer to trimmed flowline */
461 }
462 else
463 *flow = c; /* Return the entire flowline as is */
464 return (np);
465 }
466
grdspotter_set_age(struct GMT_CTRL * GMT,double * t_smt,struct GMT_GRID * A,uint64_t node,double upper_age,bool truncate)467 GMT_LOCAL bool grdspotter_set_age (struct GMT_CTRL *GMT, double *t_smt, struct GMT_GRID *A, uint64_t node, double upper_age, bool truncate) {
468 /* Returns the age of this node based on either a given seafloor age grid
469 * or the upper age, truncated if necessary */
470
471 if (!A || gmt_M_is_fnan (A->data[node])) /* Age is NaN, assign upper value */
472 *t_smt = upper_age;
473 else { /* Assign given value */
474 *t_smt = A->data[node];
475 if (*t_smt > upper_age) { /* Exceeding our limit */
476 if (truncate) /* Allowed to truncate to max age */
477 *t_smt = upper_age;
478 else { /* Consider this an error or just skip */
479 GMT_Report (GMT->parent, GMT_MSG_WARNING, "Node %" PRIu64 " has age (%g) > oldest stage (%g) (skipped)\n", node, *t_smt, upper_age);
480 return (false);
481 }
482 }
483 }
484 return (true); /* We are returning a useful age */
485 }
486
grdspotter_normalize_grid(struct GMT_CTRL * GMT,struct GMT_GRID * G,gmt_grdfloat * data)487 GMT_LOCAL void grdspotter_normalize_grid (struct GMT_CTRL *GMT, struct GMT_GRID *G, gmt_grdfloat *data) {
488 /* Determines the grid's min/max z-values and normalizes the grid
489 * so that zmax is 100% */
490 unsigned int row, col;
491 uint64_t node;
492 double CVA_scale; /* Used to normalize CVAs to percent */
493
494 G->header->z_min = +DBL_MAX;
495 G->header->z_max = -DBL_MAX;
496 gmt_M_grd_loop (GMT, G, row, col, node) { /* Loop over all output nodes */
497 if (gmt_M_is_fnan (data[node])) continue;
498 if (data[node] < G->header->z_min) G->header->z_min = data[node];
499 if (data[node] > G->header->z_max) G->header->z_max = data[node];
500 }
501 GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "CVA min/max: %g %g -> ", G->header->z_min, G->header->z_max);
502 CVA_scale = 100.0 / G->header->z_max;
503 for (node = 0; node < G->header->size; node++) data[node] *= (gmt_grdfloat)CVA_scale;
504 G->header->z_min *= CVA_scale;
505 G->header->z_max *= CVA_scale;
506 GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "%g %g\n", G->header->z_min, G->header->z_max);
507 }
508
509 #define bailout(code) {gmt_M_free_options (mode); return (code);}
510 #define Return(code) {Free_Ctrl (GMT, Ctrl); gmt_end_module (GMT, GMT_cpy); bailout (code);}
511
GMT_grdspotter(void * V_API,int mode,void * args)512 EXTERN_MSC int GMT_grdspotter (void *V_API, int mode, void *args) {
513 unsigned int n_stages; /* Number of stage rotations (poles) */
514 unsigned int try; /* Number of current bootstrap estimate */
515 unsigned int row, row2, col, col2, k_step;
516 unsigned int forth_flag; /* Holds the do_time + 10 flag passed to forthtrack */
517 bool keep_flowlines = false; /* true if Ctrl->D.active, Ctrl->PA.active, or bootstrap is true */
518 int error = GMT_NOERROR; /* nonzero when arguments are wrong */
519 int i, j; /* Signed row,col variables */
520 int *ID = NULL; /* Optional array with IDs for each node */
521
522 uint64_t ij, k, node, m, np, max_ij = 0, n_flow, n_unique_nodes = 0;
523 uint64_t n_nodes; /* Number of nodes processed */
524
525 size_t n_alloc = 0, inc_alloc = BIG_CHUNK, mem = 0;
526
527 char *processed_node = NULL; /* Pointer to array with true/false values for each grid node */
528
529 unsigned short pa = 0; /* Placeholder for PA along track */
530
531 gmt_grdfloat zz;
532
533 double sampling_int_in_km; /* Sampling interval along flowline (in km) */
534 double *x_smt = NULL; /* node longitude (input degrees, stored as radians) */
535 double *y_smt = NULL; /* node latitude (input degrees, stored as radians) */
536 double t_smt = 0.0; /* node upper age (up to age of seafloor) */
537 double *c = NULL; /* Array with one flowline */
538 double CVA_max, wesn[4], cva_contribution, yg;
539 double out[3], x_scale, y_scale, area;
540 double *lat_area = NULL; /* Area of each dx by dy note in km as function of latitude */
541 double this_wesn[4];
542 double this_pa, pa_val = 0.0, n_more_than_once = 0.0;
543 #ifdef DEBUG2
544 double *x_cva = NULL; /* x-coordinates on the CVA grid */
545 double *y_cva = NULL; /* y-coordinates on the CVA grid */
546 #endif
547 struct GMT_GRID *G = NULL; /* Grid structure for output CVA grid */
548 struct GMT_GRID *G_rad = NULL; /* Same but has radians in header (no grid) */
549 struct GMT_GRID *Z = NULL; /* Grid structure for input topo/grav grid */
550 struct GMT_GRID *A = NULL; /* Grid structure for input age grid */
551 struct GMT_GRID *L = NULL; /* Grid structure for input ID grid */
552 struct GMT_GRID *PA = NULL; /* Grid structure for output PA grid */
553 struct GMT_GRID *DI = NULL; /* Grid structure for output DI grid */
554
555 struct EULER *p = NULL; /* Array of structures with Euler stage rotations */
556
557 struct FLOWLINE *flowline = NULL; /* Array with flowline structures */
558
559 struct GRDSPOTTER_ID { /* Information regarding one chain ID */
560 double wesn[4]; /* Do not calculate flowlines outside this box */
561 bool ok; /* true if we want to calculate this CVA */
562 bool check_region; /* true if w, e, s, n is more restrictive than command line -R */
563 } *ID_info = NULL;
564 struct GMT_OPTION *ptr = NULL;
565 struct GRDSPOTTER_CTRL *Ctrl = NULL;
566 struct GMT_CTRL *GMT = NULL, *GMT_cpy = NULL;
567 struct GMT_OPTION *options = NULL;
568 struct GMTAPI_CTRL *API = gmt_get_api_ptr (V_API); /* Cast from void to GMTAPI_CTRL pointer */
569
570 /*----------------------- Standard module initialization and parsing ----------------------*/
571
572 if (API == NULL) return (GMT_NOT_A_SESSION);
573 if (mode == GMT_MODULE_PURPOSE) return (usage (API, GMT_MODULE_PURPOSE)); /* Return the purpose of program */
574 options = GMT_Create_Options (API, mode, args); if (API->error) return (API->error); /* Set or get option list */
575
576 if ((error = gmt_report_usage (API, options, 0, usage)) != GMT_NOERROR) bailout (error); /* Give usage if requested */
577
578 /* Parse the command-line arguments */
579
580 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 */
581 if ((ptr = GMT_Find_Option (API, 'f', options)) == NULL) gmt_parse_common_options (GMT, "f", 'f', "g"); /* Did not set -f, implicitly set -fg */
582 if (GMT_Parse_Common (API, THIS_MODULE_OPTIONS, options)) Return (API->error);
583 Ctrl = New_Ctrl (GMT); /* Allocate and initialize a new control structure */
584 if ((error = parse (GMT, Ctrl, options)) != 0) Return (error);
585
586 /*---------------------------- This is the grdspotter main code ----------------------------*/
587
588 /* Initialize the CVA grid and structure */
589
590 if ((G = GMT_Create_Data (API, GMT_IS_GRID, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA, NULL, NULL, NULL, \
591 GMT_GRID_DEFAULT_REG, GMT_NOTSET, NULL)) == NULL) Return (API->error);
592
593 /* ------------------- END OF PROCESSING COMMAND LINE ARGUMENTS --------------------------------------*/
594
595 /* Load in the Euler stage poles */
596
597 if ((i = spotter_init (GMT, Ctrl->E.file, &p, 1, false, Ctrl->E.mode, &Ctrl->N.t_upper)) < 0)
598 Return (-i);
599
600 n_stages = (unsigned int)i;
601
602 /* Assign grid-region variables in radians to avoid conversions inside convolution loop */
603
604 if ((G_rad = GMT_Duplicate_Data (API, GMT_IS_GRID, GMT_DUPLICATE_NONE, G)) == NULL) Return (API->error);
605 G_rad->header->inc[GMT_X] = G->header->inc[GMT_X] * D2R;
606 G_rad->header->inc[GMT_Y] = G->header->inc[GMT_Y] * D2R;
607 G_rad->header->wesn[XLO] = G->header->wesn[XLO] * D2R;
608 G_rad->header->wesn[XHI] = G->header->wesn[XHI] * D2R;
609 G_rad->header->wesn[YLO] = G->header->wesn[YLO] * D2R;
610 G_rad->header->wesn[YHI] = G->header->wesn[YHI] * D2R;
611
612 /* Get geocentric wesn region in radians */
613 gmt_M_memcpy (wesn, G_rad->header->wesn, 4, double);
614 wesn[YLO] = D2R * gmt_lat_swap (GMT, R2D * wesn[YLO], GMT_LATSWAP_G2O);
615 wesn[YHI] = D2R * gmt_lat_swap (GMT, R2D * wesn[YHI], GMT_LATSWAP_G2O);
616
617 /* Allocate T/F array */
618
619 processed_node = gmt_M_memory (GMT, NULL, G->header->size, char);
620
621 /* Set flowline sampling interval to 1/2 of the shortest distance between x-nodes */
622
623 /* sampling_int_in_km = 0.5 * G_rad->header->inc[GMT_X] * EQ_RAD * ((fabs (G_rad->header->wesn[YHI]) > fabs (G_rad->header->wesn[YLO])) ? cos (G_rad->header->wesn[YHI]) : cos (G_rad->header->wesn[YLO])); */
624 sampling_int_in_km = G_rad->header->inc[GMT_X] * EQ_RAD * ((fabs (G_rad->header->wesn[YHI]) > fabs (G_rad->header->wesn[YLO])) ?
625 cos (G_rad->header->wesn[YHI]) : cos (G_rad->header->wesn[YLO]));
626 if (Ctrl->S2.dist != 0.0) sampling_int_in_km = Ctrl->S2.dist;
627 GMT_Report (API, GMT_MSG_INFORMATION, "Flowline sampling interval = %.3f km\n", sampling_int_in_km);
628
629 if (Ctrl->T.active[TRUNC]) GMT_Report (API, GMT_MSG_INFORMATION, "Ages truncated to %g\n", Ctrl->N.t_upper);
630
631 /* Start to read input data */
632
633 if (GMT_Init_IO (API, GMT_IS_GRID, GMT_IS_SURFACE, GMT_IN, GMT_ADD_DEFAULT, 0, options) != GMT_NOERROR) { /* Establishes data input */
634 gmt_M_free (GMT, processed_node);
635 Return (API->error);
636 }
637
638 if ((Z = GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA, NULL, Ctrl->In.file, NULL)) == NULL) { /* Get data */
639 gmt_M_free (GMT, processed_node);
640 Return (API->error);
641 }
642 area = 111.195 * Z->header->inc[GMT_Y] * 111.195 * Z->header->inc[GMT_X]; /* In km^2 at Equator */
643 x_smt = gmt_M_memory (GMT, NULL, Z->header->n_columns, double);
644 for (col = 0; col < Z->header->n_columns; col++) x_smt[col] = D2R * gmt_M_grd_col_to_x (GMT, col, Z->header);
645 y_smt = gmt_M_memory (GMT, NULL, Z->header->n_rows, double);
646 for (row = 0; row < Z->header->n_rows; row++) y_smt[row] = D2R * gmt_lat_swap (GMT, gmt_M_grd_row_to_y (GMT, row, Z->header), GMT_LATSWAP_G2O); /* Convert to geocentric */
647 lat_area = gmt_M_memory (GMT, NULL, Z->header->n_rows, double);
648
649 for (row = 0; row < Z->header->n_rows; row++) lat_area[row] = area * cos (y_smt[row]);
650
651 #ifdef DEBUG2
652 x_cva = G->x;
653 y_cva = G->y;
654 #endif
655
656 if (Ctrl->A.file) {
657 if ((A = GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_ONLY, NULL, Ctrl->A.file, NULL)) == NULL) { /* Get header only */
658 Return (API->error);
659 }
660 if (!(A->header->n_columns == Z->header->n_columns && A->header->n_rows == Z->header->n_rows && A->header->wesn[XLO] == Z->header->wesn[XLO] && A->header->wesn[YLO] == Z->header->wesn[YLO])) {
661 GMT_Report (API, GMT_MSG_ERROR, "Topo grid and age grid must coregister\n");
662 Return (GMT_RUNTIME_ERROR);
663 }
664 if (GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_DATA_ONLY, NULL, Ctrl->A.file, A) == NULL) { /* Get age data */
665 Return (API->error);
666 }
667 }
668 if (Ctrl->L.file) {
669 struct GMT_GRID_HIDDEN *GH = gmt_get_G_hidden (G);
670 if ((L = GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_ONLY, NULL, Ctrl->L.file, NULL)) == NULL) { /* Get header only */
671 Return (API->error);
672 }
673 if (!(L->header->n_columns == Z->header->n_columns && L->header->n_rows == Z->header->n_rows && L->header->wesn[XLO] == Z->header->wesn[XLO] && L->header->wesn[YLO] == Z->header->wesn[YLO])) {
674 GMT_Report (API, GMT_MSG_ERROR, "Topo grid and ID grid must coregister\n");
675 Return (GMT_RUNTIME_ERROR);
676 }
677 if (GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_DATA_ONLY, NULL, Ctrl->L.file, L) == NULL) { /* Get ID data */
678 Return (API->error);
679 }
680
681 /* Store IDs in an int array instead */
682 ID = gmt_M_memory (GMT, NULL, L->header->size, int);
683 for (ij = 0; ij < L->header->size; ij++) ID[ij] = irint ((double)L->data[ij]);
684 if (GH->alloc_mode == GMT_ALLOC_INTERNALLY) gmt_M_free_aligned (GMT, L->data); /* Just free the array since we use ID; Grid struct is destroyed at end */
685 ID_info = gmt_M_memory (GMT, NULL, lrint (L->header->z_max) + 1, struct GRDSPOTTER_ID);
686 if (Ctrl->Q.mode == 1) { /* Only doing one CVA with no extra restrictions */
687 ID_info[Ctrl->Q.id].ok = true; /* Every other info in struct array is NULL or 0 */
688 }
689 else { /* Must read from file */
690 FILE *fp = NULL;
691 int Qid;
692 double wq, eq, sq, nq;
693 char line[GMT_BUFSIZ] = {""};
694
695 if ((fp = gmt_fopen (GMT, Ctrl->Q.file, "r")) == NULL) { /* Oh, oh... */
696 GMT_Report (API, GMT_MSG_ERROR, "-Q info file unreadable/nonexistent\n");
697 Return (GMT_ERROR_ON_FOPEN);
698 }
699 while (fgets (line, GMT_BUFSIZ, fp)) {
700 if (line[0] == '#') continue;
701 gmt_chop (line);
702 if (gmt_is_a_blank_line (line)) continue;
703 k = sscanf (line, "%*s %d %lf %lf %lf %lf", &Qid, &wq, &eq, &sq, &nq);
704 ID_info[Ctrl->Q.id].ok = true;
705 if (k == 5) { /* Got restricted wesn also */
706 ID_info[Qid].check_region = true;
707 sq = gmt_lat_swap (GMT, sq, GMT_LATSWAP_G2O); /* Go geocentric */
708 nq = gmt_lat_swap (GMT, nq, GMT_LATSWAP_G2O); /* Go geocentric */
709 ID_info[Qid].wesn[XLO] = wq * D2R;
710 ID_info[Qid].wesn[XHI] = eq * D2R;
711 ID_info[Qid].wesn[YLO] = sq * D2R;
712 ID_info[Qid].wesn[YHI] = nq * D2R;
713 }
714 }
715 fclose (fp);
716 }
717 Ctrl->Q.mode = true;
718 }
719
720 if (Ctrl->M.active) {
721 keep_flowlines = false; /* Do it the hard way to save memory */
722 k_step = 2; /* Reset back to 3 if needed later */
723 forth_flag = 10; /* The 10 is used to limit the flowline calculation to only resample track within the rectangular box of interest */
724 }
725 else {
726 keep_flowlines = (Ctrl->PA.active || Ctrl->D.active || Ctrl->W.active);
727 forth_flag = (Ctrl->PA.active) ? 11 : 10; /* The 10 is used to limit the flowline calculation to only resample track within the rectangular box of interest */
728 k_step = (Ctrl->PA.active) ? 3 : 2;
729 }
730 if (keep_flowlines) {
731 n_alloc = inc_alloc;
732 flowline = gmt_M_memory (GMT, NULL, n_alloc, struct FLOWLINE);
733 if (gmt_M_is_verbose (GMT, GMT_MSG_WARNING)) {
734 GMT_Report (API, GMT_MSG_WARNING, "Will attempt to keep all flowlines in memory. However, should this not be possible\n");
735 GMT_Report (API, GMT_MSG_WARNING, "then the program might crash. If so consider using the -M option.\n");
736 }
737 }
738
739 n_flow = n_nodes = 0;
740 gmt_M_grd_loop (GMT, Z, row, col, ij) { /* Loop over all input nodes */
741 /* STEP 1: Determine if z exceeds threshold and if so assign age */
742 if (gmt_M_is_fnan (Z->data[ij]) || Z->data[ij] < Ctrl->Z.min || Z->data[ij] > Ctrl->Z.max) continue; /* Skip node since it is NaN or outside the Ctrl->Z.min < z < Ctrl->Z.max range */
743 if (Ctrl->Q.mode && !ID_info[ID[ij]].ok) continue; /* Skip because of wrong ID */
744 if (!grdspotter_set_age (GMT, &t_smt, A, ij, Ctrl->N.t_upper, Ctrl->T.active[TRUNC])) continue; /* Skip if age is given and we are outside range */
745 /* STEP 2: Calculate this node's flowline */
746 if (Ctrl->Q.mode && ID_info[ID[ij]].check_region) /* Set up a box-limited flowline sampling */
747 gmt_M_memcpy (this_wesn, ID_info[ID[ij]].wesn, 4, double);
748 else
749 gmt_M_memcpy (this_wesn, wesn, 4, double);
750
751 np = grdspotter_get_flowline (GMT, x_smt[col], y_smt[row], t_smt, p, n_stages, sampling_int_in_km, k_step, forth_flag, this_wesn, &c);
752 if (np == 0) continue; /* No flowline inside this wesn */
753
754 /* STEP 3: Convolve this flowline with node shape and add to CVA grid */
755
756 /* Our convolution is approximate: We sample the flowline frequently and use
757 * one of the points on the flowline that are closest to the node. Ideally we
758 * want the nearest distance from each node to the flowline. Later versions may
759 * improve on this situation */
760
761 gmt_M_memset (processed_node, G->header->size, char); /* Fresh start for this flowline convolution */
762
763 if (keep_flowlines) {
764 flowline[n_nodes].node = gmt_M_memory (GMT, NULL, np, uint64_t);
765 if (Ctrl->PA.active) flowline[n_nodes].PA = gmt_M_memory (GMT, NULL, np, unsigned short);
766 }
767
768 /* Keep in mind that although first and last are entry/exit into the CVA grid, we must
769 * expect the flowline between those points to also exit/reentry the CVA grid; hence
770 * we must still check if points are in/out of the box. Here, we do not skip points but
771 * set the node index to UINTMAX_MAX */
772
773 cva_contribution = lat_area[row] * (Ctrl->T.active[UPPER] ? Ctrl->T.t_fix : Z->data[ij]); /* This node's contribution to the convolution */
774
775 #ifdef DEBUG2
776 printf ("> %" PRIu64 " %" PRIu64 " %d %d %g\n", n_nodes, np, row, col, cva_contribution);
777 #endif
778 for (m = 0, k = 1; m < np; m++) { /* Store nearest node indices only */
779 i = (int)gmt_M_grd_x_to_col (GMT, c[k++], G_rad->header);
780 yg = gmt_lat_swap (GMT, R2D * c[k++], GMT_LATSWAP_O2G); /* Convert back to geodetic */
781 j = (int)gmt_M_grd_y_to_row (GMT, yg, G->header);
782 if (i < 0 || (col2 = i) >= G->header->n_columns || j < 0 || (row2 = j) >= G->header->n_rows) /* Outside the CVA box, flag as outside */
783 node = UINTMAX_MAX;
784 else /* Inside the CVA box, assign node ij */
785 node = gmt_M_ijp (G->header, row2, col2);
786 if (keep_flowlines) {
787 flowline[n_nodes].node[m] = node;
788 if (Ctrl->PA.active) flowline[n_nodes].PA[m] = (unsigned short) lrint (c[k++] * T_2_PA);
789 }
790 /* If we did not keep flowlines then there is no PA to skip over (hence no k++) */
791 if (node != UINTMAX_MAX) {
792 if (!processed_node[node]) { /* Have not added to the CVA at this node yet */
793 G->data[node] += (gmt_grdfloat)cva_contribution;
794 processed_node[node] = true; /* Now we have visited this node */
795 n_unique_nodes++;
796 #ifdef DEBUG2
797 printf ("%g\t%g\n", x_cva[col2], y_cva[row2]);
798 #endif
799 }
800 n_more_than_once += 1.0;
801 }
802 }
803 if (keep_flowlines) {
804 flowline[n_nodes].n = np; /* Number of points in flowline */
805 flowline[n_nodes].ij = ij; /* Originating node in topo grid */
806 mem += sizeof (struct FLOWLINE) + np * sizeof (uint64_t) + ((Ctrl->PA.active) ? np * sizeof (unsigned short) : 0);
807 }
808
809 gmt_M_free (GMT, c); /* Free the flowline vector */
810
811 n_nodes++; /* Go to next node */
812
813 if (keep_flowlines && n_nodes == n_alloc) {
814 size_t old_n_alloc = n_alloc;
815 inc_alloc *= 2;
816 n_alloc += inc_alloc;
817 flowline = gmt_M_memory (GMT, flowline, n_alloc, struct FLOWLINE);
818 gmt_M_memset (&(flowline[old_n_alloc]), n_alloc - old_n_alloc, struct FLOWLINE); /* Set to NULL/0 */
819 }
820
821 if (!(n_nodes%100)) GMT_Report (API, GMT_MSG_INFORMATION, "Row %5ld Processed %5" PRIu64" nodes [%5ld/%.1f]\r", row, n_nodes, n_flow, mem * B_TO_MB);
822 }
823 GMT_Report (API, GMT_MSG_INFORMATION, "Row %5ld Processed %5" PRIu64 " nodes [%5ld/%.1f]\n", row, n_nodes, n_flow, mem * B_TO_MB);
824 GMT_Report (API, GMT_MSG_INFORMATION, "On average, each node was visited %g times\n", n_more_than_once / n_unique_nodes);
825
826 if (keep_flowlines && n_nodes != n_alloc) flowline = gmt_M_memory (GMT, flowline, n_nodes, struct FLOWLINE);
827
828 /* OK, Done processing, time to write out */
829
830 if (Ctrl->S.active) { /* Convert CVA values to percent of CVA maximum */
831 GMT_Report (API, GMT_MSG_INFORMATION, "Normalize CVS grid to percentages of max CVA\n");
832 grdspotter_normalize_grid (GMT, G, G->data);
833 }
834
835 GMT_Report (API, GMT_MSG_INFORMATION, "Write CVA grid %s\n", Ctrl->G.file);
836
837 if (GMT_Set_Comment (API, GMT_IS_GRID, GMT_COMMENT_IS_OPTION | GMT_COMMENT_IS_COMMAND, options, G)) {
838 gmt_M_free (GMT, y_smt); gmt_M_free (GMT, lat_area);
839 Return (API->error);
840 }
841 if (GMT_Write_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA, NULL, Ctrl->G.file, G) != GMT_NOERROR) {
842 gmt_M_free (GMT, y_smt); gmt_M_free (GMT, lat_area);
843 Return (API->error);
844 }
845
846 if (Ctrl->Z.mode) { /* Do CVA calculations for each z-slice using stored flowlines */
847 unsigned int layer, nz;
848 size_t len;
849 char file[PATH_MAX] = {""}, format[PATH_MAX] = {""};
850 double z0, z1;
851 gmt_grdfloat *CVA_inc = NULL, *old = G->data;
852
853 GMT_Report (API, GMT_MSG_INFORMATION, "Start z-slice CVA calculations\n");
854 for (len = strlen (Ctrl->G.file); len > 0 && Ctrl->G.file[len] != '.'; len--);
855 if (Ctrl->G.file[len] == '.') { /* Make a filename template from the CVA filename using the period as delimiter */
856 strncpy (format, Ctrl->G.file, len); /* Should keep the prefix from a file called prefix.ext */
857 strcat (format, "_%%d"); /* Make filenames like prefix_#.ext */
858 strncat (format, &Ctrl->G.file[len], PATH_MAX-1); /* Should add the extension from said file */
859 }
860 CVA_inc = gmt_M_memory (GMT, NULL, G->header->size, gmt_grdfloat);
861 nz = urint ((Ctrl->Z.max - Ctrl->Z.min) / Ctrl->Z.inc);
862 for (layer = 0; layer < nz; layer++) {
863 z0 = Ctrl->Z.min + layer * Ctrl->Z.inc;
864 z1 = z0 + Ctrl->Z.inc;
865 GMT_Report (API, GMT_MSG_INFORMATION, "Start z-slice %g - %g\n", z0, z1);
866 gmt_M_memset (CVA_inc, G->header->size, gmt_grdfloat); /* Fresh start for this z-slice */
867 for (m = 0; m < n_nodes; m++) { /* Loop over all active flowlines */
868 ij = flowline[m].ij;
869 if (Z->data[ij] <= z0 || Z->data[ij] > z1) continue; /* z outside current slice */
870 row = (unsigned int)gmt_M_row (Z->header, ij);
871 cva_contribution = lat_area[row] * (Ctrl->T.active[UPPER] ? Ctrl->T.t_fix : Z->data[ij]); /* This node's contribution to the convolution */
872 gmt_M_memset (processed_node, G->header->size, char); /* Fresh start for this flowline convolution */
873 for (k = 0; k < flowline[m].n; k++) { /* For each point along this flowline */
874 node = flowline[m].node[k];
875 if (node != UINTMAX_MAX && !processed_node[node]) { /* Have not added to the CVA at this node yet */
876 CVA_inc[node] += (gmt_grdfloat)cva_contribution;
877 processed_node[node] = true; /* Now we have visited this node */
878 }
879 }
880 if (!(m%10000)) GMT_Report (API, GMT_MSG_INFORMATION, "Processed %5ld flowlines\r", m);
881 }
882 GMT_Report (API, GMT_MSG_INFORMATION, "Processed %5" PRIu64 " flowlines\n", n_nodes);
883
884 /* Time to write out this z-slice grid */
885 if (Ctrl->S.active) { /* Convert CVA values to percent of CVA maximum */
886 GMT_Report (API, GMT_MSG_INFORMATION, "Normalize CVS grid to percentages of max CVA\n");
887 grdspotter_normalize_grid (GMT, G, CVA_inc);
888 }
889 snprintf (G->header->remark, GMT_GRID_REMARK_LEN160, "CVA for z-range %g - %g only", z0, z1);
890 sprintf (file, format, layer);
891 G->data = CVA_inc; /* Temporarily change the array pointer */
892 GMT_Report (API, GMT_MSG_INFORMATION, "Save z-slice CVA to file %s\n", file);
893 if (GMT_Set_Comment (API, GMT_IS_GRID, GMT_COMMENT_IS_OPTION | GMT_COMMENT_IS_COMMAND, options, G)) {
894 error = API->error;
895 gmt_M_free (GMT, CVA_inc);
896 goto END;
897 }
898 if (GMT_Write_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA, NULL, file, G) != GMT_NOERROR) {
899 error = API->error;
900 gmt_M_free (GMT, CVA_inc);
901 goto END;
902 }
903 }
904 G->data = old; /* Reset the array pointer */
905 gmt_M_free (GMT, CVA_inc);
906 }
907
908 if (Ctrl->D.active || Ctrl->PA.active) { /* Must determine max CVA along each flowline */
909 if (Ctrl->D.active) {
910 if ((DI = GMT_Duplicate_Data (API, GMT_IS_GRID, GMT_DUPLICATE_ALLOC, Z)) == NULL) {
911 error = API->error;
912 goto END;
913 }
914 }
915 if (Ctrl->PA.active) {
916 if ((PA = GMT_Duplicate_Data (API, GMT_IS_GRID, GMT_DUPLICATE_ALLOC, Z)) == NULL) {
917 error = API->error;
918 goto END;
919 }
920 }
921 GMT_Report (API, GMT_MSG_INFORMATION, "Compute DI and/or PA grids\n");
922
923 if (keep_flowlines) {
924 for (m = 0; m < n_nodes; m++) { /* Loop over all active flowlines */
925 CVA_max = 0.0; /* Fresh start for this flowline convolution */
926 for (k = 0; k < flowline[m].n; k++) { /* For each point along this flowline */
927 node = flowline[m].node[k];
928 if (node != UINTMAX_MAX && G->data[node] > CVA_max) { /* Found a higher CVA value */
929 CVA_max = G->data[node];
930 if (Ctrl->PA.active) pa = flowline[m].PA[k]; /* Keep the estimate of age at highest CVA */
931 }
932 }
933 if (Ctrl->D.active) DI->data[flowline[m].ij] = (gmt_grdfloat)CVA_max; /* Store the maximum CVA associated with this node's flowline */
934 if (Ctrl->PA.active) PA->data[flowline[m].ij] = (gmt_grdfloat) (pa * PA_2_T);
935 if (!(m%10000)) GMT_Report (API, GMT_MSG_INFORMATION, "Processed %5ld flowlines\r", m);
936 }
937 GMT_Report (API, GMT_MSG_INFORMATION, "Processed %5" PRIu64 " flowlines\n", n_nodes);
938 }
939 else { /* Must recreate flowlines */
940 k_step = 3; /* FLowlines have (x,y,t) here */
941 forth_flag = (Ctrl->PA.active) ? 11 : 10; /* The 10 is used to limit the flowline calculation to only resample track within the rectangular box of interest */
942 n_flow = n_nodes = 0;
943 gmt_M_grd_loop (GMT, Z, row, col, ij) { /* Loop over all input nodes */
944 /* STEP 1: Determine if z exceeds threshold and if so assign age */
945 if (gmt_M_is_fnan (Z->data[ij]) || Z->data[ij] < Ctrl->Z.min || Z->data[ij] > Ctrl->Z.max) continue; /* Skip node since it is NaN or outside the Ctrl->Z.min < z < Ctrl->Z.max range */
946 if (Ctrl->Q.mode && !ID_info[ID[ij]].ok) continue; /* Skip because of wrong ID */
947 if (!grdspotter_set_age (GMT, &t_smt, A, ij, Ctrl->N.t_upper, Ctrl->T.active[TRUNC])) continue; /* Skip as age is outside our range */
948 /* STEP 2: Calculate this node's flowline */
949 if (Ctrl->Q.mode && ID_info[ID[ij]].check_region) /* Set up a box-limited flowline sampling */
950 gmt_M_memcpy (this_wesn, ID_info[ID[ij]].wesn, 4, double);
951 else
952 gmt_M_memcpy (this_wesn, wesn, 4, double);
953 np = grdspotter_get_flowline (GMT, x_smt[col], y_smt[row], t_smt, p, n_stages, sampling_int_in_km, k_step, forth_flag, this_wesn, &c);
954 if (np == 0) continue; /* No flowline inside this wesn */
955 n_nodes++;
956 /* Fresh start for this flowline convolution */
957 CVA_max = 0.0;
958 this_pa = GMT->session.d_NaN;
959 for (m = 0, k = 1; m < np; m++) { /* Store nearest node indices only */
960 i = (int)gmt_M_grd_x_to_col (GMT, c[k++], G_rad->header);
961 if (i < 0 || (col = i) >= G->header->n_columns) { k += 2; continue;} /* Outside the CVA box, flag as outside */
962 yg = gmt_lat_swap (GMT, R2D * c[k++], GMT_LATSWAP_O2G); /* Convert back to geodetic */
963 j = (int)gmt_M_grd_y_to_row (GMT, yg, G->header);
964 if (j < 0 || (row = j) >= G->header->n_rows) { k++; continue;} /* Outside the CVA box, flag as outside */
965 if (Ctrl->PA.active) pa_val = c[k++];
966 node = gmt_M_ijp (G->header, row, col);
967 if (G->data[node] <= CVA_max) continue; /* Already seen higher CVA values */
968 CVA_max = G->data[node];
969 if (Ctrl->PA.active) this_pa = pa_val;
970 }
971 if (Ctrl->D.active) DI->data[ij] = (gmt_grdfloat)CVA_max; /* Store the maximum CVA associated with this node's flowline */
972 if (Ctrl->PA.active) PA->data[ij] = (gmt_grdfloat) this_pa;
973 gmt_M_free (GMT, c);
974 GMT_Report (API, GMT_MSG_INFORMATION, "Row %5ld: Processed %5" PRIu64 " flowlines\r", row, n_nodes);
975 }
976 GMT_Report (API, GMT_MSG_INFORMATION, "Row %5ld: Processed %5" PRIu64 " flowlines\n", row, n_nodes);
977 }
978
979
980 if (Ctrl->D.active) {
981 GMT_Report (API, GMT_MSG_INFORMATION, "Write DI grid %s\n", Ctrl->D.file);
982 snprintf (DI->header->remark, GMT_GRID_REMARK_LEN160, "CVA maxima along flowlines from each node");
983 if (GMT_Set_Comment (API, GMT_IS_GRID, GMT_COMMENT_IS_OPTION | GMT_COMMENT_IS_COMMAND, options, DI)) Return (API->error);
984 if (GMT_Write_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA, NULL, Ctrl->D.file, DI) != GMT_NOERROR) {
985 error = API->error;
986 goto END;
987 }
988 }
989 if (Ctrl->PA.active) {
990 GMT_Report (API, GMT_MSG_INFORMATION, "Write PA grid %s\n", Ctrl->PA.file);
991 snprintf (PA->header->remark, GMT_GRID_REMARK_LEN160, "Predicted age for each node");
992 if (GMT_Set_Comment (API, GMT_IS_GRID, GMT_COMMENT_IS_OPTION | GMT_COMMENT_IS_COMMAND, options, PA)) Return (API->error);
993 if (GMT_Write_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA, NULL, Ctrl->PA.file, PA) != GMT_NOERROR) {
994 error = API->error;
995 goto END;
996 }
997 }
998 }
999
1000 if (Ctrl->W.active) { /* Use bootstrapping to estimate confidence region for CVA maxima */
1001 struct GMT_RECORD *Out = gmt_new_record (GMT, out, NULL);
1002
1003 if (gmt_M_is_verbose (GMT, GMT_MSG_WARNING)) {
1004 GMT_Report (API, GMT_MSG_WARNING, "Preprocessed %5" PRIu64 " flowlines.\n", n_nodes);
1005 GMT_Report (API, GMT_MSG_WARNING, "%" PRIu64 " of %" PRIu64 " total flowlines entered CVA region.\n", n_nodes, n_flow);
1006 GMT_Report (API, GMT_MSG_WARNING, "Flowlines consumed %d Mb of memory.\n", lrint (mem * B_TO_MB));
1007 GMT_Report (API, GMT_MSG_WARNING, "Estimate %d CVA max locations using bootstrapping.\n", Ctrl->W.n_try);
1008 }
1009
1010 if ((error = GMT_Set_Columns (API, GMT_OUT, 3, GMT_COL_FIX_NO_TEXT)) != GMT_NOERROR) {
1011 goto END;
1012 }
1013 if (GMT_Init_IO (API, GMT_IS_DATASET, GMT_IS_POINT, GMT_OUT, GMT_ADD_DEFAULT, 0, options) != GMT_NOERROR) { /* Registers default output destination, unless already set */
1014 error = API->error;
1015 goto END;
1016 }
1017 if (GMT_Begin_IO (API, GMT_IS_DATASET, GMT_OUT, GMT_HEADER_ON) != GMT_NOERROR) { /* Enables data output and sets access mode */
1018 error = API->error;
1019 goto END;
1020 }
1021 if (GMT_Set_Geometry (API, GMT_OUT, GMT_IS_POINT) != GMT_NOERROR) { /* Sets output geometry */
1022 error = API->error;
1023 goto END;
1024 }
1025
1026 /* Now do bootstrap sampling of flowlines */
1027
1028 srand ((unsigned int)time(NULL)); /* Initialize random number generator */
1029 x_scale = (double)G->header->n_columns / (double)RAND_MAX;
1030 y_scale = (double)G->header->n_rows / (double)RAND_MAX;
1031 for (try = 1; try <= Ctrl->W.n_try; try++) {
1032 GMT_Report (API, GMT_MSG_INFORMATION, "Bootstrap try %d\r", try);
1033
1034 gmt_M_memset (G->data, G->header->size, gmt_grdfloat); /* Start with fresh grid */
1035 for (m = 0; m < n_nodes; m++) { /* Loop over all indices */
1036 row = urint (floor (gmt_rand(GMT) * y_scale)); /* Get a random integer in 0 to n_rows-1 range */
1037 col = urint (floor (gmt_rand(GMT) * x_scale)); /* Get a random integer in 0 to n_columns-1 range */
1038 ij = gmt_M_ijp (G->header, row, col); /* Get the node index */
1039 gmt_M_memset (processed_node, G->header->size, char); /* Fresh start for this flowline convolution */
1040 zz = Z->data[flowline[ij].ij];
1041 for (k = 0; k < flowline[ij].n; k++) { /* For each point along this flowline */
1042 node = flowline[ij].node[k];
1043 if (node != UINTMAX_MAX && !processed_node[node]) { /* Have not added to the CVA at this node yet */
1044 G->data[node] += zz;
1045 processed_node[node] = true; /* Now we have visited this node; flag it */
1046 }
1047 }
1048 if (!(m%10000)) GMT_Report (API, GMT_MSG_INFORMATION, "Processed %5ld flowlines\r", m);
1049 }
1050 GMT_Report (API, GMT_MSG_INFORMATION, "Processed %5" PRIu64 " flowlines\n", n_nodes);
1051
1052 /* Find max CVA location */
1053
1054 CVA_max = 0.0;
1055 for (ij = 0; ij < G->header->size; ij++) { /* Loop over all CVA nodes */
1056 if (G->data[ij] > CVA_max) { /* Update new max location */
1057 CVA_max = G->data[ij];
1058 max_ij = ij;
1059 }
1060 }
1061 col = (unsigned int)gmt_M_col (G->header, max_ij);
1062 row = (unsigned int)gmt_M_row (G->header, max_ij);
1063 out[0] = gmt_M_grd_col_to_x (GMT, col, G->header);
1064 out[1] = gmt_M_grd_row_to_y (GMT, row, G->header);
1065 out[2] = CVA_max;
1066
1067 GMT_Put_Record (API, GMT_WRITE_DATA, Out); /* Write this to output */
1068 }
1069 GMT_Report (API, GMT_MSG_INFORMATION, "Bootstrap try %d\n", Ctrl->W.n_try);
1070 gmt_M_free (GMT, Out);
1071 if (GMT_End_IO (API, GMT_OUT, 0) != GMT_NOERROR) { /* Disables further data output */
1072 error = API->error;
1073 goto END;
1074 }
1075 }
1076
1077 /* Clean up memory */
1078
1079 END:
1080 gmt_M_free (GMT, processed_node);
1081 for (m = 0; keep_flowlines && m < n_nodes; m++) {
1082 gmt_M_free (GMT, flowline[m].node);
1083 if (Ctrl->PA.active) gmt_M_free (GMT, flowline[m].PA);
1084 }
1085 gmt_M_free (GMT, p);
1086 gmt_M_free (GMT, x_smt); gmt_M_free (GMT, y_smt);
1087 gmt_M_free (GMT, lat_area);
1088 if (keep_flowlines) gmt_M_free (GMT, flowline);
1089 if (Ctrl->Q.mode) gmt_M_free (GMT, ID_info);
1090 if (GMT_Destroy_Data (API, &G_rad) != GMT_NOERROR) {
1091 GMT_Report (API, GMT_MSG_ERROR, "Failed to free G_rad\n");
1092 }
1093
1094 Return (GMT_NOERROR);
1095 }
1096