1 /*--------------------------------------------------------------------
2 *
3 * Copyright (c) 2008-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 * Spherical nearest distances - via Voronoi polygons. We read input
19 * data, assumed to be things like coastlines, and want to create a grid
20 * with distances to the nearest line. The approach here is to break
21 * the data into voronoi polygons and then visit all nodes inside each
22 * polygon and use geodesic distance calculation from each node to the
23 * unique Voronoi interior data node.
24 * Relies on STRIPACK Fortran F77 library (Renka, 1997). Reference:
25 * Renka, R, J,, 1997, Algorithm 772: STRIPACK: Delaunay Triangulation
26 * and Voronoi Diagram on the Surface of a Sphere, AMC Trans. Math.
27 * Software, 23 (3), 416-434.
28 * We translated to C using f2c -r8 and and manually edited the code
29 * so that f2c libs were not needed. For any translation errors, blame me.
30 *
31 * Author: Paul Wessel
32 * Date: 1-AUG-2011
33 * Version: 6 API
34 *
35 */
36
37 #include "gmt_dev.h"
38 #include "gmt_sph.h"
39
40 #define THIS_MODULE_CLASSIC_NAME "sphdistance"
41 #define THIS_MODULE_MODERN_NAME "sphdistance"
42 #define THIS_MODULE_LIB "core"
43 #define THIS_MODULE_PURPOSE "Create Voronoi distance, node, or natural nearest-neighbor grid on a sphere"
44 #define THIS_MODULE_KEYS "<D{,ND(,QD(,GG},Q-("
45 #define THIS_MODULE_NEEDS "R"
46 #define THIS_MODULE_OPTIONS "-:RVbdehijqr" GMT_OPT("F")
47
48 enum sphdist_modes {
49 SPHD_DIST = 0,
50 SPHD_NODES = 1,
51 SPHD_VALUES = 2};
52
53 struct SPHDISTANCE_CTRL {
54 struct SPHDISTANCE_A { /* -A[m|p|x|y|step] */
55 bool active;
56 unsigned int mode;
57 double step;
58 } A;
59 struct SPHDISTANCE_C { /* -C */
60 bool active;
61 } C;
62 struct SPHDISTANCE_D { /* -D for variable tension */
63 bool active;
64 } D;
65 struct SPHDISTANCE_E { /* -Ed|n|z[<dist>] */
66 bool active;
67 unsigned int mode;
68 double dist;
69 } E;
70 struct SPHDISTANCE_G { /* -G<maskfile> */
71 bool active;
72 char *file;
73 } G;
74 struct SPHDISTANCE_I { /* -I (for checking only) */
75 bool active;
76 } I;
77 struct SPHDISTANCE_L { /* -L<unit>] */
78 bool active;
79 char unit;
80 } L;
81 struct SPHDISTANCE_N { /* -N */
82 bool active;
83 char *file;
84 } N;
85 struct SPHDISTANCE_Q { /* -Q */
86 bool active;
87 char *file;
88 } Q;
89 };
90
sphdistance_prepare_polygon(struct GMT_CTRL * GMT,struct GMT_DATASEGMENT * P)91 GMT_LOCAL void sphdistance_prepare_polygon (struct GMT_CTRL *GMT, struct GMT_DATASEGMENT *P) {
92 /* Set the min/max extent of this polygon and determine if it
93 * is a polar cap; if so set the required metadata flags */
94 uint64_t row;
95 double lon_sum = 0.0, lat_sum = 0.0, dlon;
96 struct GMT_DATASEGMENT_HIDDEN *PH = gmt_get_DS_hidden (P);
97
98 gmt_set_seg_minmax (GMT, GMT_IS_POLY, 0, P); /* Set the domain of the segment */
99
100 /* Then loop over points to accumulate sums */
101
102 for (row = 1; row < P->n_rows; row++) { /* Start at row = 1 since (a) 0'th point is repeated at end and (b) we are doing differences */
103 gmt_M_set_delta_lon (P->data[GMT_X][row-1], P->data[GMT_X][row], dlon);
104 lon_sum += dlon;
105 lat_sum += P->data[GMT_Y][row];
106 }
107 PH->pole = 0;
108 if (gmt_M_360_range (lon_sum, 0.0)) { /* Contains a pole */
109 if (lat_sum < 0.0) { /* S */
110 PH->pole = -1;
111 PH->lat_limit = P->min[GMT_Y];
112 P->min[GMT_Y] = -90.0;
113
114 }
115 else { /* N */
116 PH->pole = +1;
117 PH->lat_limit = P->max[GMT_Y];
118 P->max[GMT_Y] = 90.0;
119 }
120 P->min[GMT_X] = 0.0; P->max[GMT_X] = 360.0;
121 }
122 }
123
New_Ctrl(struct GMT_CTRL * GMT)124 static void *New_Ctrl (struct GMT_CTRL *GMT) { /* Allocate and initialize a new control structure */
125 struct SPHDISTANCE_CTRL *C;
126
127 C = gmt_M_memory (GMT, NULL, 1, struct SPHDISTANCE_CTRL);
128 C->E.dist = 1.0; /* Default is 1 degree Voronoi edge resampling */
129 C->L.unit = 'e'; /* Default is meter distances */
130 return (C);
131 }
132
Free_Ctrl(struct GMT_CTRL * GMT,struct SPHDISTANCE_CTRL * C)133 static void Free_Ctrl (struct GMT_CTRL *GMT, struct SPHDISTANCE_CTRL *C) { /* Deallocate control structure */
134 if (!C) return;
135 gmt_M_str_free (C->G.file);
136 gmt_M_str_free (C->N.file);
137 gmt_M_str_free (C->Q.file);
138 gmt_M_free (GMT, C);
139 }
140
usage(struct GMTAPI_CTRL * API,int level)141 static int usage (struct GMTAPI_CTRL *API, int level) {
142 const char *name = gmt_show_name_and_purpose (API, THIS_MODULE_LIB, THIS_MODULE_CLASSIC_NAME, THIS_MODULE_PURPOSE);
143 if (level == GMT_MODULE_PURPOSE) return (GMT_NOERROR);
144 GMT_Usage (API, 1, "==> The hard work is done by algorithms 772 (STRIPACK) & 773 (SSRFPACK) by R. J. Renka [1997] <==\n");
145 GMT_Usage (API, 0, "usage: %s [<table>] -G%s %s [-C] [-D] [-En|z|d[<dr>]] "
146 "[-L<unit>] [-N<nodetable>] [-Q<voronoitable>] [%s] [%s] [%s] [%s] "
147 "[%s] [%s] [%s] [%s] [%s] [%s] [%s] [%s]\n",
148 name, GMT_OUTGRID, GMT_I_OPT, GMT_Rgeo_OPT, GMT_V_OPT, GMT_bi_OPT, GMT_di_OPT, GMT_e_OPT, GMT_h_OPT,
149 GMT_i_OPT, GMT_j_OPT, GMT_qi_OPT, GMT_r_OPT, GMT_colon_OPT, GMT_PAR_OPT);
150
151 if (level == GMT_SYNOPSIS) return (GMT_MODULE_SYNOPSIS);
152
153 GMT_Message (API, GMT_TIME_NONE, " REQUIRED ARGUMENTS:\n");
154 GMT_Option (API, "<");
155 gmt_outgrid_syntax (API, 'G', "Specify file name for output distance grid file");
156 GMT_Option (API, "I");
157
158 GMT_Message (API, GMT_TIME_NONE, "\n OPTIONAL ARGUMENTS:\n");
159 GMT_Usage (API, 1, "\n-C Conserve memory (Converts lon/lat <--> x/y/z when needed) [store both in memory]. Not used with -Q.");
160 GMT_Usage (API, 1, "\n-D Delete any duplicate points [Default assumes there are no duplicates].");
161 GMT_Usage (API, 1, "\n-En|z|d[<dr>]");
162 GMT_Usage (API, -2, "Specify the quantity that should be assigned to the grid nodes:");
163 GMT_Usage (API, 3, "n: The Voronoi polygon ID.");
164 GMT_Usage (API, 3, "z: The z-value of the Voronoi center node (natural nearest-neighbor gridding).");
165 GMT_Usage (API, 3, "d: The distance to the nearest data point [Default].");
166 GMT_Usage (API, -2, "Optionally append resampling interval in spherical degrees for polygon arcs [1].");
167 GMT_Usage (API, 1, "\n-L<unit>");
168 GMT_Usage (API, -2, "Set distance unit arc (d)egree, m(e)ter, (f)eet, (k)m, arc (m)inute, (M)ile, (n)autical mile, or arc (s)econd [e].");
169 GMT_Usage (API, 1, "\n-N<nodetable>");
170 GMT_Usage (API, -2, "Specify node filename for the Voronoi polygons (i.e., sphtriangulate -N output or equivalent).");
171 GMT_Usage (API, 1, "\n-Q<voronoitable>");
172 GMT_Usage (API, -2, "Specify table with Voronoi polygons in sphtriangulate -Qv format "
173 "[Default performs Voronoi construction on input data first].");
174 GMT_Option (API, "Rg");
175 if (gmt_M_showusage (API)) GMT_Usage (API, -2, "If no region is specified we default to the entire world [-Rg].");
176 GMT_Option (API, "V,bi2,di,e,h,i,j,qi,r,:,.");
177
178 return (GMT_MODULE_USAGE);
179 }
180
parse(struct GMT_CTRL * GMT,struct SPHDISTANCE_CTRL * Ctrl,struct GMT_OPTION * options)181 static int parse (struct GMT_CTRL *GMT, struct SPHDISTANCE_CTRL *Ctrl, struct GMT_OPTION *options) {
182 /* This parses the options provided to sphdistance and sets parameters in CTRL.
183 * Any GMT common options will override values set previously by other commands.
184 * It also replaces any file names specified as input or output with the data ID
185 * returned when registering these sources/destinations with the API.
186 */
187
188 unsigned int n_errors = 0, k;
189 struct GMT_OPTION *opt = NULL;
190 struct GMTAPI_CTRL *API = GMT->parent;
191
192 for (opt = options; opt; opt = opt->next) {
193 switch (opt->option) {
194
195 case '<': /* Skip input files */
196 if (GMT_Get_FilePath (API, GMT_IS_DATASET, GMT_IN, GMT_FILE_REMOTE, &(opt->arg))) n_errors++;;
197 break;
198
199 /* Processes program-specific parameters */
200
201 case 'C':
202 n_errors += gmt_M_repeated_module_option (API, Ctrl->C.active);
203 Ctrl->C.active = true;
204 break;
205 case 'D':
206 n_errors += gmt_M_repeated_module_option (API, Ctrl->D.active);
207 Ctrl->D.active = true;
208 break;
209 case 'E':
210 n_errors += gmt_M_repeated_module_option (API, Ctrl->E.active);
211 Ctrl->E.active = true;
212 switch (opt->arg[0]) { /* Select output grid mode */
213 case 'n': Ctrl->E.mode = SPHD_NODES; k = 1; break;
214 case 'z': Ctrl->E.mode = SPHD_VALUES; k = 1; break;
215 case 'd': Ctrl->E.mode = SPHD_DIST; k = 1; break;
216 default: Ctrl->E.mode = SPHD_DIST; k = 0; break;
217 }
218 if (opt->arg[k]) Ctrl->E.dist = atof (&opt->arg[k]);
219 break;
220 case 'G':
221 n_errors += gmt_M_repeated_module_option (API, Ctrl->G.active);
222 Ctrl->G.active = true;
223 if (opt->arg[0]) Ctrl->G.file = strdup (opt->arg);
224 if (GMT_Get_FilePath (API, GMT_IS_GRID, GMT_OUT, GMT_FILE_LOCAL, &(Ctrl->G.file))) n_errors++;
225 break;
226 case 'I':
227 n_errors += gmt_M_repeated_module_option (API, Ctrl->I.active);
228 Ctrl->I.active = true;
229 n_errors += gmt_parse_inc_option (GMT, 'I', opt->arg);
230 break;
231 case 'L':
232 n_errors += gmt_M_repeated_module_option (API, Ctrl->L.active);
233 Ctrl->L.active = true;
234 if (!(opt->arg && strchr (GMT_LEN_UNITS, opt->arg[0]))) {
235 GMT_Report (API, GMT_MSG_ERROR, "Expected -L%s\n", GMT_LEN_UNITS_DISPLAY);
236 n_errors++;
237 }
238 else
239 Ctrl->L.unit = opt->arg[0];
240 break;
241 case 'N':
242 n_errors += gmt_M_repeated_module_option (API, Ctrl->N.active);
243 Ctrl->N.active = true;
244 Ctrl->N.file = strdup (opt->arg);
245 break;
246 case 'Q':
247 n_errors += gmt_M_repeated_module_option (API, Ctrl->Q.active);
248 Ctrl->Q.active = true;
249 Ctrl->Q.file = strdup (opt->arg);
250 break;
251 default: /* Report bad options */
252 n_errors += gmt_default_error (GMT, opt->option);
253 break;
254 }
255 }
256
257 if (GMT->common.b.active[GMT_IN] && GMT->common.b.ncol[GMT_IN] == 0) GMT->common.b.ncol[GMT_IN] = 3;
258 n_errors += gmt_M_check_condition (GMT, GMT->common.b.active[GMT_IN] && GMT->common.b.ncol[GMT_IN] < 3, "Binary input data (-bi) must have at least 3 columns\n");
259 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");
260 n_errors += gmt_M_check_condition (GMT, Ctrl->Q.active && GMT->common.b.active[GMT_IN] && !Ctrl->N.active, "Binary input data (-bi) with -Q also requires -N.\n");
261 n_errors += gmt_M_check_condition (GMT, !Ctrl->G.file, "Option -G: Must specify output file\n");
262
263 return (n_errors ? GMT_PARSE_ERROR : GMT_NOERROR);
264 }
265
266 #define bailout(code) {gmt_M_free_options (mode); return (code);}
267 #define Return(code) {Free_Ctrl (GMT, Ctrl); gmt_end_module (GMT, GMT_cpy); bailout (code);}
268
GMT_sphdistance(void * V_API,int mode,void * args)269 EXTERN_MSC int GMT_sphdistance (void *V_API, int mode, void *args) {
270 bool first = false, periodic, duplicate_col, duplicate = false;
271 int error = 0, s_row, south_row, north_row, w_col, e_col;
272
273 unsigned int row, col, p_col, west_col, east_col, nx1, n_in = 0;
274 uint64_t n_dup = 0, n_set = 0, side, ij, node, n_new, n = 0;
275 uint64_t vertex, node_stop, node_new, vertex_new, node_last, vertex_last;
276
277 size_t n_alloc, p_alloc = 0;
278
279 double first_x = 0.0, first_y = 0.0, prev_x = 0.0, prev_y = 0.0, X[3];
280 double *grid_lon = NULL, *grid_lat = NULL, *in = NULL;
281 double *xx = NULL, *yy = NULL, *zz = NULL, *lon = NULL, *lat = NULL;
282
283 gmt_grdfloat f_val = 0.0, *z_val = NULL;
284
285 struct GMT_GRID *Grid = NULL;
286 struct GMT_RECORD *In = NULL;
287 struct SPHDISTANCE_CTRL *Ctrl = NULL;
288 struct STRIPACK T;
289 struct GMT_DATASEGMENT *P = NULL;
290 struct GMT_DATASET *Qin = NULL;
291 struct GMT_DATATABLE *Table = NULL;
292 struct STRIPACK_VORONOI *V = NULL;
293 struct GMT_CTRL *GMT = NULL, *GMT_cpy = NULL;
294 struct GMT_OPTION *options = NULL;
295 struct GMTAPI_CTRL *API = gmt_get_api_ptr (V_API); /* Cast from void to GMTAPI_CTRL pointer */
296
297 /*----------------------- Standard module initialization and parsing ----------------------*/
298
299 if (API == NULL) return (GMT_NOT_A_SESSION);
300 if (mode == GMT_MODULE_PURPOSE) return (usage (API, GMT_MODULE_PURPOSE)); /* Return the purpose of program */
301 options = GMT_Create_Options (API, mode, args); if (API->error) return (API->error); /* Set or get option list */
302
303 if ((error = gmt_report_usage (API, options, 0, usage)) != GMT_NOERROR) bailout (error); /* Give usage if requested */
304
305 /* Parse the command-line arguments */
306
307 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 */
308 gmt_parse_common_options (GMT, "f", 'f', "g"); /* Implicitly set -fg since this is spherical triangulation */
309 if (GMT_Parse_Common (API, THIS_MODULE_OPTIONS, options)) Return (API->error);
310 Ctrl = New_Ctrl (GMT); /* Allocate and initialize a new control structure */
311 if ((error = parse (GMT, Ctrl, options)) != 0) Return (error);
312
313 /*---------------------------- This is the sphdistance main code ----------------------------*/
314
315 gmt_M_memset (&T, 1, struct STRIPACK);
316
317 if (gmt_init_distaz (GMT, Ctrl->L.unit, gmt_M_sph_mode (GMT), GMT_MAP_DIST) == GMT_NOT_A_VALID_TYPE)
318 Return (GMT_NOT_A_VALID_TYPE);
319
320 if (!GMT->common.R.active[RSET]) { /* Default to a global grid */
321 GMT->common.R.wesn[XLO] = 0.0; GMT->common.R.wesn[XHI] = 360.0; GMT->common.R.wesn[YLO] = -90.0; GMT->common.R.wesn[YHI] = 90.0;
322 }
323
324 /* Now we are ready to take on some input values */
325
326 if (Ctrl->Q.active) { /* Expect a single file with Voronoi polygons */
327 GMT_Report (API, GMT_MSG_INFORMATION, "Read Volonoi polygons from %s ...", Ctrl->Q.file);
328 gmt_disable_bghio_opts (GMT); /* Do not want any -b -g -h -i -o to affect the reading from -Q files */
329 if ((Qin = GMT_Read_Data (API, GMT_IS_DATASET, GMT_IS_FILE, GMT_IS_POLY, GMT_READ_NORMAL, NULL, Ctrl->Q.file, NULL)) == NULL) {
330 Return (API->error);
331 }
332 if (Qin->n_columns < 2) {
333 GMT_Report (API, GMT_MSG_ERROR, "Input file %s has %d column(s) but at least 2 are needed\n", Ctrl->Q.file, (int)Qin->n_columns);
334 Return (GMT_DIM_TOO_SMALL);
335 }
336 gmt_reenable_bghio_opts (GMT); /* Recover settings provided by user (if -b -g -h -i were used at all) */
337 Table = Qin->table[0]; /* Only one table in a file */
338 GMT_Report (API, GMT_MSG_INFORMATION, "Found %" PRIu64 " segments\n", Table->n_segments);
339 lon = gmt_M_memory (GMT, NULL, Table->n_segments, double);
340 lat = gmt_M_memory (GMT, NULL, Table->n_segments, double);
341 if (Ctrl->N.active) { /* Must get nodes from separate file */
342 struct GMT_DATASET *Nin = NULL;
343 struct GMT_DATATABLE *NTable = NULL;
344 if ((error = GMT_Set_Columns (API, GMT_IN, 3, GMT_COL_FIX_NO_TEXT)) != GMT_NOERROR) {
345 Return (error);
346 }
347 GMT_Report (API, GMT_MSG_INFORMATION, "Read Nodes from %s ...", Ctrl->N.file);
348 gmt_disable_bghio_opts (GMT); /* Do not want any -b -g -h -i -o to affect the reading from -N files */
349 if ((Nin = GMT_Read_Data (API, GMT_IS_DATASET, GMT_IS_FILE, GMT_IS_POINT, GMT_READ_NORMAL, NULL, Ctrl->N.file, NULL)) == NULL) {
350 Return (API->error);
351 }
352 if (Nin->n_columns < 2) {
353 GMT_Report (API, GMT_MSG_ERROR, "Input file %s has %d column(s) but at least 2 are needed\n", Ctrl->N.file, (int)Nin->n_columns);
354 Return (GMT_DIM_TOO_SMALL);
355 }
356 gmt_reenable_bghio_opts (GMT); /* Recover settings provided by user (if -b -g -h -i were used at all) */
357 NTable = Nin->table[0]; /* Only one table in a file with a single segment */
358 if (NTable->n_segments != 1) {
359 GMT_Report (API, GMT_MSG_ERROR, "File %s can only have 1 segment!\n", Ctrl->N.file);
360 Return (GMT_RUNTIME_ERROR);
361 }
362 if (Table->n_segments != NTable->n_records) {
363 GMT_Report (API, GMT_MSG_ERROR, "Files %s and %s do not have same number of items!\n", Ctrl->Q.file, Ctrl->N.file);
364 Return (GMT_RUNTIME_ERROR);
365 }
366 gmt_M_memcpy (lon, NTable->segment[0]->data[GMT_X], NTable->n_records, double);
367 gmt_M_memcpy (lat, NTable->segment[0]->data[GMT_Y], NTable->n_records, double);
368 if (GMT_Destroy_Data (API, &Nin) != GMT_NOERROR) {
369 Return (API->error);
370 }
371 GMT_Report (API, GMT_MSG_INFORMATION, "Found %" PRIu64 " records\n", NTable->n_records);
372 }
373 else { /* Get extract them from the segment header */
374 for (node = 0; node < Table->n_segments; node++) {
375 if (Table->segment[node]->header == NULL) {
376 GMT_Report (API, GMT_MSG_ERROR, "No node-information found in the segment headers - must abort\n");
377 Return (GMT_RUNTIME_ERROR);
378 }
379 if (sscanf (Table->segment[node]->header, "%*s %*d %lf %lf", &lon[node], &lat[node]) != 2) {
380 GMT_Report (API, GMT_MSG_ERROR, "Could not obtain node-information from the segment headers - must abort\n");
381 Return (GMT_RUNTIME_ERROR);
382 }
383 }
384 }
385 }
386 else { /* Must process input point/line data */
387 n_in = (Ctrl->E.mode == SPHD_VALUES) ? 3 : 2;
388 if ((error = GMT_Set_Columns (API, GMT_IN, n_in, GMT_COL_FIX_NO_TEXT)) != GMT_NOERROR) {
389 Return (error);
390 }
391 if (GMT_Init_IO (API, GMT_IS_DATASET, GMT_IS_POINT, GMT_IN, GMT_ADD_DEFAULT, 0, options) != GMT_NOERROR) { /* Registers default input sources, unless already set */
392 Return (API->error);
393 }
394 if (GMT_Begin_IO (API, GMT_IS_DATASET, GMT_IN, GMT_HEADER_ON) != GMT_NOERROR) {
395 Return (API->error); /* Enables data input and sets access mode */
396 }
397
398 GMT->session.min_meminc = GMT_INITIAL_MEM_ROW_ALLOC; /* Start by allocating a 32 Mb chunk */
399
400 n_alloc = 0;
401 if (!Ctrl->C.active) gmt_M_malloc2 (GMT, lon, lat, 0, &n_alloc, double);
402 n_alloc = 0;
403 gmt_M_malloc3 (GMT, xx, yy, zz, 0, &n_alloc, double);
404 if (Ctrl->E.mode == SPHD_VALUES) z_val = gmt_M_memory (GMT, NULL, n_alloc, gmt_grdfloat);
405
406 n = 0;
407 do { /* Keep returning records until we reach EOF */
408 if ((In = GMT_Get_Record (API, GMT_READ_DATA, NULL)) == NULL) { /* Read next record, get NULL if special case */
409 if (gmt_M_rec_is_error (GMT)) { /* Bail if there are any read errors */
410 if (Ctrl->E.mode == SPHD_VALUES) gmt_M_free (GMT, z_val);
411 if (!Ctrl->C.active) {gmt_M_free (GMT, lon); gmt_M_free (GMT, lat);}
412 if (!Ctrl->Q.active) {gmt_M_free (GMT, xx); gmt_M_free (GMT, yy); gmt_M_free (GMT, zz);}
413 Return (GMT_RUNTIME_ERROR);
414 }
415 else if (gmt_M_rec_is_eof (GMT)) /* Reached end of file */
416 break;
417 else if (gmt_M_rec_is_segment_header (GMT)) /* Parse segment headers */
418 first = true;
419 continue; /* Go back and read the next record */
420 }
421
422 if (In->data == NULL) {
423 gmt_quit_bad_record (API, In);
424 Return (API->error);
425 }
426
427 /* Data record to process - avoid duplicate points with -D as gmt_stripack_lists cannot handle that */
428 in = In->data; /* Only need to process numerical part here */
429
430 if (Ctrl->D.active) { /* Check for duplicates */
431 if (first) { /* Beginning of new segment; keep track of the very first coordinate in case of duplicates */
432 first_x = prev_x = in[GMT_X]; first_y = prev_y = in[GMT_Y];
433 }
434 else { /* Look for duplicate point at end of segments that replicate start point */
435 if (in[GMT_X] == first_x && in[GMT_Y] == first_y) { /* If any point after the first matches the first */
436 n_dup++;
437 continue;
438 }
439 if (n && in[GMT_X] == prev_x && in[GMT_Y] == prev_y) { /* Identical neighbors */
440 n_dup++;
441 continue;
442 }
443 prev_x = in[GMT_X]; prev_y = in[GMT_Y];
444 }
445 }
446
447 /* Convert lon,lat in degrees to Cartesian x,y,z triplets */
448 gmt_geo_to_cart (GMT, in[GMT_Y], in[GMT_X], X, true);
449
450 xx[n] = X[GMT_X]; yy[n] = X[GMT_Y]; zz[n] = X[GMT_Z];
451 if (!Ctrl->C.active) {
452 lon[n] = in[GMT_X]; lat[n] = in[GMT_Y];
453 }
454 if (Ctrl->E.mode == SPHD_VALUES) z_val[n] = (gmt_grdfloat)in[GMT_Z];
455
456 if (++n == n_alloc) { /* Get more memory */
457 if (!Ctrl->C.active) { size_t n_tmp = n_alloc; gmt_M_malloc2 (GMT, lon, lat, n, &n_tmp, double); }
458 gmt_M_malloc3 (GMT, xx, yy, zz, n, &n_alloc, double);
459 if (Ctrl->E.mode == SPHD_VALUES) z_val = gmt_M_memory (GMT, z_val, n_alloc, gmt_grdfloat);
460 }
461 first = false;
462 } while (true);
463
464 n_alloc = n;
465 if (!Ctrl->C.active) gmt_M_malloc2 (GMT, lon, lat, 0, &n_alloc, double);
466 gmt_M_malloc3 (GMT, xx, yy, zz, 0, &n_alloc, double);
467
468 if (Ctrl->D.active) { /* Report */
469 if (n_dup)
470 GMT_Report (API, GMT_MSG_WARNING, "Skipped %" PRIu64 " duplicate points in segments\n", n_dup);
471 else
472 GMT_Report (API, GMT_MSG_INFORMATION, "No duplicate points found in the segments\n");
473 }
474 else
475 GMT_Report (API, GMT_MSG_INFORMATION, "No duplicate check performed [-D was not activated]\n");
476
477 GMT_Report (API, GMT_MSG_INFORMATION, "Do Voronoi construction using %" PRIu64 " points\n", n);
478
479 T.mode = VORONOI;
480 if (gmt_stripack_lists (GMT, n, xx, yy, zz, &T)) Return (GMT_RUNTIME_ERROR); /* Do the basic triangulation */
481 gmt_M_free (GMT, T.D.tri); /* Don't need the triangulation */
482 if (Ctrl->C.active) { /* Recompute lon,lat and set pointers */
483 gmt_n_cart_to_geo (GMT, n, xx, yy, zz, xx, yy); /* Revert to lon, lat */
484 lon = xx;
485 lat = yy;
486 }
487 gmt_M_free (GMT, zz);
488 if (GMT_End_IO (API, GMT_IN, 0) != GMT_NOERROR) { /* Disables further data input */
489 Return (API->error);
490 }
491 GMT->session.min_meminc = GMT_MIN_MEMINC; /* Reset to the default value */
492 }
493
494 /* OK, time to create and work on the distance grid */
495
496 if ((Grid = GMT_Create_Data (API, GMT_IS_GRID, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA, NULL, NULL, NULL, \
497 GMT_GRID_DEFAULT_REG, GMT_NOTSET, NULL)) == NULL) Return (API->error);
498 GMT_Report (API, GMT_MSG_INFORMATION, "Start processing distance grid\n");
499
500 grid_lon = Grid->x;
501 grid_lat = Grid->y;
502
503 nx1 = (Grid->header->registration == GMT_GRID_PIXEL_REG) ? Grid->header->n_columns : Grid->header->n_columns - 1;
504 periodic = gmt_M_360_range (GMT->common.R.wesn[XLO], GMT->common.R.wesn[XHI]);
505 duplicate_col = (periodic && Grid->header->registration == GMT_GRID_NODE_REG); /* E.g., lon = 0 column should match lon = 360 column */
506 gmt_set_inside_mode (GMT, NULL, GMT_IOO_SPHERICAL);
507
508 if (Ctrl->Q.active) { /* Pre-chewed, just get number of nodes */
509 struct GMT_DATASET_HIDDEN *QH = gmt_get_DD_hidden (Qin);
510 n = Table->n_segments;
511 duplicate = (QH->alloc_mode == GMT_ALLOC_EXTERNALLY);
512 }
513 else
514 V = &T.V;
515
516 for (node = 0; node < n; node++) {
517 GMT_Report (API, GMT_MSG_INFORMATION, "Processing polygon %7ld\r", node);
518 if (Ctrl->Q.active) { /* Just point to next polygon */
519 if (duplicate) /* Must duplicate externally allocated segment since it needs to be resampled below */
520 P = gmt_duplicate_segment (GMT, Table->segment[node]);
521 else
522 P = Table->segment[node];
523 }
524 else { /* Obtain current polygon from Voronoi listings */
525 if (P == NULL) { /* Need a single polygon structure that we reuse for each polygon */
526 P = gmt_get_segment (GMT, 2); /* Needed as pointer below */
527 P->data = gmt_M_memory (GMT, NULL, 2, double *); /* Needed as pointers below */
528 P->min = gmt_M_memory (GMT, NULL, 2, double); /* Needed to hold min lon/lat */
529 P->max = gmt_M_memory (GMT, NULL, 2, double); /* Needed to hold max lon/lat */
530 P->n_columns = 2; p_alloc = 0;
531 gmt_M_malloc2 (GMT, P->data[GMT_X], P->data[GMT_Y], GMT_TINY_CHUNK, &p_alloc, double);
532 }
533 node_new = node_stop = V->lend[node];
534 vertex_new = V->listc[node_new];
535
536 /* Each iteration of this do-loop walks along one side of the polygon,
537 considering the subtriangle NODE --> VERTEX_LAST --> VERTEX. */
538
539 vertex = 0;
540 do {
541 node_last = node_new;
542 node_new = V->lptr[node_last];
543 vertex_last = vertex_new;
544 vertex_new = V->listc[node_new];
545
546 P->data[GMT_X][vertex] = V->lon[vertex_last];
547 P->data[GMT_Y][vertex] = V->lat[vertex_last];
548 if (P->data[GMT_X][vertex] < 0.0) P->data[GMT_X][vertex] += 360.0;
549 if (P->data[GMT_X][vertex] == 360.0) P->data[GMT_X][vertex] = 0.0;
550 vertex++;
551 if (vertex == p_alloc) gmt_M_malloc2 (GMT, P->data[GMT_X], P->data[GMT_Y], vertex, &p_alloc, double);
552
553 /* When we reach the vertex where we started, we are done with this polygon */
554 } while (node_new != node_stop);
555 P->data[GMT_X][vertex] = P->data[GMT_X][0]; /* Close polygon explicitly */
556 P->data[GMT_Y][vertex] = P->data[GMT_Y][0];
557 if ((++vertex) == p_alloc) gmt_M_malloc2 (GMT, P->data[GMT_X], P->data[GMT_Y], vertex, &p_alloc, double);
558 P->n_rows = vertex;
559 switch (Ctrl->E.mode) {
560 case SPHD_NODES: f_val = (gmt_grdfloat)node; break;
561 case SPHD_VALUES: f_val = z_val[node]; break;
562 default: break; /* Must compute distances below */
563 }
564 }
565
566 /* Here we have the polygon in P */
567
568 if ((n_new = gmt_fix_up_path (GMT, &P->data[GMT_X], &P->data[GMT_Y], P->n_rows, Ctrl->E.dist, GMT_STAIRS_OFF)) == 0) {
569 gmt_M_free (GMT, P);
570 Return (GMT_RUNTIME_ERROR);
571 }
572 P->n_rows = n_new;
573 sphdistance_prepare_polygon (GMT, P); /* Determine the enclosing sector */
574
575 south_row = (int)gmt_M_grd_y_to_row (GMT, P->min[GMT_Y], Grid->header);
576 north_row = (int)gmt_M_grd_y_to_row (GMT, P->max[GMT_Y], Grid->header);
577 w_col = (int)gmt_M_grd_x_to_col (GMT, P->min[GMT_X], Grid->header);
578 while (w_col < 0) w_col += nx1;
579 west_col = w_col;
580 e_col = (int)gmt_M_grd_x_to_col (GMT, P->max[GMT_X], Grid->header);
581 while (e_col < w_col) e_col += nx1;
582 east_col = e_col;
583 /* So here, any polygon will have a positive (or 0) west_col with an east_col >= west_col */
584 for (s_row = north_row; s_row <= south_row; s_row++) { /* For each scanline intersecting this polygon */
585 if (s_row < 0) continue; /* North of region */
586 row = s_row; if (row >= Grid->header->n_rows) continue; /* South of region */
587 for (p_col = west_col; p_col <= east_col; p_col++) { /* March along the scanline using col >= 0 */
588 if (p_col >= Grid->header->n_columns) { /* Off the east end of the grid */
589 if (periodic) /* Just shuffle to the corresponding point inside the global grid */
590 col = p_col - nx1;
591 else /* Sorry, really outside the region */
592 continue;
593 }
594 else
595 col = p_col;
596 side = gmt_inonout (GMT, grid_lon[col], grid_lat[row], P);
597
598 if (side == GMT_OUTSIDE) continue; /* Outside spherical polygon */
599 ij = gmt_M_ijp (Grid->header, row, col);
600 if (Ctrl->E.mode == SPHD_DIST)
601 f_val = (gmt_grdfloat)gmt_distance (GMT, grid_lon[col], grid_lat[row], lon[node], lat[node]);
602 Grid->data[ij] = f_val;
603 n_set++;
604 if (duplicate_col) { /* Duplicate the repeating column on the other side of this one */
605 if (col == 0) Grid->data[ij+nx1] = Grid->data[ij], n_set++;
606 else if (col == nx1) Grid->data[ij-nx1] = Grid->data[ij], n_set++;
607 }
608 }
609 }
610 if (duplicate) /* Free duplicate segment */
611 gmt_free_segment (GMT, &P);
612 }
613 GMT_Report (API, GMT_MSG_INFORMATION, "Processing polygon %7ld\n", node);
614
615 if (!Ctrl->Q.active) {
616 gmt_M_free (GMT, P->data[GMT_X]);
617 gmt_M_free (GMT, P->data[GMT_Y]);
618 gmt_free_segment (GMT, &P);
619 gmt_M_free (GMT, T.V.lon);
620 gmt_M_free (GMT, T.V.lat);
621 gmt_M_free (GMT, T.V.lend);
622 gmt_M_free (GMT, T.V.listc);
623 gmt_M_free (GMT, T.V.lptr);
624 gmt_M_free (GMT, xx);
625 gmt_M_free (GMT, yy);
626 }
627 if (!Ctrl->C.active) {
628 gmt_M_free (GMT, lon);
629 gmt_M_free (GMT, lat);
630 }
631 if (Ctrl->E.mode == SPHD_VALUES) gmt_M_free (GMT, z_val);
632
633 if (GMT_Set_Comment (API, GMT_IS_GRID, GMT_COMMENT_IS_OPTION | GMT_COMMENT_IS_COMMAND, options, Grid)) Return (API->error);
634 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) {
635 Return (API->error);
636 }
637
638 if (n_set > Grid->header->nm) n_set = Grid->header->nm; /* Not confuse the public */
639 GMT_Report (API, GMT_MSG_INFORMATION, "Spherical distance calculation completed, %" PRIu64 " nodes visited (at least once)\n", n_set);
640
641 Return (GMT_NOERROR);
642 }
643