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 * grdrotater will read a grid file, apply a finite rotation to the grid
18 * coordinates, and then interpolate the old grid at the new coordinates.
19 *
20 * Author: Paul Wessel
21 * Date: 27-JAN-2006
22 * Ver: 1
23 */
24
25 #include "gmt_dev.h"
26 #include "spotter.h"
27
28 #define THIS_MODULE_CLASSIC_NAME "grdrotater"
29 #define THIS_MODULE_MODERN_NAME "grdrotater"
30 #define THIS_MODULE_LIB "spotter"
31 #define THIS_MODULE_PURPOSE "Finite rotation reconstruction of geographic grid"
32 #define THIS_MODULE_KEYS "<G{,FD(,GG},TD("
33 #define THIS_MODULE_NEEDS "g"
34 #define THIS_MODULE_OPTIONS "-:>RVbdfhno" GMT_OPT("HMmQ")
35
36 #define PAD 3 /* Used to polish up a rotated grid by checking near neighbor nodes */
37
38 struct GRDROTATER_CTRL { /* All control options for this program (except common args) */
39 /* active is true if the option has been activated */
40 struct GRDROTATER_In {
41 bool active;
42 char *file;
43 } In;
44 struct GRDROTATER_A { /* -Aw/e/s/n */
45 bool active;
46 double wesn[4];
47 } A;
48 struct GRDROTATER_D { /* -Drotpolfile or -Dtemplate */
49 bool active;
50 char *file;
51 } D;
52 struct GRDROTATER_E { /* -E[+]rotfile, -E[+]<ID1>-<ID2>, or -E<lon/lat/angle> */
53 bool active;
54 struct SPOTTER_ROT rot;
55 } E;
56 struct GRDROTATER_F { /* -Fpolfile */
57 bool active;
58 char *file;
59 } F;
60 struct GRDROTATER_G { /* -Goutfile or -Gtemplate*/
61 bool active;
62 char *file;
63 } G;
64 struct GRDROTATER_N { /* -N */
65 bool active;
66 } N;
67 struct GRDROTATER_S { /* -S */
68 bool active;
69 } S;
70 struct GRDROTATER_T { /* -T<time>, -T<start/stop/inc> or -T<tfile.txt> */
71 bool active;
72 unsigned int n_times; /* Number of reconstruction times */
73 double *value; /* Array with one or more reconstruction times */
74 } T;
75 };
76
New_Ctrl(struct GMT_CTRL * GMT)77 static void *New_Ctrl (struct GMT_CTRL *GMT) { /* Allocate and initialize a new control structure */
78 struct GRDROTATER_CTRL *C;
79
80 C = gmt_M_memory (GMT, NULL, 1, struct GRDROTATER_CTRL);
81
82 return (C);
83 }
84
Free_Ctrl(struct GMT_CTRL * GMT,struct GRDROTATER_CTRL * C)85 static void Free_Ctrl (struct GMT_CTRL *GMT, struct GRDROTATER_CTRL *C) { /* Deallocate control structure */
86 if (!C) return;
87 gmt_M_str_free (C->In.file);
88 gmt_M_str_free (C->D.file);
89 gmt_M_str_free (C->E.rot.file);
90 gmt_M_str_free (C->F.file);
91 gmt_M_str_free (C->G.file);
92 gmt_M_str_free (C->T.value);
93 gmt_M_free (GMT, C);
94 }
95
usage(struct GMTAPI_CTRL * API,int level)96 static int usage (struct GMTAPI_CTRL *API, int level) {
97 const char *name = gmt_show_name_and_purpose (API, THIS_MODULE_LIB, THIS_MODULE_CLASSIC_NAME, THIS_MODULE_PURPOSE);
98 if (level == GMT_MODULE_PURPOSE) return (GMT_NOERROR);
99 GMT_Usage (API, 0, "usage: %s %s %s -G%s [-F<polygontable>] [-A<region>] [-D<rotoutline>] [-N] "
100 "[%s] [-S] [-T<time>] [%s] [%s] [%s] [%s] [%s] [%s] [%s] [%s]\n", name, GMT_INGRID, SPOTTER_E_OPT, GMT_OUTGRID, GMT_Rgeo_OPT,
101 GMT_V_OPT, GMT_b_OPT, GMT_d_OPT, GMT_f_OPT, GMT_h_OPT, GMT_n_OPT, GMT_o_OPT, GMT_PAR_OPT);
102
103 if (level == GMT_SYNOPSIS) return (GMT_MODULE_SYNOPSIS);
104
105 GMT_Message (API, GMT_TIME_NONE, " REQUIRED ARGUMENTS:\n");
106 gmt_ingrid_syntax (API, 0, "Name of grid in geographic coordinates to be rotated");
107 gmt_outgrid_syntax (API, 'G', "Set output filename for the new, rotated grid. The boundary of the "
108 "original grid (or a subset; see -F) after rotation is written to standard output (but see -D) "
109 "unless the grid is global. If more than one reconstruction time is chosen "
110 "then -D is required unless -N is used and <outgrid> must be a filename template "
111 "containing a C-format specifier for formatting a double (for the variable time).");
112 spotter_rot_usage (API);
113 GMT_Message (API, GMT_TIME_NONE, "\n OPTIONAL ARGUMENTS:\n");
114 GMT_Usage (API, 1, "\n-A<region>");
115 GMT_Usage (API, -2, "Set the west/east/south/north bounds for the rotated grid [Default will "
116 "determine the natural extent of the rotated grid instead].");
117 GMT_Usage (API, 1, "\n-D<rotoutline>");
118 GMT_Usage (API, -2, "Write the rotated polygon or grid outline to <rotoutline> [standard output]. "
119 "Required if more than one reconstruction time is chosen and -N is not set "
120 "and must then contain a C-format specifier for formatting a double (for the variable time).");
121 GMT_Usage (API, 1, "\n-F<polygontable>");
122 GMT_Usage (API, -2, "Specify a multi-segment closed polygon table that describes the area of the grid "
123 "that should be projected [Default projects entire grid].");
124 GMT_Usage (API, 1, "\n-N Do NOT output the rotated polygon or grid outlines.");
125 GMT_Option (API, "Rg");
126 GMT_Usage (API, 1, "\n-S Do NOT rotate the grid - just produce the rotated outlines (requires -D).");
127 GMT_Usage (API, 1, "\n-T<time>");
128 GMT_Usage (API, -2, "Set the time(s) of reconstruction. Append a single time (-T<time>), "
129 "an equidistant range of times (-T<min>/<max>/<inc> or -T<min>/<max>/<npoints>+n), "
130 "or the name of a file with a list of times (-T<tfile>). If no -T is set "
131 "then the reconstruction times equal the rotation times given in -E.");
132 GMT_Option (API, "V,bi2,bo,d,f,h,n,:,.");
133
134 return (GMT_MODULE_USAGE);
135
136 }
137
parse(struct GMT_CTRL * GMT,struct GRDROTATER_CTRL * Ctrl,struct GMT_OPTION * options)138 static int parse (struct GMT_CTRL *GMT, struct GRDROTATER_CTRL *Ctrl, struct GMT_OPTION *options) {
139 /* This parses the options provided to grdrotater and sets parameters in CTRL.
140 * Any GMT common options will override values set previously by other commands.
141 * It also replaces any file names specified as input or output with the data ID
142 * returned when registering these sources/destinations with the API.
143 */
144
145 unsigned int n_errors = 0, k, n_files = 0;
146 int n_items;
147 uint64_t t = 0;
148 bool gave_e = false;
149 char txt_a[GMT_LEN256] = {""}, txt_b[GMT_LEN256] = {""}, txt_c[GMT_LEN256] = {""};
150 char txt[4][GMT_LEN64];
151 struct GMT_OPTION *opt = NULL;
152 struct GMTAPI_CTRL *API = GMT->parent;
153
154 for (opt = options; opt; opt = opt->next) if (opt->option == 'E' || opt->option == 'e') gave_e = true; /* Pre-check to see if GMT4 or newer syntax */
155
156 for (opt = options; opt; opt = opt->next) {
157 switch (opt->option) {
158
159 case '<': /* Input files */
160 if (n_files++ > 0) break;
161 Ctrl->In.active = true;
162 if (opt->arg[0]) Ctrl->In.file = strdup (opt->arg);
163 if (GMT_Get_FilePath (API, GMT_IS_GRID, GMT_IN, GMT_FILE_REMOTE, &(Ctrl->In.file))) n_errors++;
164 break;
165
166 /* Supplemental parameters */
167
168 case 'A':
169 n_errors += gmt_M_repeated_module_option (API, Ctrl->A.active);
170 Ctrl->A.active = true;
171 if (opt->arg[0] == 'g' && opt->arg[1] == '\0') { /* Got -Ag */
172 Ctrl->A.wesn[0] = 0.0; Ctrl->A.wesn[1] = 360.0; Ctrl->A.wesn[2] = -90.0; Ctrl->A.wesn[3] = 90.0;
173 break;
174 }
175 if (opt->arg[0] == 'd' && opt->arg[1] == '\0') { /* Got -Rd */
176 Ctrl->A.wesn[0] = -180.0; Ctrl->A.wesn[1] = 180.0; Ctrl->A.wesn[2] = -90.0; Ctrl->A.wesn[3] = 90.0;
177 break;
178 }
179 if (!gmt_access (GMT, opt->arg, R_OK)) { /* Gave a readable file, presumably a grid */
180 struct GMT_GRID *G = NULL;
181 if ((G = GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_ONLY, NULL, opt->arg, NULL)) == NULL) { /* Get header only */
182 return (API->error);
183 }
184 Ctrl->A.wesn[0] = G->header->wesn[XLO]; Ctrl->A.wesn[1] = G->header->wesn[XHI];
185 Ctrl->A.wesn[2] = G->header->wesn[YLO]; Ctrl->A.wesn[3] = G->header->wesn[YHI];
186 if (GMT_Destroy_Data (API, &G) != GMT_NOERROR) {
187 return (API->error);
188 }
189 break;
190 }
191 /* Only get here if the above cases did not trip */
192 n_items = sscanf (opt->arg, "%[^/]/%[^/]/%[^/]/%s", txt[0], txt[1], txt[2], txt[3]);
193 if (n_items != 4) {
194 GMT_Report (API, GMT_MSG_ERROR, "Option -A: Give g, d, <grdfile>, or <west/east/south/north>\n");
195 n_errors++;
196 }
197 n_errors += gmt_verify_expectations (GMT, gmt_M_type (GMT, GMT_IN, GMT_X), gmt_scanf_arg (GMT, txt[0], gmt_M_type (GMT, GMT_IN, GMT_X), true, &Ctrl->A.wesn[0]), txt[0]);
198 n_errors += gmt_verify_expectations (GMT, gmt_M_type (GMT, GMT_IN, GMT_X), gmt_scanf_arg (GMT, txt[1], gmt_M_type (GMT, GMT_IN, GMT_X), true, &Ctrl->A.wesn[1]), txt[1]);
199 n_errors += gmt_verify_expectations (GMT, gmt_M_type (GMT, GMT_IN, GMT_Y), gmt_scanf_arg (GMT, txt[2], gmt_M_type (GMT, GMT_IN, GMT_Y), true, &Ctrl->A.wesn[2]), txt[2]);
200 n_errors += gmt_verify_expectations (GMT, gmt_M_type (GMT, GMT_IN, GMT_Y), gmt_scanf_arg (GMT, txt[3], gmt_M_type (GMT, GMT_IN, GMT_Y), true, &Ctrl->A.wesn[3]), txt[3]);
201 break;
202 case 'C': /* Now done automatically in spotter_init */
203 if (gmt_M_compat_check (GMT, 4))
204 GMT_Report (API, GMT_MSG_COMPAT, "-C is no longer needed as total reconstruction vs stage rotation is detected automatically.\n");
205 else
206 n_errors += gmt_default_error (GMT, opt->option);
207 break;
208 case 'D':
209 n_errors += gmt_M_repeated_module_option (API, Ctrl->D.active);
210 Ctrl->D.active = true;
211 Ctrl->D.file = strdup (opt->arg);
212 break;
213 case 'e':
214 GMT_Report (API, GMT_MSG_COMPAT, "-e is deprecated and was removed in 5.3. Use -E instead.\n");
215 /* Intentionally fall through */
216 case 'E': /* File with stage poles or a single rotation pole */
217 n_errors += gmt_M_repeated_module_option (API, Ctrl->E.active);
218 Ctrl->E.active = true;
219 n_errors += spotter_parse (GMT, opt->option, opt->arg, &(Ctrl->E.rot));
220 break;
221 case 'F':
222 n_errors += gmt_M_repeated_module_option (API, Ctrl->F.active);
223 Ctrl->F.active = true;
224 if (opt->arg[0]) Ctrl->F.file = strdup (opt->arg);
225 if (GMT_Get_FilePath (API, GMT_IS_DATASET, GMT_IN, GMT_FILE_REMOTE, &(Ctrl->F.file))) n_errors++;
226 break;
227 case 'G':
228 n_errors += gmt_M_repeated_module_option (API, Ctrl->G.active);
229 Ctrl->G.active = true;
230 if (opt->arg[0]) Ctrl->G.file = strdup (opt->arg);
231 if (GMT_Get_FilePath (API, GMT_IS_GRID, GMT_OUT, GMT_FILE_LOCAL, &(Ctrl->G.file))) n_errors++;
232 break;
233 case 'N':
234 n_errors += gmt_M_repeated_module_option (API, Ctrl->N.active);
235 Ctrl->N.active = true;
236 break;
237 case 'S':
238 n_errors += gmt_M_repeated_module_option (API, Ctrl->S.active);
239 Ctrl->S.active = true;
240 break;
241 case 'T': /* New: -Tage, -Tmin/max/inc, -Tmin/max/n+, -Tfile; compat mode: -Tlon/lat/angle Finite rotation parameters */
242 n_errors += gmt_M_repeated_module_option (API, Ctrl->T.active);
243 Ctrl->T.active = true;
244 if (!gmt_access (GMT, opt->arg, R_OK)) { /* Gave a file with times in first column */
245 uint64_t seg, row;
246 struct GMT_DATASET *T = NULL;
247 struct GMT_DATASEGMENT *S = NULL;
248 if ((T = GMT_Read_Data (API, GMT_IS_DATASET, GMT_IS_FILE, GMT_IS_NONE, GMT_READ_NORMAL, NULL, opt->arg, NULL)) == NULL) {
249 GMT_Report (API, GMT_MSG_ERROR, "Failure while reading file %s\n", opt->arg);
250 n_errors++;
251 continue;
252 }
253 /* Single table, build t array */
254 Ctrl->T.n_times = (unsigned int)T->n_records;
255 Ctrl->T.value = gmt_M_memory (GMT, NULL, Ctrl->T.n_times, double);
256 for (seg = t = 0; seg < T->table[0]->n_segments; seg++) {
257 S = T->table[0]->segment[seg]; /* Shorthand to current segment */
258 for (row = 0; row < S->n_rows; seg++, t++)
259 Ctrl->T.value[t] = S->data[GMT_X][row];
260 }
261 if (GMT_Destroy_Data (API, &T) != GMT_NOERROR)
262 n_errors++;
263 break;
264 }
265 /* Not a file */
266 k = sscanf (opt->arg, "%[^/]/%[^/]/%s", txt_a, txt_b, txt_c);
267 if (k == 3) { /* Gave -Ttstart/tstop/tinc[+n] or possibly ancient -Tlon/lat/angle */
268 if (gmt_M_compat_check (GMT, 4) && !gave_e) { /* No -E|e so likely ancient syntax */
269 GMT_Report (API, GMT_MSG_COMPAT, "-T<lon>/<lat>/<angle> is deprecated; use -E<lon>/<lat>/<angle> instead.\n");
270 Ctrl->E.rot.single = Ctrl->E.active = true;
271 Ctrl->T.active = false;
272 Ctrl->E.rot.w = atof (txt_c);
273 Ctrl->E.rot.age = GMT->session.d_NaN;
274 n_errors += gmt_verify_expectations (GMT, gmt_M_type (GMT, GMT_IN, GMT_X), gmt_scanf_arg (GMT, txt_a, gmt_M_type (GMT, GMT_IN, GMT_X), false, &Ctrl->E.rot.lon), txt_a);
275 n_errors += gmt_verify_expectations (GMT, gmt_M_type (GMT, GMT_IN, GMT_Y), gmt_scanf_arg (GMT, txt_b, gmt_M_type (GMT, GMT_IN, GMT_Y), false, &Ctrl->E.rot.lat), txt_b);
276 }
277 else { /* Must be Ttstart/tstop/tinc */
278 double min, max, inc;
279 min = atof (txt_a); max = atof (txt_b); inc = atof (txt_c);
280 if (strstr(opt->arg, "+n") || opt->arg[strlen(opt->arg)-1] == '+') /* Gave number of points instead; calculate inc */
281 inc = (max - min) / (inc - 1.0);
282 if (inc <= 0.0) {
283 GMT_Report (API, GMT_MSG_ERROR, "Option -T: Age increment must be positive\n");
284 n_errors++;
285 }
286 else {
287 Ctrl->T.n_times = lrint ((max - min) / inc) + 1;
288 Ctrl->T.value = gmt_M_memory (GMT, NULL, Ctrl->T.n_times, double);
289 for (t = 0; t < Ctrl->T.n_times; t++) Ctrl->T.value[t] = (t == (Ctrl->T.n_times-1)) ? max: min + t * inc;
290 }
291 }
292 }
293 else { /* Got a single time */
294 Ctrl->T.n_times = 1;
295 Ctrl->T.value = gmt_M_memory (GMT, NULL, Ctrl->T.n_times, double);
296 Ctrl->T.value[0] = atof (txt_a);
297 }
298 break;
299
300 default: /* Report bad options */
301 n_errors += gmt_default_error (GMT, opt->option);
302 break;
303 }
304 }
305
306 if (GMT->common.b.active[GMT_IN] && GMT->common.b.ncol[GMT_IN] == 0) GMT->common.b.ncol[GMT_IN] = 2;
307 n_errors += gmt_M_check_condition (GMT, Ctrl->S.active && Ctrl->G.active, "No output grid file allowed with -S\n");
308 n_errors += gmt_M_check_condition (GMT, Ctrl->S.active && Ctrl->N.active, "Cannot use -N with -S\n");
309 n_errors += gmt_M_check_condition (GMT, !Ctrl->S.active && !Ctrl->In.file, "Must specify input file\n");
310 n_errors += gmt_M_check_condition (GMT, !Ctrl->S.active && !Ctrl->G.file, "Option -G: Must specify output file\n");
311 n_errors += gmt_M_check_condition (GMT, Ctrl->S.active && Ctrl->N.active, "-N and -S cannot both be given\n");
312 n_errors += gmt_M_check_condition (GMT, GMT->common.b.active[GMT_IN] && GMT->common.b.ncol[GMT_IN] < 3,
313 "Binary input data (-bi) must have at least 2 columns\n");
314 n_errors += gmt_M_check_condition (GMT, Ctrl->D.active && Ctrl->N.active, "-N and -D cannot both be given\n");
315 n_errors += gmt_M_check_condition (GMT, !Ctrl->E.active, "Option -E is required\n");
316
317 return (n_errors ? GMT_PARSE_ERROR : GMT_NOERROR);
318 }
319
grdrotater_get_grid_path(struct GMT_CTRL * GMT,struct GMT_GRID_HEADER * h)320 GMT_LOCAL struct GMT_DATASET * grdrotater_get_grid_path (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *h) {
321 /* Return a single polygon that encloses this geographic grid exactly.
322 * It is used in the case when no particular clip polygon has been given.
323 * Note that the path is the same for pixel or grid-registered grids.
324 */
325
326 unsigned int np = 0, add, col, row;
327 uint64_t dim[GMT_DIM_SIZE] = {1, 1, 0, 2};
328 struct GMT_DATASET *D = NULL;
329 struct GMT_DATASEGMENT *S = NULL;
330 struct GMT_DATASEGMENT_HIDDEN *SH = NULL;
331
332 if ((D = GMT_Create_Data (GMT->parent, GMT_IS_DATASET, GMT_IS_POLY, 0, dim, NULL, NULL, 0, 0, NULL)) == NULL)
333 return (NULL); /* An empty table with one segment, two cols */
334
335 S = D->table[0]->segment[0]; /* Short hand */
336
337 /* Add south border w->e */
338 if (h->wesn[YLO] == -90.0) { /* If at the S pole we just add it twice for end longitudes */
339 add = 2;
340 S->data[GMT_X] = gmt_M_memory (GMT, NULL, add, double);
341 S->data[GMT_Y] = gmt_M_memory (GMT, NULL, add, double);
342 S->data[GMT_X][0] = h->wesn[XLO]; S->data[GMT_X][1] = h->wesn[XHI];
343 S->data[GMT_Y][0] = S->data[GMT_Y][1] = h->wesn[YLO];
344 }
345 else { /* Loop along south border from west to east */
346 add = h->n_columns - !h->registration;
347 S->data[GMT_X] = gmt_M_memory (GMT, NULL, add, double);
348 S->data[GMT_Y] = gmt_M_memory (GMT, NULL, add, double);
349 for (col = 0; col < add; col++) {
350 S->data[GMT_X][col] = gmt_M_col_to_x (GMT, col, h->wesn[XLO], h->wesn[XHI], h->inc[GMT_X], 0.0, h->n_columns);
351 S->data[GMT_Y][col] = h->wesn[YLO];
352 }
353 }
354 np += add;
355 /* Add east border s->n */
356 add = h->n_rows - !h->registration;
357 S->data[GMT_X] = gmt_M_memory (GMT, S->data[GMT_X], add + np, double);
358 S->data[GMT_Y] = gmt_M_memory (GMT, S->data[GMT_Y], add + np, double);
359 for (row = 0; row < add; row++) { /* Loop along east border from south to north */
360 S->data[GMT_X][np+row] = h->wesn[XHI];
361 S->data[GMT_Y][np+row] = gmt_M_row_to_y (GMT, h->n_rows - 1 - row, h->wesn[YLO], h->wesn[YHI], h->inc[GMT_Y], 0.0, h->n_rows);
362 }
363 np += add;
364 /* Add north border e->w */
365 if (h->wesn[YHI] == 90.0) { /* If at the N pole we just add it twice for end longitudes */
366 add = 2;
367 S->data[GMT_X] = gmt_M_memory (GMT, S->data[GMT_X], add + np, double);
368 S->data[GMT_Y] = gmt_M_memory (GMT, S->data[GMT_Y], add + np, double);
369 S->data[GMT_X][np] = h->wesn[XHI]; S->data[GMT_X][np+1] = h->wesn[XLO];
370 S->data[GMT_Y][np] = S->data[GMT_Y][np+1] = h->wesn[YHI];
371 }
372 else { /* Loop along north border from east to west */
373 add = h->n_columns - !h->registration;
374 S->data[GMT_X] = gmt_M_memory (GMT, S->data[GMT_X], add + np, double);
375 S->data[GMT_Y] = gmt_M_memory (GMT, S->data[GMT_Y], add + np, double);
376 for (col = 0; col < add; col++) {
377 S->data[GMT_X][np+col] = gmt_M_col_to_x (GMT, h->n_columns - 1 - col, h->wesn[XLO], h->wesn[XHI], h->inc[GMT_X], 0.0, h->n_columns);
378 S->data[GMT_Y][np+col] = h->wesn[YHI];
379 }
380 }
381 np += add;
382 /* Add west border n->s */
383 add = h->n_rows - !h->registration;
384 S->data[GMT_X] = gmt_M_memory (GMT, S->data[GMT_X], add + np + 1, double);
385 S->data[GMT_Y] = gmt_M_memory (GMT, S->data[GMT_Y], add + np + 1, double);
386 for (row = 0; row < add; row++) { /* Loop along west border from north to south */
387 S->data[GMT_X][np+row] = h->wesn[XLO];
388 S->data[GMT_Y][np+row] = gmt_M_row_to_y (GMT, row, h->wesn[YLO], h->wesn[YHI], h->inc[GMT_Y], 0.0, h->n_rows);
389 }
390 np += add;
391 S->data[GMT_X][np] = S->data[GMT_X][0]; /* Close polygon explicitly */
392 S->data[GMT_Y][np] = S->data[GMT_Y][0];
393 np++;
394 S->n_rows = np;
395 S->n_columns = 2;
396 S->min[GMT_X] = h->wesn[XLO]; S->max[GMT_X] = h->wesn[XHI];
397 S->min[GMT_Y] = h->wesn[YLO]; S->max[GMT_Y] = h->wesn[YHI];
398 SH = gmt_get_DS_hidden (S);
399 SH->pole = 0;
400
401 return (D);
402 }
403
grdrotater_skip_if_outside(struct GMT_CTRL * GMT,struct GMT_DATATABLE * P,double lon,double lat)404 GMT_LOCAL bool grdrotater_skip_if_outside (struct GMT_CTRL *GMT, struct GMT_DATATABLE *P, double lon, double lat) {
405 /* Returns true if the selected point is outside the polygon */
406 uint64_t seg;
407 unsigned int inside = 0;
408 for (seg = 0; seg < P->n_segments && !inside; seg++) { /* Use degrees since function expects it */
409 if (gmt_polygon_is_hole (GMT, P->segment[seg])) continue; /* Holes are handled within gmt_inonout */
410 inside = (gmt_inonout (GMT, lon, lat, P->segment[seg]) > GMT_OUTSIDE);
411 }
412 return ((inside) ? false : true); /* true if outside */
413 }
414
415 #define bailout(code) {gmt_M_free_options (mode); return (code);}
416 #define Return(code) {Free_Ctrl (GMT, Ctrl); gmt_end_module (GMT, GMT_cpy); bailout (code);}
417
GMT_grdrotater(void * V_API,int mode,void * args)418 EXTERN_MSC int GMT_grdrotater (void *V_API, int mode, void *args) {
419 int scol, srow, error = 0; /* Signed row, col */
420 int n_stages;
421 bool not_global, global = false;
422 unsigned int col, row, col_o, row_o, start_row, stop_row, start_col, stop_col;
423 char gfile[PATH_MAX] = {""};
424
425 uint64_t ij, ij_rot, seg, rec, t;
426
427 double xx, yy, lon, P_original[3], P_rotated[3], R[3][3], plon, plat, pw;
428 double *grd_x = NULL, *grd_y = NULL, *grd_yc = NULL;
429
430 struct EULER *p = NULL; /* Pointer to array of stage poles */
431 struct GMT_DATASET *D = NULL, *Dr = NULL;
432 struct GMT_DATATABLE *pol = NULL, *polr = NULL;
433 struct GMT_DATASEGMENT *S = NULL, *Sr = NULL;
434 struct GMT_OPTION *ptr = NULL;
435 struct GMT_GRID *G = NULL, *G_rot = NULL;
436 struct GMT_GRID_HEADER_HIDDEN *HH = NULL;
437 struct GRDROTATER_CTRL *Ctrl = NULL;
438 struct GMT_CTRL *GMT = NULL, *GMT_cpy = NULL;
439 struct GMT_OPTION *options = NULL;
440 struct GMTAPI_CTRL *API = gmt_get_api_ptr (V_API); /* Cast from void to GMTAPI_CTRL pointer */
441
442 /*----------------------- Standard module initialization and parsing ----------------------*/
443
444 if (API == NULL) return (GMT_NOT_A_SESSION);
445 if (mode == GMT_MODULE_PURPOSE) return (usage (API, GMT_MODULE_PURPOSE)); /* Return the purpose of program */
446 options = GMT_Create_Options (API, mode, args); if (API->error) return (API->error); /* Set or get option list */
447
448 if ((error = gmt_report_usage (API, options, 0, usage)) != GMT_NOERROR) bailout (error); /* Give usage if requested */
449
450 /* Parse the command-line arguments */
451
452 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 */
453 if ((ptr = GMT_Find_Option (API, 'f', options)) == NULL) gmt_parse_common_options (GMT, "f", 'f', "g"); /* Did not set -f, implicitly set -fg */
454 if (GMT_Parse_Common (API, THIS_MODULE_OPTIONS, options)) Return (API->error);
455 Ctrl = New_Ctrl (GMT); /* Allocate and initialize a new control structure */
456 if ((error = parse (GMT, Ctrl, options)) != 0) Return (error);
457
458 /*---------------------------- This is the grdrotater main code ----------------------------*/
459
460 gmt_set_pad (GMT, 2U); /* Ensure space for BCs in case an API passed pad == 0 */
461
462 /* Check limits and get data file */
463
464 if (Ctrl->In.file) { /* Provided an input grid */
465 if ((G = 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 */
466 Return (API->error);
467 }
468 HH = gmt_get_H_hidden (G->header);
469 if (!GMT->common.R.active[RSET]) gmt_M_memcpy (GMT->common.R.wesn, G->header->wesn, 4, double); /* -R was not set so we use the grid domain */
470
471 /* Determine the wesn to be used to read the Ctrl->In.file; or exit if file is outside -R */
472
473 if (!gmt_grd_setregion (GMT, G->header, GMT->common.R.wesn, BCR_BILINEAR)) {
474 GMT_Report (API, GMT_MSG_ERROR, "No grid values inside selected region - aborting\n");
475 Return (GMT_RUNTIME_ERROR);
476 }
477 global = (doubleAlmostEqual (GMT->common.R.wesn[XHI] - GMT->common.R.wesn[XLO], 360.0)
478 && doubleAlmostEqual (GMT->common.R.wesn[YHI] - GMT->common.R.wesn[YLO], 180.0));
479 if (!GMT->common.R.active[RSET]) global = gmt_grd_is_global (GMT, G->header);
480 }
481 not_global = !global;
482
483 if (!Ctrl->S.active) { /* Read the input grid */
484 GMT_Report (API, GMT_MSG_INFORMATION, "Allocates memory and read grid file\n");
485 if (GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_DATA_ONLY, GMT->common.R.wesn, Ctrl->In.file, G) == NULL) {
486 Return (API->error);
487 }
488 }
489
490 if (Ctrl->F.active) { /* Read the user's polygon file */
491 if ((D = GMT_Read_Data (API, GMT_IS_DATASET, GMT_IS_FILE, GMT_IS_POLY, GMT_READ_NORMAL, NULL, Ctrl->F.file, NULL)) == NULL) {
492 Return (API->error);
493 }
494 pol = D->table[0]; /* Since we know it is a single file */
495 }
496 else if (not_global) { /* Make a single grid-outline polygon */
497 if (!G) {
498 GMT_Report (API, GMT_MSG_ERROR, "No grid give so cannot determine grid outline path\n");
499 Return (API->error);
500 }
501 if ((D = grdrotater_get_grid_path (GMT, G->header)) == NULL) Return (API->error);
502 pol = D->table[0]; /* Since we know it is a single file */
503 }
504 if (pol) {
505 Dr = GMT_Duplicate_Data (API, GMT_IS_DATASET, GMT_DUPLICATE_ALLOC, D); /* Same table length as D */
506 polr = Dr->table[0]; /* Since we know it is also a single file */
507 }
508 gmt_set_inside_mode (GMT, NULL, GMT_IOO_SPHERICAL);
509
510 if (Ctrl->E.rot.single) { /* Got a single rotation, create a rotation table with one entry */
511 Ctrl->T.n_times = n_stages = 1;
512 p = gmt_M_memory (GMT, NULL, n_stages, struct EULER);
513 p[0].lon = Ctrl->E.rot.lon; p[0].lat = Ctrl->E.rot.lat; p[0].omega = Ctrl->E.rot.w;
514 if (gmt_M_is_dnan (Ctrl->E.rot.age)) { /* No age, use fake age = 1 everywhere */
515 Ctrl->T.value = gmt_M_memory (GMT, NULL, Ctrl->T.n_times, double);
516 Ctrl->T.value[0] = p[0].t_start = 1.0;
517 }
518 spotter_setrot (GMT, &(p[0]));
519 }
520 else { /* Got a rotation file with multiple rotations in total reconstruction format */
521 double t_max = 0.0;
522 if ((n_stages = spotter_init (GMT, Ctrl->E.rot.file, &p, 0, true, Ctrl->E.rot.invert, &t_max)) < 0)
523 Return (-n_stages);
524 gmt_set_segmentheader (GMT, GMT_OUT, true);
525 }
526
527 if (!Ctrl->T.active && !Ctrl->E.rot.single) { /* Gave no time to go with the rotations, use rotation times */
528 GMT_Report (API, GMT_MSG_INFORMATION, "No reconstruction times specified; using %d reconstruction times from rotation table\n", n_stages);
529 Ctrl->T.n_times = n_stages;
530 Ctrl->T.value = gmt_M_memory (GMT, NULL, Ctrl->T.n_times, double);
531 for (t = 0; t < Ctrl->T.n_times; t++) Ctrl->T.value[t] = p[t].t_start;
532 }
533
534 if (Ctrl->T.n_times > 1) { /* Requires that template names be given */
535 if (!Ctrl->N.active) { /* Did not give -N so require -D template */
536 if (!Ctrl->D.file || !strchr (Ctrl->D.file, '%')) { /* No file given or filename without C-format specifiers */
537 GMT_Report (API, GMT_MSG_WARNING, "Multiple output times requires a template name via -D (unless -N is set)\n");
538 Return (API->error);
539 }
540 }
541 if (!Ctrl->S.active && !strchr (Ctrl->G.file, '%')) { /* Grid filename without C-format specifiers */
542 GMT_Report (API, GMT_MSG_WARNING, "Multiple output times requires a template gridfile name via -G\n");
543 Return (API->error);
544 }
545 }
546
547 for (t = 0; t < Ctrl->T.n_times; t++) { /* For each reconstruction time */
548 if (Ctrl->E.rot.single) {
549 plon = Ctrl->E.rot.lon; plat = Ctrl->E.rot.lat; pw = Ctrl->E.rot.w;
550 GMT_Report (API, GMT_MSG_INFORMATION, "Using rotation (%g, %g, %g)\n", plon, plat, pw);
551 }
552 else { /* Extract rotation for given time */
553 if (Ctrl->T.value[t] < p[0].t_stop || Ctrl->T.value[t] > p[n_stages-1].t_start) {
554 GMT_Report (API, GMT_MSG_ERROR, "Requested a reconstruction time outside range of rotation table - skipped\n");
555 continue;
556 }
557 spotter_get_rotation (GMT, p, n_stages, Ctrl->T.value[t], &plon, &plat, &pw);
558 GMT_Report (API, GMT_MSG_INFORMATION, "Time %g Ma: Using rotation (%g, %g, %g)\n", Ctrl->T.value[t], plon, plat, pw);
559 }
560 gmt_make_rot_matrix (GMT, plon, plat, pw, R); /* Make rotation matrix from rotation parameters */
561
562 if (Ctrl->E.rot.single)
563 GMT_Report (API, GMT_MSG_INFORMATION, "Reconstruct polygon outline\n");
564 else
565 GMT_Report (API, GMT_MSG_INFORMATION, "Reconstruct polygon outline for time %g\n", Ctrl->T.value[t]);
566
567 /* First reconstruct the polygon outline */
568
569 for (seg = 0; pol && seg < pol->n_segments; seg++) {
570 S = pol->segment[seg]; /* Shorthand for current original segment */
571 Sr = polr->segment[seg]; /* Shorthand for current rotated segment */
572 for (rec = 0; rec < pol->segment[seg]->n_rows; rec++) {
573 Sr->data[GMT_X][rec] = S->data[GMT_X][rec];
574 Sr->data[GMT_Y][rec] = gmt_lat_swap (GMT, S->data[GMT_Y][rec], GMT_LATSWAP_G2O); /* Convert to geocentric */
575 gmt_geo_to_cart (GMT, Sr->data[GMT_Y][rec], Sr->data[GMT_X][rec], P_original, true); /* Convert to a Cartesian x,y,z vector; true since we have degrees */
576 gmt_matrix_vect_mult (GMT, 3U, R, P_original, P_rotated); /* Rotate the vector */
577 gmt_cart_to_geo (GMT, &Sr->data[GMT_Y][rec], &Sr->data[GMT_X][rec], P_rotated, true); /* Recover lon lat representation; true to get degrees */
578 Sr->data[GMT_Y][rec] = gmt_lat_swap (GMT, Sr->data[GMT_Y][rec], GMT_LATSWAP_O2G); /* Convert back to geodetic */
579 }
580 gmt_set_seg_polar (GMT, Sr); /* Determine if it is a polar cap */
581 }
582 gmt_set_tbl_minmax (GMT, GMT_IS_POLY, polr); /* Update table domain */
583 if (!Ctrl->N.active && not_global) {
584 char dfile[PATH_MAX] = {""}, *file = NULL;
585 if (Ctrl->T.n_times > 1) {
586 sprintf (dfile, Ctrl->D.file, Ctrl->T.value[t]);
587 file = dfile;
588 }
589 else
590 file = Ctrl->D.file;
591 if (!Ctrl->E.rot.single) { /* Add a segment header with the age via -Z */
592 char txt[PATH_MAX] = {""};
593 sprintf (txt, "-Z%g", Ctrl->T.value[t]);
594 for (seg = 0; seg < polr->n_segments; seg++) {
595 Sr = polr->segment[seg]; /* Shorthand for current rotated segment */
596 gmt_M_str_free (Sr->header);
597 Sr->header = strdup (txt);
598 }
599 }
600 if (GMT_Write_Data (API, GMT_IS_DATASET, GMT_IS_FILE, GMT_IS_POLY, GMT_WRITE_SET, NULL, file, Dr) != GMT_NOERROR) {
601 Return (API->error);
602 }
603 }
604 if (Ctrl->S.active) /* No grids will be rotated */
605 continue;
606
607 /* Then, find min/max of reconstructed outline */
608
609 if (global)
610 gmt_M_memcpy (GMT->common.R.wesn, G->header->wesn, 4, double);
611 else if (Ctrl->A.active) {
612 gmt_M_memcpy (GMT->common.R.wesn, Ctrl->A.wesn, 4, double);
613 /* Adjust longitude range, as indicated by FORMAT_GEO_OUT */
614 gmt_lon_range_adjust (GMT->current.io.geo.range, &GMT->common.R.wesn[XLO]);
615 gmt_lon_range_adjust (GMT->current.io.geo.range, &GMT->common.R.wesn[XHI]);
616 }
617 else {
618 GMT->common.R.wesn[XLO] = floor (polr->min[GMT_X] * HH->r_inc[GMT_X]) * G->header->inc[GMT_X];
619 GMT->common.R.wesn[XHI] = ceil (polr->max[GMT_X] * HH->r_inc[GMT_X]) * G->header->inc[GMT_X];
620 GMT->common.R.wesn[YLO] = floor (polr->min[GMT_Y] * HH->r_inc[GMT_Y]) * G->header->inc[GMT_Y];
621 GMT->common.R.wesn[YHI] = ceil (polr->max[GMT_Y] * HH->r_inc[GMT_Y]) * G->header->inc[GMT_Y];
622 /* Adjust longitude range, as indicated by FORMAT_GEO_OUT */
623 gmt_lon_range_adjust (GMT->current.io.geo.range, &GMT->common.R.wesn[XLO]);
624 gmt_lon_range_adjust (GMT->current.io.geo.range, &GMT->common.R.wesn[XHI]);
625 if (GMT->common.R.wesn[XLO] >= GMT->common.R.wesn[XHI]) GMT->common.R.wesn[XHI] += 360.0;
626 }
627 GMT->common.R.active[RSET] = true;
628
629 if ((G_rot = GMT_Create_Data (API, GMT_IS_GRID, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA, NULL, NULL, G->header->inc, \
630 GMT_GRID_DEFAULT_REG, GMT_NOTSET, NULL)) == NULL) Return (API->error);
631
632 /* Precalculate node coordinates in both degrees and radians */
633 grd_x = G_rot->x;
634 grd_y = G_rot->y;
635 grd_yc = gmt_M_memory (GMT, NULL, G_rot->header->n_rows, double);
636 for (row = 0; row < G_rot->header->n_rows; row++) grd_yc[row] = gmt_lat_swap (GMT, grd_y[row], GMT_LATSWAP_G2O);
637
638 /* Loop over all nodes in the new rotated grid and find those inside the reconstructed polygon */
639
640 if (Ctrl->E.rot.single)
641 GMT_Report (API, GMT_MSG_INFORMATION, "Interpolate reconstructed grid\n");
642 else
643 GMT_Report (API, GMT_MSG_INFORMATION, "Interpolate reconstructed grid for time %g\n", Ctrl->T.value[t]);
644
645 gmt_make_rot_matrix (GMT, plon, plat, -pw, R); /* Make inverse rotation using negative angle */
646
647 gmt_M_grd_loop (GMT, G_rot, row, col, ij_rot) {
648 G_rot->data[ij_rot] = GMT->session.f_NaN;
649 if (not_global && grdrotater_skip_if_outside (GMT, polr, grd_x[col], grd_y[row])) continue; /* Outside rotated polygon */
650
651 /* Here we are inside; get the coordinates and rotate back to original grid coordinates */
652
653 gmt_geo_to_cart (GMT, grd_yc[row], grd_x[col], P_rotated, true); /* Convert degree lon,lat to a Cartesian x,y,z vector */
654 gmt_matrix_vect_mult (GMT, 3U, R, P_rotated, P_original); /* Rotate the vector */
655 gmt_cart_to_geo (GMT, &yy, &xx, P_original, true); /* Recover degree lon lat representation */
656 yy = gmt_lat_swap (GMT, yy, GMT_LATSWAP_O2G); /* Convert back to geodetic */
657 xx -= 360.0;
658 while (xx < G->header->wesn[XLO]) xx += 360.0; /* Make sure we deal with 360 issues */
659 G_rot->data[ij_rot] = (gmt_grdfloat)gmt_bcr_get_z (GMT, G, xx, yy);
660 }
661
662 /* Also loop over original node locations to make sure the nearest nodes are set */
663
664 for (seg = 0; not_global && seg < pol->n_segments; seg++) {
665 for (rec = 0; rec < pol->segment[seg]->n_rows; rec++) {
666 lon = pol->segment[seg]->data[GMT_X][rec];
667 while (lon < G_rot->header->wesn[XLO]) lon += 360.0;
668 scol = (int)gmt_M_grd_x_to_col (GMT, lon, G_rot->header);
669 srow = (int)gmt_M_grd_y_to_row (GMT, pol->segment[seg]->data[GMT_Y][rec], G_rot->header);
670 /* Visit the PAD * PAD number of cells centered on col, row and make sure they have been set */
671 start_row = (srow > PAD) ? srow - PAD : 0;
672 stop_row = ((srow + PAD) >= 0) ? srow + PAD : 0;
673 start_col = (scol > PAD) ? scol - PAD : 0;
674 stop_col = ((scol + PAD) >= 0) ? scol + PAD : 0;
675 for (row = start_row; row <= stop_row; row++) {
676 if (row >= G_rot->header->n_rows) continue;
677 for (col = start_col; col <= stop_col; col++) {
678 if (col >= G_rot->header->n_columns) continue;
679 ij_rot = gmt_M_ijp (G_rot->header, row, col);
680 if (!gmt_M_is_fnan (G_rot->data[ij_rot])) continue; /* Already done this */
681 if (not_global && grdrotater_skip_if_outside (GMT, pol, grd_x[col], grd_yc[row])) continue; /* Outside input polygon */
682 gmt_geo_to_cart (GMT, grd_yc[row], grd_x[col], P_rotated, true); /* Convert degree lon,lat to a Cartesian x,y,z vector */
683 gmt_matrix_vect_mult (GMT, 3U, R, P_rotated, P_original); /* Rotate the vector */
684 gmt_cart_to_geo (GMT, &xx, &yy, P_original, true); /* Recover degree lon lat representation */
685 yy = gmt_lat_swap (GMT, yy, GMT_LATSWAP_O2G); /* Convert back to geodetic */
686 scol = (int)gmt_M_grd_x_to_col (GMT, xx, G->header);
687 if (scol < 0) continue;
688 col_o = scol; if (col_o >= G->header->n_columns) continue;
689 srow = (int)gmt_M_grd_y_to_row (GMT, yy, G->header);
690 if (srow < 0) continue;
691 row_o = srow; if (row_o >= G->header->n_rows) continue;
692 ij = gmt_M_ijp (G->header, row_o, col_o);
693 G_rot->data[ij_rot] = G->data[ij];
694 }
695 }
696 }
697 }
698
699 /* Now write rotated grid */
700
701 if (Ctrl->E.rot.single)
702 GMT_Report (API, GMT_MSG_INFORMATION, "Write reconstructed grid\n");
703 else
704 GMT_Report (API, GMT_MSG_INFORMATION, "Write reconstructed grid for time %g\n", Ctrl->T.value[t]);
705
706 gmt_set_pad (GMT, API->pad); /* Reset to session default pad before output */
707
708 if (Ctrl->E.rot.single)
709 snprintf (G_rot->header->remark, GMT_GRID_REMARK_LEN160, "Grid reconstructed using R[lon lat omega] = %g %g %g", plon, plat, pw);
710 else
711 snprintf (G_rot->header->remark, GMT_GRID_REMARK_LEN160, "Grid reconstructed using R[lon lat omega] = %g %g %g for time %g", plon, plat, pw, Ctrl->T.value[t]);
712 if (GMT_Set_Comment (API, GMT_IS_GRID, GMT_COMMENT_IS_OPTION | GMT_COMMENT_IS_COMMAND, options, G_rot)) Return (API->error);
713 if (Ctrl->T.n_times > 1) /* Use template to create name */
714 sprintf (gfile, Ctrl->G.file, Ctrl->T.value[t]);
715 else
716 strncpy (gfile, Ctrl->G.file, PATH_MAX-1);
717 if (GMT_Write_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA, NULL, gfile, G_rot) != GMT_NOERROR) {
718 Return (API->error);
719 }
720 if (G_rot && GMT_Destroy_Data (API, &G_rot) != GMT_NOERROR)
721 Return (API->error);
722
723 gmt_M_free (GMT, grd_yc);
724 } /* End of loop over reconstruction times */
725
726 gmt_M_free (GMT, p);
727
728 if (Ctrl->S.active) {
729 if (Ctrl->F.active && GMT_Destroy_Data (API, &D) != GMT_NOERROR) {
730 Return (API->error);
731 }
732 else if (not_global)
733 gmt_free_dataset (GMT, &D);
734 }
735 gmt_M_free (GMT, Ctrl->T.value);
736 if (D && GMT_Destroy_Data (API, &D) != GMT_NOERROR)
737 Return (API->error);
738 if (Dr && GMT_Destroy_Data (API, &Dr) != GMT_NOERROR)
739 Return (API->error);
740
741 Return (GMT_NOERROR);
742 }
743