1 /*--------------------------------------------------------------------
2 *
3 * Copyright (c) 1991-2021 by the GMT Team (https://www.generic-mapping-tools.org/team.html)
4 * See LICENSE.TXT file for copying and redistribution conditions.
5 *
6 * This program is free software; you can redistribute it and/or modify
7 * it under the terms of the GNU Lesser General Public License as published by
8 * the Free Software Foundation; version 3 or any later version.
9 *
10 * This program is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 * GNU Lesser General Public License for more details.
14 *
15 * Contact info: www.generic-mapping-tools.org
16 *--------------------------------------------------------------------*/
17 /*
18 * Author: Paul Wessel
19 * Date: 1-NOV-2021
20 * Version: 6 API
21 *
22 * Brief synopsis: grdselect reads one or more grid/cube files and determines the
23 * union or intersection of the regions, modulated by other options. Alternatively
24 * it just writes the names of the data sources that pass the specified tests.
25 *
26 */
27
28 #include "gmt_dev.h"
29
30 #define THIS_MODULE_CLASSIC_NAME "grdselect"
31 #define THIS_MODULE_MODERN_NAME "grdselect"
32 #define THIS_MODULE_LIB "core"
33 #define THIS_MODULE_PURPOSE "Make selections or determine common regions from 2-D grids, images or 3-D cubes"
34 #define THIS_MODULE_KEYS "<?{+,>D}"
35 #define THIS_MODULE_NEEDS ""
36 #define THIS_MODULE_OPTIONS "->RVfhor"
37
38 /* Control structure for grdselect */
39
40 enum Opt_A_modes {
41 GRDSELECT_INTERSECTION = 0, /* -Ai */
42 GRDSELECT_UNION = 1, /* -Au */
43 GRDSELECT_NO_INC = 0,
44 GRDSELECT_MIN_INC = 1, /* +il */
45 GRDSELECT_MAX_INC = 2, /* +ih */
46 GRDSELECT_SET_INC = 3}; /* +i<incs> */
47
48 enum Opt_N_modes {
49 GRDSELECT_LESS_NANS = 1, /* -Nl */
50 GRDSELECT_MORE_NANS = 2}; /* -Nh */
51
52 enum GRDSELECT { /* Indices for the various tests */
53 GRD_SELECT_C = 0,
54 GRD_SELECT_D,
55 GRD_SELECT_F,
56 GRD_SELECT_L,
57 GRD_SELECT_N,
58 GRD_SELECT_R,
59 GRD_SELECT_W,
60 GRD_SELECT_Z,
61 GRD_SELECT_r,
62 GRD_SELECT_N_TESTS /* Number of specific tests available */
63 };
64
65 #define GRDSELECT_STRING "CDFLNRWZr" /* These are all the possible items to invert */
66
67 #define WLO 6/* Internal array indices */
68 #define WHI 7
69
70 struct GRDSELECT_CTRL {
71 unsigned int n_files; /* How many data sources given */
72 struct GRDSELECT_A { /* -A[i|u][+il|h|arg>]*/
73 bool active;
74 bool round; /* True if +i is used to set rounding */
75 unsigned int i_mode; /* 0 = no increment rounding, 1 is use smallest inc, 2 = use largest inc, 3 = use appended inc */
76 unsigned int mode; /* 0 = intersection [i], 1 = union [r] */
77 double inc[2]; /* Optional increments */
78 } A;
79 struct GRDSELECT_C { /* -C<pointfile> */
80 bool active;
81 char *file;
82 } C;
83 struct GRDSELECT_D { /* -D<dx[/dy[/dz]] */
84 bool active;
85 double inc[3];
86 } D;
87 struct GRDSELECT_E { /* -E[b] */
88 bool active;
89 bool box;
90 } E;
91 struct GRDSELECT_F { /* -F<polyfile>[+i|o] */
92 bool active;
93 unsigned int mode; /* 0 = outside, 1 = crossing, 2 = inside */
94 char *file;
95 } F;
96 struct GRDSELECT_G { /* -G */
97 bool active;
98 } G;
99 struct GRDSELECT_I { /* -ICDFLNRWZr */
100 bool active;
101 bool pass[GRD_SELECT_N_TESTS]; /* One flag for each setting */
102 } I;
103 struct GRDSELECT_L { /* -F<linefile> */
104 bool active;
105 char *file;
106 } L;
107 struct GRDSELECT_M { /* -M */
108 bool active;
109 double margin[4]; /* Optional margins [none] */
110 } M;
111 struct GRDSELECT_N { /* -Nl|h<nans> */
112 bool active;
113 unsigned int mode;
114 int64_t n_nans; /* Limit on nans [0] */
115 } N;
116 struct GRDSELECT_Q { /* -Q */
117 bool active;
118 } Q;
119 struct GRDSELECT_W { /* -Wwmin/wmax */
120 bool active;
121 double w_min, w_max;
122 } W;
123 struct GRDSELECT_Z { /* -Zzmin/zmax */
124 bool active;
125 double z_min, z_max;
126 } Z;
127 };
128
New_Ctrl(struct GMT_CTRL * GMT)129 static void *New_Ctrl (struct GMT_CTRL *GMT) { /* Allocate and initialize a new control structure */
130 struct GRDSELECT_CTRL *C;
131
132 C = gmt_M_memory (GMT, NULL, 1, struct GRDSELECT_CTRL);
133
134 /* Initialize values whose defaults are not 0/false/NULL */
135 C->F.mode = GMT_ONEDGE; /* Partial overlap passes the test */
136 for (unsigned int k = 0; k < GRD_SELECT_N_TESTS; k++) C->I.pass[k] = true; /* Default is to include the grid if we pass the test */
137 return (C);
138 }
139
Free_Ctrl(struct GMT_CTRL * GMT,struct GRDSELECT_CTRL * C)140 static void Free_Ctrl (struct GMT_CTRL *GMT, struct GRDSELECT_CTRL *C) { /* Deallocate control structure */
141 if (!C) return;
142 gmt_M_str_free (C->C.file);
143 gmt_M_str_free (C->F.file);
144 gmt_M_str_free (C->L.file);
145 gmt_M_free (GMT, C);
146 }
147
usage(struct GMTAPI_CTRL * API,int level)148 static int usage (struct GMTAPI_CTRL *API, int level) {
149 const char *name = gmt_show_name_and_purpose (API, THIS_MODULE_LIB, THIS_MODULE_CLASSIC_NAME, THIS_MODULE_PURPOSE);
150 if (level == GMT_MODULE_PURPOSE) return (GMT_NOERROR);
151 GMT_Usage (API, 0, "usage: %s source1 source2 ... [-Ai|u[+il|h|<inc>]] [-C<pointtable>] [-D<inc>] [-E[b]] [-F<polygontable>[+i|o]] [-G] [-I%s] "
152 " [-L<linetable>] [-M<margins>] [-Nl|h[<n>]] [-Q] [%s] [-W[<min>]/[<max>]] [-Z[<min>]/[<max>]] [%s] [%s] [%s] [%s] [%s]\n",
153 name, GRDSELECT_STRING, GMT_Rgeo_OPT, GMT_V_OPT, GMT_f_OPT, GMT_ho_OPT, GMT_o_OPT, GMT_PAR_OPT);
154
155 if (level == GMT_SYNOPSIS) return (GMT_MODULE_SYNOPSIS);
156
157 GMT_Message (API, GMT_TIME_NONE, " REQUIRED ARGUMENTS:\n");
158 gmt_ingrid_syntax (API, 0, "Name of one or more grid, image or cube files");
159 GMT_Usage (API, -2, "Note: Options -N and -W are ignored for images.");
160 GMT_Message (API, GMT_TIME_NONE, "\n OPTIONAL ARGUMENTS:\n");
161 GMT_Usage (API, 1, "\n-Ai|u[+il|h|<inc>]");
162 GMT_Usage (API, -2, "Report a unifying data region for all the data sources. Append:");
163 GMT_Usage (API, 3, "i: Determine the intersection of all passed data source regions.");
164 GMT_Usage (API, 3, "u: Determine the union of all passed data source regions.");
165 GMT_Usage (API, -2, "Optionally, append +i to round region using the (l)owest, (h)ighest, or specified increments [no rounding]."
166 "Note: Without -A we simply report the list of data sources passing any imposed tests.");
167 GMT_Usage (API, 1, "\n-C<pointtable>");
168 GMT_Usage (API, -2, "Only consider data sources that contain at least one of these points [Consider all input data sources].");
169 GMT_Usage (API, 1, "\n-D<inc>");
170 GMT_Usage (API, -2, "Only consider data sources that match the given increments [Consider all input data sources].");
171 GMT_Usage (API, 1, "\n-E Report information in form of a single data record using the format "
172 "<w e s n {b t} w0 w1> [Default reports a -R<w/e/s/n>[/b/t] string]. If -Ai is in effect then we recompute <w0 w1> for that common region. "
173 "Alternatively, append b to output the determined regions bounding box polygon instead.");
174 GMT_Usage (API, 1, "\n-F<polygontable>");
175 GMT_Usage (API, -2, "Only consider data sources that partially or fully overlap with these polygons in map view [Consider all input data sources].");
176 GMT_Usage (API, 1, "\n-G Force possible download of all tiles for a remote <grid> if given as input [no report for tiled grids].");
177 GMT_Usage (API, 3, "+i: Only pass data sources that are fully inside the polygon.");
178 GMT_Usage (API, 3, "+o: Only pass data sources that are fully outside the polygon.");
179 GMT_Usage (API, 1, "\n-I%s", GRDSELECT_STRING);
180 GMT_Usage (API, -2, "Reverse the tests, i.e., pass data sources when the test fails. "
181 "Supply any combination of %s where each flag means:", GRDSELECT_STRING);
182 GMT_Usage (API, 3, "C: Pass data sources that do not contain any of the points set in -C.");
183 GMT_Usage (API, 3, "D: Pass data sources that do not match the increment set in -D.");
184 GMT_Usage (API, 3, "F: Pass data sources that do not overlap with polygons set in -F.");
185 GMT_Usage (API, 3, "L: Pass data sources that are not traversed by lines set in -L.");
186 GMT_Usage (API, 3, "N: Pass data sources that fail the NaN-test in -N.");
187 GMT_Usage (API, 3, "R: Pass data sources that do not overlap with region in -R.");
188 GMT_Usage (API, 3, "W: Pass data sources outside the data range given in -W.");
189 GMT_Usage (API, 3, "Z: Pass cubes outside the z-coordinate range given in -Z (requires -Q).");
190 GMT_Usage (API, 3, "r: Pass data sources with the opposite registration than given in -r.");
191 GMT_Usage (API, -2, "Note: If no argument is given then we default to -I%s.", GRDSELECT_STRING);
192 GMT_Usage (API, 1, "\n-L<linetable>");
193 GMT_Usage (API, -2, "Only consider data sources that are traversed by these lines in map view [Consider all input data sources].");
194 GMT_Usage (API, 1, "\n-M<margins>");
195 GMT_Usage (API, -2, "Add padding around the final (rounded) region. Append a uniform <margin>, separate <xmargin>/<ymargin>, "
196 "or individual <wmargin>/<emargin>/<smargin>/<nmargin> for each side [no padding].");
197 GMT_Usage (API, 1, "\n-Nl|h[<n>]");
198 GMT_Usage (API, -2, "Only consider data sources that satisfy a NaN-condition [Consider all input data sources]:");
199 GMT_Usage (API, 3, "l: Only data sources with lower than <n> NaNs will pass [0].");
200 GMT_Usage (API, 3, "h: Only data sources with higher than <n> NaNs will pass [0].");
201 GMT_Usage (API, 1, "\n-Q Input file(s) is 3-D data cube(s), not grid(s) [2-D grids].");
202 GMT_Option (API, "R");
203 GMT_Usage (API, 1, "\n-W[<min>]/[<max>]");
204 GMT_Usage (API, -2, "Only consider data sources that have data-values in the given range [Consider all input data sources]. "
205 "At least one of <min> or <max> must be specified, as well as the slash [-infinity/+infinity].");
206 GMT_Usage (API, 1, "\n-Z[<min>]/[<max>]");
207 GMT_Usage (API, -2, "Only consider cubes that have z-coordinates in the given range [Consider all input cubes]. "
208 "At least one of <min> or <max> must be specified, as well as the slash [-infinity/+infinity]. Requires -Q.");
209 GMT_Usage (API, 1, "\n-r[g|p]");
210 GMT_Usage (API, -2, "Only consider data sources that have the specified registration [Consider all input data sources].");
211 GMT_Option (API, "V,f,h,o,r,.");
212
213 return (GMT_MODULE_USAGE);
214 }
215
parse(struct GMT_CTRL * GMT,struct GRDSELECT_CTRL * Ctrl,struct GMT_OPTION * options)216 static int parse (struct GMT_CTRL *GMT, struct GRDSELECT_CTRL *Ctrl, struct GMT_OPTION *options) {
217
218 /* This parses the options provided to grdselect and sets parameters in CTRL.
219 * Any GMT common options will override values set previously by other commands.
220 * It also replaces any file names specified as input or output with the data ID
221 * returned when registering these sources/destinations with the API.
222 */
223
224 unsigned int n_errors = 0, k = 0;
225 char string[GMT_LEN128] = {""}, *c = NULL;
226 struct GMT_OPTION *opt = NULL;
227 struct GMTAPI_CTRL *API = GMT->parent;
228
229 for (opt = options; opt; opt = opt->next) { /* Process all the options given */
230
231 switch (opt->option) {
232 /* Common parameters */
233
234 case '<': /* Input files */
235 if (GMT_Get_FilePath (API, GMT_IS_GRID, GMT_IN, GMT_FILE_REMOTE, &(opt->arg)))
236 n_errors++;
237 else
238 Ctrl->n_files++;
239 break;
240
241 /* Processes program-specific parameters */
242
243 case 'A': /* Area comparisons */
244 n_errors += gmt_M_repeated_module_option (API, Ctrl->A.active);
245 Ctrl->A.active = true;
246 switch (opt->arg[0]) {
247 case 'u': Ctrl->A.mode = GRDSELECT_UNION; k = 1; break;
248 case 'i': Ctrl->A.mode = GRDSELECT_INTERSECTION; k = 1; break;
249 default: k = 0; break; /* No directive, default to intersection */
250 }
251 if (gmt_get_modifier (opt->arg, 'i', string)) { /* Want to control rounding */
252 Ctrl->A.round = true;
253 switch (string[0]) {
254 case 'l': Ctrl->A.i_mode = GRDSELECT_MIN_INC; break;
255 case 'h': Ctrl->A.i_mode = GRDSELECT_MAX_INC; break;
256 default:
257 if (gmt_getinc (GMT, string, Ctrl->A.inc)) n_errors++;
258 Ctrl->A.i_mode = GRDSELECT_SET_INC;
259 break;
260 }
261 }
262 break;
263 case 'C': /* Point file */
264 n_errors += gmt_M_repeated_module_option (API, Ctrl->C.active);
265 Ctrl->C.active = true;
266 if (opt->arg[0]) Ctrl->C.file = strdup (opt->arg);
267 if (GMT_Get_FilePath (API, GMT_IS_DATASET, GMT_IN, GMT_FILE_REMOTE, &(Ctrl->C.file))) n_errors++;
268 break;
269 case 'D': /* Specified grid increments */
270 n_errors += gmt_M_repeated_module_option (API, Ctrl->D.active);
271 Ctrl->D.active = true;
272 if (opt->arg[0] && gmt_getinc (GMT, opt->arg, Ctrl->D.inc)) {
273 gmt_inc_syntax (GMT, 'D', 1);
274 n_errors++;
275 }
276 break;
277 case 'E': /* Column format */
278 n_errors += gmt_M_repeated_module_option (API, Ctrl->E.active);
279 Ctrl->E.active = true;
280 if (opt->arg[0] == 'b') Ctrl->E.box = true;
281 break;
282 case 'F': /* Polygon file */
283 n_errors += gmt_M_repeated_module_option (API, Ctrl->F.active);
284 Ctrl->F.active = true;
285 if ((c = strstr (opt->arg, "+i")))
286 Ctrl->F.mode = GMT_INSIDE;
287 else if ((c = strstr (opt->arg, "+o")))
288 Ctrl->F.mode = GMT_OUTSIDE;
289 if (c) c[0] = '\0'; /* Temporarily chop off modifier */
290 if (opt->arg[0]) Ctrl->F.file = strdup (opt->arg);
291 if (GMT_Get_FilePath (API, GMT_IS_DATASET, GMT_IN, GMT_FILE_REMOTE, &(Ctrl->F.file))) n_errors++;
292 if (c) c[0] = '\0'; /* Restore modifier */
293 break;
294 case 'I': /* Invert selected tests */
295 n_errors += gmt_M_repeated_module_option (API, Ctrl->I.active);
296 Ctrl->I.active = true;
297 if (opt->arg[0] == '\0') { /* If -I is given then all tests are reversed */
298 for (k = 0; k < GRD_SELECT_N_TESTS; k++) Ctrl->I.pass[k] = false;
299 }
300 else { /* Reverse those requested only */
301 for (k = 0; opt->arg[k]; k++) {
302 switch (opt->arg[k]) {
303 case 'C': Ctrl->I.pass[GRD_SELECT_C] = false; break;
304 case 'D': Ctrl->I.pass[GRD_SELECT_D] = false; break;
305 case 'F': Ctrl->I.pass[GRD_SELECT_F] = false; break;
306 case 'L': Ctrl->I.pass[GRD_SELECT_L] = false; break;
307 case 'N': Ctrl->I.pass[GRD_SELECT_N] = false; break;
308 case 'R': Ctrl->I.pass[GRD_SELECT_R] = false; break;
309 case 'W': Ctrl->I.pass[GRD_SELECT_W] = false; break;
310 case 'Z': Ctrl->I.pass[GRD_SELECT_Z] = false; break;
311 case 'r': Ctrl->I.pass[GRD_SELECT_r] = false; break;
312 default:
313 GMT_Report (API, GMT_MSG_ERROR, "Option -I: Expects any combination of %s (%c is not valid)\n", GRDSELECT_STRING, opt->arg[k]);
314 n_errors++;
315 break;
316 }
317 }
318 }
319 break;
320 case 'L': /* Line file */
321 n_errors += gmt_M_repeated_module_option (API, Ctrl->L.active);
322 Ctrl->L.active = true;
323 if (opt->arg[0]) Ctrl->L.file = strdup (opt->arg);
324 if (GMT_Get_FilePath (API, GMT_IS_DATASET, GMT_IN, GMT_FILE_REMOTE, &(Ctrl->L.file))) n_errors++;
325 break;
326 case 'M': /* Extend the region */
327 n_errors += gmt_M_repeated_module_option (API, Ctrl->M.active);
328 Ctrl->M.active = true;
329 k = GMT_Get_Values (GMT->parent, opt->arg, Ctrl->M.margin, 4);
330 if (k == 1) /* Same increments in all directions */
331 Ctrl->M.margin[XHI] = Ctrl->M.margin[YLO] = Ctrl->M.margin[YHI] = Ctrl->M.margin[XLO];
332 else if (k == 2) { /* Separate increments in x and y */
333 Ctrl->M.margin[YLO] = Ctrl->M.margin[YHI] = Ctrl->M.margin[XHI];
334 Ctrl->M.margin[XHI] = Ctrl->M.margin[XLO];
335 }
336 else if (k != 4) { /* The only other option is 4 but somehow we failed */
337 GMT_Report (GMT->parent, GMT_MSG_ERROR, "Option -M: Bad number of increments given (%s)\n", opt->arg);
338 n_errors++;
339 }
340 break;
341 case 'N': /* NaN condition */
342 n_errors += gmt_M_repeated_module_option (API, Ctrl->N.active);
343 Ctrl->N.active = true;
344 switch (opt->arg[0]) {
345 case 'l': Ctrl->N.mode = GRDSELECT_LESS_NANS; break;
346 case 'h': Ctrl->N.mode = GRDSELECT_MORE_NANS; break;
347 default:
348 GMT_Report (GMT->parent, GMT_MSG_ERROR, "Option -M: Bad directive, must be -Ml or -Mh\n");
349 n_errors++;
350 }
351 if (opt->arg[1]) Ctrl->N.n_nans = atoi (&opt->arg[1]);
352 break;
353 case 'Q': /* Expect cubes */
354 n_errors += gmt_M_repeated_module_option (API, Ctrl->Q.active);
355 Ctrl->Q.active = true;
356 break;
357 case 'Z': /* z range */
358 n_errors += gmt_M_repeated_module_option (API, Ctrl->Z.active);
359 Ctrl->Z.active = true;
360 n_errors += gmt_get_limits (GMT, 'Z', opt->arg, 1, &Ctrl->Z.z_min, &Ctrl->Z.z_max);
361 break;
362
363 case 'W': /* w range */
364 n_errors += gmt_M_repeated_module_option (API, Ctrl->W.active);
365 Ctrl->W.active = true;
366 n_errors += gmt_get_limits (GMT, 'W', opt->arg, 1, &Ctrl->W.w_min, &Ctrl->W.w_max);
367 break;
368
369 default: /* Report bad options */
370 n_errors += gmt_default_error (GMT, opt->option);
371 break;
372 }
373 }
374
375 n_errors += gmt_M_check_condition (GMT, Ctrl->n_files == 0,
376 "Must specify one or more input files\n");
377 n_errors += gmt_M_check_condition (GMT, Ctrl->A.i_mode == GRDSELECT_SET_INC && (Ctrl->A.inc[GMT_X] <= 0.0 || Ctrl->A.inc[GMT_Y] <= 0.0),
378 "Option -A: Must specify a positive increment(s) via +i<inc>\n");
379 n_errors += gmt_M_check_condition (GMT, Ctrl->D.active && (Ctrl->D.inc[GMT_X] <= 0.0 || Ctrl->D.inc[GMT_Y] <= 0.0),
380 "Option -D: Must specify a positive increment(s)\n");
381 n_errors += gmt_M_check_condition (GMT, GMT->common.o.active && Ctrl->E.active,
382 "The -o option requires -C\n");
383 n_errors += gmt_M_check_condition (GMT, Ctrl->N.active && Ctrl->N.n_nans < 0,
384 "The -N option argument cannot be negative\n");
385 n_errors += gmt_M_check_condition (GMT, Ctrl->Z.active && !Ctrl->Q.active,
386 "The -Z option requires -Q since the limits applies to the 3-D z-dimension (see -W for data range limits)\n");
387
388 return (n_errors ? GMT_PARSE_ERROR : GMT_NOERROR);
389 }
390
grdselect_line_overlaps(struct GMT_DATASEGMENT * S,uint64_t k1,struct GMT_GRID_HEADER * header,double west,double east)391 GMT_LOCAL bool grdselect_line_overlaps (struct GMT_DATASEGMENT *S, uint64_t k1, struct GMT_GRID_HEADER *header, double west, double east) {
392 /* Examine if the line defined by point k1-1 and k1 intersects or are inside the header->wesn region */
393 uint64_t k, k0 = k1 - 1;
394 int x_status[2], y_status[2];
395 double a, cross;
396 for (unsigned int s = 0, k = k0; k <= k1; s++, k++) {
397 x_status[s] = (S->data[GMT_X][k] < west) ? -1 : ((S->data[GMT_X][k] > east) ? +1 : 0);
398 y_status[s] = (S->data[GMT_Y][k] < header->wesn[YLO]) ? -1 : ((S->data[GMT_Y][k] > header->wesn[YHI]) ? +1 : 0);
399 }
400 if ((x_status[0] * x_status[1]) == 1) return false; /* Overlap not possible since both x-coordinates are outside either on west or east side */
401 if ((y_status[0] * y_status[1]) == 1) return false; /* Overlap not possible since both y-coordinates are outside either on south or north side */
402 /* Simple cases ruled out. Here the line may intersect the region. More checking needed */
403 if (x_status[0] == 0 && x_status[1] == 0) return true; /* With both x-values in range, at least one y value is inside or they are outside on opposite sides */
404 if (y_status[0] == 0 && y_status[1] == 0) return true; /* With both y-values in range, at least one x value is inside or they are outside on opposite sides */
405 /* Here we have remaining crossing checks */
406 a = (S->data[GMT_Y][k1] - S->data[GMT_Y][k0]) / (S->data[GMT_X][k1] - S->data[GMT_X][k0]); /* Line slope guaranteed not to be infinity */
407 for (k = XLO; k <= XHI; k++) { /* Try the west and east boundaries for intersections */
408 cross = a * (header->wesn[k] - S->data[GMT_X][k0]) + S->data[GMT_Y][k0]; /* y-coordinate of intersection */
409 if (cross >= header->wesn[YLO] && cross <= header->wesn[YHI]) return true; /* Yes found valid intersection */
410 }
411 /* If we are still here then we failed on the w/e boundaries. Now try n/s */
412 for (k = YLO; k <= YHI; k++) { /* Try the south and north boundaries for intersections */
413 cross = (header->wesn[k] - S->data[GMT_Y][k0]) / a + S->data[GMT_X][k0]; /* x-coordinate of intersection */
414 if (cross >= header->wesn[XLO] && cross <= header->wesn[XHI]) return true; /* Yes found valid intersection */
415 }
416 return false; /* Nope */
417 }
418
419 #define bailout(code) {gmt_M_free_options (mode); return (code);}
420 #define Return(code) {Free_Ctrl (GMT, Ctrl); gmt_set_pad (GMT, GMT_PAD_DEFAULT); gmt_M_free (GMT, use); gmt_M_free (GMT, Out); gmt_end_module (GMT, GMT_cpy); bailout (code);}
421
GMT_grdselect(void * V_API,int mode,void * args)422 EXTERN_MSC int GMT_grdselect (void *V_API, int mode, void *args) {
423 bool first_r = true, first = true, is_cube, pass, *use = NULL, one_layer;
424 int f = -1, error = 0, k_data, k_tile_id;
425 unsigned int n_cols = 0, cmode = GMT_COL_FIX, geometry = GMT_IS_TEXT;
426 uint64_t seg, tbl, row;
427
428 double wesn[8], out[8], last_inc[2], *subset = NULL;
429
430 char record[GMT_BUFSIZ] = {""};
431 static char *type[2] = {"grid", "cube"};
432
433 struct GRDSELECT_CTRL *Ctrl = NULL;
434 struct GMT_GRID *G = NULL;
435 struct GMT_GRID_HEADER *header = NULL;
436 struct GMT_CUBE *U = NULL;
437 struct GMT_DATASET *Df = NULL, *Dl = NULL, *Dc = NULL;
438 struct GMT_DATASEGMENT *S = NULL;
439 struct GMT_RECORD *Out = NULL;
440 struct GMT_OPTION *opt = NULL;
441 struct GMT_CTRL *GMT = NULL, *GMT_cpy = NULL;
442 struct GMT_OPTION *options = NULL;
443 struct GMTAPI_CTRL *API = gmt_get_api_ptr (V_API); /* Cast from void to GMTAPI_CTRL pointer */
444
445 /*----------------------- Standard module initialization and parsing ----------------------*/
446
447 if (API == NULL) return (GMT_NOT_A_SESSION);
448 if (mode == GMT_MODULE_PURPOSE) return (usage (API, GMT_MODULE_PURPOSE)); /* Return the purpose of program */
449 options = GMT_Create_Options (API, mode, args); if (API->error) return (API->error); /* Set or get option list */
450
451 if ((error = gmt_report_usage (API, options, 0, usage)) != GMT_NOERROR) bailout (error); /* Give usage if requested */
452
453 /* Parse the command-line arguments */
454
455 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 */
456 if (GMT_Parse_Common (API, THIS_MODULE_OPTIONS, options)) Return (API->error);
457 Ctrl = New_Ctrl (GMT); /* Allocate and initialize a new control structure */
458 if ((error = parse (GMT, Ctrl, options)) != 0) Return (error);
459
460 /*---------------------------- This is the grdselect main code ----------------------------*/
461
462 /* OK, done parsing, now process all input data sources in a loop */
463
464 gmt_M_memset (wesn, 8, double); /* Initialize */
465 gmt_M_memset (out, 8, double); /* Initialize */
466 use = gmt_M_memory (GMT, NULL, Ctrl->n_files, bool);
467
468 if (Ctrl->E.active) {
469 if (Ctrl->E.box) { /* Write polygon */
470 n_cols = 2;
471 cmode = GMT_COL_FIX_NO_TEXT;
472 geometry = GMT_IS_POLYGON;
473 }
474 else {
475 n_cols = (Ctrl->Q.active) ? 8 : 6; /* w e s n [z0 z1] [w0 w1] */
476 cmode = GMT_COL_FIX_NO_TEXT;
477 geometry = GMT_IS_NONE;
478 }
479 }
480
481 if (Ctrl->C.active) { /* Read the user's data point file */
482 if ((Dc = GMT_Read_Data (API, GMT_IS_DATASET, GMT_IS_FILE, GMT_IS_POINT, GMT_READ_NORMAL, NULL, Ctrl->C.file, NULL)) == NULL) {
483 Return (API->error);
484 }
485 }
486 if (Ctrl->F.active) { /* Read the user's polygon file */
487 if ((Df = GMT_Read_Data (API, GMT_IS_DATASET, GMT_IS_FILE, GMT_IS_POLY, GMT_READ_NORMAL, NULL, Ctrl->F.file, NULL)) == NULL) {
488 Return (API->error);
489 }
490 }
491 if (Ctrl->L.active) { /* Read the user's line file */
492 if ((Dl = GMT_Read_Data (API, GMT_IS_DATASET, GMT_IS_FILE, GMT_IS_LINE, GMT_READ_NORMAL, NULL, Ctrl->L.file, NULL)) == NULL) {
493 Return (API->error);
494 }
495 }
496
497 GMT_Set_Columns (GMT->parent, GMT_OUT, n_cols, cmode);
498
499 if (GMT_Init_IO (API, GMT_IS_DATASET, geometry, GMT_OUT, GMT_ADD_DEFAULT, 0, options) != GMT_NOERROR) { /* Registers default output destination, unless already set */
500 Return (API->error);
501 }
502 if (GMT_Begin_IO (API, GMT_IS_DATASET, GMT_OUT, GMT_HEADER_OFF) != GMT_NOERROR) { /* Enables data output and sets access mode */
503 Return (API->error);
504 }
505
506 Out = gmt_new_record (GMT, (n_cols) ? out : NULL, (cmode) == GMT_COL_FIX ? record : NULL); /* the two items */
507
508 gmt_set_pad (API->GMT, 0); /* Not read pads to simplify operations */
509
510 if (GMT->common.R.active[RSET]) { /* Gave a sub region so set the first region limit */
511 gmt_M_memcpy (wesn, GMT->common.R.wesn, 6U, double);
512 first_r = false;
513 subset = GMT->common.R.wesn; /* Read subsets */
514 }
515 wesn[ZLO] = wesn[WLO] = DBL_MAX; wesn[ZHI] = wesn[WHI] = -DBL_MAX;
516
517 for (opt = options; opt; opt = opt->next) { /* Loop over arguments, skip options */
518
519 if (opt->option != '<') continue; /* We are only processing filenames here */
520
521 f++; /* Here we are at data source number f */
522 gmt_set_cartesian (GMT, GMT_IN); /* Reset since we may get a bunch of files, some geo, some not */
523 one_layer = true; /* Assume we are reading grids */
524 k_tile_id = k_data = GMT_NOTSET;
525 if ((k_data = gmt_remote_dataset_id (API, opt->arg)) != GMT_NOTSET || (k_tile_id = gmt_get_tile_id (API, opt->arg)) != GMT_NOTSET) {
526 if (k_tile_id != GMT_NOTSET && !Ctrl->G.active) {
527 GMT_Report (API, GMT_MSG_WARNING, "Information on tiled remote global grids requires -G since download or all tiles may be required\n");
528 continue;
529 }
530 gmt_set_geographic (GMT, GMT_IN); /* Since this will be returned as a memory grid */
531 }
532
533 is_cube = gmt_nc_is_cube (API, opt->arg); /* Determine if this file is a cube or not */
534 if (Ctrl->Q.active != is_cube) {
535 if (is_cube)
536 GMT_Report (API, GMT_MSG_WARNING, "Detected a data cube (%s) but -Q not set - skipping\n", opt->arg);
537 else
538 GMT_Report (API, GMT_MSG_WARNING, "Detected a data grid (%s) but -Q is set - skipping\n", opt->arg);
539 continue;
540 }
541 if (is_cube) {
542 if ((U = GMT_Read_Data (API, GMT_IS_CUBE, GMT_IS_FILE, GMT_IS_VOLUME, GMT_CONTAINER_ONLY, NULL, opt->arg, NULL)) == NULL) {
543 Return (API->error);
544 }
545 header = U->header;
546 }
547 else {
548 if ((G = GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_ONLY, NULL, opt->arg, NULL)) == NULL) {
549 Return (API->error);
550 }
551 header = G->header;
552 one_layer = (header->n_bands == 1);
553 }
554
555 if (Ctrl->D.active) { /* Only pass data sources with the same increments */
556 pass = (doubleAlmostEqual (Ctrl->D.inc[GMT_X], header->inc[GMT_X]) && doubleAlmostEqual (Ctrl->D.inc[GMT_Y], header->inc[GMT_Y]));
557 if (pass != Ctrl->I.pass[GRD_SELECT_D]) /* Skip if grid spacing is different (or the same via -ID) */
558 continue;
559 }
560 if (Ctrl->W.active) { /* Skip data sources outside the imposed w-range */
561 if (!one_layer) {
562 GMT_Report (API, GMT_MSG_WARNING, "Cannot use -W with images - skipping that test\n", opt->arg);
563 }
564 else {
565 pass = !(header->z_max < Ctrl->W.w_min || header->z_min > Ctrl->W.w_max);
566 if (pass != Ctrl->I.pass[GRD_SELECT_W]) /* Skip if outside (or inside via -IW) the range */
567 continue;
568 }
569 }
570 if (Ctrl->Z.active) { /* Skip cubes outside the imposed z-range */
571 pass = !(U->z_range[1] < Ctrl->Z.z_min || U->z_range[0] > Ctrl->Z.z_max);
572 if (pass != Ctrl->I.pass[GRD_SELECT_Z]) /* Skip if outside (or inside via -IZ) the range */
573 continue;
574 }
575 if (GMT->common.R.active[GSET]) { /* Specified -r to only keep those data sources with the same registration */
576 pass = (header->registration == GMT->common.R.registration);
577 if (pass != Ctrl->I.pass[GRD_SELECT_r]) /* Skip if wrong (or right via -Ir) registration */
578 continue;
579 }
580 if (GMT->common.R.active[RSET]) { /* Specified -R to check for overlap */
581 bool pass_x, pass_y = !(header->wesn[YLO] > GMT->common.R.wesn[YHI] || header->wesn[YHI] < GMT->common.R.wesn[YLO]); /* true if there is y-overlap */
582 pass = false;
583 if (pass_y) { /* Still alive, now check for x-overlap */
584 if (gmt_M_is_cartesian (GMT, GMT_IN)) /* Easy, peasy */
585 pass_x = !(header->wesn[XLO] > GMT->common.R.wesn[XHI] || header->wesn[XHI] < GMT->common.R.wesn[XLO]); /* true if there is x-overlap */
586 else { /* Must worry about wrapping */
587 double w = header->wesn[XLO] - 720.0, e = header->wesn[XHI] - 720.0; /* Ensure we are far west */
588 while (e < GMT->common.R.wesn[XLO]) w += 360.0, e += 360.0; /* Wind to the west */
589 pass_x = (w < GMT->common.R.wesn[XHI]);
590 }
591 pass = pass_x && pass_y; /* Only true if we have both x and y overlap */
592 }
593 if (pass != Ctrl->I.pass[GRD_SELECT_R]) /* Skip if wrong (or right via -IR) registration */
594 continue;
595 }
596 if (Ctrl->C.active) { /* Specified -P to only keep those data sources that contain one or more of these points */
597 bool overlap, crossed;
598 double w = header->wesn[XLO], e = header->wesn[XHI]; /* This is set once if Cartesian; otherwise per segment */
599 uint64_t row;
600 pass = true;
601 for (tbl = 0; pass && tbl < Dc->n_tables; tbl++) {
602 for (seg = 0; pass && seg < Dc->table[tbl]->n_segments; seg++) {
603 S = Dc->table[tbl]->segment[seg]; /* Shorthand to current segment */
604 for (row = 0; pass && row < S->n_rows; row++) {
605 if (is_cube && S->n_columns > GMT_Y) { /* Check if point is outside cube z-range */
606 if (S->data[GMT_Z][row] < U->z_range[0] || S->data[GMT_Z][row] > U->z_range[1])
607 continue; /* Point outside cube layer */
608 }
609 if (S->data[GMT_Y][row] < header->wesn[YLO] || S->data[GMT_Y][row] > header->wesn[YHI])
610 continue; /* Point outside grid */
611 if (gmt_M_is_cartesian (GMT, GMT_IN)) /* Easy, peasy */
612 pass = (S->data[GMT_X][row] < w || S->data[GMT_X][row] > e); /* true if there is no x-overlap */
613 else { /* Must worry about wrapping */
614 w = header->wesn[XLO] - 720.0; e = header->wesn[XHI] - 720.0; /* Ensure we are far west */
615 while (e < S->data[GMT_X][row]) w += 360.0, e += 360.0; /* Wind to the west */
616 pass = (w > S->data[GMT_X][row]);
617 }
618 }
619 }
620 }
621 if (pass == Ctrl->I.pass[GRD_SELECT_C]) /* Skip if not overlap (or overlap via -IC) */
622 continue;
623 }
624 if (Ctrl->L.active) { /* Specified -L to only keep those data sources that are traversed by these lines in map view */
625 bool overlap_x, overlap_y, overlap, crossed;
626 double w = header->wesn[XLO], e = header->wesn[XHI]; /* This is set once if Cartesian; otherwise per segment */
627 pass = true;
628 for (tbl = 0; pass && tbl < Dl->n_tables; tbl++) {
629 for (seg = 0; pass && seg < Dl->table[tbl]->n_segments; seg++) {
630 S = Dl->table[tbl]->segment[seg]; /* Shorthand to current segment */
631 overlap_y = !(header->wesn[YLO] > S->max[GMT_Y] || header->wesn[YHI] < S->min[GMT_Y]); /* true if there is potential y-overlap */
632 if (overlap_y) { /* Still alive, now check for x-overlap */
633 if (gmt_M_is_cartesian (GMT, GMT_IN)) { /* Easy, peasy */
634 overlap_x = !(w > S->max[GMT_X] || e < S->min[GMT_X]); /* true if there is x-overlap */
635 }
636 else { /* Must worry about wrapping */
637 w = header->wesn[XLO] - 720.0; e = header->wesn[XHI] - 720.0; /* Ensure we are far west */
638 while (e < S->min[GMT_X]) w += 360.0, e += 360.0; /* Wind to the west */
639 overlap_x = (w < S->max[GMT_X]);
640 }
641 overlap = overlap_x && overlap_y; /* Only true if we have both x and y overlap */
642 if (!overlap) /* No overlap possible with this segment so move to next */
643 continue;
644 /* Here we have potential overlap and must not do the harder work of checking if there are crossings between polygon and BB */
645 for (row = 1, crossed = false; !crossed && row < S->n_rows; row++) {
646 if (grdselect_line_overlaps (S, row, header, w, e))
647 crossed = true;
648 }
649 if (crossed) /* Cannot pass since we have overlap */
650 pass = false;
651 }
652 }
653 }
654 if (pass == Ctrl->I.pass[GRD_SELECT_L]) /* Skip if not overlap (or overlap via -IL) */
655 continue;
656 }
657 if (Ctrl->F.active) { /* Specified -F to only keep those data sources that overlap fully or partially with these polygons in map view */
658 /* There are two checks to make here:
659 * 1. If center point of grid is inside polygon then we know we have overlap.
660 * 2. If not, then try all points in polygon and if any is inside grid boundary then we have overlap
661 * If we fail both then we have no overlap. These tests must be made on all segments in the polygon file
662 * and if any segment yields overlap then we have overlap and can jump to decision time
663 */
664 bool overlap_x, overlap_y, overlap, crossed;
665 unsigned int io_status = (Ctrl->F.mode != GMT_ONEDGE) ? Ctrl->F.mode : GMT_INSIDE;
666 double x0 = 0.5 * (header->wesn[XLO] + header->wesn[XHI]); /* Mod point of grid */
667 double y0 = 0.5 * (header->wesn[YLO] + header->wesn[YHI]); /* Mod point of grid */
668 double w = header->wesn[XLO], e = header->wesn[XHI]; /* This is set once if Cartesian; otherwise per segment */
669 pass = true;
670 for (tbl = 0; pass && tbl < Df->n_tables; tbl++) {
671 for (seg = 0; pass && seg < Df->table[tbl]->n_segments; seg++) {
672 S = Df->table[tbl]->segment[seg]; /* Shorthand to current segment */
673 overlap_y = !(header->wesn[YLO] > S->max[GMT_Y] || header->wesn[YHI] < S->min[GMT_Y]); /* true if there is potential y-overlap */
674 if (overlap_y) { /* Still alive, now check for x-overlap */
675 if (gmt_M_is_cartesian (GMT, GMT_IN)) { /* Easy, peasy */
676 overlap_x = !(w > S->max[GMT_X] || e < S->min[GMT_X]); /* true if there is x-overlap */
677 }
678 else { /* Must worry about wrapping */
679 w = header->wesn[XLO] - 720.0; e = header->wesn[XHI] - 720.0; /* Ensure we are far west */
680 while (e < S->min[GMT_X]) w += 360.0, e += 360.0; /* Wind to the west */
681 overlap_x = (w < S->max[GMT_X]);
682 }
683 overlap = overlap_x && overlap_y; /* Only true if we have both x and y overlap */
684 if (!overlap) /* No overlap possible with this segment so move to next */
685 continue;
686 /* Here we have potential overlap and must not do the harder work of checking if there are crossings between polygon and BB */
687 for (row = 1, crossed = false; !crossed && row < S->n_rows; row++) {
688 if (grdselect_line_overlaps (S, row, header, w, e))
689 crossed = true;
690 }
691 if (crossed) /* Cannot pass since we have overlap, but depends on +i or not */
692 pass = (Ctrl->F.mode != GMT_ONEDGE);
693 else /* Either polygon swallowed grid or grid swallowed polygon - which is it? */
694 pass = (gmt_inonout (GMT, x0, y0, S) != io_status);
695 }
696 }
697 }
698 if (pass == Ctrl->I.pass[GRD_SELECT_F]) /* Skip if not overlap (or overlap via -IF) */
699 continue;
700 }
701 if (Ctrl->N.active && !one_layer)
702 GMT_Report (API, GMT_MSG_WARNING, "Cannot use -N with images - skipping that test\n", opt->arg);
703 else if (Ctrl->N.active) { /* Must read the data to know how many NaNs, then skip data sources that fail the test (or pass if -In) */
704 uint64_t level, here = 0, ij, n_nan = 0;
705 gmt_grdfloat *data = NULL;
706 if (is_cube) {
707 if (GMT_Read_Data (API, GMT_IS_CUBE, GMT_IS_FILE, GMT_IS_VOLUME, GMT_DATA_ONLY, subset, opt->arg, U) == NULL) {
708 Return (API->error);
709 }
710 }
711 else {
712 if (GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_DATA_ONLY, subset, opt->arg, G) == NULL) {
713 Return (API->error);
714 }
715 }
716 data = (is_cube) ? U->data : G->data;
717
718 for (level = 0; level < header->n_bands; level++) {
719 for (ij = 0; ij < header->size; ij++) {
720 if (gmt_M_is_fnan (data[ij + here])) n_nan++;
721 }
722 here += header->size;
723 }
724 if (is_cube && GMT_Destroy_Data (API, &U) != GMT_NOERROR) {
725 Return (API->error);
726 }
727 else if (!is_cube && GMT_Destroy_Data (API, &G) != GMT_NOERROR) {
728 Return (API->error);
729 }
730 pass = true;
731 if (Ctrl->N.mode == GRDSELECT_LESS_NANS && n_nan > Ctrl->N.n_nans) pass = false; /* Skip this item */
732 else if (Ctrl->N.mode == GRDSELECT_MORE_NANS && n_nan < Ctrl->N.n_nans) pass = false; /* Skip this item */
733 if (pass != Ctrl->I.pass[GRD_SELECT_N]) /* Skip if outside (or inside via -IN) the range */
734 continue;
735 }
736
737 /* OK, here the grid/cube passed any other obstacles */
738 if (first_r) { /* No w/e/s/n initialized via -R, use the first header to initialize wesn */
739 gmt_M_memcpy (wesn, header->wesn, 4U, double); /* The first grid needs to set the extent of the region for now */
740 first_r = false;
741 }
742 if (first) { /* This is the first grid/cube */
743 if (Ctrl->Q.active) { /* Cube: Keep track of z-dimension range */
744 wesn[ZLO] = U->z[0];
745 wesn[ZHI] = U->z[header->n_bands-1];
746 }
747 /* Set the initial data range */
748 wesn[WLO] = header->z_min;
749 wesn[WHI] = header->z_max;
750 if (Ctrl->A.i_mode == GRDSELECT_MIN_INC || Ctrl->A.i_mode == GRDSELECT_MAX_INC) /* Copy the first grid increments */
751 gmt_M_memcpy (Ctrl->A.inc, header->inc, 2U, double);
752 gmt_M_memcpy (last_inc, header->inc, 2U, double);
753 first = false;
754 }
755 if (!(doubleAlmostEqual (header->inc[GMT_X], last_inc[GMT_X]) && doubleAlmostEqual (header->inc[GMT_Y], last_inc[GMT_Y]))) {
756 /* The grids have different increments */
757 if (Ctrl->E.active && !Ctrl->A.round && !Ctrl->E.box) {
758 GMT_Report (API, GMT_MSG_ERROR, "Using -C with mixed-increment grids requires a suitable -A+i rounding selection\n");
759 Return (GMT_RUNTIME_ERROR);
760 }
761 }
762 use[f] = true; /* May need to read a subset from this grid */
763 if (Ctrl->A.active) { /* Want region information */
764 if (Ctrl->A.mode == GRDSELECT_INTERSECTION) { /* Retain only region common to all data sources (i.e., their intersection) */
765 wesn[YLO] = MAX (wesn[YLO], header->wesn[YLO]);
766 wesn[YHI] = MIN (wesn[YHI], header->wesn[YHI]);
767 if (wesn[YHI] <= wesn[YLO]) { /* No common region can be found - no point to continue */
768 GMT_Report (API, GMT_MSG_WARNING, "No common area for all %ss is possible.\n", type[is_cube]);
769 goto nothing;
770 }
771 if (is_cube) {
772 wesn[ZLO] = MAX (wesn[ZLO], U->z[0]);
773 wesn[ZHI] = MIN (wesn[ZHI], U->z[header->n_bands-1]);
774 if (wesn[ZHI] <= wesn[ZLO]) { /* No common cube region can be found - no point to continue */
775 GMT_Report (API, GMT_MSG_WARNING, "No common area for all cubes is possible.\n");
776 goto nothing;
777 }
778 }
779 if (gmt_M_is_cartesian (GMT, GMT_IN)) { /* Easy, peasy */
780 wesn[XLO] = MAX (wesn[XLO], header->wesn[XLO]);
781 wesn[XHI] = MIN (wesn[XHI], header->wesn[XHI]);
782 }
783 else { /* Must worry about 360 wrapping */
784 double w = header->wesn[XLO] - 720.0, e = header->wesn[XHI] - 720.0; /* Ensure we are far west */
785 while (e < wesn[XLO]) w += 360.0, e += 360.0; /* Wind to the west */
786 wesn[XLO] = MAX (wesn[XLO], w);
787 wesn[XHI] = MIN (wesn[XHI], e);
788 }
789 if (wesn[XHI] <= wesn[XLO]) { /* No common region can be found - no point to continue */
790 GMT_Report (API, GMT_MSG_WARNING, "No common area for all %ss is possible.\n", type[is_cube]);
791 goto nothing;
792 }
793 }
794 else { /* Find the union of all regions */
795 wesn[YLO] = MIN (wesn[YLO], header->wesn[YLO]);
796 wesn[YHI] = MAX (wesn[YHI], header->wesn[YHI]);
797 if (gmt_M_is_cartesian (GMT, GMT_IN)) { /* Easy, peasy */
798 wesn[XLO] = MIN (wesn[XLO], header->wesn[XLO]);
799 wesn[XHI] = MAX (wesn[XHI], header->wesn[XHI]);
800 }
801 else { /* Must worry about 360 wrapping */
802 double w = header->wesn[XLO] - 720.0, e = header->wesn[XHI] - 720.0; /* Ensure we are far west */
803 while (e < wesn[XLO]) w += 360.0, e += 360.0; /* Wind to the west */
804 wesn[XLO] = MIN (wesn[XLO], w);
805 wesn[XHI] = MAX (wesn[XHI], e);
806 }
807 if (is_cube) {
808 wesn[ZLO] = MIN (wesn[ZLO], U->z[0]);
809 wesn[ZHI] = MAX (wesn[ZHI], U->z[header->n_bands-1]);
810 }
811 }
812 wesn[WLO] = MIN (wesn[WLO], header->z_min);
813 wesn[WHI] = MAX (wesn[WHI], header->z_max);
814 if (Ctrl->A.i_mode == GRDSELECT_MIN_INC) { /* Update the smallest increments found */
815 if (header->inc[GMT_X] < Ctrl->A.inc[GMT_X]) Ctrl->A.inc[GMT_X] = header->inc[GMT_X]; /* Update the smallest x-increments found */
816 if (header->inc[GMT_Y] < Ctrl->A.inc[GMT_Y]) Ctrl->A.inc[GMT_Y] = header->inc[GMT_Y]; /* Update the smallest y-increments found */
817 }
818 if (Ctrl->A.i_mode == GRDSELECT_MAX_INC) { /* Update the largest increments found */
819 if (header->inc[GMT_X] > Ctrl->A.inc[GMT_X]) Ctrl->A.inc[GMT_X] = header->inc[GMT_X]; /* Update the largest x-increments found */
820 if (header->inc[GMT_Y] > Ctrl->A.inc[GMT_Y]) Ctrl->A.inc[GMT_Y] = header->inc[GMT_Y]; /* Update the largest y-increments found */
821 }
822 }
823 else { /* Only wanted to list the grid files that passed the tests */
824 sprintf (record, "%s", opt->arg);
825 GMT_Put_Record (API, GMT_WRITE_DATA, Out);
826 }
827 }
828
829 if (Ctrl->A.active) { /* Finalize region via rounding/padding, then report it */
830 if (Ctrl->A.round) { /* Must round the region via the increments */
831 wesn[XLO] = floor (wesn[XLO] / Ctrl->A.inc[GMT_X]) * Ctrl->A.inc[GMT_X];
832 wesn[XHI] = ceil (wesn[XHI] / Ctrl->A.inc[GMT_X]) * Ctrl->A.inc[GMT_X];
833 wesn[YLO] = floor (wesn[YLO] / Ctrl->A.inc[GMT_Y]) * Ctrl->A.inc[GMT_Y];
834 wesn[YHI] = ceil (wesn[YHI] / Ctrl->A.inc[GMT_Y]) * Ctrl->A.inc[GMT_Y];
835 }
836 if (Ctrl->M.active) { /* Extend region by given margins */
837 wesn[XLO] -= Ctrl->M.margin[XLO];
838 wesn[XHI] += Ctrl->M.margin[XHI];
839 wesn[YLO] -= Ctrl->M.margin[YLO];
840 wesn[YHI] += Ctrl->M.margin[YHI];
841 }
842 if (gmt_M_is_geographic (GMT, GMT_IN)) { /* Trim to fit the sphere */
843 if ((wesn[XHI] - wesn[XLO] > 360.0)) { /* Trim to 360 range */
844 if (wesn[XLO] >= 0.0 || wesn[XHI] > 300.0) /* Mostly positive longitudes so use 0/360 */
845 wesn[XLO] = 0.0, wesn[XHI] = 360.0;
846 else /* Use dateline as the jump */
847 wesn[XLO] = -180.0, wesn[XHI] = 18.0;
848 }
849 /* Ensure latitudes are not out of bound */
850 if (wesn[YLO] < -90.0) wesn[YLO] = -90.0;
851 if (wesn[YHI] > +90.0) wesn[YHI] = +90.0;
852 }
853 if (Ctrl->E.box) { /* Output closed polygon */
854 out[GMT_X] = wesn[XLO]; out[GMT_Y] = wesn[YLO]; GMT_Put_Record (API, GMT_WRITE_DATA, Out);
855 out[GMT_X] = wesn[XHI]; GMT_Put_Record (API, GMT_WRITE_DATA, Out);
856 out[GMT_Y] = wesn[YHI]; GMT_Put_Record (API, GMT_WRITE_DATA, Out);
857 out[GMT_X] = wesn[XLO]; GMT_Put_Record (API, GMT_WRITE_DATA, Out);
858 out[GMT_Y] = wesn[YLO]; GMT_Put_Record (API, GMT_WRITE_DATA, Out);
859 }
860 else if (Ctrl->E.active) { /* Want numerical output */
861 if (Ctrl->A.mode == GRDSELECT_INTERSECTION) { /* Must revisit the grid subsets to learn about actual w-range */
862 wesn[WLO] = DBL_MAX; wesn[WHI] = -DBL_MAX; /* Need to redo these within the intersection region only */
863 f = -1; /* Start index over again */
864 for (opt = options; opt; opt = opt->next) { /* Loop over arguments, skip options */
865
866 if (opt->option != '<') continue; /* We are only processing filenames here */
867
868 f++; /* Here we are at data source number f */
869 if (!use[f]) continue; /* This one failed some test earlier */
870 if (is_cube) {
871 if (GMT_Read_Data (API, GMT_IS_CUBE, GMT_IS_FILE, GMT_IS_VOLUME, GMT_CONTAINER_AND_DATA, wesn, opt->arg, U) == NULL) {
872 Return (API->error);
873 }
874 }
875 else {
876 if (GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA, wesn, opt->arg, G) == NULL) {
877 Return (API->error);
878 }
879 }
880 /* Update the common data range */
881 wesn[WLO] = MIN (wesn[WLO], header->z_min);
882 wesn[WHI] = MAX (wesn[WHI], header->z_max);
883 if (is_cube && GMT_Destroy_Data (API, &U) != GMT_NOERROR) {
884 Return (API->error);
885 }
886 else if (!is_cube && GMT_Destroy_Data (API, &G) != GMT_NOERROR) {
887 Return (API->error);
888 }
889 }
890 }
891 if (is_cube) { /* Need w e s n b t w0 w1 */
892 gmt_M_memcpy (out, wesn, 6, double);
893 gmt_M_memcpy (&out[6], &wesn[WLO], 2, double);
894 }
895 else { /* Need w e s n w0 w1 */
896 gmt_M_memcpy (out, wesn, 4, double);
897 gmt_M_memcpy (&out[4], &wesn[WLO], 2, double);
898 }
899 GMT_Put_Record (API, GMT_WRITE_DATA, Out);
900 }
901 else { /* Give the string -Rw/e/s/n[/b/t] only */
902 char text[GMT_LEN64] = {""};
903 sprintf (record, "-R");
904 gmt_ascii_format_col (GMT, text, wesn[XLO], GMT_OUT, GMT_X); strcat (record, text); strcat (record, "/");
905 gmt_ascii_format_col (GMT, text, wesn[XHI], GMT_OUT, GMT_X); strcat (record, text); strcat (record, "/");
906 gmt_ascii_format_col (GMT, text, wesn[YLO], GMT_OUT, GMT_Y); strcat (record, text); strcat (record, "/");
907 gmt_ascii_format_col (GMT, text, wesn[YHI], GMT_OUT, GMT_Y); strcat (record, text);
908 if (is_cube) { /* Must also report the z-range */
909 gmt_ascii_format_col (GMT, text, wesn[ZLO], GMT_OUT, GMT_Z); strcat (record, text); strcat (record, "/");
910 gmt_ascii_format_col (GMT, text, wesn[ZHI], GMT_OUT, GMT_Z); strcat (record, text);
911 }
912 GMT_Put_Record (API, GMT_WRITE_DATA, Out);
913 }
914 }
915
916 nothing:
917
918 gmt_M_free (GMT, Out);
919 gmt_M_free (GMT, use);
920
921 if (GMT_End_IO (API, GMT_OUT, 0) != GMT_NOERROR) { /* Disables further data output */
922 Return (API->error);
923 }
924
925 Return (GMT_NOERROR);
926 }
927