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 * grdpmodeler will read an age grid file and a plate motion model and
18 * evaluates chosen attributes such as speed, direction, lon/lat origin
19 * of the grid location.
20 *
21 * Author: Paul Wessel
22 * Date: 1-NOV-2010
23 * Ver: 5
24 */
25
26 #include "gmt_dev.h"
27 #include "spotter.h"
28
29 #define THIS_MODULE_CLASSIC_NAME "grdpmodeler"
30 #define THIS_MODULE_MODERN_NAME "grdpmodeler"
31 #define THIS_MODULE_LIB "spotter"
32 #define THIS_MODULE_PURPOSE "Evaluate a plate motion model on a geographic grid"
33 #define THIS_MODULE_KEYS "<G{,FD(,GG),>DG"
34 #define THIS_MODULE_NEEDS "g"
35 #define THIS_MODULE_OPTIONS "-:RVbdhor"
36
37 #define N_PM_ITEMS 9
38 #define PM_AZIM 0
39 #define PM_DIST 1
40 #define PM_STAGE 2
41 #define PM_VEL 3
42 #define PM_OMEGA 4
43 #define PM_DLON 5
44 #define PM_DLAT 6
45 #define PM_LON 7
46 #define PM_LAT 8
47
48 struct GRDROTATER_CTRL { /* All control options for this program (except common args) */
49 /* active is true if the option has been activated */
50 struct GRDPMODELER_In {
51 bool active;
52 char *file;
53 } In;
54 struct GRDPMODELER_E { /* -E[+]rotfile, -E[+]<ID1>-<ID2>, or -E<lon/lat/angle> */
55 bool active;
56 struct SPOTTER_ROT rot;
57 } E;
58 struct GRDPMODELER_F { /* -Fpolfile */
59 bool active;
60 char *file;
61 } F;
62 struct GRDPMODELER_G { /* -Goutfile */
63 bool active;
64 char *file;
65 } G;
66 struct GRDPMODELER_I { /* -I (for checking only) */
67 bool active;
68 } I;
69 struct GRDPMODELER_N { /* -N */
70 bool active;
71 double t_upper;
72 } N;
73 struct GRDPMODELER_S { /* -Sa|d|s|v|w|x|y|X|Y */
74 bool active, center;
75 unsigned int mode[N_PM_ITEMS];
76 unsigned int n_items;
77 } S;
78 struct GRDPMODELER_T { /* -T<fixtime> */
79 bool active;
80 double value;
81 } T;
82 };
83
New_Ctrl(struct GMT_CTRL * GMT)84 static void *New_Ctrl (struct GMT_CTRL *GMT) { /* Allocate and initialize a new control structure */
85 struct GRDROTATER_CTRL *C;
86
87 C = gmt_M_memory (GMT, NULL, 1, struct GRDROTATER_CTRL);
88
89 /* Initialize values whose defaults are not 0/false/NULL */
90
91 return (C);
92 }
93
Free_Ctrl(struct GMT_CTRL * GMT,struct GRDROTATER_CTRL * C)94 static void Free_Ctrl (struct GMT_CTRL *GMT, struct GRDROTATER_CTRL *C) { /* Deallocate control structure */
95 if (!C) return;
96 gmt_M_str_free (C->In.file);
97 gmt_M_str_free (C->E.rot.file);
98 gmt_M_str_free (C->F.file);
99 gmt_M_str_free (C->G.file);
100 gmt_M_free (GMT, C);
101 }
102
usage(struct GMTAPI_CTRL * API,int level)103 static int usage (struct GMTAPI_CTRL *API, int level) {
104 const char *name = gmt_show_name_and_purpose (API, THIS_MODULE_LIB, THIS_MODULE_CLASSIC_NAME, THIS_MODULE_PURPOSE);
105 if (level == GMT_MODULE_PURPOSE) return (GMT_NOERROR);
106 GMT_Usage (API, 0, "usage: %s %s %s [-F<polygontable>] [-G%s] [%s] [%s] [-N<upper_age>] "
107 "[-SadrswxyXY] [-T<time>] [%s] [%s] [%s] [%s] [%s] [%s]\n",
108 name, GMT_INGRID, SPOTTER_E_OPT, GMT_OUTGRID, GMT_Id_OPT, GMT_Rgeo_OPT, GMT_V_OPT, GMT_b_OPT, GMT_d_OPT, GMT_h_OPT, GMT_r_OPT, GMT_PAR_OPT);
109
110 if (level == GMT_SYNOPSIS) return (GMT_MODULE_SYNOPSIS);
111
112 GMT_Message (API, GMT_TIME_NONE, " REQUIRED ARGUMENTS:\n");
113 gmt_ingrid_syntax (API, 0, "Name of input grid in geographic coordinates with crustal ages");
114 spotter_rot_usage (API);
115 GMT_Message (API, GMT_TIME_NONE, "\n OPTIONAL ARGUMENTS:\n");
116 GMT_Usage (API, 1, "\n-F<polygontable>");
117 GMT_Usage (API, -2, "Specify a multi-segment closed polygon file that describes the area "
118 "of the grid to work on [Default works on the entire grid].");
119 gmt_outgrid_syntax (API, 'G', "Set output filename with the model predictions. "
120 "Must contain %%s if more than one item is specified in -S. "
121 "Default writes x,y,<predictions> to standard output.");
122 GMT_Usage (API, 1, "\n%s", GMT_Id_OPT);
123 GMT_Usage (API, -2, "Specify grid interval(s); Append m [or s] to <dx> and/or <dy> for minutes [or seconds].");
124 GMT_Usage (API, 1, "\n-N<upper_age>");
125 GMT_Usage (API, -2, "Extend earliest stage pole back to <upper_age> [no extension].");
126 GMT_Option (API, "Rg");
127 GMT_Usage (API, 1, "\n-SadrswxyXY");
128 GMT_Usage (API, -2, "Select one or more model predictions as a function of crustal age. Choose from:");
129 GMT_Usage (API, 3, "a: Plate spreading azimuth.");
130 GMT_Usage (API, 3, "d: Distance to origin of crust in km.");
131 GMT_Usage (API, 3, "r: Plate motion rate in mm/yr or km/Myr.");
132 GMT_Usage (API, 3, "s: Plate motion stage ID (1 is youngest).");
133 GMT_Usage (API, 3, "w: Rotation rate in degrees/Myr.");
134 GMT_Usage (API, 3, "x: Change in longitude since formation.");
135 GMT_Usage (API, 3, "y: Change in latitude since formation.");
136 GMT_Usage (API, 3, "X: Longitude at origin of crust.");
137 GMT_Usage (API, 3, "Y: Latitude at origin of crust.");
138 GMT_Usage (API, -2, "Default writes separate grids for adrswxyXY.");
139 GMT_Usage (API, 1, "\n-T<time>");
140 GMT_Usage (API, -2, "Set fixed time of reconstruction to override age grid.");
141 GMT_Option (API, "V,bi2,d,h,r,.");
142
143 return (GMT_MODULE_USAGE);
144 }
145
parse(struct GMT_CTRL * GMT,struct GRDROTATER_CTRL * Ctrl,struct GMT_OPTION * options)146 static int parse (struct GMT_CTRL *GMT, struct GRDROTATER_CTRL *Ctrl, struct GMT_OPTION *options) {
147 /* This parses the options provided to grdpmodeler and sets parameters in CTRL.
148 * Any GMT common options will override values set previously by other commands.
149 * It also replaces any file names specified as input or output with the data ID
150 * returned when registering these sources/destinations with the API.
151 */
152
153 unsigned int n_errors = 0, n_files = 0, k;
154 struct GMT_OPTION *opt = NULL;
155 struct GMTAPI_CTRL *API = GMT->parent;
156
157 for (opt = options; opt; opt = opt->next) {
158 switch (opt->option) {
159
160 case '<': /* Input files */
161 if (n_files++ > 0) {n_errors++; continue; }
162 Ctrl->In.active = true;
163 if (opt->arg[0]) Ctrl->In.file = strdup (opt->arg);
164 if (GMT_Get_FilePath (API, GMT_IS_DATASET, GMT_IN, GMT_FILE_REMOTE, &(Ctrl->In.file))) n_errors++;
165 break;
166
167 /* Supplemental parameters */
168
169 case 'E': /* File with stage poles */
170 n_errors += gmt_M_repeated_module_option (API, Ctrl->E.active);
171 Ctrl->E.active = true;
172 n_errors += spotter_parse (GMT, opt->option, opt->arg, &(Ctrl->E.rot));
173 break;
174 case 'F':
175 n_errors += gmt_M_repeated_module_option (API, Ctrl->F.active);
176 Ctrl->F.active = true;
177 if (opt->arg[0]) Ctrl->F.file = strdup (opt->arg);
178 if (GMT_Get_FilePath (API, GMT_IS_DATASET, GMT_IN, GMT_FILE_REMOTE, &(Ctrl->F.file))) n_errors++;
179 break;
180 case 'G':
181 n_errors += gmt_M_repeated_module_option (API, Ctrl->G.active);
182 Ctrl->G.active = true;
183 if (opt->arg[0]) Ctrl->G.file = strdup (opt->arg);
184 if (GMT_Get_FilePath (API, GMT_IS_GRID, GMT_OUT, GMT_FILE_LOCAL, &(Ctrl->G.file))) n_errors++;
185 break;
186 case 'I':
187 n_errors += gmt_M_repeated_module_option (API, Ctrl->I.active);
188 Ctrl->I.active = true;
189 n_errors += gmt_parse_inc_option (GMT, 'I', opt->arg);
190 break;
191 case 'N': /* Extend oldest stage back to this time [no extension] */
192 n_errors += gmt_M_repeated_module_option (API, Ctrl->N.active);
193 Ctrl->N.active = true;
194 Ctrl->N.t_upper = atof (opt->arg);
195 break;
196 case 'S':
197 n_errors += gmt_M_repeated_module_option (API, Ctrl->S.active);
198 Ctrl->S.active = true;
199 while (opt->arg[Ctrl->S.n_items]) {
200 switch (opt->arg[Ctrl->S.n_items]) {
201 case 'a': /* Plate spreading azimuth */
202 Ctrl->S.mode[Ctrl->S.n_items] = PM_AZIM; break;
203 case 'd': /* Distance from point to origin at ridge */
204 Ctrl->S.mode[Ctrl->S.n_items] = PM_DIST; break;
205 case 's': /* Plate motion stage ID */
206 Ctrl->S.mode[Ctrl->S.n_items] = PM_STAGE; break;
207 case 'v': case 'r': /* Plate spreading rate [r is backwards compatible] */
208 Ctrl->S.mode[Ctrl->S.n_items] = PM_VEL; break;
209 case 'w': /* Plate rotation omega */
210 Ctrl->S.mode[Ctrl->S.n_items] = PM_OMEGA; break;
211 case 'x': /* Change in longitude since origin */
212 Ctrl->S.mode[Ctrl->S.n_items] = PM_DLON; Ctrl->S.center = true; break;
213 case 'y': /* Change in latitude since origin */
214 Ctrl->S.mode[Ctrl->S.n_items] = PM_DLAT; break;
215 case 'X': /* Plate longitude at crust origin */
216 Ctrl->S.mode[Ctrl->S.n_items] = PM_LON; break;
217 case 'Y': /* Plate latitude at crust origin */
218 Ctrl->S.mode[Ctrl->S.n_items] = PM_LAT; break;
219 default:
220 n_errors++; break;
221 }
222 Ctrl->S.n_items++;
223 }
224 break;
225 case 'T':
226 n_errors += gmt_M_repeated_module_option (API, Ctrl->T.active);
227 Ctrl->T.active = true;
228 Ctrl->T.value = atof (opt->arg);
229 break;
230
231 default: /* Report bad options */
232 n_errors += gmt_default_error (GMT, opt->option);
233 break;
234 }
235 }
236
237 if (Ctrl->S.n_items == 0) { /* Set default which are all the items */
238 Ctrl->S.active = true;
239 Ctrl->S.n_items = N_PM_ITEMS;
240 for (k = 0; k < Ctrl->S.n_items; k++) Ctrl->S.mode[k] = k;
241 }
242
243 if (!Ctrl->In.file) { /* Must have -R -I [-r] */
244 n_errors += gmt_M_check_condition (GMT, !GMT->common.R.active[RSET] && !GMT->common.R.active[ISET], "Must specify input file or -R -I [-r]\n");
245 n_errors += gmt_M_check_condition (GMT, !Ctrl->T.active, "Must specify -T if no age grid is given.\n");
246 }
247 else { /* Must not have -I -r */
248 n_errors += gmt_M_check_condition (GMT, GMT->common.R.active[ISET] || GMT->common.R.active[GSET], "Cannot specify input file AND -R -r\n");
249 }
250 if (Ctrl->G.active) { /* Specified output grid(s) */
251 n_errors += gmt_M_check_condition (GMT, !Ctrl->G.file, "Option -G: Must specify output file\n");
252 n_errors += gmt_M_check_condition (GMT, Ctrl->S.n_items > 1 && !strstr (Ctrl->G.file, "%s"), "Option -G: File name must be a template containing \"%s\"\n");
253 }
254 else
255 n_errors += gmt_M_check_condition (GMT, !Ctrl->In.file, "Must specify input file when no output grids are created\n");
256 n_errors += gmt_M_check_condition (GMT, !Ctrl->E.active, "Must specify -E\n");
257 n_errors += gmt_M_check_condition (GMT, !Ctrl->S.active, "Must specify -S\n");
258 n_errors += gmt_M_check_condition (GMT, Ctrl->S.n_items == 0, "Must specify one or more fields with -S\n");
259 n_errors += gmt_M_check_condition (GMT, Ctrl->T.value < 0.0, "Option -T: Must specify positive age.\n");
260
261 return (n_errors ? GMT_PARSE_ERROR : GMT_NOERROR);
262 }
263
264 #define bailout(code) {gmt_M_free_options (mode); return (code);}
265 #define Return(code) {gmt_M_free (GMT, p); Free_Ctrl (GMT, Ctrl); gmt_end_module (GMT, GMT_cpy); bailout (code);}
266
GMT_grdpmodeler(void * V_API,int mode,void * args)267 EXTERN_MSC int GMT_grdpmodeler (void *V_API, int mode, void *args) {
268 unsigned int col, row, inside, stage, n_stages, registration, k;
269 int retval, error = 0;
270
271 bool skip, spotted;
272
273 uint64_t node, seg, n_old = 0, n_outside = 0, n_NaN = 0;
274
275 double lon = 0, lat = 0, d, value = 0.0, age, wesn[4], inc[2], *grd_x = NULL, *grd_y = NULL, *grd_yc = NULL, *out = NULL;
276
277 char *quantity[N_PM_ITEMS] = { "azimuth", "distance displacement", "stage", "velocity", "rotation rate", "longitude displacement", \
278 "latitude displacement", "reconstructed longitude", "reconstructed latitude"};
279 char *tag[N_PM_ITEMS] = { "az", "dist", "stage", "vel", "omega", "dlon", "dlat", "lon", "lat" };
280 struct GMT_DATASET *D = NULL;
281 struct GMT_DATATABLE *pol = NULL;
282 struct EULER *p = NULL; /* Pointer to array of stage poles */
283 struct GMT_OPTION *ptr = NULL;
284 struct GMT_GRID *G_age = NULL, **G_mod = NULL, *G = NULL;
285 struct GMT_RECORD *Out = NULL;
286 struct GRDROTATER_CTRL *Ctrl = NULL;
287 struct GMT_CTRL *GMT = NULL, *GMT_cpy = NULL;
288 struct GMT_OPTION *options = NULL;
289 struct GMTAPI_CTRL *API = gmt_get_api_ptr (V_API); /* Cast from void to GMTAPI_CTRL pointer */
290
291 /*----------------------- Standard module initialization and parsing ----------------------*/
292
293 if (API == NULL) return (GMT_NOT_A_SESSION);
294 if (mode == GMT_MODULE_PURPOSE) return (usage (API, GMT_MODULE_PURPOSE)); /* Return the purpose of program */
295 options = GMT_Create_Options (API, mode, args); if (API->error) return (API->error); /* Set or get option list */
296
297 if ((error = gmt_report_usage (API, options, 0, usage)) != GMT_NOERROR) bailout (error); /* Give usage if requested */
298
299 /* Parse the command-line arguments */
300
301 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 */
302 if ((ptr = GMT_Find_Option (API, 'f', options)) == NULL) gmt_parse_common_options (GMT, "f", 'f', "g"); /* Did not set -f, implicitly set -fg */
303 if (GMT_Parse_Common (API, THIS_MODULE_OPTIONS, options)) Return (API->error);
304 Ctrl = New_Ctrl (GMT); /* Allocate and initialize a new control structure */
305 if ((error = parse (GMT, Ctrl, options)) != 0) Return (error);
306
307 /*---------------------------- This is the grdpmodeler main code ----------------------------*/
308
309 /* Check limits and get data file */
310
311 if (Ctrl->In.file) { /* Gave an age grid */
312 if ((G_age = GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_ONLY, NULL, Ctrl->In.file, NULL)) == NULL) { /* Get header only */
313 Return (API->error);
314 }
315 gmt_M_memcpy (wesn, (GMT->common.R.active[RSET] ? GMT->common.R.wesn : G_age->header->wesn), 4, double);
316 if (GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_DATA_ONLY, wesn, Ctrl->In.file, G_age) == NULL) {
317 Return (API->error); /* Get header only */
318 }
319 gmt_M_memcpy (inc, G_age->header->inc, 2, double); /* Use same increment for output grid */
320 registration = G_age->header->registration;
321 gmt_M_memcpy (GMT->common.R.wesn, G_age->header->wesn, 4, double);
322 GMT->common.R.active[RSET] = true;
323 }
324 else { /* Use the input options of -R -I [and -r] */
325 gmt_M_memcpy (inc, GMT->common.R.inc, 2, double);
326 registration = GMT->common.R.registration;
327 }
328
329 if (Ctrl->F.active) { /* Read the user's clip polygon file */
330 if ((D = GMT_Read_Data (API, GMT_IS_DATASET, GMT_IS_FILE, GMT_IS_POLY, GMT_READ_NORMAL, NULL, Ctrl->F.file, NULL)) == NULL) {
331 Return (API->error);
332 }
333 pol = D->table[0]; /* Since it is a single file */
334 gmt_set_inside_mode (GMT, D, GMT_IOO_UNKNOWN);
335 GMT_Report (API, GMT_MSG_INFORMATION, "Restrict evaluation to within polygons in file %s\n", Ctrl->F.file);
336 }
337
338 if (Ctrl->E.rot.single) { /* Got a single rotation, no time, create a rotation table with one entry */
339 n_stages = 1;
340 p = gmt_M_memory (GMT, NULL, n_stages, struct EULER);
341 p[0].lon = Ctrl->E.rot.lon; p[0].lat = Ctrl->E.rot.lat; p[0].omega = Ctrl->E.rot.w;
342 if (gmt_M_is_dnan (Ctrl->E.rot.age)) { /* No age, use fake age = 1 everywhere */
343 Ctrl->T.active = true;
344 Ctrl->T.value = Ctrl->N.t_upper = p[0].t_start = 1.0;
345 }
346 spotter_setrot (GMT, &(p[0]));
347 }
348 else { /* Got a file or Gplates plate pair */
349 if ((retval = spotter_init (GMT, Ctrl->E.rot.file, &p, 0, false, Ctrl->E.rot.invert, &Ctrl->N.t_upper)) < 0)
350 Return (-retval);
351 n_stages = (unsigned int)retval;
352 }
353 for (stage = 0; stage < n_stages; stage++) {
354 if (p[stage].omega < 0.0) { /* Ensure all stages have positive rotation angles */
355 p[stage].omega = -p[stage].omega;
356 p[stage].lat = -p[stage].lat;
357 p[stage].lon += 180.0;
358 if (p[stage].lon > 360.0) p[stage].lon -= 360.0;
359 }
360 }
361 if (Ctrl->T.active && Ctrl->T.value > Ctrl->N.t_upper) {
362 GMT_Report (API, GMT_MSG_ERROR, "Requested a fixed reconstruction time outside range of rotation table\n");
363 Return (GMT_RUNTIME_ERROR);
364 }
365
366 if (Ctrl->G.active) { /* Need one or more output grids */
367 G_mod = gmt_M_memory (GMT, NULL, Ctrl->S.n_items, struct GMT_GRID *);
368 for (k = 0; k < Ctrl->S.n_items; k++) {
369 if ((G_mod[k] = GMT_Create_Data (API, GMT_IS_GRID, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA, NULL, NULL, inc, \
370 registration, GMT_NOTSET, NULL)) == NULL) { /* Free previous grids and bail */
371 unsigned int kk;
372 for (kk = 0; kk < k; kk++)
373 (void)GMT_Destroy_Data (API, &G_mod[kk]);
374 gmt_M_free (GMT, G_mod);
375 Return (API->error);
376 }
377
378 switch (Ctrl->S.mode[k]) {
379 case PM_AZIM: /* Compute plate motion direction at this point in time/space */
380 strcpy (G_mod[k]->header->z_units, "degree"); break;
381 case PM_DIST: /* Distance to origin in km */
382 strcpy (G_mod[k]->header->z_units, "km"); break;
383 case PM_STAGE: /* Compute plate motion stage at this point in time/space */
384 strcpy (G_mod[k]->header->z_units, "integer"); break;
385 case PM_VEL: /* Compute plate motion speed at this point in time/space */
386 strcpy (G_mod[k]->header->z_units, "mm/yr"); break;
387 case PM_OMEGA: /* Compute plate rotation rate omega */
388 strcpy (G_mod[k]->header->z_units, "degree/Myr"); break;
389 case PM_DLAT: /* Difference in latitude relative to where this point was formed in the model */
390 case PM_LAT: /* Latitude where this point was formed in the model */
391 strcpy (G_mod[k]->header->z_units, "degrees_north"); break;
392 case PM_DLON: /* Difference in longitude relative to where this point was formed in the model */
393 case PM_LON: /* Longitude where this point was formed in the model */
394 strcpy (G_mod[k]->header->z_units, "degrees_east"); break;
395 }
396 }
397 /* Just need on common set of x/y arrays; select G_mod[0] as our template */
398 G = G_mod[0];
399 GMT_Report (API, GMT_MSG_INFORMATION, "Evaluate %d model prediction grids based on %s\n", Ctrl->S.n_items, Ctrl->E.rot.file);
400 }
401 else { /* No output grids, must have input age grid to rely on */
402 G = G_age;
403 if (GMT_Init_IO (API, GMT_IS_DATASET, GMT_IS_POINT, GMT_OUT, GMT_ADD_DEFAULT, 0, options) != GMT_NOERROR) { /* Establishes data output */
404 Return (API->error);
405 }
406 if ((error = GMT_Set_Columns (API, GMT_OUT, Ctrl->S.n_items + 3, GMT_COL_FIX_NO_TEXT)) != GMT_NOERROR) {
407 Return (error);
408 }
409 if (GMT_Begin_IO (API, GMT_IS_DATASET, GMT_OUT, GMT_HEADER_ON) != GMT_NOERROR) { /* Enables data output and sets access mode */
410 Return (API->error);
411 }
412 if (GMT_Set_Geometry (API, GMT_OUT, GMT_IS_POINT) != GMT_NOERROR) { /* Sets output geometry */
413 Return (API->error);
414 }
415 GMT_Report (API, GMT_MSG_INFORMATION, "Evaluate %d model predictions based on %s\n", Ctrl->S.n_items, Ctrl->E.rot.file);
416 out = gmt_M_memory (GMT, NULL, Ctrl->S.n_items + 3, double);
417 Out = gmt_new_record (GMT, out, NULL); /* Since we only need to worry about numerics in this module */
418 }
419
420 grd_x = gmt_M_memory (GMT, NULL, G->header->n_columns, double);
421 grd_y = gmt_M_memory (GMT, NULL, G->header->n_rows, double);
422 grd_yc = gmt_M_memory (GMT, NULL, G->header->n_rows, double);
423 /* Precalculate node coordinates in both degrees and radians */
424 for (row = 0; row < G->header->n_rows; row++) {
425 grd_y[row] = gmt_M_grd_row_to_y (GMT, row, G->header);
426 grd_yc[row] = gmt_lat_swap (GMT, grd_y[row], GMT_LATSWAP_G2O);
427 }
428 for (col = 0; col < G->header->n_columns; col++) grd_x[col] = gmt_M_grd_col_to_x (GMT, col, G->header);
429
430 /* Loop over all nodes in the new rotated grid and find those inside the reconstructed polygon */
431
432
433 if (gmt_init_distaz (GMT, 'd', GMT_GREATCIRCLE, GMT_MAP_DIST) == GMT_NOT_A_VALID_TYPE) /* Great circle distances in degrees */
434 Return (GMT_NOT_A_VALID_TYPE);
435 if (Ctrl->S.center) GMT->current.io.geo.range = GMT_IS_M180_TO_P180_RANGE; /* Need +- around 0 here */
436
437 gmt_M_grd_loop (GMT, G, row, col, node) {
438 skip = false;
439 if (Ctrl->F.active) { /* Use the bounding polygon */
440 for (seg = inside = 0; seg < pol->n_segments && !inside; seg++) { /* Use degrees since function expects it */
441 if (gmt_polygon_is_hole (GMT, pol->segment[seg])) continue; /* Holes are handled within gmt_inonout */
442 inside = (gmt_inonout (GMT, grd_x[col], grd_y[row], pol->segment[seg]) > GMT_OUTSIDE);
443 }
444 if (!inside) skip = true, n_outside++; /* Outside the polygon(s); set all output grids to NaN for this node */
445 }
446 /* Get age for this node (or the constant age) */
447 if (Ctrl->T.active)
448 age = Ctrl->T.value;
449 else {
450 age = G_age->data[node];
451 if (gmt_M_is_dnan (age)) skip = true, n_NaN++; /* No crustal age */
452 else if (age > Ctrl->N.t_upper) skip = true, n_old++; /* Outside of model range */
453 }
454 if (skip) {
455 if (Ctrl->G.active) for (k = 0; k < Ctrl->S.n_items; k++) G_mod[k]->data[node] = GMT->session.f_NaN;
456 continue;
457 }
458 /* Here we are inside; get the coordinates and rotate back to original grid coordinates */
459 if ((retval = spotter_stage (GMT, age, p, n_stages)) < 0) continue; /* Outside valid stage rotation range */
460 stage = retval; /* Current rotation stage */
461 spotted = false; /* Not yet called spotter_backtrack at this node */
462 if (!Ctrl->G.active) { /* Need x,y,t for rec-by-rec output */
463 out[GMT_X] = grd_x[col]; out[GMT_Y] = grd_y[row]; out[GMT_Z] = age;
464 }
465 for (k = 0; k < Ctrl->S.n_items; k++) {
466 switch (Ctrl->S.mode[k]) {
467 case PM_AZIM: /* Compute plate motion direction at this point in time/space */
468 value = gmt_az_backaz (GMT, grd_x[col], grd_yc[row], p[stage].lon, p[stage].lat, false) - 90.0;
469 gmt_lon_range_adjust (GMT->current.io.geo.range, &value);
470 break;
471 case PM_DIST: /* Compute great-circle distance between node and point of origin at ridge */
472 if (!spotted) {
473 lon = grd_x[col] * D2R; lat = grd_yc[row] * D2R;
474 if (spotter_backtrack (GMT, &lon, &lat, &age, 1U, p, n_stages, 0.0, 0.0, 0, NULL, NULL) < 0)
475 Return (GMT_RUNTIME_ERROR);
476 spotted = true;
477 }
478 value = GMT->current.proj.DIST_KM_PR_DEG * gmt_distance (GMT, grd_x[col], grd_yc[row], lon * R2D, lat * R2D);
479 break;
480 case PM_STAGE: /* Compute plate rotation stage */
481 value = stage;
482 break;
483 case PM_VEL: /* Compute plate motion speed at this point in time/space */
484 d = gmt_distance (GMT, grd_x[col], grd_yc[row], p[stage].lon, p[stage].lat);
485 value = sind (d) * p[stage].omega * GMT->current.proj.DIST_KM_PR_DEG; /* km/Myr or mm/yr */
486 break;
487 case PM_OMEGA: /* Compute plate rotation rate omega */
488 value = p[stage].omega; /* degree/Myr */
489 break;
490 case PM_DLON: /* Compute latitude where this point was formed in the model */
491 if (!spotted) {
492 lon = grd_x[col] * D2R; lat = grd_yc[row] * D2R;
493 if (spotter_backtrack (GMT, &lon, &lat, &age, 1U, p, n_stages, 0.0, 0.0, 0, NULL, NULL) < 0)
494 Return (GMT_RUNTIME_ERROR);
495 spotted = true;
496 }
497 value = grd_x[col] - lon * R2D;
498 if (fabs (value) > 180.0) value = copysign (360.0 - fabs (value), -value);
499 break;
500 case PM_DLAT: /* Compute latitude where this point was formed in the model */
501 if (!spotted) {
502 lon = grd_x[col] * D2R; lat = grd_yc[row] * D2R;
503 if (spotter_backtrack (GMT, &lon, &lat, &age, 1U, p, n_stages, 0.0, 0.0, 0, NULL, NULL) < 0)
504 Return (GMT_RUNTIME_ERROR);
505 spotted = true;
506 }
507 value = grd_y[row] - gmt_lat_swap (GMT, lat * R2D, GMT_LATSWAP_O2G); /* Convert back to geodetic */
508 break;
509 case PM_LON: /* Compute latitude where this point was formed in the model */
510 if (!spotted) {
511 lon = grd_x[col] * D2R; lat = grd_yc[row] * D2R;
512 if (spotter_backtrack (GMT, &lon, &lat, &age, 1U, p, n_stages, 0.0, 0.0, 0, NULL, NULL) < 0)
513 Return (GMT_RUNTIME_ERROR);
514 spotted = true;
515 }
516 value = lon * R2D;
517 break;
518 case PM_LAT: /* Compute latitude where this point was formed in the model */
519 if (!spotted) {
520 lon = grd_x[col] * D2R; lat = grd_yc[row] * D2R;
521 if (spotter_backtrack (GMT, &lon, &lat, &age, 1U, p, n_stages, 0.0, 0.0, 0, NULL, NULL) < 0)
522 Return (GMT_RUNTIME_ERROR);
523 spotted = true;
524 }
525 value = gmt_lat_swap (GMT, lat * R2D, GMT_LATSWAP_O2G); /* Convert back to geodetic */
526 break;
527 }
528 if (Ctrl->G.active)
529 G_mod[k]->data[node] = (gmt_grdfloat)value;
530 else
531 out[k+3] = value;
532 }
533 if (!Ctrl->G.active) GMT_Put_Record (API, GMT_WRITE_DATA, Out);
534 }
535
536 if (n_outside) GMT_Report (API, GMT_MSG_INFORMATION, "%" PRIu64 " points fell outside the polygonal boundary\n", n_outside);
537 if (n_old) GMT_Report (API, GMT_MSG_INFORMATION, "%" PRIu64 " points had ages that exceeded the limit of the rotation model\n", n_old);
538 if (n_NaN) GMT_Report (API, GMT_MSG_INFORMATION, "%" PRIu64 " points had ages that were NaN\n", n_NaN);
539 if (Ctrl->G.active) { /* Need one or more output grids */
540 /* Now write model prediction grid(s) */
541 char file[PATH_MAX] = {""};
542 for (k = 0; k < Ctrl->S.n_items; k++) {
543 sprintf (file, Ctrl->G.file, tag[Ctrl->S.mode[k]]);
544 GMT_Report (API, GMT_MSG_INFORMATION, "Write model prediction grid for %s (%s) to file %s\n", quantity[Ctrl->S.mode[k]], G_mod[k]->header->z_units, file);
545 strcpy (G_mod[k]->header->x_units, "degrees_east");
546 strcpy (G_mod[k]->header->y_units, "degrees_north");
547 snprintf (G_mod[k]->header->remark, GMT_GRID_REMARK_LEN160, "Plate Model predictions of %s for model %s", quantity[Ctrl->S.mode[k]], Ctrl->E.rot.file);
548 if (GMT_Set_Comment (API, GMT_IS_GRID, GMT_COMMENT_IS_OPTION | GMT_COMMENT_IS_COMMAND, options, G_mod[k])) Return (API->error);
549 if (GMT_Write_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA, NULL, file, G_mod[k]) != GMT_NOERROR) {
550 Return (API->error);
551 }
552 }
553 gmt_M_free (GMT, G_mod);
554 }
555 else {
556 if (GMT_End_IO (API, GMT_OUT, 0) != GMT_NOERROR) { /* Disables further data output */
557 Return (API->error);
558 }
559 gmt_M_free (GMT, out);
560 gmt_M_free (GMT, Out);
561 }
562
563 gmt_M_free (GMT, grd_x);
564 gmt_M_free (GMT, grd_y);
565 gmt_M_free (GMT, grd_yc);
566
567 Return (GMT_NOERROR);
568 }
569