1 /*--------------------------------------------------------------------
2 *
3 * Copyright (c) 1991-2021 by the GMT Team (https://www.generic-mapping-tools.org/team.html)
4 * See LICENSE.TXT file for copying and redistribution conditions.
5 *
6 * This program is free software; you can redistribute it and/or modify
7 * it under the terms of the GNU Lesser General Public License as published by
8 * the Free Software Foundation; version 3 or any later version.
9 *
10 * This program is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 * GNU Lesser General Public License for more details.
14 *
15 * Contact info: www.generic-mapping-tools.org
16 *--------------------------------------------------------------------*/
17 /*
18 * Brief synopsis: grdlandmask defines a grid based on region and xinc/yinc values,
19 * reads a shoreline data base, and sets the grid nodes inside, on the
20 * boundary, and outside of the polygons to the user-defined values
21 * <in>, <on>, and <out>. These may be any number, including NaN.
22 *
23 * Author: P. Wessel
24 * Date: 1-JAN-2010
25 * Version: 6 API
26 */
27
28 #include "gmt_dev.h"
29
30 #define THIS_MODULE_CLASSIC_NAME "grdlandmask"
31 #define THIS_MODULE_MODERN_NAME "grdlandmask"
32 #define THIS_MODULE_LIB "core"
33 #define THIS_MODULE_PURPOSE "Create a \"wet-dry\" mask grid from shoreline data base"
34 #define THIS_MODULE_KEYS "GG}"
35 #define THIS_MODULE_NEEDS "R"
36 #define THIS_MODULE_OPTIONS "-VRr" GMT_ADD_x_OPT GMT_OPT("F")
37
38 #define GRDLANDMASK_N_CLASSES (2*GSHHS_MAX_LEVEL + 1) /* Number of bands separated by the levels */
39 /* There are 4 GSHHG levels, but we may find points exactly on the border between levels.
40 * To handle these cases we must double the number of levels. The N.mask[] array will hold the
41 * data values we wish to assign to each of these 9 potential cases. As we determine where a
42 * node fits into this hierarchy we assign values 0-9 to the grid. Then, at the end we replace
43 * these values with the corresponding entries in N.mask[].
44 */
45
46 struct GRDLANDMASK_CTRL { /* All control options for this program (except common args) */
47 /* ctive is true if the option has been activated */
48 struct GRDLANDMASK_A { /* -A<min_area>[/<min_level>/<max_level>] */
49 bool active;
50 struct GMT_SHORE_SELECT info;
51 } A;
52 struct GRDLANDMASK_D { /* -D<resolution>[+f] */
53 bool active;
54 bool force; /* if true, select next highest level if current set is not available */
55 char set; /* One of f, h, i, l, c, or a for auto */
56 } D;
57 struct GRDLANDMASK_E { /* -E[<border> | <cborder>/<lborder>/<iborder>/<pborder>] */
58 bool active;
59 bool linetrace; /* True if -E was given arguments as we must trace the lines through the grid */
60 unsigned int inside; /* if 2, then a point exactly on a polygon boundary is considered OUTSIDE, else 1 */
61 } E;
62 struct GRDLANDMASK_G { /* -G<maskfile> */
63 bool active;
64 char *file;
65 } G;
66 struct GRDLANDMASK_I { /* -I (for checking only) */
67 bool active;
68 } I;
69 struct GRDLANDMASK_N { /* -N<maskvalues> */
70 bool active;
71 unsigned int wetdry; /* 1 if dry/wet only, 0 if 5 mask levels */
72 gmt_grdfloat mask[GRDLANDMASK_N_CLASSES]; /* values for each level */
73 } N;
74 };
75
76 struct GRDLANDMASK_BINCROSS {
77 double x, y, d;
78 };
79
grdlandmask_comp_bincross(const void * p1,const void * p2)80 GMT_LOCAL int grdlandmask_comp_bincross (const void *p1, const void *p2) {
81 const struct GRDLANDMASK_BINCROSS *a = p1, *b = p2;
82
83 if (a->d < b->d) return (-1);
84 if (a->d > b->d) return (+1);
85 return (0);
86 }
87
New_Ctrl(struct GMT_CTRL * GMT)88 static void *New_Ctrl (struct GMT_CTRL *GMT) { /* Allocate and initialize a new control structure */
89 struct GRDLANDMASK_CTRL *C;
90
91 C = gmt_M_memory (GMT, NULL, 1, struct GRDLANDMASK_CTRL);
92
93 /* Initialize values whose defaults are not 0/false/NULL */
94
95 C->A.info.high = GSHHS_MAX_LEVEL; /* Include all GSHHS levels */
96 C->D.set = 'l'; /* Low-resolution coastline data */
97 C->E.inside = GMT_ONEDGE; /* Default is that points on a boundary are inside */
98 gmt_M_memset (C->N.mask, GRDLANDMASK_N_CLASSES, gmt_grdfloat); /* Default "wet" value = 0 */
99 C->N.mask[2] = C->N.mask[6] = 1.0f; /* Default for "dry" areas = 1 (inside) */
100
101 return (C);
102 }
103
Free_Ctrl(struct GMT_CTRL * GMT,struct GRDLANDMASK_CTRL * C)104 static void Free_Ctrl (struct GMT_CTRL *GMT, struct GRDLANDMASK_CTRL *C) { /* Deallocate control structure */
105 if (!C) return;
106 gmt_M_str_free (C->G.file);
107 gmt_M_free (GMT, C);
108 }
109
usage(struct GMTAPI_CTRL * API,int level)110 static int usage (struct GMTAPI_CTRL *API, int level) {
111 const char *name = gmt_show_name_and_purpose (API, THIS_MODULE_LIB, THIS_MODULE_CLASSIC_NAME, THIS_MODULE_PURPOSE);
112 if (level == GMT_MODULE_PURPOSE) return (GMT_NOERROR);
113 GMT_Usage (API, 0, "usage: %s -G%s %s %s [%s] [-D<resolution>[+f]] [-E[<bordervalues>]] [-N<maskvalues>] "
114 "[%s] [%s] %s [%s]\n", name, GMT_OUTGRID, GMT_I_OPT, GMT_Rgeo_OPT, GMT_A_OPT, GMT_V_OPT, GMT_r_OPT, GMT_x_OPT, GMT_PAR_OPT);
115
116 if (level == GMT_SYNOPSIS) return (GMT_MODULE_SYNOPSIS);
117
118 GMT_Message (API, GMT_TIME_NONE, " REQUIRED ARGUMENTS:\n");
119 gmt_outgrid_syntax (API, 'G', "Specify file name for output mask grid");
120 GMT_Option (API, "I,R");
121 GMT_Message (API, GMT_TIME_NONE, "\n OPTIONAL ARGUMENTS:\n");
122 gmt_GSHHG_syntax (API->GMT, 'A');
123 gmt_GSHHG_resolution_syntax (API->GMT, 'D', "Alternatively, choose (a)uto to automatically select the best "
124 "resolution given the chosen region.");
125 GMT_Usage (API, 1, "\n-E[<bordervalues>]");
126 GMT_Usage (API, -2, "Indicate that nodes exactly on a polygon boundary are outside [inside]. Optionally append "
127 "<border> or <cborder>/<lborder>/<iborder>/<pborder>. We will then trace lines through the grid and reset the "
128 "cells crossed by the lines to the indicated values [Default is no line tracing].");
129 GMT_Usage (API, 1, "\n-N<maskvalues>");
130 GMT_Usage (API, -2, "Give values to use if a node is outside or inside a feature. Specify this information using 1 "
131 "of 2 formats:");
132 GMT_Usage (API, 3, "%s <wet>/<dry>.", GMT_LINE_BULLET);
133 GMT_Usage (API, 3, "%s <ocean>/<land>/<lake>/<island>/<pond>.", GMT_LINE_BULLET);
134 GMT_Usage (API, -2, "NaN is a valid entry. [Default values are 0/1/0/1/0 (i.e., 0/1)].");
135 GMT_Option (API, "V,r,x,.");
136
137 return (GMT_MODULE_USAGE);
138 }
139
parse(struct GMT_CTRL * GMT,struct GRDLANDMASK_CTRL * Ctrl,struct GMT_OPTION * options)140 static int parse (struct GMT_CTRL *GMT, struct GRDLANDMASK_CTRL *Ctrl, struct GMT_OPTION *options) {
141 /* This parses the options provided to grdlandmask and sets parameters in CTRL.
142 * Any GMT common options will override values set previously by other commands.
143 * It also replaces any file names specified as input or output with the data ID
144 * returned when registering these sources/destinations with the API.
145 */
146
147 unsigned int n_errors = 0, j, pos, n_files = 0;
148 char line[GMT_LEN256] = {""}, ptr[GMT_BUFSIZ] = {""};
149 struct GMT_OPTION *opt = NULL;
150 struct GMTAPI_CTRL *API = GMT->parent;
151
152 for (opt = options; opt; opt = opt->next) {
153 switch (opt->option) {
154
155 case '<': /* Input files */
156 n_files++;
157 break;
158
159 /* Processes program-specific parameters */
160
161 case 'A': /* Restrict GSHHS features */
162 n_errors += gmt_M_repeated_module_option (API, Ctrl->A.active);
163 Ctrl->A.active = true;
164 n_errors += gmt_set_levels (GMT, opt->arg, &Ctrl->A.info);
165 break;
166 case 'D': /* Set GSHHS resolution */
167 n_errors += gmt_M_repeated_module_option (API, Ctrl->D.active);
168 Ctrl->D.active = true;
169 Ctrl->D.set = opt->arg[0];
170 Ctrl->D.force = (opt->arg[1] == '+' && (opt->arg[2] == 'f' || opt->arg[2] == '\0'));
171 break;
172 case 'E': /* On-boundary setting */
173 n_errors += gmt_M_repeated_module_option (API, Ctrl->E.active);
174 Ctrl->E.active = true;
175 if (opt->arg[0]) { /* Trace lines through grid */
176 Ctrl->E.linetrace = true;
177 j = pos = 0;
178 strncpy (line, opt->arg, GMT_LEN256-1);
179 while (j < 4 && (gmt_strtok (line, "/", &pos, ptr))) {
180 Ctrl->N.mask[2*j+1] = (ptr[0] == 'N' || ptr[0] == 'n') ? GMT->session.f_NaN : (gmt_grdfloat)atof (ptr);
181 j++;
182 }
183 if (!(j == 1 || j == 4)) {
184 GMT_Report (API, GMT_MSG_ERROR, "Option -E: Specify 1 or 4 border values\n");
185 n_errors++;
186 }
187 if (j == 1) Ctrl->N.mask[3] = Ctrl->N.mask[5] = Ctrl->N.mask[7] = Ctrl->N.mask[1]; /* Duplicate border values */
188 }
189 else
190 Ctrl->E.inside = GMT_INSIDE;
191 break;
192 case 'G': /* Output filename */
193 n_errors += gmt_M_repeated_module_option (API, Ctrl->G.active);
194 Ctrl->G.active = true;
195 if (opt->arg[0]) Ctrl->G.file = strdup (opt->arg);
196 if (GMT_Get_FilePath (API, GMT_IS_GRID, GMT_OUT, GMT_FILE_LOCAL, &(Ctrl->G.file))) n_errors++;
197 break;
198 case 'I': /* Grid spacings */
199 n_errors += gmt_M_repeated_module_option (API, Ctrl->I.active);
200 Ctrl->I.active = true;
201 n_errors += gmt_parse_inc_option (GMT, 'I', opt->arg);
202 break;
203 case 'N': /* Mask values */
204 Ctrl->N.active = true;
205 strncpy (line, opt->arg, GMT_LEN256);
206 if (line[strlen(line)-1] == 'o' && gmt_M_compat_check (GMT, 4)) { /* Edge is considered outside */
207 GMT_Report (API, GMT_MSG_COMPAT, "Option -N...o is deprecated; use -E instead\n");
208 Ctrl->E.active = true;
209 Ctrl->E.inside = GMT_INSIDE;
210 line[strlen(line)-1] = 0;
211 }
212 j = pos = 0;
213 while (j < 5 && (gmt_strtok (line, "/", &pos, ptr))) {
214 Ctrl->N.mask[2*j] = (ptr[0] == 'N' || ptr[0] == 'n') ? GMT->session.f_NaN : (gmt_grdfloat)atof (ptr);
215 j++;
216 }
217 if (!(j == 2 || j == 5)) {
218 GMT_Report (API, GMT_MSG_ERROR, "Option -N: Specify 2 or 5 mask values\n");
219 n_errors++;
220 }
221 Ctrl->N.wetdry = (j == 2);
222 break;
223 default: /* Report bad options */
224 n_errors += gmt_default_error (GMT, opt->option);
225 break;
226 }
227 }
228
229 n_errors += gmt_M_check_condition (GMT, !GMT->common.R.active[RSET], "Must specify -R option\n");
230 n_errors += gmt_M_check_condition (GMT, GMT->common.R.inc[GMT_X] <= 0.0 || GMT->common.R.inc[GMT_Y] <= 0.0, "Option -I: Must specify positive increment(s)\n");
231 n_errors += gmt_M_check_condition (GMT, !Ctrl->G.file, "Option -G: Must specify an output file\n");
232 n_errors += gmt_M_check_condition (GMT, n_files, "No input files allowed.\n");
233
234 return (n_errors ? GMT_PARSE_ERROR : GMT_NOERROR);
235 }
236
grdlandmask_inside(struct GMT_CTRL * GMT,struct GMT_GRID_HEADER * C,double x,double y)237 GMT_LOCAL bool grdlandmask_inside (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *C, double x, double y) {
238 /* Return true if this point is inside the grid */
239 int row, col;
240 gmt_M_unused(GMT);
241 row = gmt_M_grd_y_to_row (GMT, y, C);
242 if (row < 0 || row > (int)C->n_rows) return false;
243 col = gmt_M_grd_x_to_col (GMT, x, C);
244 if (col < 0 || col > (int)C->n_columns) return false;
245 return true;
246 }
247
grdlandmask_assign_node(struct GMT_CTRL * GMT,struct GMT_GRID * G,struct GMT_GRID_HEADER * C,double x,double y,gmt_grdfloat f_level,int64_t * ij)248 GMT_LOCAL void grdlandmask_assign_node (struct GMT_CTRL *GMT, struct GMT_GRID *G, struct GMT_GRID_HEADER *C, double x, double y, gmt_grdfloat f_level, int64_t *ij) {
249 /* Set node value and return node index if this point is inside the grid */
250 int row, col;
251 gmt_M_unused(GMT);
252 row = gmt_M_grd_y_to_row (GMT, y, C);
253 if (row < 0 || row > (int)C->n_rows) return;
254 col = gmt_M_grd_x_to_col (GMT, x, C);
255 if (col < 0 || col > (int)C->n_columns) return;
256 *ij = gmt_M_ijp (G->header, row, col);
257 G->data[*ij] = f_level;
258 }
259
260 #define bailout(code) {gmt_M_free_options (mode); return (code);}
261 #define Return(code) {Free_Ctrl (GMT, Ctrl); gmt_end_module (GMT, GMT_cpy); bailout (code);}
262
GMT_grdlandmask(void * V_API,int mode,void * args)263 EXTERN_MSC int GMT_grdlandmask (void *V_API, int mode, void *args) {
264 bool temp_shift = false, wrap, used_polygons, double_dip;
265 unsigned int base = 3, k, bin, np, side, np_new;
266 int row, row_min, row_max, ii, col, col_min, col_max, i, direction, err, ind, nx1, ny1, error = 0;
267
268 int64_t ij;
269 uint64_t count[GRDLANDMASK_N_CLASSES];
270
271 size_t nx_alloc = GMT_SMALL_CHUNK;
272
273 char line[GMT_LEN256] = {""};
274 char *shore_resolution[5] = {"full", "high", "intermediate", "low", "crude"};
275
276 double xmin, xmax, ymin, ymax, west_border, east_border, i_dx_inch, i_dy_inch, inc_inch[2];
277 double dummy, *x = NULL, *y = NULL;
278
279 gmt_grdfloat f_level = 0.0f;
280
281 struct GMT_SHORE c;
282 struct GMT_GRID *Grid = NULL, *Cart = NULL;
283 struct GMT_GRID_HEADER *C = NULL;
284 struct GMT_GRID_HEADER_HIDDEN *HH = NULL;
285 struct GMT_GSHHS_POL *p = NULL;
286 struct GRDLANDMASK_BINCROSS *X = NULL;
287 struct GRDLANDMASK_CTRL *Ctrl = NULL;
288 struct GMT_CTRL *GMT = NULL, *GMT_cpy = NULL;
289 struct GMT_OPTION *options = NULL;
290 struct GMTAPI_CTRL *API = gmt_get_api_ptr (V_API); /* Cast from void to GMTAPI_CTRL pointer */
291
292 /*----------------------- Standard module initialization and parsing ----------------------*/
293
294 if (API == NULL) return (GMT_NOT_A_SESSION);
295 if (mode == GMT_MODULE_PURPOSE) return (usage (API, GMT_MODULE_PURPOSE)); /* Return the purpose of program */
296 options = GMT_Create_Options (API, mode, args); if (API->error) return (API->error); /* Set or get option list */
297
298 if ((error = gmt_report_usage (API, options, 0, usage)) != GMT_NOERROR) bailout (error); /* Give usage if requested */
299
300 /* Parse the command-line arguments */
301
302 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 */
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 grdlandmask main code ----------------------------*/
308
309 /* We know coastline data are geographic so we hardwire this here so parsing of -R works: */
310 gmt_set_geographic (GMT, GMT_IN);
311 gmt_enable_threads (GMT); /* Set number of active threads, if supported */
312
313 /* Create the empty grid and allocate space */
314 if ((Grid = GMT_Create_Data (API, GMT_IS_GRID, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA, NULL, NULL, NULL, \
315 GMT_GRID_DEFAULT_REG, GMT_NOTSET, NULL)) == NULL) Return (API->error);
316
317 if (Grid->header->wesn[XLO] < 0.0 && Grid->header->wesn[XHI] < 0.0) { /* Shift longitudes */
318 temp_shift = true;
319 Grid->header->wesn[XLO] += 360.0;
320 Grid->header->wesn[XHI] += 360.0;
321 }
322 HH = gmt_get_H_hidden (Grid->header);
323
324 if (Ctrl->D.force) Ctrl->D.set = gmt_shore_adjust_res (GMT, Ctrl->D.set, true);
325 base = gmt_set_resolution (GMT, &Ctrl->D.set, 'D');
326 gmt_M_memset (count, GRDLANDMASK_N_CLASSES, uint64_t); /* Counts of each level */
327
328 if (Ctrl->N.wetdry) { /* Must duplicate wet/dry settings */
329 Ctrl->N.mask[6] = Ctrl->N.mask[2];
330 Ctrl->N.mask[4] = Ctrl->N.mask[8] = Ctrl->N.mask[0];
331 }
332
333 if ((err = gmt_init_shore (GMT, Ctrl->D.set, &c, Grid->header->wesn, &Ctrl->A.info))) {
334 GMT_Report (API, GMT_MSG_ERROR, "%s [GSHHG %s resolution shorelines]\n", GMT_strerror(err), shore_resolution[base]);
335 Return (GMT_RUNTIME_ERROR);
336 }
337 if (gmt_M_is_verbose (GMT, GMT_MSG_INFORMATION)) {
338 GMT_Report (API, GMT_MSG_INFORMATION, "GSHHG version %s\n", c.version);
339 GMT_Report (API, GMT_MSG_INFORMATION, "%s\n", c.title);
340 GMT_Report (API, GMT_MSG_INFORMATION, "%s\n", c.source);
341
342 sprintf (line, "%s\n", GMT->current.setting.format_float_out);
343 if (Ctrl->N.wetdry) {
344 GMT_Report (API, GMT_MSG_INFORMATION, "Nodes in water will be set to ");
345 (gmt_M_is_fnan (Ctrl->N.mask[0])) ? GMT_Message (API, GMT_TIME_NONE, "NaN\n") : GMT_Message (API, GMT_TIME_NONE, line, Ctrl->N.mask[0]);
346 GMT_Report (API, GMT_MSG_INFORMATION, "Nodes on land will be set to ");
347 (gmt_M_is_fnan (Ctrl->N.mask[2])) ? GMT_Message (API, GMT_TIME_NONE, "NaN\n") : GMT_Message (API, GMT_TIME_NONE, line, Ctrl->N.mask[2]);
348 }
349 else {
350 GMT_Report (API, GMT_MSG_INFORMATION, "Nodes in the oceans will be set to ");
351 (gmt_M_is_fnan (Ctrl->N.mask[0])) ? GMT_Message (API, GMT_TIME_NONE, "NaN\n") : GMT_Message (API, GMT_TIME_NONE, line, Ctrl->N.mask[0]);
352 GMT_Report (API, GMT_MSG_INFORMATION, "Nodes on land will be set to ");
353 (gmt_M_is_fnan (Ctrl->N.mask[2])) ? GMT_Message (API, GMT_TIME_NONE, "NaN\n") : GMT_Message (API, GMT_TIME_NONE, line, Ctrl->N.mask[2]);
354 GMT_Report (API, GMT_MSG_INFORMATION, "Nodes in lakes will be set to ");
355 (gmt_M_is_fnan (Ctrl->N.mask[4])) ? GMT_Message (API, GMT_TIME_NONE, "NaN\n") : GMT_Message (API, GMT_TIME_NONE, line, Ctrl->N.mask[4]);
356 GMT_Report (API, GMT_MSG_INFORMATION, "Nodes in islands will be set to ");
357 (gmt_M_is_fnan (Ctrl->N.mask[6])) ? GMT_Message (API, GMT_TIME_NONE, "NaN\n") : GMT_Message (API, GMT_TIME_NONE, line, Ctrl->N.mask[6]);
358 GMT_Report (API, GMT_MSG_INFORMATION, "Nodes in ponds will be set to ");
359 (gmt_M_is_fnan (Ctrl->N.mask[8])) ? GMT_Message (API, GMT_TIME_NONE, "NaN\n") : GMT_Message (API, GMT_TIME_NONE, line, Ctrl->N.mask[8]);
360 }
361 if (Ctrl->E.linetrace) {
362 GMT_Report (API, GMT_MSG_INFORMATION, "Nodes near shoreline will be set to ");
363 (gmt_M_is_fnan (Ctrl->N.mask[1])) ? GMT_Message (API, GMT_TIME_NONE, "NaN\n") : GMT_Message (API, GMT_TIME_NONE, line, Ctrl->N.mask[1]);
364 GMT_Report (API, GMT_MSG_INFORMATION, "Nodes near lakeline will be set to ");
365 (gmt_M_is_fnan (Ctrl->N.mask[3])) ? GMT_Message (API, GMT_TIME_NONE, "NaN\n") : GMT_Message (API, GMT_TIME_NONE, line, Ctrl->N.mask[3]);
366 GMT_Report (API, GMT_MSG_INFORMATION, "Nodes near islandline will be set to ");
367 (gmt_M_is_fnan (Ctrl->N.mask[5])) ? GMT_Message (API, GMT_TIME_NONE, "NaN\n") : GMT_Message (API, GMT_TIME_NONE, line, Ctrl->N.mask[5]);
368 GMT_Report (API, GMT_MSG_INFORMATION, "Nodes near pondline will be set to ");
369 (gmt_M_is_fnan (Ctrl->N.mask[7])) ? GMT_Message (API, GMT_TIME_NONE, "NaN\n") : GMT_Message (API, GMT_TIME_NONE, line, Ctrl->N.mask[7]);
370 }
371 }
372
373 GMT_Report (API, GMT_MSG_INFORMATION, "Apply linear scale to avoid wrap-around: -Jx100id; this temporarily changes the domain\n");
374 gmt_parse_common_options (GMT, "J", 'J', "x100id"); /* Fake linear projection so the shore machinery will work. Used bo be -Jx1id but had trouble with roundoff */
375 if (gmt_M_err_pass (GMT, gmt_proj_setup (GMT, Grid->header->wesn), "")) Return (GMT_PROJECTION_ERROR);
376 GMT->current.map.parallel_straight = GMT->current.map.meridian_straight = 2; /* No resampling along bin boundaries */
377 wrap = GMT->current.map.is_world = gmt_grd_is_global (GMT, Grid->header);
378 double_dip = (wrap && Grid->header->registration == GMT_GRID_NODE_REG); /* Must duplicate the west nodes to east */
379 /* Using -Jx1d means output is Cartesian but we want to force geographic */
380 gmt_set_geographic (GMT, GMT_OUT);
381 /* All data nodes are thus initialized to 0 */
382 x = gmt_M_memory (GMT, NULL, Grid->header->n_columns, double);
383 y = gmt_M_memory (GMT, NULL, Grid->header->n_rows, double);
384 if (Ctrl->E.linetrace)
385 X = gmt_M_memory (GMT, X, nx_alloc, struct GRDLANDMASK_BINCROSS);
386
387 nx1 = Grid->header->n_columns - 1; ny1 = Grid->header->n_rows - 1;
388
389 /* Fill out gridnode coordinates and apply the implicit linear projection */
390
391 for (col = 0; col <= nx1; col++) gmt_geo_to_xy (GMT, gmt_M_grd_col_to_x (GMT, col, Grid->header), 0.0, &x[col], &dummy);
392 for (row = 0; row <= ny1; row++) gmt_geo_to_xy (GMT, 0.0, gmt_M_grd_row_to_y (GMT, row, Grid->header), &dummy, &y[row]);
393 inc_inch[GMT_X] = fabs (x[1] - x[0]);
394 inc_inch[GMT_Y] = fabs (y[1] - y[0]);
395 i_dx_inch = 1.0 / inc_inch[GMT_X];
396 i_dy_inch = 1.0 / inc_inch[GMT_Y];
397
398 if ((Cart = GMT_Create_Data (API, GMT_IS_GRID, GMT_IS_SURFACE, GMT_CONTAINER_ONLY, NULL, GMT->current.proj.rect, inc_inch, \
399 Grid->header->registration, GMT_NOTSET, NULL)) == NULL) Return (API->error);
400 C = Cart->header;
401
402 west_border = floor (GMT->common.R.wesn[XLO] / c.bsize) * c.bsize;
403 east_border = ceil (GMT->common.R.wesn[XHI] / c.bsize) * c.bsize;
404
405 for (ind = 0; ind < c.nb; ind++) { /* Loop over necessary bins only */
406
407 bin = c.bins[ind];
408 GMT_Report (API, GMT_MSG_DEBUG, "Working on block # %5ld\n", bin);
409
410 if ((err = gmt_get_shore_bin (GMT, ind, &c))) {
411 GMT_Report (API, GMT_MSG_ERROR, "%s [%s resolution shoreline]\n", GMT_strerror(err), shore_resolution[base]);
412 Return (GMT_RUNTIME_ERROR);
413 }
414
415 /* Use polygons, if any. Go in both directions to cover both land and sea */
416
417 used_polygons = false;
418
419 for (direction = -1; c.ns > 0 && direction < 2; direction += 2) {
420
421 /* Assemble one or more segments into polygons */
422
423 np = gmt_assemble_shore (GMT, &c, direction, true, west_border, east_border, &p);
424
425 /* Get clipped polygons in x,y inches that can be processed */
426
427 np_new = gmt_prep_shore_polygons (GMT, &p, np, false, 0.0, -1);
428
429 for (k = 0; k < np_new; k++) {
430
431 if (p[k].n == 0) continue;
432
433 used_polygons = true; /* At least some points made it to here */
434
435 /* Find min/max of polygon in inches */
436
437 xmin = xmax = p[k].lon[0];
438 ymin = ymax = p[k].lat[0];
439 for (i = 1; i < p[k].n; i++) {
440 if (p[k].lon[i] < xmin) xmin = p[k].lon[i];
441 if (p[k].lon[i] > xmax) xmax = p[k].lon[i];
442 if (p[k].lat[i] < ymin) ymin = p[k].lat[i];
443 if (p[k].lat[i] > ymax) ymax = p[k].lat[i];
444 }
445 col_min = MAX (0, irint (ceil (xmin * i_dx_inch - Grid->header->xy_off - GMT_CONV8_LIMIT)));
446 if (col_min > nx1) col_min = 0;
447 /* So col_min is in range [0,nx1] */
448 col_max = MIN (nx1, irint (floor (xmax * i_dx_inch - Grid->header->xy_off + GMT_CONV8_LIMIT)));
449 if (col_max <= 0 || col_max < col_min) col_max = nx1;
450 /* So col_max is in range [1,nx1] */
451 row_min = MAX (0, irint (ceil ((GMT->current.proj.rect[YHI] - ymax) * i_dy_inch - Grid->header->xy_off - GMT_CONV8_LIMIT)));
452 /* So row_min is in range [0,?] */
453 row_max = MIN (ny1, irint (floor ((GMT->current.proj.rect[YHI] - ymin) * i_dy_inch - Grid->header->xy_off + GMT_CONV8_LIMIT)));
454 /* So row_max is in range [?,ny1] */
455 f_level = (gmt_grdfloat)(2*p[k].level);
456
457 #ifdef _OPENMP
458 #pragma omp parallel for private(row,col,side,ij) shared(row_min,row_max,col_min,col_max,GMT,x,y,p,k,Ctrl,Grid,f_level)
459 #endif
460 for (row = row_min; row <= row_max; row++) {
461 assert (row >= 0); /* Just in case we have a logic bug somewhere */
462 for (col = col_min; col <= col_max; col++) {
463
464 if ((side = gmt_non_zero_winding (GMT, x[col], y[row], p[k].lon, p[k].lat, p[k].n)) < Ctrl->E.inside) continue; /* Outside */
465
466 /* Here, point is inside, we must assign value */
467
468 ij = gmt_M_ijp (Grid->header, row, col);
469 if (f_level > Grid->data[ij]) Grid->data[ij] = f_level;
470 }
471 }
472 }
473
474 gmt_free_shore_polygons (GMT, p, np_new);
475 gmt_M_free (GMT, p);
476
477 if (Ctrl->E.linetrace) { /* Deal with tracing the lines and resetting gridnodes */
478 int64_t last_ij;
479 uint64_t nx;
480 unsigned int curr_x_pt, prev_x_pt, pt;
481 int last_col, last_row, start_col, end_col, brow, bcol;
482 double dx, del_x, del_y, xc, yc, xb, yb;
483 bool last_not_set;
484
485 if ((np = gmt_assemble_shore (GMT, &c, 1, false, west_border, east_border, &p)) == 0) { /* Just get segments */
486 gmt_free_shore (GMT, &c);
487 continue;
488 }
489 for (k = 0; k < np; k++) {
490 if (p[k].n == 0) continue;
491 for (pt = 0; pt < (unsigned int)p[k].n; pt++) { /* Convert geo coordinates to inches */
492 gmt_geo_to_xy (GMT, p[k].lon[pt], p[k].lat[pt], &xc, &yc);
493 p[k].lon[pt] = xc; p[k].lat[pt] = yc;
494 }
495 f_level = (gmt_grdfloat)(2*p[k].level-1); /* Border index is half-way between the two levels */
496 last_ij = UINT_MAX;
497 last_not_set = true;
498 last_col = last_row = -1;
499
500 /* To handle lines that exit the grid we pursue the entire line even if outside.
501 * We only check if (row,col) is inside when filling in between points and assigning nodes */
502
503 for (pt = 0; pt < (unsigned int)p[k].n; pt++) {
504 /* Get (row,col) and index to nearest node for this point */
505 row = gmt_M_grd_y_to_row (GMT, p[k].lat[pt], C);
506 col = gmt_M_grd_x_to_col (GMT, p[k].lon[pt], C);
507 ij = gmt_M_ijp (Grid->header, row, col);
508
509 if (last_not_set) { /* Initialize last node to this node the first time */
510 last_ij = ij;
511 last_not_set = false;
512 }
513
514 if (ij != last_ij) { /* Crossed into another cell */
515 start_col = MIN (last_col, col) + 1;
516 end_col = MAX (last_col, col);
517 dx = p[k].lon[pt] - p[k].lon[pt-1];
518
519 /* Find all the bin-line intersections - modified from x2sys_binlist */
520
521 /* Add the start and stop coordinates to the xc/yc arrays so we can get mid-points of each interval */
522
523 nx = 0;
524 if (grdlandmask_inside (GMT, C, last_row, last_col)) { /* Previous point is inside, add it to our list with zero distance */
525 X[nx].x = p[k].lon[pt-1]; X[nx].y = p[k].lat[pt-1]; X[nx++].d = 0.0;
526 }
527 if (grdlandmask_inside (GMT, C, row, col)) { /* This point is inside, add it to our list with max distance */
528 X[nx].x = p[k].lon[pt]; X[nx].y = p[k].lat[pt]; X[nx++].d = hypot (dx, p[k].lat[pt] - p[k].lat[pt-1]);
529 }
530
531 /* Now add all crossings between this line segment and the gridlines outlining the cells centered on the nodes */
532
533 for (brow = MIN (last_row, row) + 1; brow <= MAX (last_row, row); brow++) { /* If we go in here we know dy is non-zero */
534 if (brow < 0 || brow > (int)C->n_rows) continue; /* Outside grid */
535 /* Determine intersections between the line segment and parallels */
536 yb = gmt_M_grd_row_to_y (GMT, brow, C) + 0.5 * inc_inch[GMT_Y];
537 del_y = yb - p[k].lat[pt-1];
538 del_x = del_y * dx / (p[k].lat[pt] - p[k].lat[pt-1]);
539 xc = p[k].lon[pt-1] + del_x;
540 X[nx].x = xc; X[nx].y = yb; X[nx].d = hypot (del_x , del_y);
541 nx++;
542 if (nx == nx_alloc) {
543 nx_alloc <<= 1;
544 X = gmt_M_memory (GMT, X, nx_alloc, struct GRDLANDMASK_BINCROSS);
545 }
546 }
547 for (bcol = start_col; bcol <= end_col; bcol++) { /* If we go in here we think dx is non-zero (we do a last-ditch dx check just in case) */
548 if (bcol < 0 || bcol > (int)C->n_columns) continue;
549 /* Determine intersections between the line segment and meridians */
550 xb = gmt_M_grd_col_to_x (GMT, bcol, C) - 0.5 * inc_inch[GMT_X];
551 del_x = xb - p[k].lon[pt-1];
552 del_y = (dx == 0.0) ? 0.5 * (p[k].lat[pt] - p[k].lat[pt-1]) : del_x * (p[k].lat[pt] - p[k].lat[pt-1]) / dx;
553 yc = p[k].lat[pt-1] + del_y;
554 X[nx].x = xb; X[nx].y = yc; X[nx].d = hypot (del_x, del_y);
555 nx++;
556 if (nx == nx_alloc) {
557 nx_alloc <<= 1;
558 X = gmt_M_memory (GMT, X, nx_alloc, struct GRDLANDMASK_BINCROSS);
559 }
560 }
561
562 if (nx) { /* Here we have intersections */
563 qsort (X, nx, sizeof (struct GRDLANDMASK_BINCROSS), grdlandmask_comp_bincross); /* Sort on distances along the line segment */
564 grdlandmask_assign_node (GMT, Grid, C, X[0].x, X[0].y, f_level, &ij); /* Possibly assign f_level to nearest node */
565 for (curr_x_pt = 1, prev_x_pt = 0; curr_x_pt < nx; curr_x_pt++, prev_x_pt++) { /* Process the intervals, getting mid-points and using that to get bin */
566 grdlandmask_assign_node (GMT, Grid, C, X[curr_x_pt].x, X[curr_x_pt].y, f_level, &ij); /* Possibly assign f_level to nearest node */
567 /* Add mid-point between the cuts */
568 xc = 0.5 * (X[curr_x_pt].x + X[prev_x_pt].x);
569 yc = 0.5 * (X[curr_x_pt].y + X[prev_x_pt].y);
570 grdlandmask_assign_node (GMT, Grid, C, xc, yc, f_level, &ij); /* Possibly assign f_level to nearest node */
571 }
572 }
573 }
574 /* Update last row/col and index */
575 last_ij = ij;
576 last_col = col;
577 last_row = row;
578 }
579
580 }
581 gmt_free_shore_polygons (GMT, p, np);
582 gmt_M_free (GMT, p);
583 }
584 }
585
586 if (!used_polygons) { /* Lack of polygons or clipping etc resulted in no polygons after all, must deal with background */
587 k = INT_MAX; /* Initialize to outside range of levels (4 is highest) */
588 /* Visit each of the 4 nodes, test if it is inside -R, and if so update lowest level found so far */
589
590 if (!gmt_map_outside (GMT, c.lon_sw, c.lat_sw)) k = MIN (k, c.node_level[0]); /* SW */
591 if (!gmt_map_outside (GMT, c.lon_sw + c.bsize, c.lat_sw)) k = MIN (k, c.node_level[1]); /* SE */
592 if (!gmt_map_outside (GMT, c.lon_sw + c.bsize, c.lat_sw - c.bsize)) k = MIN (k, c.node_level[2]); /* NE */
593 if (!gmt_map_outside (GMT, c.lon_sw, c.lat_sw - c.bsize)) k = MIN (k, c.node_level[3]); /* NW */
594
595 /* If k is still INT_MAX we must assume this patch should have the min level of the bin */
596
597 if (k == INT_MAX) k = MIN (MIN (c.node_level[0], c.node_level[1]) , MIN (c.node_level[2], c.node_level[3]));
598 f_level = (gmt_grdfloat)(2*k);
599
600 /* Determine nodes to initialize */
601
602 row_min = MAX (0, irint (ceil ((Grid->header->wesn[YHI] - c.lat_sw - c.bsize) * HH->r_inc[GMT_Y] - Grid->header->xy_off)));
603 row_max = MIN (ny1, irint (floor ((Grid->header->wesn[YHI] - c.lat_sw) * HH->r_inc[GMT_Y] - Grid->header->xy_off)));
604 if (wrap) { /* Handle jumps */
605 col_min = irint (ceil (fmod (c.lon_sw - Grid->header->wesn[XLO], 360.0) * HH->r_inc[GMT_X] - Grid->header->xy_off));
606 col_max = irint (floor (fmod (c.lon_sw + c.bsize - Grid->header->wesn[XLO], 360.0) * HH->r_inc[GMT_X] - Grid->header->xy_off));
607 if (col_max < col_min) col_max += Grid->header->n_columns;
608 }
609 else { /* Make sure we are inside our grid */
610 double lon_w, lon_e;
611 lon_w = c.lon_sw - Grid->header->wesn[XLO]; lon_e = c.lon_sw + c.bsize - Grid->header->wesn[XLO];
612 if (lon_w < Grid->header->wesn[XLO] && (lon_w+360.0) < Grid->header->wesn[XHI]) {
613 lon_w += 360.0; lon_e += 360.0;
614 }
615 else if (lon_e > Grid->header->wesn[XHI] && (lon_e-360.0) > Grid->header->wesn[XLO]) {
616 lon_w -= 360.0; lon_e -= 360.0;
617 }
618 col_min = irint (ceil (lon_w * HH->r_inc[GMT_X] - Grid->header->xy_off));
619 col_max = irint (floor (lon_e * HH->r_inc[GMT_X] - Grid->header->xy_off));
620 if (col_min < 0) col_min = 0;
621 if (col_max > nx1) col_max = nx1;
622 }
623 #ifdef _OPENMP
624 #pragma omp parallel for private(row,col,ii,ij) shared(row_min,row_max,col_min,col_max,wrap,nx1,Grid,f_level)
625 #endif
626 for (row = row_min; row <= row_max; row++) {
627 for (col = col_min; col <= col_max; col++) {
628 ii = (wrap) ? col % (int)Grid->header->n_columns : col;
629 if (ii < 0 || ii > nx1) continue;
630 ij = gmt_M_ijp (Grid->header, row, ii);
631 Grid->data[ij] = f_level;
632 }
633 }
634 }
635
636 gmt_free_shore (GMT, &c);
637 }
638
639 gmt_shore_cleanup (GMT, &c);
640 gmt_M_free (GMT, x);
641 gmt_M_free (GMT, y);
642 if (Ctrl->E.linetrace) gmt_M_free (GMT, X);
643 if (GMT_Destroy_Data (API, &Cart) != GMT_NOERROR) {
644 Return (API->error);
645 }
646
647 #ifdef _OPENMP
648 #pragma omp parallel for private(row,col,k,ij) shared(GMT,Grid,Ctrl)
649 #endif
650
651 gmt_M_grd_loop (GMT, Grid, row, col, ij) { /* Turn levels into mask values */
652 k = urint (Grid->data[ij]);
653 Grid->data[ij] = Ctrl->N.mask[k];
654 count[k]++;
655 if (col == 0 && double_dip) count[k]++; /* Count these guys twice */
656 }
657
658 if (double_dip) { /* Copy over values to the repeating right column */
659 unsigned int row_l;
660 for (row_l = 0, ij = gmt_M_ijp (Grid->header, row_l, 0); row_l < Grid->header->n_rows; row_l++, ij += Grid->header->mx) Grid->data[ij+nx1] = Grid->data[ij];
661 }
662
663 if (temp_shift) {
664 Grid->header->wesn[XLO] -= 360.0;
665 Grid->header->wesn[XHI] -= 360.0;
666 }
667
668 sprintf (line, "Derived from the %s resolution shorelines", shore_resolution[base]);
669 if (GMT_Set_Comment (API, GMT_IS_GRID, GMT_COMMENT_IS_REMARK, line, Grid)) return (API->error);
670 if (GMT_Set_Comment (API, GMT_IS_GRID, GMT_COMMENT_IS_OPTION | GMT_COMMENT_IS_COMMAND, options, Grid)) Return (API->error);
671 if (GMT_Write_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA, NULL, Ctrl->G.file, Grid) != GMT_NOERROR) {
672 Return (API->error);
673 }
674
675 if (gmt_M_is_verbose (GMT, GMT_MSG_WARNING)) {
676 for (k = 0; k < GRDLANDMASK_N_CLASSES; k++) {
677 if (count[k] == 0) continue;
678 if (k%2 == 0)
679 GMT_Report (API, GMT_MSG_INFORMATION, "Level %d set for %" PRIu64 " nodes\n", k/2, count[k]);
680 else
681 GMT_Report (API, GMT_MSG_INFORMATION, "Border between Levels %d-%d set for %" PRIu64 " nodes\n", (k-1)/2, (k+1)/2, count[k]);
682 }
683 }
684
685 Return (GMT_NOERROR);
686 }
687