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 * This program is free software; you can redistribute it and/or modify
6 * it under the terms of the GNU Lesser General Public License as published by
7 * the Free Software Foundation; version 3 or any later version.
8 *
9 * This program is distributed in the hope that it will be useful,
10 * but WITHOUT ANY WARRANTY; without even the implied warranty of
11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 * GNU Lesser General Public License for more details.
13 *
14 * Contact info: www.generic-mapping-tools.org
15 *--------------------------------------------------------------------*/
16 /*
17 * gmtspatial performs miscellaneous geospatial operations on polygons, such
18 * as truncating them against a clipping polygon, calculate areas, find
19 * crossings with other polygons, etc.
20 *
21 * Author: Paul Wessel
22 * Date: 10-Jun-2009
23 * Version: 6 API
24 */
25
26 #include "gmt_dev.h"
27
28 #define THIS_MODULE_CLASSIC_NAME "gmtspatial"
29 #define THIS_MODULE_MODERN_NAME "gmtspatial"
30 #define THIS_MODULE_LIB "core"
31 #define THIS_MODULE_PURPOSE "Geospatial operations on points, lines and polygons"
32 #define THIS_MODULE_KEYS "<D{,DD(=f,ND(=,TD(,>D}"
33 #define THIS_MODULE_NEEDS ""
34 #define THIS_MODULE_OPTIONS "-:RVabdefghijoqs" GMT_OPT("HMm")
35
36 #define GMT_W 3
37
38 #define POL_UNION 1
39 #define POL_INTERSECTION 2
40 #define POL_CLIP 3
41 #define POL_SPLIT 4
42 #define POL_JOIN 5
43 #define POL_HOLE 6
44 #define POL_BUFFER 7
45 #define POL_CENTROID 8
46
47 #define MIN_AREA_DIFF 0.01; /* If two polygons have areas that differ more than 1 % of each other then they are not the same feature */
48 #define MIN_SEPARATION 0 /* If the two closest points for two features are > 0 units apart then they are not the same feature */
49 #define MIN_CLOSENESS 0.01 /* If two close segments has an mean separation exceeding 1% of segment length, then they are not the same feature */
50 #define MIN_SUBSET 2.0 /* If two close segments deemed approximate fits has lengths that differ by this factor then they are sub/super sets of each other */
51
52 #ifdef HAVE_GEOS
53 #include <geos_c.h>
54 int geos_methods(struct GMT_CTRL *GMT, struct GMT_DATASET *D, char *fname, double buf_dist, char *method);
55 int geos_method_polygon(struct GMT_CTRL *GMT, struct GMT_DATASET *Din, struct GMT_DATASET *Dout, char *method);
56 int geos_method_linestring(struct GMT_CTRL *GMT, struct GMT_DATASET *Din, struct GMT_DATASET *Dout, double buf_dist, char *method);
57 #endif
58
59 struct GMTSPATIAL_DUP { /* Holds information on which single segment is closest to the current test segment */
60 uint64_t point;
61 uint64_t segment;
62 unsigned int table;
63 int mode;
64 bool inside;
65 double distance;
66 double mean_distance;
67 double closeness;
68 double setratio;
69 double a_threshold;
70 double d_threshold;
71 double c_threshold;
72 double s_threshold;
73 };
74
75 struct DUP_INFO {
76 uint64_t point;
77 int mode;
78 double distance;
79 double closeness;
80 double setratio;
81 };
82
83 struct GMTSPATIAL_CTRL {
84 struct GMTSPATIAL_Out { /* -> */
85 bool active;
86 char *file;
87 } Out;
88 struct GMTSPATIAL_A { /* -Aa<min_dist>, -A */
89 bool active;
90 unsigned int mode;
91 int smode;
92 double min_dist;
93 char unit;
94 } A;
95 struct GMTSPATIAL_C { /* -C */
96 bool active;
97 } C;
98 struct GMTSPATIAL_D { /* -D[pol] */
99 bool active;
100 int mode;
101 char unit;
102 char *file;
103 struct GMTSPATIAL_DUP I;
104 } D;
105 struct GMTSPATIAL_E { /* -E+n|p */
106 bool active;
107 unsigned int mode;
108 } E;
109 struct GMTSPATIAL_F { /* -F */
110 bool active;
111 unsigned int geometry;
112 } F;
113 struct GMTSPATIAL_I { /* -I[i|e] */
114 bool active;
115 unsigned int mode;
116 } I;
117 struct GMTSPATIAL_L { /* -L */
118 bool active;
119 char unit;
120 double s_cutoff, path_noise, box_offset;
121 } L;
122 struct GMTSPATIAL_N { /* -N<file>[+a][+p>ID>][+r][+z] */
123 bool active;
124 bool all; /* All points in lines and polygons must be inside a polygon for us to report ID */
125 unsigned int mode; /* 0 for reporting ID in -Z<ID> header, 1 via data column, 2 just as a report */
126 unsigned int ID; /* If 1 we use running numbers */
127 char *file;
128 } N;
129 struct GMTSPATIAL_Q { /* -Q[+c<min>/<max>][+h][+l][+p][+s[a|d]] */
130 bool active;
131 bool header; /* Place dimension and centroid in segment headers */
132 bool area; /* Apply range test on dimension */
133 bool sort; /* Sort segments based on dimension */
134 unsigned int mode; /* 0 use input as is, 1 force to line, 2 force to polygon */
135 unsigned int dmode; /* for geo data: 1 = flat earth, 2 = great circle, 3 = geodesic (for distances) */
136 int dir; /* For segment sorting: -1 is descending, +1 is ascending [Default] */
137 double limit[2]; /* Min and max area or length for output segments */
138 char unit;
139 } Q;
140 struct GMTSPATIAL_S { /* -S[u|i|c|j|h] */
141 bool active;
142 unsigned int mode;
143 double width;
144 } S;
145 struct GMTSPATIAL_T { /* -T[pol] */
146 bool active;
147 char *file;
148 } T;
149 };
150
151 struct GMTSPATIAL_PAIR {
152 double node;
153 uint64_t pos;
154 };
155
156 #ifdef __APPLE__
157 /* macOX has it built in, so ensure we define this flag */
158 #define HAVE_MERGESORT
159 #endif
160
161 #ifndef HAVE_MERGESORT
162 #include "mergesort.c"
163 #endif
164
New_Ctrl(struct GMT_CTRL * GMT)165 static void *New_Ctrl (struct GMT_CTRL *GMT) { /* Allocate and initialize a new control structure */
166 struct GMTSPATIAL_CTRL *C;
167
168 C = gmt_M_memory (GMT, NULL, 1, struct GMTSPATIAL_CTRL);
169
170 /* Initialize values whose defaults are not 0/false/NULL */
171 C->A.unit = 'X'; /* Cartesian units as default */
172 C->L.box_offset = 1.0e-10; /* Minimum significant amplitude */
173 C->L.path_noise = 1.0e-10; /* Minimum significant amplitude */
174 C->D.unit = 'X'; /* Cartesian units as default */
175 C->D.I.a_threshold = MIN_AREA_DIFF;
176 C->D.I.d_threshold = MIN_SEPARATION;
177 C->D.I.c_threshold = MIN_CLOSENESS;
178 C->D.I.s_threshold = MIN_SUBSET;
179 C->Q.mode = GMT_IS_POINT; /* Undecided on line vs poly */
180 C->Q.dmode = 2; /* Great-circle distance if not specified */
181 return (C);
182 }
183
Free_Ctrl(struct GMT_CTRL * GMT,struct GMTSPATIAL_CTRL * C)184 static void Free_Ctrl (struct GMT_CTRL *GMT, struct GMTSPATIAL_CTRL *C) { /* Deallocate control structure */
185 if (!C) return;
186 gmt_M_str_free (C->Out.file);
187 gmt_M_str_free (C->D.file);
188 gmt_M_str_free (C->N.file);
189 gmt_M_str_free (C->T.file);
190 gmt_M_free (GMT, C);
191 }
192
gmtspatial_area_size(struct GMT_CTRL * GMT,double x[],double y[],uint64_t n,double * out,int geo)193 GMT_LOCAL unsigned int gmtspatial_area_size (struct GMT_CTRL *GMT, double x[], double y[], uint64_t n, double *out, int geo) {
194 double size = gmt_centroid_area (GMT, x, y, n, geo, out);
195 out[GMT_Z] = fabs (size);
196 return ((size < 0.0) ? GMT_POL_IS_CCW : GMT_POL_IS_CW);
197 }
198
199 #if 0
200 GMT_LOCAL unsigned int gmtspatial_gmtspatial_area_size_old (struct GMT_CTRL *GMT, double x[], double y[], uint64_t n, double *out, int geo) {
201 uint64_t i;
202 double wesn[4], xx, yy, size, ix, iy;
203 double *xp = NULL, *yp = NULL;
204
205 wesn[XLO] = wesn[YLO] = DBL_MAX;
206 wesn[XHI] = wesn[YHI] = -DBL_MAX;
207
208 gmt_centroid (GMT, x, y, n, out, geo); /* Get mean location */
209
210 for (i = 0; i < n; i++) {
211 wesn[XLO] = MIN (wesn[XLO], x[i]);
212 wesn[XHI] = MAX (wesn[XHI], x[i]);
213 wesn[YLO] = MIN (wesn[YLO], y[i]);
214 wesn[YHI] = MAX (wesn[YHI], y[i]);
215 }
216 xp = gmt_M_memory (GMT, NULL, n, double); yp = gmt_M_memory (GMT, NULL, n, double);
217 if (geo == 1) { /* Initializes GMT projection parameters to the -JA settings */
218 GMT->current.proj.projection_GMT = GMT->current.proj.projection = GMT_LAMB_AZ_EQ;
219 GMT->current.proj.unit = 1.0;
220 GMT->current.proj.pars[3] = 39.3700787401574814;
221 GMT->common.R.oblique = false;
222 GMT->common.J.active = true;
223 GMT->current.setting.map_line_step = 1.0e7; /* To avoid nlon/nlat being huge */
224 GMT->current.proj.pars[0] = out[GMT_X];
225 GMT->current.proj.pars[1] = out[GMT_Y];
226 if (gmt_M_err_pass (GMT, gmt_proj_setup (GMT, wesn), "")) {
227 gmt_M_free (GMT, xp); gmt_M_free (GMT, yp);
228 return (0);
229 }
230
231 ix = 1.0 / GMT->current.proj.scale[GMT_X];
232 iy = 1.0 / GMT->current.proj.scale[GMT_Y];
233
234 for (i = 0; i < n; i++) {
235 gmt_geo_to_xy (GMT, x[i], y[i], &xx, &yy);
236 xp[i] = (xx - GMT->current.proj.origin[GMT_X]) * ix;
237 yp[i] = (yy - GMT->current.proj.origin[GMT_Y]) * iy;
238 }
239 }
240 else { /* Just take out mean coordinates */
241 for (i = 0; i < n; i++) {
242 xp[i] = x[i] - out[GMT_X];
243 yp[i] = y[i] - out[GMT_Y];
244 }
245 }
246
247 size = gmt_pol_area (xp, yp, n);
248 gmt_M_free (GMT, xp);
249 gmt_M_free (GMT, yp);
250 if (geo) size *= (GMT->current.map.dist[GMT_MAP_DIST].scale * GMT->current.map.dist[GMT_MAP_DIST].scale);
251 out[GMT_Z] = fabs (size);
252 return ((size < 0.0) ? GMT_POL_IS_CCW : GMT_POL_IS_CW);
253 }
254 #endif
255
gmtspatial_length_size(struct GMT_CTRL * GMT,double x[],double y[],uint64_t n,double * out)256 GMT_LOCAL void gmtspatial_length_size (struct GMT_CTRL *GMT, double x[], double y[], uint64_t n, double *out) {
257 uint64_t i;
258 double length = 0.0, mid, f, *s = NULL;
259
260 assert (n > 0);
261
262 /* Estimate 'average' position */
263
264 s = gmt_M_memory (GMT, NULL, n, double);
265 for (i = 1; i < n; i++) {
266 length += gmt_distance (GMT, x[i-1], y[i-1], x[i], y[i]);
267 s[i] = length;
268 }
269 mid = 0.5 * length;
270 i = 0;
271 while (s[i] <= mid) i++; /* Find mid point length-wise */
272 f = (mid - s[i-1]) / (s[i] - s[i-1]);
273 out[GMT_X] = x[i-1] + f * (x[i] - x[i-1]);
274 out[GMT_Y] = y[i-1] + f * (y[i] - y[i-1]);
275 out[GMT_Z] = length;
276 gmt_M_free (GMT, s);
277 }
278
gmtspatial_comp_pairs(const void * a,const void * b)279 GMT_LOCAL int gmtspatial_comp_pairs (const void *a, const void *b) {
280 const struct GMTSPATIAL_PAIR *xa = a, *xb = b;
281 /* Sort on node value */
282 if (xa->node < xb->node) return (-1);
283 if (xa->node > xb->node) return (+1);
284 return (0);
285 }
286
gmtspatial_write_record(struct GMT_CTRL * GMT,double ** R,uint64_t n,uint64_t p)287 GMT_LOCAL void gmtspatial_write_record (struct GMT_CTRL *GMT, double **R, uint64_t n, uint64_t p) {
288 uint64_t c;
289 double out[GMT_MAX_COLUMNS];
290 struct GMT_RECORD Out;
291 Out.data = out; Out.text = NULL;
292 for (c = 0; c < n; c++) out[c] = R[c][p];
293 GMT_Put_Record (GMT->parent, GMT_WRITE_DATA, &Out);
294 }
295
gmtspatial_is_duplicate(struct GMT_CTRL * GMT,struct GMT_DATASEGMENT * S,struct GMT_DATASET * D,struct GMTSPATIAL_DUP * I,struct DUP_INFO ** L)296 GMT_LOCAL int gmtspatial_is_duplicate (struct GMT_CTRL *GMT, struct GMT_DATASEGMENT *S, struct GMT_DATASET *D, struct GMTSPATIAL_DUP *I, struct DUP_INFO **L) {
297 /* Given single line segment S and a dataset of many line segments in D, determine the closest neighbor
298 * to S in D (call it S'), and if "really close" it might be a duplicate or slight revision to S.
299 * There might be several features S' in D close to S so we return how many near or exact matches we
300 * found and pass the details via L. We return 0 if no duplicates are found. If some are found then
301 * these are the possible types of result per match:
302 * +1 : S has identical duplicate in D (S == S')
303 * -1 : S has identical duplicate in D (S == S'), but is reversed
304 * +2 : S is very close to S', probably a slight revised version
305 * -2 : S is very close to S', probably a slight revised version, but is reversed
306 * +3 : Same as +2 but S is likely a subset of S'
307 * -3 : Same as -2 but S is likely a reversed subset of S'
308 * +4 : Same as +2 but S is likely a superset of S'
309 * -4 : Same as -2 but S is likely a reversed superset of S'
310 * 5 : Two lines split exactly at the Dateline
311 * The algorithm first finds the smallest separation from any point on S to the line Sp.
312 * We then reverse the situation and check the smallest separation from any point on any
313 * of the Sp candidates to the line S.
314 * If the smallest distance found exceeds I->d_threshold then we move on (not close enough).
315 * We next compute the mean (or median, with +Ccmax) closeness, defined as the ratio between
316 * the mean (or median) separation between points on S and the line Sp (or vice versa) and
317 * the length of S (or Sp). We compute it both ways since the segments may differ in lengths
318 * and degree of overlap (i.e., subsets or supersets).
319 * NOTE: Lines that intersect obviously has a zero min distance but since we only compute
320 * distances from the nodes on one line to the nearest point along the other line we will
321 * not detect such crossings. For most situations this should not matter much (?).
322 */
323
324 bool status;
325 unsigned int k, tbl, n_close = 0, n_dup = 0, mode1, mode3;
326 uint64_t row, seg, pt, np, sno, *n_sum = NULL;
327 int k_signed;
328 double dist, f_seg, f_pt, d1, d2, closest, length[2], separation[2], close[2], *d_mean = NULL;
329 double med_separation[2], med_close[2], high = 0, low = 0, use_length, use_sep, use_close, *sep = NULL;
330 struct GMT_DATASEGMENT *Sp = NULL;
331 struct GMT_DATASEGMENT_HIDDEN *SH = NULL;
332
333 GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Determine the segments in D closest to our segment\n");
334 I->distance = I->mean_distance = DBL_MAX;
335 mode3 = 3 + 10 * I->inside; /* Set gmt_near_lines modes */
336 mode1 = 1 + 10 * I->inside; /* Set gmt_near_lines modes */
337
338 d_mean = gmt_M_memory (GMT, NULL, D->n_segments, double); /* Mean distances from points along S to other lines */
339 n_sum = gmt_M_memory (GMT, NULL, D->n_segments, uint64_t); /* Number of distances from points along S to other lines */
340
341 /* We first want to obtain the mean distance from one segment to all others. We do this by computing the
342 * nearest distance from each point along our segment S to the other segments and compute the average of
343 * all those distances. Then we reverse the search and continue to accumulate averages */
344
345 /* Process each point along the trace in S and find sum of nearest distances for each segment in table */
346 for (row = 0; row < S->n_rows; row++) {
347 for (sno = tbl = 0; tbl < D->n_tables; tbl++) { /* For each table to compare it to */
348 for (seg = 0; seg < D->table[tbl]->n_segments; seg++, sno++) { /* For each segment in current table */
349 if (D->table[tbl]->segment[seg]->n_rows == 0) continue; /* Skip segments with no records (may be itself) */
350 dist = DBL_MAX; /* Reset for each line to find distance to that line */
351 (void) gmt_near_a_line (GMT, S->data[GMT_X][row], S->data[GMT_Y][row], seg, D->table[tbl]->segment[seg], mode3, &dist, &f_seg, &f_pt);
352 d_mean[sno] += dist; /* Sum of distances to this segment */
353 n_sum[sno]++; /* Number of such distances */
354 }
355 }
356 }
357
358 /* Must also do the reverse test: for each point for each line in the table, compute distance to S; it might be shorter */
359 for (sno = tbl = 0; tbl < D->n_tables; tbl++) { /* For each table to compare it to */
360 for (seg = 0; seg < D->table[tbl]->n_segments; seg++, sno++) { /* For each segment in current table */
361 Sp = D->table[tbl]->segment[seg]; /* This is S', one of the segments that is close to S */
362 for (row = 0; row < Sp->n_rows; row++) { /* Process each point along the trace in S and find nearest distance for each segment in table */
363 dist = DBL_MAX; /* Reset for each line to find distance to that line */
364 (void) gmt_near_a_line (GMT, Sp->data[GMT_X][row], Sp->data[GMT_Y][row], seg, S, mode3, &dist, &f_seg, &f_pt);
365 d_mean[sno] += dist; /* Sum of distances to this segment */
366 n_sum[sno]++; /* Number of such distances */
367 }
368 }
369 }
370 /* Compute the average distances */
371 for (sno = 0; sno < D->n_segments; sno++) {
372 d_mean[sno] = (n_sum[sno] > 0) ? d_mean[sno] / n_sum[sno] : DBL_MAX;
373 }
374 gmt_M_free (GMT, n_sum);
375
376 /* Process each point along the trace in S and find nearest distance for each segment in table */
377 for (row = 0; row < S->n_rows; row++) {
378 for (sno = tbl = 0; tbl < D->n_tables; tbl++) { /* For each table to compare it to */
379 for (seg = 0; seg < D->table[tbl]->n_segments; seg++, sno++) { /* For each segment in current table */
380 SH = gmt_get_DS_hidden (D->table[tbl]->segment[seg]);
381 dist = DBL_MAX; /* Reset for each line to find distance to that line */
382 status = gmt_near_a_line (GMT, S->data[GMT_X][row], S->data[GMT_Y][row], seg, D->table[tbl]->segment[seg], mode3, &dist, &f_seg, &f_pt);
383 if (!status && I->inside) continue; /* Only consider points that project perpendicularly within the line segment */
384 pt = lrint (f_pt); /* We know f_seg == seg so no point assigning that */
385 if (dist < I->distance && d_mean[sno] < I->mean_distance) { /* Keep track of the single closest feature */
386 I->point = pt;
387 I->segment = seg;
388 I->distance = dist;
389 I->mean_distance = d_mean[sno];
390 I->table = tbl;
391 }
392 if (dist > I->d_threshold) continue; /* Not close enough for duplicate consideration */
393 if (SH->mode == 0) {
394 n_close++;
395 L[tbl][seg].distance = DBL_MAX;
396 }
397 SH->mode = 1; /* Use mode to flag segments that are close enough */
398 if (dist < L[tbl][seg].distance) { /* Keep track of the closest feature */
399 L[tbl][seg].point = pt;
400 L[tbl][seg].distance = dist;
401 }
402 }
403 }
404 }
405
406 /* Must also do the reverse test: for each point for each line in the table, compute distance to S; it might be shorter */
407
408 for (sno = tbl = 0; tbl < D->n_tables; tbl++) { /* For each table to compare it to */
409 for (seg = 0; seg < D->table[tbl]->n_segments; seg++, sno++) { /* For each segment in current table */
410 Sp = D->table[tbl]->segment[seg]; /* This is S', one of the segments that is close to S */
411 SH = gmt_get_DS_hidden (Sp);
412 for (row = 0; row < Sp->n_rows; row++) { /* Process each point along the trace in S and find nearest distance for each segment in table */
413 dist = DBL_MAX; /* Reset for each line to find distance to that line */
414 status = gmt_near_a_line (GMT, Sp->data[GMT_X][row], Sp->data[GMT_Y][row], seg, S, mode3, &dist, &f_seg, &f_pt);
415 if (!status && I->inside) continue; /* Only consider points that project perpendicularly within the line segment */
416 pt = lrint (f_pt);
417 if (dist < I->distance && d_mean[sno] < I->mean_distance) { /* Keep track of the single closest feature */
418 I->point = pt;
419 I->segment = seg;
420 I->distance = dist;
421 I->mean_distance = d_mean[sno];
422 I->table = tbl;
423 }
424 if (dist > I->d_threshold) continue; /* Not close enough for duplicate consideration */
425 if (SH->mode == 0) {
426 n_close++;
427 L[tbl][seg].distance = DBL_MAX;
428 }
429 SH->mode = 1; /* Use mode to flag segments that are close enough */
430 if (dist < L[tbl][seg].distance) { /* Keep track of the closest feature */
431 L[tbl][seg].point = pt;
432 L[tbl][seg].distance = dist;
433 }
434 }
435 }
436 }
437 gmt_M_free (GMT, d_mean);
438
439 if (n_close == 0)
440 GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "No other segment found within dmax [probably due to +p requirement]\n");
441 else
442 GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Closest segment (Table %d, segment %" PRIu64 ") is %.3f km away; %d segments found within dmax\n", I->table, I->segment, I->distance, n_close);
443
444 /* Here we have found the shortest distance from a point on S to another segment; S' is segment number seg */
445
446 if (I->distance > I->d_threshold) return (0); /* We found no duplicates or slightly different versions of S */
447
448 /* Here we need to compare one or more S' candidates and S a bit more. */
449
450 for (row = 1, length[0] = 0.0; row < S->n_rows; row++) { /* Compute length of S once */
451 length[0] += gmt_distance (GMT, S->data[GMT_X][row], S->data[GMT_Y][row], S->data[GMT_X][row-1], S->data[GMT_Y][row-1]);
452 }
453
454 for (tbl = 0; tbl < D->n_tables; tbl++) { /* For each table to compare it to */
455 for (seg = 0; seg < D->table[tbl]->n_segments; seg++) { /* For each segment in current table */
456 SH = gmt_get_DS_hidden (D->table[tbl]->segment[seg]);
457 if (SH->mode != 1) continue; /* Not one of the close segments we flagged earlier */
458 SH->mode = 0; /* Remove this temporary flag */
459
460 Sp = D->table[tbl]->segment[seg]; /* This is S', one of the segments that is close to S */
461 if (S->n_rows == Sp->n_rows) { /* Exactly the same number of data points; check for identical duplicate (possibly reversed) */
462 for (row = 0, d1 = d2 = 0.0; row < S->n_rows; row++) { /* Compare each point along the trace in S and Sp and S and Sp-reversed */
463 d1 += gmt_distance (GMT, S->data[GMT_X][row], S->data[GMT_Y][row], Sp->data[GMT_X][row], Sp->data[GMT_Y][row]);
464 d2 += gmt_distance (GMT, S->data[GMT_X][row], S->data[GMT_Y][row], Sp->data[GMT_X][S->n_rows-row-1], Sp->data[GMT_Y][S->n_rows-row-1]);
465 }
466 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "S to Sp of same length and gave d1 = %g and d2 = %g\n", d1, d2);
467 if (gmt_M_is_zero (d1)) { /* Exact duplicate */
468 L[tbl][seg].mode = +1;
469 n_dup++;
470 continue;
471 }
472 else if (gmt_M_is_zero (d2)) { /* Exact duplicate but reverse order */
473 L[tbl][seg].mode = -1;
474 n_dup++;
475 continue;
476 }
477 }
478
479 /* We get here when S' is not an exact duplicate of S, but approximate.
480 * We compute the mean [and possibly median] separation between S' and S
481 * and the closeness ratio (separation/length). */
482
483 separation[0] = 0.0;
484 if (I->mode) { /* Various items needed to compute median separation */
485 sep = gmt_M_memory (GMT, NULL, MAX (S->n_rows, Sp->n_rows), double);
486 low = DBL_MAX;
487 high = -DBL_MAX;
488 }
489 for (row = np = 0; row < S->n_rows; row++) { /* Process each point along the trace in S */
490 dist = DBL_MAX;
491 status = gmt_near_a_line (GMT, S->data[GMT_X][row], S->data[GMT_Y][row], seg, Sp, mode1, &dist, NULL, NULL);
492 if (!status && I->inside) continue; /* Only consider points that project perpendicularly within the Sp segment */
493 separation[0] += dist;
494 if (I->mode) { /* Special processing for median calculation */
495 sep[np] = dist;
496 if (dist < low) low = dist;
497 if (dist > high) high = dist;
498 }
499 np++; /* Number of points within the overlap zone */
500 }
501 separation[0] = (np > 1) ? separation[0] / np : DBL_MAX; /* Mean distance between S and S' */
502 use_length = (np) ? length[0] * np / S->n_rows : length[0]; /* ~reduce length to overlap section assuming equal point spacing */
503 close[0] = (np > 1) ? separation[0] / use_length : DBL_MAX; /* Closeness as viewed from S */
504 use_sep = (separation[0] == DBL_MAX) ? GMT->session.d_NaN : separation[0];
505 use_close = (close[0] == DBL_MAX) ? GMT->session.d_NaN : close[0];
506 GMT_Report (GMT->parent, GMT_MSG_INFORMATION,
507 "S has length %.3f km, has mean separation to Sp of %.3f km, and a closeness ratio of %g [n = %" PRIu64 "/%" PRIu64 "]\n",
508 length[0], use_sep, use_close, np, S->n_rows);
509 if (I->mode) {
510 if (np > 1) {
511 gmt_median (GMT, sep, np, low, high, separation[0], &med_separation[0]);
512 med_close[0] = med_separation[0] / use_length;
513 }
514 else med_close[0] = DBL_MAX;
515 use_sep = (np == 0 || med_separation[0] == DBL_MAX) ? GMT->session.d_NaN : med_separation[0];
516 use_close = (np == 0 || med_close[0] == DBL_MAX) ? GMT->session.d_NaN : med_close[0];
517 GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "S has median separation to Sp of %.3f km, and a robust closeness ratio of %g\n",
518 use_sep, use_close);
519 }
520
521 /* Must now compare the other way */
522
523 separation[1] = length[1] = 0.0;
524 if (I->mode) { /* Reset for median calculation */
525 low = +DBL_MAX;
526 high = -DBL_MAX;
527 }
528 for (row = np = 0; row < Sp->n_rows; row++) { /* Process each point along the trace in S' */
529 dist = DBL_MAX; /* Reset for each line to find distance to that line */
530 status = gmt_near_a_line (GMT, Sp->data[GMT_X][row], Sp->data[GMT_Y][row], seg, S, mode1, &dist, NULL, NULL);
531 if (row) length[1] += gmt_distance (GMT, Sp->data[GMT_X][row], Sp->data[GMT_Y][row], Sp->data[GMT_X][row-1], Sp->data[GMT_Y][row-1]);
532 if (!status && I->inside) continue; /* Only consider points that project perpendicularly within the Sp segment */
533 separation[1] += dist;
534 if (I->mode) { /* Special processing for median calculation */
535 sep[np] = dist;
536 if (dist < low) low = dist;
537 if (dist > high) high = dist;
538 }
539 np++; /* Number of points within the overlap zone */
540 }
541 separation[1] = (np > 1) ? separation[1] / np : DBL_MAX; /* Mean distance between S' and S */
542 use_length = (np) ? length[1] * np / Sp->n_rows : length[1]; /* ~reduce length to overlap section assuming equal point spacing */
543 close[1] = (np > 1) ? separation[1] / use_length : DBL_MAX; /* Closeness as viewed from S' */
544 use_sep = (separation[1] == DBL_MAX) ? GMT->session.d_NaN : separation[1];
545 use_close = (close[1] == DBL_MAX) ? GMT->session.d_NaN : close[1];
546 GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Sp has length %.3f km, has mean separation to S of %.3f km, and a closeness ratio of %g [n = %" PRIu64 "/%" PRIu64 "]\n",
547 length[1], use_sep, use_close, np, Sp->n_rows);
548 if (I->mode) {
549 if (np > 1) {
550 gmt_median (GMT, sep, np, low, high, separation[1], &med_separation[1]);
551 med_close[1] = med_separation[1] / use_length;
552 }
553 else med_close[1] = DBL_MAX;
554 gmt_M_free (GMT, sep);
555 use_sep = (med_separation[1] == DBL_MAX) ? GMT->session.d_NaN : med_separation[1];
556 use_close = (med_close[1] == DBL_MAX) ? GMT->session.d_NaN : med_close[1];
557 GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Sp has median separation to S of %.3f km, and a robust closeness ratio of %g\n",
558 use_sep, use_close);
559 k = (med_close[0] <= med_close[1]) ? 0 : 1; /* Pick the setup with the smallest robust closeness */
560 closest = med_close[k]; /* The longer the segment and the closer they are, the smaller the closeness */
561 }
562 else {
563 k = (close[0] <= close[1]) ? 0 : 1; /* Pick the setup with the smallest robust closeness */
564 closest = close[k]; /* The longer the segment and the closer they are, the smaller the closeness */
565 }
566
567 if (closest > I->c_threshold) continue; /* Not a duplicate or slightly different version of S */
568 L[tbl][seg].closeness = closest; /* The longer the segment and the closer they are the smaller the I->closeness */
569
570 /* Compute distances from first point in S to first (d1) and last (d2) in S' and use this info to see if S' is reversed */
571 d1 = gmt_distance (GMT, S->data[GMT_X][0], S->data[GMT_Y][0], Sp->data[GMT_X][0], Sp->data[GMT_Y][0]);
572 d2 = gmt_distance (GMT, S->data[GMT_X][0], S->data[GMT_Y][0], Sp->data[GMT_X][Sp->n_rows-1], Sp->data[GMT_Y][Sp->n_rows-1]);
573
574 L[tbl][seg].setratio = (length[0] > length[1]) ? length[0] / length[1] : length[1] / length[0]; /* Ratio of length is used to detect subset (or superset) */
575
576 if (length[0] < (length[1] / I->s_threshold))
577 k_signed = 3; /* S is a subset of Sp */
578 else if (length[1] < (length[0] / I->s_threshold))
579 k_signed = 4; /* S is a superset of Sp */
580 else
581 k_signed = 2; /* S and S' are practically the same feature */
582 L[tbl][seg].mode = (d1 < d2) ? k_signed : -k_signed; /* Negative means Sp is reversed w.r.t. S */
583 /* Last check: If two lines are split at 180 the above will find them to be approximate duplicates even
584 * though they just share one point (at 180 longitude) */
585 if (L[tbl][seg].closeness < GMT_CONV8_LIMIT && L[tbl][seg].distance < GMT_CONV8_LIMIT) {
586 if (fabs (S->data[GMT_X][0]) == 180.0 && (fabs (Sp->data[GMT_X][0]) == 180.0 || fabs (Sp->data[GMT_X][Sp->n_rows-1]) == 180.0))
587 L[tbl][seg].mode = 5;
588 else if (fabs (S->data[GMT_X][S->n_rows-1]) == 180.0 && (fabs (Sp->data[GMT_X][0]) == 180.0 || fabs (Sp->data[GMT_X][Sp->n_rows-1]) == 180.0))
589 L[tbl][seg].mode = 5;
590 }
591 n_dup++;
592 }
593 }
594 return (n_dup);
595 }
596
597 struct NN_DIST {
598 double data[4]; /* Up to x,y,z,weight */
599 double distance; /* Distance to nearest neighbor */
600 int64_t ID; /* Input ID # of this point (original input record number 0,1,...)*/
601 int64_t neighbor; /* Input ID # of this point's neighbor */
602 };
603
604 struct NN_INFO {
605 int64_t sort_rec; /* Input ID # of this point's neighbor */
606 int64_t orig_rec; /* Rec # of this point */
607 };
608
gmtspatial_compare_nn_points(const void * point_1v,const void * point_2v)609 GMT_LOCAL int gmtspatial_compare_nn_points (const void *point_1v, const void *point_2v) {
610 /* Routine for qsort to sort NN data structure on distance.
611 */
612 const struct NN_DIST *point_1 = point_1v, *point_2 = point_2v;
613
614 if (gmt_M_is_dnan (point_1->distance)) return (+1);
615 if (gmt_M_is_dnan (point_2->distance)) return (-1);
616 if (point_1->distance < point_2->distance) return (-1);
617 if (point_1->distance > point_2->distance) return (+1);
618 return (0);
619 }
620
gmtspatial_NNA_update_dist(struct GMT_CTRL * GMT,struct NN_DIST * P,uint64_t * n_points)621 GMT_LOCAL struct NN_DIST *gmtspatial_NNA_update_dist (struct GMT_CTRL *GMT, struct NN_DIST *P, uint64_t *n_points) {
622 /* Return array of NN results sorted on smallest distances */
623 int64_t k, k2, np;
624 double *distance = gmt_M_memory (GMT, NULL, *n_points, double);
625 #ifdef DEBUG
626 static int iteration = 0;
627 #endif
628 np = *n_points;
629 for (k = 0; k < (np-1); k++) {
630 if (gmt_M_is_dnan (P[k].distance)) continue; /* Skip deleted point */
631 P[k].distance = DBL_MAX;
632 /* We split the loop over calculation of distance from the loop over assignments since
633 * if OpenMP is used then we cannot interchange k and k2 as there may be overprinting.
634 */
635 #ifdef _OPENMP
636 #pragma omp parallel for private(k2) shared(k,np,GMT,P,distance)
637 #endif
638 for (k2 = k + 1; k2 < np; k2++) {
639 if (gmt_M_is_dnan (P[k2].distance)) continue; /* Skip deleted point */
640 distance[k2] = gmt_distance (GMT, P[k].data[GMT_X], P[k].data[GMT_Y], P[k2].data[GMT_X], P[k2].data[GMT_Y]);
641 }
642 for (k2 = k + 1; k2 < np; k2++) {
643 if (gmt_M_is_dnan (P[k2].distance)) continue; /* Skip deleted point */
644 if (distance[k2] < P[k].distance) {
645 P[k].distance = distance[k2];
646 P[k].neighbor = P[k2].ID;
647 }
648 if (distance[k2] < P[k2].distance) {
649 P[k2].distance = distance[k2];
650 P[k2].neighbor = P[k].ID;
651 }
652 }
653 }
654 gmt_M_free (GMT, distance);
655
656 /* Prefer mergesort since qsort is not stable for equalities */
657 mergesort (P, np, sizeof (struct NN_DIST), gmtspatial_compare_nn_points); /* Sort on small to large distances */
658
659 for (k = np; k > 0 && gmt_M_is_dnan (P[k-1].distance); k--); /* Skip the NaN distances that were placed at end */
660 *n_points = k; /* Update point count */
661 #ifdef DEBUG
662 if (gmt_M_is_verbose (GMT, GMT_MSG_DEBUG)) {
663 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "===> Iteration = %d\n", iteration);
664 for (k = 0; k < (int64_t)(*n_points); k++)
665 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "%6d\tID=%6d\tNeighbor=%6d\tDistance = %.12g\n", (int)k, (int)P[k].ID, P[k].neighbor, P[k].distance);
666 }
667 iteration++;
668 #endif
669 return (P);
670 }
671
gmtspatial_NNA_init_dist(struct GMT_CTRL * GMT,struct GMT_DATASET * D,uint64_t * n_points)672 GMT_LOCAL struct NN_DIST *gmtspatial_NNA_init_dist (struct GMT_CTRL *GMT, struct GMT_DATASET *D, uint64_t *n_points) {
673 /* Return array of NN results sorted on smallest distances */
674 uint64_t tbl, seg, row, col, n_cols;
675 int64_t k, np = 0; /* Must be signed due to Win OpenMP retardedness */
676 double *distance = NULL;
677 struct GMT_DATASEGMENT *S = NULL;
678 struct NN_DIST *P = gmt_M_memory (GMT, NULL, D->n_records, struct NN_DIST);
679
680 n_cols = MIN (D->n_columns, 4); /* Expects lon,lat and makes room for at least z, w and other columns */
681 np = (int64_t)(D->n_records * (D->n_records - 1)) / 2;
682 distance = gmt_M_memory (GMT, NULL, np, double);
683 for (tbl = 0, np = 0; tbl < D->n_tables; tbl++) {
684 for (seg = 0; seg < D->table[tbl]->n_segments; seg++) {
685 S = D->table[tbl]->segment[seg];
686 for (row = 0; row < S->n_rows; row++) {
687 for (col = 0; col < n_cols; col++) P[np].data[col] = S->data[col][row]; /* Duplicate coordinates */
688 if (n_cols < 4) P[np].data[GMT_W] = 1.0; /* No weight provided, set to unity */
689 P[np].ID = (uint64_t)np; /* Assign ID based on input record # from 0 */
690 P[np].distance = DBL_MAX;
691 /* We split the loop over calculation of distance from the loop over assignments since
692 * if OpenMP is used then we cannot interchange k and np as there may be overprinting.
693 */
694 #ifdef _OPENMP
695 #pragma omp parallel for private(k) shared(np,GMT,P,distance)
696 #endif
697 for (k = 0; k < np; k++) { /* Get distances to other points */
698 distance[k] = gmt_distance (GMT, P[k].data[GMT_X], P[k].data[GMT_Y], P[np].data[GMT_X], P[np].data[GMT_Y]);
699 }
700 for (k = 0; k < np; k++) { /* Compare this point to all previous points */
701 if (distance[k] < P[k].distance) { /* Update shortest distance so far, and with which neighbor */
702 P[k].distance = distance[k];
703 P[k].neighbor = np;
704 }
705 if (distance[k] < P[np].distance) { /* Update shortest distance so far, and with which neighbor */
706 P[np].distance = distance[k];
707 P[np].neighbor = k;
708 }
709 }
710 np++;
711 }
712 }
713 }
714 gmt_M_free (GMT, distance);
715
716 /* Prefer mergesort since qsort is not stable for equalities */
717 mergesort (P, np, sizeof (struct NN_DIST), gmtspatial_compare_nn_points);
718
719 *n_points = (uint64_t)np;
720 #ifdef DEBUG
721 if (gmt_M_is_verbose (GMT, GMT_MSG_DEBUG)) {
722 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "===> Initialization\n");
723 for (k = 0; k < (int64_t)(*n_points); k++)
724 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "%6d\tID=%6d\tNeighbor=%6d\tDistance = %.12g\n", (int)k, (int)P[k].ID, P[k].neighbor, P[k].distance);
725 }
726 #endif
727 return (P);
728 }
729
gmtspatial_compare_nn_info(const void * point_1v,const void * point_2v)730 GMT_LOCAL int gmtspatial_compare_nn_info (const void *point_1v, const void *point_2v) {
731 /* Routine for qsort to sort NN rec numbers structure on original record order.
732 */
733 const struct NN_INFO *point_1 = point_1v, *point_2 = point_2v;
734
735 if (point_1->orig_rec < point_2->orig_rec) return (-1);
736 if (point_1->orig_rec > point_2->orig_rec) return (+1);
737 return (0);
738 }
739
gmtspatial_NNA_update_info(struct GMT_CTRL * GMT,struct NN_INFO * I,struct NN_DIST * NN_dist,uint64_t n_points)740 GMT_LOCAL struct NN_INFO *gmtspatial_NNA_update_info (struct GMT_CTRL *GMT, struct NN_INFO * I, struct NN_DIST *NN_dist, uint64_t n_points) {
741 /* Return revised array of NN ID lookups via sorting on neighbor IDs */
742 uint64_t k;
743 struct NN_INFO *info = (I) ? I : gmt_M_memory (GMT, NULL, n_points, struct NN_INFO);
744 for (k = 0; k < n_points; k++) {
745 info[k].sort_rec = k;
746 info[k].orig_rec = int64_abs (NN_dist[k].ID);
747 }
748
749 /* Prefer mergesort since qsort is not stable for equalities */
750 mergesort (info, n_points, sizeof (struct NN_INFO), gmtspatial_compare_nn_info);
751
752 /* Now, I[k].sort_rec will take the original record # k and return the corresponding record in the sorted array */
753 return (info);
754 }
755
usage(struct GMTAPI_CTRL * API,int level)756 static int usage (struct GMTAPI_CTRL *API, int level) {
757 const char *name = gmt_show_name_and_purpose (API, THIS_MODULE_LIB, THIS_MODULE_CLASSIC_NAME, THIS_MODULE_PURPOSE);
758 if (level == GMT_MODULE_PURPOSE) return (GMT_NOERROR);
759 GMT_Usage (API, 0, "usage: %s [<table>] [-A[a<min_dist>]] [-C] [-D[+a<amax>][+c|C<cmax>][+d<dmax>][+f<file>][+p][+s<sfact>]] [-E+n|p] "
760 "[-F[l]] [-I[i|e]] [-L%s/<noise>/<offset>] [-N<pfile>[+a][+p<ID>][+r][+z]] [-Q[+c<min>[/<max>]][+h][+l][+p][+s[a|d]]] [%s] "
761 "[-Sb<width>|h|i|j|s|u] [-T[<cpol>]] [%s] [%s] [%s] [%s] [%s] [%s] [%s] [%s] [%s] [%s] [%s] [%s] [%s] [%s] [%s]\n", name, GMT_DIST_OPT, GMT_Rgeo_OPT,
762 GMT_V_OPT, GMT_a_OPT, GMT_b_OPT, GMT_d_OPT, GMT_e_OPT, GMT_f_OPT, GMT_g_OPT, GMT_h_OPT, GMT_i_OPT, GMT_j_OPT, GMT_o_OPT, GMT_q_OPT, GMT_s_OPT, GMT_colon_OPT, GMT_PAR_OPT);
763
764 if (level == GMT_SYNOPSIS) return (GMT_MODULE_SYNOPSIS);
765
766 GMT_Message (API, GMT_TIME_NONE, " OPTIONAL ARGUMENTS:\n");
767 GMT_Option (API, "<");
768 GMT_Usage (API, 1, "\n-A[a<min_dist>]");
769 GMT_Usage (API, -2, "Nearest Neighbor (NN) Analysis. Compute minimum distances between NN point pairs. "
770 "Append unit used for NN distance calculation. Returns minimum distances and point IDs for pairs. "
771 "Use -Aa to replace close neighbor pairs with their weighted average location until "
772 "no point pair has a NN distance less than the specified <min_dist> distance [0]. "
773 "Considers 3rd column as z (if present) and 4th as w, if present [weight = 1].");
774 GMT_Usage (API, 1, "\n-C Clip polygons to the given region box (requires -R), possibly yielding new closed polygons. "
775 "For truncation instead (possibly yielding open polygons, i.e., lines), see -T.");
776 GMT_Usage (API, 1, "\n-D[+a<amax>][+c|C<cmax>][+d<dmax>][+f<file>][+p][+s<sfact>]");
777 GMT_Usage (API, -2, "Look for (near-)duplicate lines or polygons. Duplicate lines have a minimum point separation less than <dmax> and a closeness "
778 "ratio (mean separation/length) less than <cmax>. "
779 "If near-duplicates have lengths that differ by <sfact> or more then they are subsets or supersets. "
780 "To flag duplicate polygons, the fractional difference in areas must be less than <amax>. "
781 "By default, we consider all points when comparing two lines. Change these settings via modifiers:");
782 GMT_Usage (API, 3, "+a Set limit of fractional difference of polygon areas [0.01].");
783 GMT_Usage (API, 3, "+c Set closeness ratio (mean separation/length) [0.01]. "
784 "Use +C to use median separation instead.");
785 GMT_Usage (API, 3, "+d Set minimum mean point separation [0].");
786 GMT_Usage (API, 3, "+f Compare <table> against <file> instead of itself.");
787 GMT_Usage (API, 3, "+p limit comparison to points that project perpendicularly on to the other line [all points].");
788 GMT_Usage (API, 3, "+s Set length ratio difference threshold [2].");
789 GMT_Usage (API, 1, "\n-E+n|p");
790 GMT_Usage (API, -2, "Orient all polygons to have the same handedness:");
791 GMT_Usage (API, 3, "+n Impose clockwise (negative) handedness.");
792 GMT_Usage (API, 3, "+p Impose counter-clockwise (positive) handedness.");
793 GMT_Usage (API, 1, "\n-F[l]");
794 GMT_Usage (API, -2, "Force all input segments to become closed polygons on output by adding repeated point if needed. "
795 "Use -Fl instead to ensure input lines are not treated as polygons.");
796 GMT_Usage (API, 1, "\n-I[i|e]");
797 GMT_Usage (API, -2, "Compute Intersection locations between input polygon(s). Optionally append a directive:");
798 GMT_Usage (API, 3, "e: Compute external crossings only [Default is both].");
799 GMT_Usage (API, 3, "i: Compute internal crossings only [Default is both].");
800 GMT_Usage (API, -2, "Use uppercase E or I to consider all segments within the same table as one entity [separate].");
801 GMT_Usage (API, 1, "\n-L%s/<noise>/<offset>", GMT_DIST_OPT);
802 GMT_Usage (API, -2, "Remove tile Lines. These are superfluous lines along the -R border. "
803 "Append <dist> (in m) [0], coordinate noise [1e-10], and max offset from gridline [1e-10].");
804 GMT_Usage (API, 1, "\n-N<pfile>[+a][+p<ID>][+r][+z]");
805 GMT_Usage (API, -2, "Determine ID of polygon (in <pfile>) enclosing each input feature. The ID is set as follows:");
806 GMT_Usage (API, 3, "%s If OGR/GMT polygons, get polygon ID via -a for Z column, else", GMT_LINE_BULLET);
807 GMT_Usage (API, 3, "%s Interpret segment labels (-Z<value>) as polygon IDs, else", GMT_LINE_BULLET);
808 GMT_Usage (API, 3, "%s Interpret segment labels (-L<label>) as polygon IDs, else", GMT_LINE_BULLET);
809 GMT_Usage (API, 3, "%s Append +p<ID> to set origin for auto-incrementing polygon IDs [0].", GMT_LINE_BULLET);
810 GMT_Usage (API, -2, "Additional modifiers are available:");
811 GMT_Usage (API, 3, "+a All points of a feature (line, polygon) must be inside the ID polygon [mid point].");
812 GMT_Usage (API, 3, "+r No table output; just reports which polygon a feature is inside.");
813 GMT_Usage (API, 3, "+z Append the ID as a new output data column [Default adds -Z<ID> to segment header].");
814 GMT_Usage (API, 1, "\n-Q[+c<min>[/<max>]][+h][+l][+p][+s[a|d]]");
815 GMT_Usage (API, -2, "Measure area and handedness of polygon(s) or length of line segments. If -fg is used "
816 "you may append unit %s [k]; otherwise it will be based on the input Cartesian data unit. "
817 "We also compute polygon centroid or line mid-point. See documentation for more information. Optional modifiers:", GMT_LEN_UNITS_DISPLAY);
818 GMT_Usage (API, 3, "+c Limit output segments to those with area or length within specified range [output all segments]. "
819 "If <max> is not given then it defaults to infinity.");
820 GMT_Usage (API, 3, "+h Place the (area, handedness) or length result in the segment header on output.");
821 GMT_Usage (API, 3, "+p Consider all input as polygons and close them if necessary "
822 "[only closed polygons are considered polygons].");
823 GMT_Usage (API, 3, "+l Consider all input as lines even if closed [closed polygons are considered polygons].");
824 GMT_Usage (API, 3, "+s Sort segments based on area or length; append a for ascending or d for descending [ascending].");
825 GMT_Usage (API, -2, "[Default only reports results to standard output].\n");
826 GMT_Option (API, "R");
827 GMT_Usage (API, 1, "\n-Sb<width>|h|i|j|s|u");
828 GMT_Usage (API, -2, "Spatial manipulation of polygons; choose a directive:");
829 #ifdef HAVE_GEOS
830 GMT_Usage (API, 3, "b: Compute buffer polygon around line/polygon. Append <width> of buffer zone. "
831 "Note: this is a purely Cartesian operation so <width> must be in data units.");
832 #endif
833 GMT_Usage (API, 3, "h: Detect holes and reverse them relative to perimeters.");
834 GMT_Usage (API, 3, "i: Find intersection [Not implemented yet].");
835 GMT_Usage (API, 3, "j: Join polygons that were split by the Dateline [Not implemented yet].");
836 GMT_Usage (API, 3, "s: Split polygons that straddle the Dateline.");
837 GMT_Usage (API, 3, "u: Find union [Not implemented yet].");
838 GMT_Usage (API, 1, "\n-T[<cpol>]");
839 GMT_Usage (API, -2, "Truncate polygons against the clip polygon <cpol>; if <cpol> is not given we require -R "
840 "and clip against a polygon derived from the region border.");
841 GMT_Option (API, "V,a,bi2,bo,d,e,f,g,h,i,j,o,q,s,:,.");
842
843 return (GMT_MODULE_USAGE);
844 }
845
parse(struct GMT_CTRL * GMT,struct GMTSPATIAL_CTRL * Ctrl,struct GMT_OPTION * options)846 static int parse (struct GMT_CTRL *GMT, struct GMTSPATIAL_CTRL *Ctrl, struct GMT_OPTION *options) {
847
848 /* This parses the options provided to grdsample and sets parameters in CTRL.
849 * Any GMT common options will override values set previously by other commands.
850 * It also replaces any file names specified as input or output with the data ID
851 * returned when registering these sources/destinations with the API.
852 */
853
854 unsigned int n_files[2] = {0, 0}, pos, n_errors = 0;
855 int n;
856 char txt_a[GMT_LEN64] = {""}, txt_b[GMT_LEN64] = {""}, txt_c[GMT_LEN64] = {""}, p[GMT_LEN256] = {""}, *s = NULL;
857 struct GMT_OPTION *opt = NULL;
858 struct GMTAPI_CTRL *API = GMT->parent;
859
860 if (gmt_M_is_geographic (GMT, GMT_IN)) Ctrl->Q.unit = 'k'; /* Default geographic distance unit is km */
861
862 for (opt = options; opt; opt = opt->next) {
863 switch (opt->option) {
864
865 case '<': /* Skip input files */
866 if (GMT_Get_FilePath (API, GMT_IS_DATASET, GMT_IN, GMT_FILE_REMOTE, &(opt->arg))) n_errors++;;
867 n_files[GMT_IN]++;
868 break;
869 case '>': /* Got named output file */
870 if (n_files[GMT_OUT]++ == 0) Ctrl->Out.file = strdup (opt->arg);
871 break;
872
873 /* Processes program-specific parameters */
874
875 case 'A': /* Do nearest neighbor analysis */
876 n_errors += gmt_M_repeated_module_option (API, Ctrl->A.active);
877 Ctrl->A.active = true;
878 if (opt->arg[0] == 'a' || opt->arg[0] == 'A') { /* Spatially average points until minimum NN distance is less than given distance */
879 Ctrl->A.mode = (opt->arg[0] == 'A') ? 2 : 1; /* Slow mode is an undocumented test mode */
880 Ctrl->A.smode = gmt_get_distance (GMT, &opt->arg[1], &(Ctrl->A.min_dist), &(Ctrl->A.unit));
881 }
882 else if (((opt->arg[0] == '-' || opt->arg[0] == '+') && strchr (GMT_LEN_UNITS, opt->arg[1])) || (opt->arg[0] && strchr (GMT_LEN_UNITS, opt->arg[0]))) {
883 /* Just compute NN distances and return the unique ones */
884 if (opt->arg[0] == '-') /* Flat earth calculation */
885 Ctrl->A.smode = 1;
886 else if (opt->arg[0] == '+') /* Geodetic */
887 Ctrl->A.smode = 3;
888 else
889 Ctrl->A.smode = 2; /* Great circle */
890 Ctrl->A.unit = (Ctrl->A.smode == 2) ? opt->arg[0] : opt->arg[1];
891 if (gmt_M_is_cartesian (GMT, GMT_IN)) { /* Data was not geographic, revert to Cartesian settings */
892 Ctrl->A.smode = 0;
893 Ctrl->A.unit = 'X';
894 }
895 }
896 else if (opt->arg[0]) {
897 GMT_Report (API, GMT_MSG_ERROR, "Bad modifier in %c in -A option\n", opt->arg[0]);
898 n_errors++;
899 }
900 break;
901 case 'C': /* Clip to given region */
902 n_errors += gmt_M_repeated_module_option (API, Ctrl->C.active);
903 Ctrl->C.active = true;
904 break;
905 case 'D': /* Look for duplications */
906 n_errors += gmt_M_repeated_module_option (API, Ctrl->D.active);
907 Ctrl->D.active = true;
908 pos = 0;
909 while (gmt_strtok (opt->arg, "+", &pos, p)) {
910 switch (p[0]) {
911 case 'a': /* Gave a new +a<dmax> value */
912 GMT_Report (API, GMT_MSG_ERROR, "+a not implemented yet\n");
913 Ctrl->D.I.a_threshold = atof (&p[1]);
914 break;
915 case 'd': /* Gave a new +d<dmax> value */
916 Ctrl->D.mode = gmt_get_distance (GMT, &p[1], &(Ctrl->D.I.d_threshold), &(Ctrl->D.unit));
917 break;
918 case 'C': /* Gave a new +C<cmax> value */
919 Ctrl->D.I.mode = 1; /* Median instead of mean */
920 /* Intentionally fall through */
921 case 'c': /* Gave a new +c<cmax> value */
922 if (p[1]) Ctrl->D.I.c_threshold = atof (&p[1]); /* This allows +C by itself just to change to median */
923 break;
924 case 's': /* Gave a new +s<fact> value */
925 Ctrl->D.I.s_threshold = atof (&p[1]);
926 break;
927 case 'p': /* Consider only inside projections */
928 Ctrl->D.I.inside = 1;
929 break;
930 case 'f': /* Gave a file name */
931 Ctrl->D.file = strdup (&p[1]);
932 if (GMT_Get_FilePath (API, GMT_IS_DATASET, GMT_IN, GMT_FILE_REMOTE, &(Ctrl->D.file))) n_errors++;
933 break;
934 }
935 }
936 break;
937 case 'E': /* Orient polygons -E+n|p (old -E-|+) */
938 n_errors += gmt_M_repeated_module_option (API, Ctrl->E.active);
939 Ctrl->E.active = true;
940 if (opt->arg[0] == '-' || strstr (opt->arg, "+n"))
941 Ctrl->E.mode = GMT_POL_IS_CW;
942 else if (opt->arg[0] == '+' || strstr (opt->arg, "+p"))
943 Ctrl->E.mode = GMT_POL_IS_CCW;
944 else
945 n_errors++;
946 break;
947 case 'F': /* Force polygon or line mode */
948 n_errors += gmt_M_repeated_module_option (API, Ctrl->F.active);
949 Ctrl->F.active = true;
950 Ctrl->F.geometry = (opt->arg[0] == 'l') ? GMT_IS_LINE : GMT_IS_POLY;
951 break;
952 case 'I': /* Compute intersections between polygons */
953 n_errors += gmt_M_repeated_module_option (API, Ctrl->I.active);
954 Ctrl->I.active = true;
955 if (opt->arg[0] == 'i') Ctrl->I.mode = 1;
956 if (opt->arg[0] == 'I') Ctrl->I.mode = 2;
957 if (opt->arg[0] == 'e') Ctrl->I.mode = 4;
958 if (opt->arg[0] == 'E') Ctrl->I.mode = 8;
959 break;
960 case 'L': /* Remove tile lines */
961 n_errors += gmt_M_repeated_module_option (API, Ctrl->L.active);
962 Ctrl->L.active = true;
963 n = sscanf (opt->arg, "%[^/]/%[^/]/%s", txt_a, txt_b, txt_c);
964 if (n >= 1) Ctrl->L.s_cutoff = atof (txt_a);
965 if (n >= 2) Ctrl->L.path_noise = atof (txt_b);
966 if (n == 3) Ctrl->L.box_offset = atof (txt_c);
967 break;
968 case 'N': /* Determine containing polygons for features */
969 n_errors += gmt_M_repeated_module_option (API, Ctrl->N.active);
970 Ctrl->N.active = true;
971 if ((s = strchr (opt->arg, '+')) == NULL) { /* No modifiers */
972 Ctrl->N.file = strdup (opt->arg);
973 continue;
974 }
975 s[0] = '\0'; Ctrl->N.file = strdup (opt->arg); s[0] = '+';
976 pos = 0;
977 while (gmt_strtok (s, "+", &pos, p)) {
978 switch (p[0]) {
979 case 'a': /* All points must be inside polygon */
980 Ctrl->N.all = true;
981 break;
982 case 'p': /* Set start of running numbers [0] */
983 Ctrl->N.ID = (p[1]) ? atoi (&p[1]) : 1;
984 break;
985 case 'r': /* Just give a report */
986 Ctrl->N.mode = 1;
987 break;
988 case 'z': /* Gave a new +s<fact> value */
989 Ctrl->N.mode = 2;
990 break;
991 }
992 }
993 break;
994 case 'Q': /* Measure area/length and handedness of polygons */
995 n_errors += gmt_M_repeated_module_option (API, Ctrl->Q.active);
996 Ctrl->Q.active = true;
997 s = opt->arg;
998 /* Handle +sa|d versus +s for deprecated ellipsoidal arc seconds */
999 if (s[0] && !strcmp (s, "+s") && (s[2] == '\0' || !strchr ("ad", s[2]))) { /* Since [-|+] is deprecated as of GMT 6 we must assume this is +s[a] */
1000 GMT_Report (API, GMT_MSG_WARNING, "-Q+s is assumed to mean -A+sa. If you meant ellipsoidal arc second distances then use -Qs -je\n");
1001 }
1002 else if (s[0] && strchr ("-+", s[0]) && strchr (GMT_LEN_UNITS, s[1])) { /* Since [-|+] is deprecated as of GMT 6 */
1003 if (gmt_M_compat_check (GMT, 6))
1004 GMT_Report (API, GMT_MSG_COMPAT, "Leading -|+ with unit to set flat Earth or ellipsoidal mode is deprecated; use -j<mode> instead\n");
1005 else {
1006 GMT_Report (API, GMT_MSG_ERROR, "Signed unit is not allowed - ignored\n");
1007 n_errors++;
1008 }
1009
1010 }
1011 if (s[0] == '-' && strchr (GMT_LEN_UNITS, s[1])) { /* Flat earth distances [deprecated; use -j instead] */
1012 Ctrl->Q.dmode = 1; Ctrl->Q.unit = s[1]; s += 2;
1013 }
1014 else if (s[0] == '+' && s[1] == 's' && (s[2] == '\0' || !strchr ("ad", s[2]))) { /* Geodesic distances using arc sec [deprecated; use -j instead] */
1015 Ctrl->Q.dmode = 3; Ctrl->Q.unit = s[1]; s += 2;
1016 }
1017 else if (s[0] == '+' && strchr ("dmefkMnu", s[1])) { /* Geodesic distances (except arc second) [deprecated; use -j instead] */
1018 Ctrl->Q.dmode = 3; Ctrl->Q.unit = s[1]; s += 2;
1019 }
1020 else if (s[0] && strchr (GMT_LEN_UNITS, s[0])) { /* Great circle distances */
1021 Ctrl->Q.dmode = 2; Ctrl->Q.unit = s[0]; s++;
1022 }
1023 pos = 0;
1024 while (gmt_strtok (s, "+", &pos, p)) {
1025 switch (p[0]) {
1026 case 'c': /* Set limits on output segments based on lengths or area */
1027 Ctrl->Q.area = true;
1028 n = sscanf (&p[1], "%[^/]/%s", txt_a, txt_b);
1029 Ctrl->Q.limit[0] = atof (txt_a);
1030 if (n == 1) /* Just got the minimum cutoff */
1031 Ctrl->Q.limit[1] = DBL_MAX;
1032 else
1033 Ctrl->Q.limit[1] = atof (txt_b);
1034 break;
1035 case 'l': /* Consider input as lines, even if closed */
1036 Ctrl->Q.mode = GMT_IS_LINE;
1037 break;
1038 case 'p': /* Consider input as polygones, close if necessary */
1039 Ctrl->Q.mode = GMT_IS_POLY;
1040 break;
1041 case 'h': /* Place result in output header */
1042 Ctrl->Q.header = true;
1043 break;
1044 case 's': /* Sort segments on output based on dimension */
1045 Ctrl->Q.sort = true;
1046 switch (p[1]) { /* Ascending or descending sort */
1047 case 'a': case '\0': Ctrl->Q.dir = +1; break; /* Ascending [Default] */
1048 case 'd': Ctrl->Q.dir = -1; break; /* Descending */
1049 default:
1050 GMT_Report (API, GMT_MSG_ERROR, "Option -Q: Unrecognized direction %c given to modifier +s\n", p[1]);
1051 n_errors++;
1052 break;
1053 }
1054 break;
1055 default:
1056 GMT_Report (API, GMT_MSG_ERROR, "Option -Q: Unrecognized modifier +%c\n", p[0]);
1057 n_errors++;
1058 break;
1059 }
1060 }
1061 if (strstr (s, "++") || (s[0] && s[strlen(s)-1] == '+')) { /* Deal with the old-style single "+" to mean header */
1062 Ctrl->Q.header = true;
1063 GMT_Report (API, GMT_MSG_WARNING, "Option -Q+ is interpreted as -Q+h\n");
1064 }
1065 break;
1066 case 'S': /* Spatial polygon operations */
1067 n_errors += gmt_M_repeated_module_option (API, Ctrl->S.active);
1068 Ctrl->S.active = true;
1069 if (opt->arg[0] == 'u') {
1070 Ctrl->S.mode = POL_UNION;
1071 GMT_Report (API, GMT_MSG_ERROR, "Option -Su not implemented yet\n");
1072 }
1073 #ifdef HAVE_GEOS
1074 else if (opt->arg[0] == 'b') {
1075 Ctrl->S.mode = POL_BUFFER;
1076 Ctrl->S.width = atof (&opt->arg[1]);
1077 if (isnan(Ctrl->S.width) || Ctrl->S.width <= 0) {
1078 GMT_Report (API, GMT_MSG_ERROR, "Option -Sb<val> must provide a width > 0\n");
1079 n_errors++;
1080 }
1081 }
1082 else if (opt->arg[0] == 'c')
1083 Ctrl->S.mode = POL_CENTROID;
1084 #endif
1085 else if (opt->arg[0] == 'i') {
1086 Ctrl->S.mode = POL_INTERSECTION;
1087 GMT_Report (API, GMT_MSG_ERROR, "Option -Si not implemented yet\n");
1088 }
1089 else if (opt->arg[0] == 's')
1090 Ctrl->S.mode = POL_SPLIT;
1091 else if (opt->arg[0] == 'h')
1092 Ctrl->S.mode = POL_HOLE;
1093 else if (opt->arg[0] == 'j') {
1094 Ctrl->S.mode = POL_JOIN;
1095 GMT_Report (API, GMT_MSG_ERROR, "Option -Sj not implemented yet\n");
1096 }
1097 else
1098 n_errors++;
1099 break;
1100 case 'T': /* Truncate against polygon */
1101 n_errors += gmt_M_repeated_module_option (API, Ctrl->T.active);
1102 Ctrl->T.active = true;
1103 Ctrl->C.active = Ctrl->S.active = true;
1104 Ctrl->S.mode = POL_CLIP;
1105 if (opt->arg[0]) Ctrl->T.file = strdup (opt->arg);
1106 break;
1107 default: /* Report bad options */
1108 n_errors += gmt_default_error (GMT, opt->option);
1109 break;
1110 }
1111 }
1112
1113 if (Ctrl->E.active) Ctrl->Q.active = true;
1114
1115 if (GMT->common.b.active[GMT_IN] && GMT->common.b.ncol[GMT_IN] == 0) GMT->common.b.ncol[GMT_IN] = 2;
1116 n_errors += gmt_M_check_condition (GMT, GMT->common.b.active[GMT_IN] && GMT->common.b.ncol[GMT_IN] < 2, "Binary input data (-bi) must have at least %d columns\n", 2);
1117 n_errors += gmt_M_check_condition (GMT, Ctrl->S.mode == POL_CLIP && !Ctrl->T.file && !GMT->common.R.active[RSET], "Option -T without a polygon requires -R\n");
1118 n_errors += gmt_M_check_condition (GMT, Ctrl->C.active && !Ctrl->T.active && !GMT->common.R.active[RSET], "Option -C requires -R\n");
1119 n_errors += gmt_M_check_condition (GMT, Ctrl->L.active && !GMT->common.R.active[RSET], "Option -L requires -R\n");
1120 n_errors += gmt_M_check_condition (GMT, Ctrl->L.active && Ctrl->L.s_cutoff < 0.0, "Option -L requires a positive cutoff in meters\n");
1121 n_errors += gmt_M_check_condition (GMT, Ctrl->D.active && Ctrl->D.file && gmt_access (GMT, Ctrl->D.file, R_OK), "Option -D: Cannot read file %s!\n", Ctrl->D.file);
1122 n_errors += gmt_M_check_condition (GMT, Ctrl->T.active && Ctrl->T.file && gmt_access (GMT, Ctrl->T.file, R_OK), "Option -T: Cannot read file %s!\n", Ctrl->T.file);
1123 n_errors += gmt_M_check_condition (GMT, n_files[GMT_OUT] > 1, "Only one output destination can be specified\n");
1124
1125 return (n_errors ? GMT_PARSE_ERROR : GMT_NOERROR);
1126 }
1127
1128 #define bailout(code) {gmt_M_free_options (mode); return (code);}
1129 #define Return(code) {Free_Ctrl (GMT, Ctrl); gmt_end_module (GMT, GMT_cpy); bailout (code);}
1130
GMT_gmtspatial(void * V_API,int mode,void * args)1131 EXTERN_MSC int GMT_gmtspatial (void *V_API, int mode, void *args) {
1132 int error = 0;
1133 unsigned int geometry = GMT_IS_POLY, internal = 0, external = 0, smode = GMT_NO_STRINGS;
1134
1135 static char *kind[2] = {"CCW", "CW"};
1136
1137 double out[GMT_MAX_COLUMNS];
1138
1139 struct GMT_DATASET *D = NULL;
1140 struct GMT_DATASEGMENT *S = NULL;
1141 struct GMT_RECORD Out;
1142 struct GMTSPATIAL_CTRL *Ctrl = NULL;
1143 struct GMT_CTRL *GMT = NULL, *GMT_cpy = NULL;
1144 struct GMT_OPTION *options = NULL;
1145 struct GMTAPI_CTRL *API = gmt_get_api_ptr (V_API); /* Cast from void to GMTAPI_CTRL pointer */
1146
1147 /*----------------------- Standard module initialization and parsing ----------------------*/
1148
1149 if (API == NULL) return (GMT_NOT_A_SESSION);
1150 if (mode == GMT_MODULE_PURPOSE) return (usage (API, GMT_MODULE_PURPOSE)); /* Return the purpose of program */
1151 options = GMT_Create_Options (API, mode, args); if (API->error) return (API->error); /* Set or get option list */
1152
1153 if ((error = gmt_report_usage (API, options, 0, usage)) != GMT_NOERROR) bailout (error); /* Give usage if requested */
1154
1155 /* Parse the command-line arguments */
1156
1157 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 */
1158 if (GMT_Parse_Common (API, THIS_MODULE_OPTIONS, options)) Return (API->error);
1159 Ctrl = New_Ctrl (GMT); /* Allocate and initialize a new control structure */
1160 if ((error = parse (GMT, Ctrl, options)) != 0) Return (error);
1161
1162 /*---------------------------- This is the gmtspatial main code ----------------------------*/
1163
1164 if (Ctrl->I.active) {
1165 switch (Ctrl->I.mode) {
1166 case 0:
1167 internal = external = 1;
1168 break;
1169 case 1:
1170 internal = 1;
1171 break;
1172 case 2:
1173 internal = 2;
1174 break;
1175 case 4:
1176 external = 1;
1177 break;
1178 case 8:
1179 external = 2;
1180 break;
1181 }
1182 }
1183
1184 /* Read input data set */
1185
1186 if (Ctrl->D.active) geometry = GMT_IS_LINE|GMT_IS_POLY; /* May be lines, may be polygons... */
1187 else if (Ctrl->Q.active) geometry = Ctrl->Q.mode; /* May be lines, may be polygons... */
1188 else if (Ctrl->A.active) geometry = GMT_IS_POINT; /* NN analysis involves points */
1189 else if (Ctrl->F.active) geometry = Ctrl->F.geometry; /* Forcing polygon or line mode */
1190 else if (Ctrl->S.active) geometry = GMT_IS_POLY; /* Forcing polygon mode */
1191 if (GMT_Init_IO (API, GMT_IS_DATASET, geometry, GMT_IN, GMT_ADD_DEFAULT, 0, options) != GMT_NOERROR) { /* Registers default input sources, unless already set */
1192 Return (API->error);
1193 }
1194 GMT_Report (API, GMT_MSG_INFORMATION, "Processing input table data\n");
1195 if ((D = GMT_Read_Data (API, GMT_IS_DATASET, GMT_IS_FILE, 0, GMT_READ_NORMAL, NULL, NULL, NULL)) == NULL) {
1196 Return (API->error);
1197 }
1198 if (D->n_columns < 2) {
1199 GMT_Report (API, GMT_MSG_ERROR, "Input data have %d column(s) but at least 2 are needed\n", (int)D->n_columns);
1200 Return (GMT_DIM_TOO_SMALL);
1201 }
1202
1203 /* Allocate memory and read in all the files; each file can have many lines (-m) */
1204
1205 if (D->n_records == 0) { /* Empty files, nothing to do */
1206 GMT_Report (API, GMT_MSG_WARNING, "No data records found.\n");
1207 if (GMT_Destroy_Data (API, &D) != GMT_NOERROR) {
1208 Return (API->error);
1209 }
1210 Return (GMT_NOERROR);
1211 }
1212
1213 if (Ctrl->S.active && !(Ctrl->S.mode == POL_SPLIT || Ctrl->S.mode == POL_HOLE || Ctrl->S.mode == POL_BUFFER || Ctrl->S.mode == POL_CENTROID))
1214 external = 1;
1215
1216 if (gmt_init_distaz (GMT, 'X', 0, GMT_MAP_DIST) == GMT_NOT_A_VALID_TYPE) /* Use Cartesian calculations and user units */
1217 Return (GMT_NOT_A_VALID_TYPE);
1218
1219 gmt_set_inside_mode (GMT, NULL, (gmt_M_is_geographic (GMT, GMT_IN)) ? GMT_IOO_SPHERICAL : GMT_IOO_CARTESIAN);
1220
1221 /* OK, with data in hand we can do some damage */
1222
1223 if (Ctrl->A.active) { /* Nearest neighbor analysis. We compute distances between all point pairs and sort on minimum distance */
1224 uint64_t n_points, k, a, b, n, col, n_pairs;
1225 double A[3], B[3], w, iw, d_bar, out[7];
1226 struct NN_DIST *NN_dist = NULL;
1227 struct NN_INFO *NN_info = NULL;
1228
1229 if (gmt_init_distaz (GMT, Ctrl->A.unit, Ctrl->A.smode, GMT_MAP_DIST) == GMT_NOT_A_VALID_TYPE) /* Set the unit and distance calculation we requested */
1230 Return (GMT_NOT_A_VALID_TYPE);
1231
1232 NN_dist = gmtspatial_NNA_init_dist (GMT, D, &n_points); /* Return array of NN results sorted on smallest distances */
1233 NN_info = gmtspatial_NNA_update_info (GMT, NN_info, NN_dist, n_points); /* Return array of NN ID record look-ups */
1234 if (GMT_Destroy_Data (API, &D) != GMT_NOERROR) { /* All data now in NN_dist so free original dataset */
1235 gmt_M_free (GMT, NN_dist); gmt_M_free (GMT, NN_info);
1236 Return (API->error);
1237 }
1238 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 */
1239 gmt_M_free (GMT, NN_dist); gmt_M_free (GMT, NN_info);
1240 Return (API->error);
1241 }
1242 if (GMT_Begin_IO (API, GMT_IS_DATASET, GMT_OUT, GMT_HEADER_ON) != GMT_NOERROR) { /* Enables data output and sets access mode */
1243 gmt_M_free (GMT, NN_dist); gmt_M_free (GMT, NN_info);
1244 Return (API->error);
1245 }
1246 if (GMT_Set_Geometry (API, GMT_OUT, GMT_IS_POINT) != GMT_NOERROR) { /* Sets output geometry */
1247 gmt_M_free (GMT, NN_dist); gmt_M_free (GMT, NN_info);
1248 Return (API->error);
1249 }
1250 Out.data = out; Out.text = NULL;
1251 if (Ctrl->A.mode) { /* Need to combine close neighbors until minimum distance >= min_dist, then output revised dataset */
1252 GMT_Report (API, GMT_MSG_INFORMATION, "NNA using min separation of %g%c\n", Ctrl->A.min_dist, Ctrl->A.unit);
1253 n = 0;
1254 while (n < n_points && NN_dist[n].distance < Ctrl->A.min_dist) n++; /* Find # of pairs that are too close together */
1255 while (n) { /* Must do more combining since n pairs exceed threshold distance */
1256 if (Ctrl->A.mode == 2) {
1257 GMT_Report (API, GMT_MSG_INFORMATION, "Slow mode: Replace the single closest pair with its weighted average, then redo NNA\n");
1258 n = 1;
1259 }
1260 for (k = n_pairs = 0; k < n; k++) { /* Loop over pairs that are too close */
1261 if (gmt_M_is_dnan (NN_dist[k].distance)) continue; /* Already processed */
1262 a = k; /* The current point */
1263 b = NN_info[int64_abs(NN_dist[a].neighbor)].sort_rec; /* a's neighbor location in the sorted NN_dist array */
1264 GMT_Report (API, GMT_MSG_DEBUG, "Replace pair %" PRIu64 " and %" PRIu64 " with its weighted average location\n", a, b);
1265 w = NN_dist[a].data[GMT_W] + NN_dist[b].data[GMT_W]; /* Weight sum */
1266 iw = 1.0 / w; /* Inverse weight for scaling */
1267 /* Compute weighted average z */
1268 NN_dist[a].data[GMT_Z] = iw * (NN_dist[a].data[GMT_Z] * NN_dist[a].data[GMT_W] + NN_dist[b].data[GMT_Z] * NN_dist[b].data[GMT_W]);
1269 if (gmt_M_is_geographic (GMT, GMT_IN)) { /* Must do vector averaging */
1270 gmt_geo_to_cart (GMT, NN_dist[a].data[GMT_Y], NN_dist[a].data[GMT_X], A, true);
1271 gmt_geo_to_cart (GMT, NN_dist[b].data[GMT_Y], NN_dist[b].data[GMT_X], B, true);
1272 for (col = 0; col < 3; col++) { /* Get weighted average vector */
1273 A[col] = iw * (A[col] * NN_dist[a].data[GMT_W] + B[col] * NN_dist[b].data[GMT_W]);
1274 }
1275 gmt_normalize3v (GMT, A); /* Get unit vector */
1276 gmt_cart_to_geo (GMT, &NN_dist[a].data[GMT_Y], &NN_dist[a].data[GMT_X], A, true); /* Get lon, lat of the average point */
1277 }
1278 else { /* Cartesian weighted averaging */
1279 NN_dist[a].data[GMT_X] = iw * (NN_dist[a].data[GMT_X] * NN_dist[a].data[GMT_W] + NN_dist[b].data[GMT_X] * NN_dist[b].data[GMT_W]);
1280 NN_dist[a].data[GMT_Y] = iw * (NN_dist[a].data[GMT_Y] * NN_dist[a].data[GMT_W] + NN_dist[b].data[GMT_Y] * NN_dist[b].data[GMT_W]);
1281 }
1282 NN_dist[a].data[GMT_W] = 0.5 * w; /* Replace with the average weight */
1283 NN_dist[a].ID = -int64_abs (NN_dist[a].ID); /* Negative means it was averaged with other points */
1284 NN_dist[b].distance = GMT->session.d_NaN; /* Flag this point as used. gmtspatial_NNA_update_dist will sort it and place all NaNs at the end */
1285 n_pairs++;
1286 }
1287 GMT_Report (API, GMT_MSG_INFORMATION, "NNA Found %" PRIu64 " points, %" PRIu64 " pairs were too close and were replaced by their weighted average\n", n_points, n_pairs);
1288 NN_dist = gmtspatial_NNA_update_dist (GMT, NN_dist, &n_points); /* Return recomputed array of NN NN_dist sorted on smallest distances */
1289 NN_info = gmtspatial_NNA_update_info (GMT, NN_info, NN_dist, n_points); /* Return resorted array of NN ID lookups */
1290 n = 0;
1291 while (n < n_points && NN_dist[n].distance < Ctrl->A.min_dist) n++; /* Any more pairs with distances less than the threshold? */
1292 }
1293 if ((error = GMT_Set_Columns (API, GMT_OUT, 7, GMT_COL_FIX_NO_TEXT)) != 0) {
1294 gmt_M_free (GMT, NN_dist); gmt_M_free (GMT, NN_info);
1295 Return (error);
1296 }
1297 if (GMT->common.h.add_colnames) {
1298 char header[GMT_LEN256] = {""}, *name[2][2] = {{"x", "y"}, {"lon", "lat"}};
1299 k = (gmt_M_is_geographic (GMT, GMT_IN)) ? 1 : 0;
1300 sprintf (header, "%s[0]\t%s[1]\tz[2]\tweight[3]\tNN_dist[4]\tID[5]\tNN_ID[6]", name[k][GMT_X], name[k][GMT_Y]);
1301 GMT_Put_Record (API, GMT_WRITE_TABLE_HEADER, header);
1302 }
1303 }
1304 else { /* Just output the NN analysis results */
1305 gmt_set_cartesian (GMT, GMT_OUT); /* Since we are not writing coordinates */
1306 if ((error = GMT_Set_Columns (API, GMT_OUT, 3, GMT_COL_FIX_NO_TEXT)) != 0) {
1307 gmt_M_free (GMT, NN_dist); gmt_M_free (GMT, NN_info);
1308 Return (error);
1309 }
1310 if (GMT->common.h.add_colnames) {
1311 char header[GMT_LEN64];
1312 sprintf (header, "NN_dist[0]\tID[1]\tNN_ID[2]");
1313 GMT_Put_Record (API, GMT_WRITE_TABLE_HEADER, header);
1314 }
1315 }
1316 d_bar = 0.0;
1317 for (k = 0; k < n_points; k++) {
1318 d_bar += NN_dist[k].distance;
1319 if (Ctrl->A.mode) { /* Want all points out, including their coordinates */
1320 gmt_M_memcpy (out, NN_dist[k].data, 4, double);
1321 col = 4;
1322 }
1323 else { /* Just get the unique NN distances and the points involved */
1324 if (k && NN_dist[k].ID == NN_dist[k-1].neighbor && NN_dist[k].neighbor == NN_dist[k-1].ID) continue; /* Skip duplicate pairs */
1325 col = 0;
1326 }
1327 out[col++] = NN_dist[k].distance;
1328 if (NN_dist[k].ID < NN_dist[k].neighbor) {
1329 out[col++] = (double)NN_dist[k].ID;
1330 out[col++] = (double)NN_dist[k].neighbor;
1331 }
1332 else {
1333 out[col++] = (double)NN_dist[k].neighbor;
1334 out[col++] = (double)NN_dist[k].ID;
1335 }
1336 GMT_Put_Record (API, GMT_WRITE_DATA, &Out); /* Write points of NN info to stdout */
1337 }
1338 if (gmt_M_is_verbose (GMT, GMT_MSG_INFORMATION)) {
1339 d_bar /= n_points;
1340 if (GMT->common.R.active[RSET]) {
1341 int geo = gmt_M_is_geographic (GMT, GMT_IN) ? 1 : 0;
1342 double info[3], d_expect, R_index;
1343 struct GMT_DATASEGMENT *S = GMT_Alloc_Segment (GMT->parent, GMT_NO_STRINGS, 5, 2, NULL, NULL);
1344 S->data[GMT_X][0] = S->data[GMT_X][3] = S->data[GMT_X][4] = GMT->common.R.wesn[XLO];
1345 S->data[GMT_X][1] = S->data[GMT_X][2] = GMT->common.R.wesn[XHI];
1346 S->data[GMT_Y][0] = S->data[GMT_Y][1] = S->data[GMT_Y][4] = GMT->common.R.wesn[YLO];
1347 S->data[GMT_Y][2] = S->data[GMT_Y][3] = GMT->common.R.wesn[YHI];
1348 (void)gmtspatial_area_size (GMT, S->data[GMT_X], S->data[GMT_Y], S->n_rows, info, geo);
1349 gmt_free_segment (GMT, &S);
1350 d_expect = 0.5 * sqrt (info[GMT_Z]/n_points);
1351 R_index = d_bar / d_expect;
1352 GMT_Report (API, GMT_MSG_INFORMATION, "NNA Found %" PRIu64 " points, D_bar = %g, D_expect = %g, Spatial index = %g\n", n_points, d_bar, d_expect, R_index);
1353 }
1354 else
1355 GMT_Report (API, GMT_MSG_INFORMATION, "NNA Found %" PRIu64 " points, D_bar = %g\n", n_points, d_bar);
1356
1357 }
1358 if (GMT_End_IO (API, GMT_OUT, 0) != GMT_NOERROR) { /* Disables further data output */
1359 gmt_M_free (GMT, NN_dist);
1360 gmt_M_free (GMT, NN_info);
1361 Return (API->error);
1362 }
1363 gmt_M_free (GMT, NN_dist);
1364 gmt_M_free (GMT, NN_info);
1365 Return (GMT_NOERROR);
1366 }
1367 if (Ctrl->L.active) { /* Remove tile lines only */
1368 int gap = 0, first, prev_OK;
1369 uint64_t row, seg, tbl;
1370 double dx, dy, DX, DY, dist;
1371
1372 if (gmt_init_distaz (GMT, GMT_MAP_DIST_UNIT, 2, GMT_MAP_DIST) == GMT_NOT_A_VALID_TYPE) /* Default is m using great-circle distances */
1373 Return (GMT_NOT_A_VALID_TYPE);
1374
1375 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 */
1376 Return (API->error);
1377 }
1378 if (GMT_Begin_IO (API, GMT_IS_DATASET, GMT_OUT, GMT_HEADER_ON) != GMT_NOERROR) { /* Enables data output and sets access mode */
1379 Return (API->error);
1380 }
1381 if (GMT_Set_Geometry (API, GMT_OUT, geometry) != GMT_NOERROR) { /* Sets output geometry */
1382 Return (API->error);
1383 }
1384 for (tbl = 0; tbl < D->n_tables; tbl++) {
1385 for (seg = 0; seg < D->table[tbl]->n_segments; seg++) {
1386 S = D->table[tbl]->segment[seg];
1387 if (S->n_rows == 0) continue;
1388 for (row = 1, first = true, prev_OK = false; row < S->n_rows; row++) {
1389 dx = S->data[GMT_X][row] - S->data[GMT_X][row-1];
1390 dy = S->data[GMT_Y][row] - S->data[GMT_Y][row-1];
1391 DX = MIN (fabs (S->data[GMT_X][row] - GMT->common.R.wesn[XLO]), fabs (S->data[GMT_X][row] - GMT->common.R.wesn[XHI]));
1392 DY = MIN (fabs (S->data[GMT_Y][row] - GMT->common.R.wesn[YLO]), fabs (S->data[GMT_Y][row] - GMT->common.R.wesn[YHI]));
1393 gap = false;
1394 if ((fabs (dx) < Ctrl->L.path_noise && DX < Ctrl->L.box_offset) || (fabs (dy) < Ctrl->L.path_noise && DY < Ctrl->L.box_offset)) { /* Going along a tile boundary */
1395 dist = gmt_distance (GMT, S->data[GMT_X][row], S->data[GMT_Y][row], S->data[GMT_X][row-1], S->data[GMT_Y][row-1]);
1396 if (dist > Ctrl->L.s_cutoff) gap = true;
1397 }
1398 if (gap) { /* Distance exceed threshold, start new segment */
1399 first = true;
1400 if (prev_OK) gmtspatial_write_record (GMT, S->data, S->n_columns, row-1);
1401 prev_OK = false;
1402 }
1403 else {
1404 if (first && S->header) {
1405 strncpy (GMT->current.io.segment_header, S->header, GMT_BUFSIZ-1);
1406 GMT_Put_Record (API, GMT_WRITE_SEGMENT_HEADER, NULL);
1407 }
1408 gmtspatial_write_record (GMT, S->data, S->n_columns, row-1);
1409 first = false;
1410 prev_OK = true;
1411 }
1412 }
1413 if (!gap) gmtspatial_write_record (GMT, S->data, S->n_columns, row-1);
1414 }
1415
1416 }
1417 if (GMT_End_IO (API, GMT_OUT, 0) != GMT_NOERROR) { /* Disables further data output */
1418 Return (API->error);
1419 }
1420 if (GMT_Destroy_Data (API, &D) != GMT_NOERROR) {
1421 Return (API->error);
1422 }
1423 Return (GMT_NOERROR);
1424 }
1425
1426 if (Ctrl->Q.active) { /* Calculate centroid and polygon areas or line lengths and place in segment headers */
1427 double out[3];
1428 static char *type[2] = {"length", "area"}, upper[GMT_LEN32] = {"infinity"};
1429 bool new_data = (Ctrl->Q.header || Ctrl->Q.sort || Ctrl->E.active);
1430 uint64_t seg, row_f, row_l, tbl, col, n_seg = 0, n_alloc_seg = 0;
1431 unsigned int handedness = 0;
1432 int qmode, poly = 0, geo;
1433 struct GMT_DATASET *Dout = NULL;
1434 struct GMT_DATASEGMENT *Sout = NULL;
1435 struct GMT_ORDER *Q = NULL;
1436
1437 char line[GMT_LEN128] = {""};
1438
1439 if (gmt_M_is_cartesian (GMT, GMT_IN) && Ctrl->Q.unit && strchr (GMT_LEN_UNITS, Ctrl->Q.unit)) {
1440 gmt_parse_common_options (GMT, "f", 'f', "g"); /* Set -fg if -Q uses unit */
1441 }
1442 geo = gmt_M_is_geographic (GMT, GMT_IN);
1443 if (geo) {
1444 if (gmt_init_distaz (GMT, Ctrl->Q.unit, Ctrl->Q.dmode, GMT_MAP_DIST) == GMT_NOT_A_VALID_TYPE) /* Default is m using great-circle distances */
1445 Return (GMT_NOT_A_VALID_TYPE);
1446 }
1447
1448 if (Ctrl->Q.header) { /* Add line length or polygon area stuff to segment header */
1449 qmode = Ctrl->Q.mode; /* Don't know if line or polygon but passing GMT_IS_POLY would close any open polygon, which we want with +p */
1450 GMT->current.io.multi_segments[GMT_OUT] = true; /* To ensure we can write headers */
1451 }
1452 else {
1453 qmode = GMT_IS_POINT;
1454 if ((error = GMT_Set_Columns (API, GMT_OUT, 3, GMT_COL_FIX_NO_TEXT)) != 0) Return (error);
1455 }
1456 if (new_data) { /* Must create an output dataset */
1457 enum GMT_enum_geometry gmtry;
1458 uint64_t dim[GMT_DIM_SIZE] = {1, 0, 0, 0}; /* One table, no rows yet to avoid allocations */
1459 dim[GMT_COL] = D->n_columns; /* Same number of columns as the input */
1460 switch (Ctrl->Q.mode) { /* Set geometry */
1461 case GMT_IS_LINE: gmtry = GMT_IS_LINE; break;
1462 case GMT_IS_POLY: gmtry = GMT_IS_POLY; break;
1463 default: gmtry = D->geometry; break;
1464 }
1465 if ((Dout = GMT_Create_Data (API, GMT_IS_DATASET, gmtry, 0, dim, NULL, NULL, 0, 0, NULL)) == NULL)
1466 Return (API->error);
1467 if (Ctrl->Q.area && Ctrl->Q.limit[1] < DBL_MAX) sprintf (upper, "%.12g", Ctrl->Q.limit[1]);
1468 if (Ctrl->Q.sort) Q = gmt_M_memory (GMT, NULL, D->n_segments, struct GMT_ORDER);
1469 }
1470 else { /* Just write results as one line per segment to stdout */
1471 if (GMT_Init_IO (API, GMT_IS_DATASET, qmode, GMT_OUT, GMT_ADD_DEFAULT, 0, options) != GMT_NOERROR) { /* Registers default output destination, unless already set */
1472 Return (API->error);
1473 }
1474 if (GMT_Begin_IO (API, GMT_IS_DATASET, GMT_OUT, GMT_HEADER_ON) != GMT_NOERROR) { /* Enables data output and sets access mode */
1475 Return (API->error);
1476 }
1477 if (GMT_Set_Geometry (API, GMT_OUT, GMT_IS_NONE) != GMT_NOERROR) { /* Sets output geometry */
1478 Return (API->error);
1479 }
1480 }
1481 Out.data = out; Out.text = NULL;
1482 for (tbl = 0; tbl < D->n_tables; tbl++) {
1483 for (seg = 0; seg < D->table[tbl]->n_segments; seg++) {
1484 S = D->table[tbl]->segment[seg];
1485 if (S->n_rows == 0) continue;
1486 switch (Ctrl->Q.mode) { /* Set or determine line vs polygon type */
1487 case GMT_IS_LINE: poly = 0; break;
1488 case GMT_IS_POLY: poly = 1; break;
1489 default:
1490 poly = (S->n_rows > 2 && !gmt_polygon_is_open (GMT, S->data[GMT_X], S->data[GMT_Y], S->n_rows)); /* Line or polygon */
1491 break;
1492 }
1493 if (poly) /* Polygon */
1494 handedness = gmtspatial_area_size (GMT, S->data[GMT_X], S->data[GMT_Y], S->n_rows, out, geo);
1495 else /* Line */
1496 gmtspatial_length_size (GMT, S->data[GMT_X], S->data[GMT_Y], S->n_rows, out);
1497 /* Must determine if this segment passes our dimension test */
1498 if (Ctrl->Q.area && (out[GMT_Z] < Ctrl->Q.limit[0] || out[GMT_Z] > Ctrl->Q.limit[1])) {
1499 GMT_Report (API, GMT_MSG_INFORMATION, "Input segment %s %g is outside the chosen range %g to %s\n", type[poly], out[GMT_Z], Ctrl->Q.limit[0], upper);
1500 continue;
1501 }
1502 if (Ctrl->Q.header) { /* Add the information to the segment header */
1503 if (S->header) {
1504 if (poly)
1505 snprintf (line, GMT_LEN128, "%s -Z%.12g (area) %.12g/%.12g (centroid) %s", S->header, out[GMT_Z], out[GMT_X], out[GMT_Y], kind[handedness]);
1506 else
1507 snprintf (line, GMT_LEN128, "%s -Z%.12g (length) %.12g/%.12g (midpoint)", S->header, out[GMT_Z], out[GMT_X], out[GMT_Y]);
1508 }
1509 else {
1510 if (poly)
1511 snprintf (line, GMT_LEN128, "-Z%.12g (area) %.12g/%.12g (centroid) %s", out[GMT_Z], out[GMT_X], out[GMT_Y], kind[handedness]);
1512 else
1513 snprintf (line, GMT_LEN128, "-Z%.12g (length) %.12g/%.12g (midpoint)", out[GMT_Z], out[GMT_X], out[GMT_Y]);
1514 }
1515 }
1516 if (new_data) { /* Must create the output segment and duplicate */
1517 if (n_alloc_seg <= n_seg) {
1518 n_alloc_seg = (n_seg == 0) ? D->n_segments : (n_alloc_seg < 1);
1519 Dout->table[0]->segment = gmt_M_memory (GMT, Dout->table[0]->segment, n_alloc_seg, struct GMT_DATASEGMENT *);
1520 }
1521 smode = (S->text) ? GMT_WITH_STRINGS : GMT_NO_STRINGS;
1522 Sout = GMT_Alloc_Segment (API, smode, S->n_rows, S->n_columns, line, NULL);
1523 for (col = 0; col < S->n_columns; col++) gmt_M_memcpy (Sout->data[col], S->data[col], S->n_rows, double);
1524 Dout->table[0]->segment[n_seg] = Sout;
1525 if (Ctrl->Q.sort) Q[n_seg].value = out[GMT_Z], Q[n_seg].order = n_seg;
1526 n_seg++;
1527 }
1528 if (poly && Ctrl->E.active && handedness != Ctrl->E.mode) { /* Must reverse line */
1529 for (row_f = 0, row_l = Sout->n_rows - 1; row_f < Sout->n_rows/2; row_f++, row_l--) {
1530 for (col = 0; col < Sout->n_columns; col++) gmt_M_double_swap (Sout->data[col][row_f], Sout->data[col][row_l]);
1531 }
1532 handedness = Ctrl->E.mode;
1533 }
1534 if (!new_data) /* Just write centroid and area|length to output */
1535 GMT_Put_Record (API, GMT_WRITE_DATA, &Out);
1536 }
1537 }
1538 if (new_data) { /* Must write out a revised dataset */
1539 Dout->n_segments = Dout->table[0]->n_segments = n_seg;
1540 Dout->table[0]->segment = gmt_M_memory (GMT, Dout->table[0]->segment, n_seg, struct GMT_DATASEGMENT *);
1541 if (Ctrl->Q.sort) { /* Sort on area or length and shuffle the order of segments before output */
1542 struct GMT_DATASEGMENT **tmp = gmt_M_memory (GMT, NULL, n_seg, struct GMT_DATASEGMENT *);
1543 gmt_sort_order (GMT, Q, n_seg, Ctrl->Q.dir);
1544 for (seg = 0; seg < n_seg; seg++) tmp[seg] = Dout->table[0]->segment[Q[seg].order];
1545 gmt_M_memcpy (Dout->table[0]->segment, tmp, n_seg, struct GMT_DATASEGMENT *);
1546 gmt_M_free (GMT, tmp);
1547 gmt_M_free (GMT, Q);
1548 }
1549 if (Dout->n_segments > 1) gmt_set_segmentheader (GMT, GMT_OUT, true); /* Turn on "-mo" */
1550 if (GMT_Write_Data (API, GMT_IS_DATASET, GMT_IS_FILE, GMT_IS_POLY, GMT_WRITE_SET, NULL, Ctrl->Out.file, Dout) != GMT_NOERROR) {
1551 Return (API->error);
1552 }
1553 n_seg = D->n_segments - Dout->n_segments; /* Lost segments */
1554 if (n_seg) GMT_Report (API, GMT_MSG_INFORMATION, "%" PRIu64 " segments were outside and %" PRIu64 " were inside the chosen %s range of %g to %s\n",
1555 n_seg, Dout->n_segments, type[poly], Ctrl->Q.limit[0], upper);
1556 }
1557 if (!new_data && GMT_End_IO (API, GMT_OUT, 0) != GMT_NOERROR) { /* Disables further data output */
1558 Return (API->error);
1559 }
1560 if (GMT_Destroy_Data (API, &D) != GMT_NOERROR) {
1561 Return (API->error);
1562 }
1563 Return (GMT_NOERROR);
1564 }
1565
1566 if (Ctrl->I.active || external) { /* Crossovers between polygons */
1567 bool same_feature;
1568 unsigned int in, wtype, n_columns;
1569 uint64_t tbl1, tbl2, col, nx, row, seg1, seg2;
1570 struct GMT_XSEGMENT *ylist1 = NULL, *ylist2 = NULL;
1571 struct GMT_XOVER XC;
1572 char record[GMT_BUFSIZ] = {""}, fmt[GMT_BUFSIZ] = {""};
1573 struct GMT_DATASET *C = NULL;
1574 struct GMT_DATASEGMENT *S1 = NULL, *S2 = NULL;
1575
1576 if (Ctrl->S.mode == POL_CLIP) { /* Need to set up a separate table with the clip polygon */
1577 if (Ctrl->T.file) {
1578 gmt_disable_bghio_opts (GMT); /* Do not want any -b -g -h -i -o to affect the reading from -C,-F,-L files */
1579 if ((C = GMT_Read_Data (API, GMT_IS_DATASET, GMT_IS_FILE, GMT_IS_POLY, GMT_READ_NORMAL, NULL, Ctrl->T.file, NULL)) == NULL) {
1580 Return (API->error);
1581 }
1582 if (C->n_columns < 2) {
1583 GMT_Report (API, GMT_MSG_ERROR, "Input data have %d column(s) but at least 2 are needed\n", (int)C->n_columns);
1584 Return (GMT_DIM_TOO_SMALL);
1585 }
1586 gmt_reenable_bghio_opts (GMT); /* Recover settings provided by user (if -b -g -h -i were used at all) */
1587 }
1588 else { /* Design a table based on -Rw/e/s/n */
1589 uint64_t dim[GMT_DIM_SIZE] = {1, 1, 5, 2};
1590 if ((C = GMT_Create_Data (API, GMT_IS_DATASET, GMT_IS_POLY, 0, dim, NULL, NULL, 0, 0, NULL)) == NULL) Return (API->error);
1591 S1 = C->table[0]->segment[0];
1592 S1->data[GMT_X][0] = S1->data[GMT_X][3] = S1->data[GMT_X][4] = GMT->common.R.wesn[XLO];
1593 S1->data[GMT_X][1] = S1->data[GMT_X][2] = GMT->common.R.wesn[XHI];
1594 S1->data[GMT_Y][0] = S1->data[GMT_Y][1] = S1->data[GMT_Y][4] = GMT->common.R.wesn[YLO];
1595 S1->data[GMT_Y][2] = S1->data[GMT_Y][3] = GMT->common.R.wesn[YHI];
1596 S1->n_rows = 5;
1597 }
1598 }
1599 else
1600 C = D; /* Compare with itself */
1601
1602 geometry = (Ctrl->S.active) ? GMT_IS_PLP : GMT_IS_NONE;
1603 n_columns = (Ctrl->S.active) ? (unsigned int)C->n_columns : 4;
1604 wtype = (Ctrl->S.active) ? GMT_COL_FIX_NO_TEXT : GMT_COL_FIX;
1605 if ((error = GMT_Set_Columns (API, GMT_OUT, n_columns, wtype)) != GMT_NOERROR) {
1606 Return (error);
1607 }
1608 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 */
1609 Return (API->error);
1610 }
1611 if (GMT_Begin_IO (API, GMT_IS_DATASET, GMT_OUT, GMT_HEADER_ON) != GMT_NOERROR) { /* Enables data output and sets access mode */
1612 Return (API->error);
1613 }
1614 if (GMT_Set_Geometry (API, GMT_OUT, geometry) != GMT_NOERROR) { /* Sets output geometry */
1615 Return (API->error);
1616 }
1617
1618 sprintf (fmt, "%s%s%s%s%s%s%s%s%%s%s%%s\n", GMT->current.setting.format_float_out, GMT->current.setting.io_col_separator, GMT->current.setting.format_float_out, \
1619 GMT->current.setting.io_col_separator, GMT->current.setting.format_float_out, GMT->current.setting.io_col_separator, GMT->current.setting.format_float_out, \
1620 GMT->current.setting.io_col_separator, GMT->current.setting.io_col_separator);
1621 Out.data = out; Out.text = NULL;
1622 for (tbl1 = 0; tbl1 < C->n_tables; tbl1++) {
1623 for (seg1 = 0; seg1 < C->table[tbl1]->n_segments; seg1++) {
1624 S1 = C->table[tbl1]->segment[seg1];
1625 if (S1->n_rows == 0) continue;
1626 gmt_init_track (GMT, S1->data[GMT_Y], S1->n_rows, &ylist1);
1627 for (tbl2 = (Ctrl->S.mode == POL_CLIP) ? 0 : tbl1; tbl2 < D->n_tables; tbl2++) {
1628 for (seg2 = 0; seg2 < D->table[tbl2]->n_segments; seg2++) {
1629 S2 = D->table[tbl2]->segment[seg2];
1630 if (S2->n_rows == 0) continue;
1631 if (Ctrl->S.mode != POL_CLIP) { /* So there is only one dataset being compared with itself */
1632 same_feature = (external == 2 || internal == 2) ? (tbl1 == tbl2) : (tbl1 == tbl2 && seg1 == seg2); /* What constitutes the same feature */
1633 if (!internal && same_feature) continue; /* Do not do internal crossings */
1634 if (!external && !same_feature) continue; /* Do not do external crossings */
1635 }
1636 gmt_init_track (GMT, S2->data[GMT_Y], S2->n_rows, &ylist2);
1637 nx = gmt_crossover (GMT, S1->data[GMT_X], S1->data[GMT_Y], NULL, ylist1, S1->n_rows, S2->data[GMT_X], S2->data[GMT_Y], NULL, ylist2, S2->n_rows, (internal > 0), gmt_M_is_geographic (GMT, GMT_IN), &XC);
1638 if (nx) { /* Polygon pair generated crossings */
1639 uint64_t px;
1640 if (Ctrl->S.active) { /* Do the spatial clip operation */
1641 uint64_t row0;
1642 bool go, first;
1643 double *xx = NULL, *yy = NULL, *kk = NULL;
1644 struct GMTSPATIAL_PAIR *pair = NULL;
1645
1646 pair = gmt_M_memory (GMT, NULL, nx, struct GMTSPATIAL_PAIR);
1647 xx = gmt_M_memory (GMT, NULL, nx, double);
1648 yy = gmt_M_memory (GMT, NULL, nx, double);
1649 kk = gmt_M_memory (GMT, NULL, nx, double);
1650 for (px = 0; px < nx; px++) pair[px].node = XC.xnode[1][px], pair[px].pos = px;
1651 qsort (pair, nx, sizeof (struct GMTSPATIAL_PAIR), gmtspatial_comp_pairs);
1652 for (px = 0; px < nx; px++) {
1653 xx[px] = XC.x[pair[px].pos];
1654 yy[px] = XC.y[pair[px].pos];
1655 kk[px] = XC.xnode[1][pair[px].pos];
1656 }
1657 gmt_M_free (GMT, pair);
1658 in = gmt_non_zero_winding (GMT, S2->data[GMT_X][0], S2->data[GMT_Y][0], S1->data[GMT_X], S1->data[GMT_Y], S1->n_rows);
1659 go = first = true;
1660 row0 = px = 0;
1661 while (go) {
1662 for (row = row0; row < S2->n_rows && row < kk[px]; row++) {
1663 if (!in) continue;
1664 for (col = 0; col < S2->n_columns; col++) out[col] = S2->data[col][row];
1665 if (first && GMT->current.io.multi_segments[GMT_OUT]) { /* Must find unique edges to output only once */
1666 if (S2->header)
1667 strncpy (GMT->current.io.segment_header, S2->header, GMT_BUFSIZ-1);
1668 else
1669 sprintf (GMT->current.io.segment_header, "New segment");
1670 GMT_Put_Record (API, GMT_WRITE_SEGMENT_HEADER, NULL);
1671 first = false;
1672 }
1673 GMT_Put_Record (API, GMT_WRITE_DATA, &Out); /* Write this to output */
1674 }
1675 /* Always output crossover point */
1676 if (first && GMT->current.io.multi_segments[GMT_OUT]) { /* Must find unique edges to output only once */
1677 if (S2->header)
1678 strncpy (GMT->current.io.segment_header, S2->header, GMT_BUFSIZ-1);
1679 else
1680 sprintf (GMT->current.io.segment_header, "New segment");
1681 GMT_Put_Record (API, GMT_WRITE_SEGMENT_HEADER, NULL);
1682 first = false;
1683 }
1684 for (col = 2; col < S2->n_columns; col++) out[col] = 0.0;
1685 out[GMT_X] = xx[px]; out[GMT_Y] = yy[px];
1686 GMT_Put_Record (API, GMT_WRITE_DATA, &Out); /* Write this to output */
1687 px++;
1688 in = !in; /* Go from out to in or vice versa */
1689 if (!in) first = true; /* Since we went outside */
1690 row0 = row;
1691 if (px == nx) {
1692 for (row = row0; row < S2->n_rows; row++) {
1693 if (!in) continue;
1694 for (col = 0; col < S2->n_columns; col++) out[col] = S2->data[col][row];
1695 if (first && GMT->current.io.multi_segments[GMT_OUT]) { /* Must find unique edges to output only once */
1696 if (S2->header)
1697 strncpy (GMT->current.io.segment_header, S2->header, GMT_BUFSIZ-1);
1698 else
1699 sprintf (GMT->current.io.segment_header, "New segment");
1700 GMT_Put_Record (API, GMT_WRITE_SEGMENT_HEADER, NULL);
1701 first = false;
1702 }
1703 GMT_Put_Record (API, GMT_WRITE_DATA, &Out); /* Write this to output */
1704 }
1705 go = false;
1706 }
1707 }
1708 gmt_M_free (GMT, xx);
1709 gmt_M_free (GMT, yy);
1710 gmt_M_free (GMT, kk);
1711 }
1712 else { /* Just report */
1713 struct GMT_DATATABLE_HIDDEN *TH1 = gmt_get_DT_hidden (C->table[tbl1]), *TH2 = gmt_get_DT_hidden (C->table[tbl2]);
1714 char F1[PATH_MAX] = {""}, F2[PATH_MAX] = {""};
1715 if (TH1->file[GMT_IN] == NULL) snprintf (F1, PATH_MAX, "%" PRIu64 "-%" PRIu64, tbl1, seg1); else snprintf (F1, PATH_MAX, "%s-%" PRIu64, TH1->file[GMT_IN], seg1);
1716 if (TH2->file[GMT_IN] == NULL) snprintf (F2, PATH_MAX, "%" PRIu64 "-%" PRIu64, tbl2, seg2); else snprintf (F2, PATH_MAX, "%s-%" PRIu64, TH2->file[GMT_IN], seg2);
1717 sprintf (record, "%s%s%s", F1, GMT->current.setting.io_col_separator, F2);
1718 Out.text = record;
1719 for (px = 0; px < nx; px++) { /* Write these to output */
1720 out[GMT_X] = XC.x[px]; out[GMT_Y] = XC.y[px];
1721 out[2] = XC.xnode[0][px]; out[3] = XC.xnode[1][px];
1722 GMT_Put_Record (API, GMT_WRITE_DATA, &Out);
1723 }
1724 }
1725 gmt_x_free (GMT, &XC);
1726 }
1727 else if (Ctrl->S.mode == POL_CLIP) { /* No crossings; see if it is inside or outside C */
1728 if ((in = gmt_non_zero_winding (GMT, S2->data[GMT_X][0], S2->data[GMT_Y][0], S1->data[GMT_X], S1->data[GMT_Y], S1->n_rows)) != 0) {
1729 /* Inside, copy out the entire polygon */
1730 if (GMT->current.io.multi_segments[GMT_OUT]) { /* Must find unique edges to output only once */
1731 if (S2->header)
1732 strncpy (GMT->current.io.segment_header, S2->header, GMT_BUFSIZ-1);
1733 else
1734 sprintf (GMT->current.io.segment_header, "New segment");
1735 GMT_Put_Record (API, GMT_WRITE_SEGMENT_HEADER, NULL);
1736 }
1737 for (row = 0; row < S2->n_rows; row++) {
1738 for (col = 0; col < S2->n_columns; col++) out[col] = S2->data[col][row];
1739 GMT_Put_Record (API, GMT_WRITE_DATA, &Out); /* Write this to output */
1740 }
1741 }
1742 }
1743 gmt_M_free (GMT, ylist2);
1744 if (Ctrl->S.mode == POL_UNION) {
1745 GMT_Report (API, GMT_MSG_ERROR, "Computing polygon union not implemented yet\n");
1746 }
1747 if (Ctrl->S.mode == POL_INTERSECTION) {
1748 GMT_Report (API, GMT_MSG_ERROR, "Computing polygon intersection not implemented yet\n");
1749 }
1750 }
1751 }
1752 gmt_M_free (GMT, ylist1);
1753 }
1754 }
1755 if (GMT_End_IO (API, GMT_OUT, 0) != GMT_NOERROR) { /* Disables further data output */
1756 Return (API->error);
1757 }
1758 if (GMT_Destroy_Data (API, &D) != GMT_NOERROR) {
1759 Return (API->error);
1760 }
1761 if (Ctrl->S.mode == POL_CLIP) {
1762 if (GMT_Destroy_Data (API, &C) != GMT_NOERROR) {
1763 Return (API->error);
1764 }
1765 }
1766 Return (GMT_NOERROR);
1767 }
1768
1769 if (Ctrl->D.active) { /* Look for duplicates of lines or polygons */
1770 unsigned int n_dup, poly_D, poly_S2;
1771 uint64_t tbl, tbl2, col, seg, seg2;
1772 bool same_feature = false;
1773 char *kind[10] = {"approximate-reversed-superset", "approximate-reversed-subset", "approximate-reversed", "exact-reversed" , "", "exact", "approximate", "approximate-subset", "approximate-superset", "Dateline-split"};
1774 char record[GMT_BUFSIZ] = {""}, format[GMT_BUFSIZ] = {""}, src[GMT_BUFSIZ] = {""}, dup[GMT_BUFSIZ] = {""}, *feature[2] = {"polygon", "line"}, *from = NULL;
1775 char *in = "the same data set", *verdict = "NY~-+/"; /* No, Yes, Approximate, Subsection, Supersection */
1776 struct GMT_DATASET *C = NULL;
1777 struct GMT_DATASEGMENT *S1 = NULL, *S2 = NULL;
1778 struct DUP_INFO **Info = NULL, *I = NULL;
1779
1780 if (Ctrl->D.file) { /* Get trial features via a file */
1781 gmt_disable_bghio_opts (GMT); /* Do not want any -b -g -h -i -o to affect the reading from -D files */
1782 if ((C = GMT_Read_Data (API, GMT_IS_DATASET, GMT_IS_FILE, GMT_IS_LINE|GMT_IS_POLY, GMT_READ_NORMAL, NULL, Ctrl->D.file, NULL)) == NULL) {
1783 Return (API->error);
1784 }
1785 if (C->n_columns < 2) {
1786 GMT_Report (API, GMT_MSG_ERROR, "Input data have %d column(s) but at least 2 are needed\n", (int)C->n_columns);
1787 Return (GMT_DIM_TOO_SMALL);
1788 }
1789 gmt_reenable_bghio_opts (GMT); /* Recover settings provided by user (if -b -g -h -i were used at all) */
1790 from = Ctrl->D.file;
1791 }
1792 else {
1793 C = D; /* Compare with itself */
1794 same_feature = true;
1795 from = in;
1796 smode = (C->table[0]->segment[0]->text) ? GMT_WITH_STRINGS : GMT_NO_STRINGS;
1797 S2 = GMT_Alloc_Segment (GMT->parent, smode, 0, C->n_columns, NULL, NULL);
1798 }
1799 if (GMT_Init_IO (API, GMT_IS_DATASET, C->geometry, GMT_OUT, GMT_ADD_DEFAULT, 0, options) != GMT_NOERROR) {
1800 gmt_free_segment (GMT, &S2);
1801 Return (API->error); /* Registers default output destination, unless already set */
1802 }
1803 if (GMT_Begin_IO (API, GMT_IS_DATASET, GMT_OUT, GMT_HEADER_ON) != GMT_NOERROR) {
1804 gmt_free_segment (GMT, &S2);
1805 Return (API->error); /* Enables data output and sets access mode */
1806 }
1807 if (GMT_Set_Geometry (API, GMT_OUT, C->geometry) != GMT_NOERROR) { /* Sets output geometry */
1808 Return (API->error);
1809 }
1810
1811 Info = gmt_M_memory (GMT, NULL, C->n_tables, struct DUP_INFO *);
1812 for (tbl = 0; tbl < C->n_tables; tbl++) Info[tbl] = gmt_M_memory (GMT, NULL, C->table[tbl]->n_segments, struct DUP_INFO);
1813
1814 if (gmt_init_distaz (GMT, Ctrl->D.unit, Ctrl->D.mode, GMT_MAP_DIST) == GMT_NOT_A_VALID_TYPE)
1815 Return (GMT_NOT_A_VALID_TYPE);
1816
1817 sprintf (format, "%%c : Input %%s %%s is an %%s duplicate of a %%s %%s in %%s, with d = %s c = %%.6g s = %%.4g",
1818 GMT->current.setting.format_float_out);
1819
1820 Out.text = record;
1821 for (tbl = 0; tbl < D->n_tables; tbl++) {
1822 for (seg = 0; seg < D->table[tbl]->n_segments; seg++) {
1823 S1 = D->table[tbl]->segment[seg];
1824 if (S1->n_rows == 0) continue;
1825
1826 GMT_Report (API, GMT_MSG_INFORMATION, "Check if segment %" PRIu64 " from Table %d has duplicates:\n", seg, tbl);
1827 if (same_feature) { /* We must exclude this segment from the comparison otherwise we end up finding itself as a duplicate */
1828 S2->n_rows = S1->n_rows;
1829 for (col = 0; col < S1->n_columns; col++) S2->data[col] = S1->data[col];
1830 S1->n_rows = 0; /* This means it will be skipped by gmtspatial_is_duplicate */
1831 }
1832 else
1833 S2 = S1;
1834 poly_S2 = (S2->n_rows > 2 && gmt_polygon_is_open (GMT, S2->data[GMT_X], S2->data[GMT_Y], S2->n_rows)) ? 1 : 0;
1835 for (tbl2 = 0; tbl2 < C->n_tables; tbl2++) gmt_M_memset (Info[tbl2], C->table[tbl2]->n_segments, struct DUP_INFO);
1836 n_dup = gmtspatial_is_duplicate (GMT, S2, C, &(Ctrl->D.I), Info); /* Returns -3, -2, -1, 0, +1, +2, or +3 */
1837 if (same_feature) {
1838 S1->n_rows = S2->n_rows; /* Reset the count */
1839 if (Ctrl->D.I.table < tbl || (Ctrl->D.I.table == tbl && Ctrl->D.I.segment < seg)) n_dup = 0; /* To avoid reporting the same pair twice */
1840 }
1841 if (n_dup == 0) { /* No duplicate found for this segment */
1842 if (!same_feature) {
1843 (D->n_tables == 1) ? sprintf (src, "[ segment %" PRIu64 " ]", seg) : sprintf (src, "[ table %" PRIu64 " segment %" PRIu64 " ]", tbl, seg);
1844 poly_D = (D->table[tbl]->segment[seg]->n_rows > 2 && gmt_polygon_is_open (GMT, D->table[tbl]->segment[seg]->data[GMT_X], D->table[tbl]->segment[seg]->data[GMT_Y], D->table[tbl]->segment[seg]->n_rows)) ? 1 : 0;
1845 sprintf (record, "N : Input %s %s not present in %s", feature[poly_D], src, from);
1846 GMT_Put_Record (API, GMT_WRITE_DATA, &Out);
1847 }
1848 continue;
1849 }
1850 for (tbl2 = 0; tbl2 < C->n_tables; tbl2++) {
1851 for (seg2 = 0; seg2 < C->table[tbl2]->n_segments; seg2++) {
1852 I = &(Info[tbl2][seg2]);
1853 if (I->mode == 0) continue;
1854 /* Report on all the close/exact matches */
1855 poly_D = (C->table[tbl2]->segment[seg2]->n_rows > 2 && gmt_polygon_is_open (GMT, C->table[tbl2]->segment[seg2]->data[GMT_X],
1856 C->table[tbl2]->segment[seg2]->data[GMT_Y], C->table[tbl2]->segment[seg2]->n_rows)) ? 1 : 0;
1857 (D->n_tables == 1) ? sprintf (src, "[ segment %" PRIu64 " ]", seg) : sprintf (src, "[ table %" PRIu64 " segment %" PRIu64 " ]", tbl, seg);
1858 (C->n_tables == 1) ? sprintf (dup, "[ segment %" PRIu64 " ]", seg2) : sprintf (dup, "[ table %" PRIu64 " segment %" PRIu64 " ]", tbl2, seg2);
1859 if (I->mode == 5)
1860 sprintf (record, "| : Input %s %s was separated at the Dateline from %s %s in %s", feature[poly_D], src, feature[poly_S2], dup, from);
1861 else
1862 sprintf (record, format, verdict[abs(I->mode)], feature[poly_D], src, kind[I->mode+4], feature[poly_S2],
1863 dup, from, I->distance, I->closeness, I->setratio);
1864 GMT_Put_Record (API, GMT_WRITE_DATA, &Out);
1865 }
1866 }
1867 }
1868 }
1869 if (same_feature) {
1870 for (col = 0; col < S2->n_columns; col++) S2->data[col] = NULL; /* Since they were not allocated */
1871 gmt_free_segment (GMT, &S2);
1872 }
1873 for (tbl = 0; tbl < C->n_tables; tbl++) gmt_M_free (GMT, Info[tbl]);
1874 gmt_M_free (GMT, Info);
1875 if (GMT_End_IO (API, GMT_OUT, 0) != GMT_NOERROR) { /* Disables further data output */
1876 Return (API->error);
1877 }
1878 if (GMT_Destroy_Data (API, &D) != GMT_NOERROR) {
1879 Return (API->error);
1880 }
1881 if (Ctrl->D.file && GMT_Destroy_Data (API, &C) != GMT_NOERROR) {
1882 Return (API->error);
1883 }
1884 Return (GMT_NOERROR);
1885 }
1886
1887 if (Ctrl->C.active) { /* Clip polygon to bounding box */
1888 uint64_t np, p, nx, tbl, seg, col;
1889 double *cp[2] = {NULL, NULL};
1890 struct GMT_DATASEGMENT_HIDDEN *SH = NULL;
1891 if (!GMT->common.J.active) { /* -J not specified, set one implicitly */
1892 /* Supply dummy linear proj for Cartesian or geographic data */
1893 if (gmt_M_is_geographic (GMT, GMT_IN))
1894 gmt_parse_common_options (GMT, "J", 'J', "x1d"); /* Fake linear degree projection */
1895 else
1896 gmt_parse_common_options (GMT, "J", 'J', "x1"); /* Fake linear Cartesian projection */
1897 }
1898 if (gmt_M_err_pass (GMT, gmt_proj_setup (GMT, GMT->common.R.wesn), "")) Return (GMT_RUNTIME_ERROR);
1899 for (tbl = 0; tbl < D->n_tables; tbl++) {
1900 for (seg = 0; seg < D->table[tbl]->n_segments; seg++) {
1901 S = D->table[tbl]->segment[seg];
1902 SH = gmt_get_DS_hidden (S);
1903 if ((np = (*GMT->current.map.clip) (GMT, S->data[GMT_X], S->data[GMT_Y], S->n_rows, &cp[GMT_X], &cp[GMT_Y], &nx)) == 0) {
1904 /* Everything is outside, let go */
1905 SH->mode = GMT_WRITE_SKIP;
1906 continue;
1907 }
1908 smode = (S->text) ? GMT_WITH_STRINGS : GMT_NO_STRINGS;
1909 if (np > S->n_rows) S = GMT_Alloc_Segment (GMT->parent, smode, np, S->n_columns, NULL, S);
1910 for (p = 0; p < np; p++) gmt_xy_to_geo (GMT, &cp[GMT_X][p], &cp[GMT_Y][p], cp[GMT_X][p], cp[GMT_Y][p]);
1911 for (col = 0; col < 2; col++) {
1912 gmt_M_memcpy (S->data[col], cp[col], np, double);
1913 gmt_M_free (GMT, cp[col]);
1914 }
1915 S->n_rows = np;
1916 }
1917 }
1918 if (GMT_Write_Data (API, GMT_IS_DATASET, GMT_IS_FILE, GMT_IS_POLY, GMT_WRITE_SET, NULL, Ctrl->Out.file, D) != GMT_NOERROR) {
1919 Return (API->error);
1920 }
1921 }
1922
1923 if (Ctrl->N.active) { /* Report the polygons that contain the given features */
1924 uint64_t tbl, row, first, last, n, p, np, seg, seg2, n_inside;
1925 unsigned int *count = NULL, nmode;
1926 int ID = -1;
1927 char seg_label[GMT_LEN64] = {""}, record[GMT_BUFSIZ] = {""}, *kind[2] = {"Middle point", "All points"};
1928 struct GMT_DATASET *C = NULL;
1929 struct GMT_DATATABLE *T = NULL;
1930 struct GMT_DATASEGMENT *S = NULL, *S2 = NULL;
1931 struct GMT_DATASEGMENT_HIDDEN *SH = NULL;
1932
1933 gmt_disable_bghio_opts (GMT); /* Do not want any -b -g -h -i -o to affect the reading from -CN files */
1934 if ((C = GMT_Read_Data (API, GMT_IS_DATASET, GMT_IS_FILE, GMT_IS_POLY, GMT_READ_NORMAL, NULL, Ctrl->N.file, NULL)) == NULL) {
1935 Return (API->error);
1936 }
1937 if (C->n_columns < 2) {
1938 GMT_Report (API, GMT_MSG_ERROR, "Input data have %d column(s) but at least 2 are needed\n", (int)C->n_columns);
1939 Return (GMT_DIM_TOO_SMALL);
1940 }
1941 gmt_reenable_bghio_opts (GMT); /* Recover settings provided by user (if -b -g -h -i were used at all) */
1942 nmode = (Ctrl->N.mode == 1) ? GMT_IS_NONE : GMT_IS_LINE;
1943 if (GMT_Init_IO (API, GMT_IS_DATASET, nmode, GMT_OUT, GMT_ADD_DEFAULT, 0, options) != GMT_NOERROR) { /* Registers default output destination, unless already set */
1944 Return (API->error);
1945 }
1946 if (GMT_Begin_IO (API, GMT_IS_DATASET, GMT_OUT, GMT_HEADER_ON) != GMT_NOERROR) { /* Enables data output and sets access mode */
1947 Return (API->error);
1948 }
1949 if (GMT_Set_Geometry (API, GMT_OUT, nmode) != GMT_NOERROR) { /* Sets output geometry */
1950 Return (API->error);
1951 }
1952
1953 if (Ctrl->N.mode == 2) gmt_adjust_dataset (GMT, D, D->n_columns + 1); /* Add one more output column */
1954
1955 T = C->table[0]; /* Only one input file so only one table */
1956 count = gmt_M_memory (GMT, NULL, D->n_segments, unsigned int);
1957 Out.text = record;
1958 for (seg2 = 0; seg2 < T->n_segments; seg2++) { /* For all polygons */
1959 S2 = T->segment[seg2];
1960 SH = gmt_get_DS_hidden (S2);
1961 if (gmt_polygon_is_hole (GMT, S2)) continue; /* Holes are handled in gmt_inonout */
1962 GMT_Report (API, GMT_MSG_INFORMATION, "Look for points/features inside polygon segment %" PRIu64 " :\n", seg2);
1963 if (Ctrl->N.ID == 0) { /* Look for polygon IDs in the data headers */
1964 if (SH->ogr) /* OGR data */
1965 ID = irint (gmt_get_aspatial_value (GMT, GMT_IS_Z, S2));
1966 else if (gmt_parse_segment_item (GMT, S2->header, "-Z", seg_label)) /* Look for segment header ID */
1967 ID = atoi (seg_label);
1968 else if (gmt_parse_segment_item (GMT, S2->header, "-L", seg_label)) /* Look for segment header ID */
1969 ID = atoi (seg_label);
1970 else
1971 GMT_Report (API, GMT_MSG_ERROR, "No polygon ID found; ID set to NaN\n");
1972 }
1973 else /* Increment running polygon ID */
1974 ID++;
1975
1976 for (tbl = p = 0; tbl < D->n_tables; tbl++) {
1977 for (seg = 0; seg < D->table[tbl]->n_segments; seg++, p++) {
1978 S = D->table[tbl]->segment[seg];
1979 if (S->n_rows == 0) continue;
1980 if (Ctrl->N.all) { first = 0; last = S->n_rows - 1; np = S->n_rows; } else { first = last = S->n_rows / 2; np = 1; }
1981 for (row = first, n = 0; row <= last; row++) { /* Check one or all points if they are inside */
1982 n += (gmt_inonout (GMT, S->data[GMT_X][row], S->data[GMT_Y][row], S2) == GMT_INSIDE);
1983 }
1984 if (n < np) continue; /* Not inside this polygon */
1985 if (count[p]) {
1986 GMT_Report (API, GMT_MSG_ERROR, "Segment %" PRIu64 "-%" PRIu64 " already inside another polygon; skipped\n", tbl, seg);
1987 continue;
1988 }
1989 count[p]++;
1990 /* Here we are inside */
1991 if (Ctrl->N.mode == 1) { /* Just report on which polygon contains each feature */
1992 sprintf (record, "%s from table %" PRIu64 " segment %" PRIu64 " is inside polygon # %d", kind[Ctrl->N.all], tbl, seg, ID);
1993 GMT_Put_Record (API, GMT_WRITE_DATA, &Out);
1994 }
1995 else if (Ctrl->N.mode == 2) { /* Add ID as last data column */
1996 for (row = 0, n = S->n_columns-1; row < S->n_rows; row++) S->data[n][row] = (double)ID;
1997 GMT_Report (API, GMT_MSG_INFORMATION, "%s from table %" PRIu64 " segment %" PRIu64 " is inside polygon # %d\n", kind[Ctrl->N.all], tbl, seg, ID);
1998 }
1999 else { /* Add ID via the segment header -Z */
2000 if (gmt_parse_segment_item (GMT, S->header, "-Z", NULL))
2001 GMT_Report (API, GMT_MSG_ERROR, "Segment header %d-%" PRIu64 " already has a -Z flag, skipped\n", tbl, seg);
2002 else { /* Add -Z<ID< to the segment header */
2003 char buffer[GMT_BUFSIZ] = {""}, txt[GMT_LEN64] = {""};
2004 buffer[0] = txt[0] = 0;
2005 if (S->header) { strncpy (buffer, S->header, GMT_BUFSIZ-1); gmt_M_str_free (S->header); }
2006 sprintf (txt, " -Z%d", ID);
2007 strcat (buffer, txt);
2008 S->header = strdup (buffer);
2009 GMT_Report (API, GMT_MSG_INFORMATION, "%s from table %" PRIu64 " segment %" PRIu64 " is inside polygon # %d\n", kind[Ctrl->N.all], tbl, seg, ID);
2010 }
2011 }
2012 }
2013 }
2014 }
2015 for (p = n_inside = 0; p < D->n_segments; p++) if (count[p]) n_inside++;
2016 if (Ctrl->N.mode != 1) { /* Write out results */
2017 if (GMT_Write_Data (API, GMT_IS_DATASET, GMT_IS_FILE, GMT_IS_POLY, GMT_WRITE_SET, NULL, Ctrl->Out.file, D) != GMT_NOERROR) {
2018 Return (API->error);
2019 }
2020 }
2021 if (GMT_End_IO (API, GMT_OUT, 0) != GMT_NOERROR) { /* Disables further data output */
2022 Return (API->error);
2023 }
2024 GMT_Report (API, GMT_MSG_INFORMATION, "%" PRIu64 " segments found to be inside polygons, %" PRIu64 " were outside and skipped\n", n_inside, D->n_segments - n_inside);
2025 if (GMT_Destroy_Data (API, &D) != GMT_NOERROR) {
2026 Return (API->error);
2027 }
2028 if (GMT_Destroy_Data (API, &C) != GMT_NOERROR) {
2029 Return (API->error);
2030 }
2031 gmt_M_free (GMT, count);
2032 Return (GMT_NOERROR);
2033 }
2034 if (Ctrl->S.active && Ctrl->S.mode == POL_SPLIT) { /* Split polygons at dateline */
2035 bool crossing;
2036 uint64_t n_split = 0, tbl, seg_out, seg, n_segs, kseg, n_split_tot = 0;
2037 uint64_t dim[GMT_DIM_SIZE] = {0, 0, 0, 0};
2038 struct GMT_DATASET *Dout = NULL;
2039 struct GMT_DATATABLE *T = NULL;
2040 struct GMT_DATASEGMENT **L = NULL;
2041
2042 dim[GMT_TBL] = D->n_tables; dim[GMT_COL] = D->n_columns; /* Same number of tables and columns as the input */
2043 if ((Dout = GMT_Create_Data (API, GMT_IS_DATASET, GMT_IS_POLY, 0, dim, NULL, NULL, 0, 0, NULL)) == NULL) Return (API->error);
2044 /* Dout has no allocated segments yet */
2045 Dout->n_segments = 0;
2046 for (tbl = 0; tbl < D->n_tables; tbl++) {
2047 T = Dout->table[tbl];
2048 n_segs = D->table[tbl]->n_segments;
2049 Dout->table[tbl]->n_segments = 0;
2050 T->segment = gmt_M_memory (GMT, NULL, n_segs, struct GMT_DATASEGMENT *); /* Need at least this many segments */
2051 for (seg = seg_out = 0; seg < D->table[tbl]->n_segments; seg++) {
2052 S = D->table[tbl]->segment[seg]; /* Current input segment */
2053 if (S->n_rows == 0) continue; /* Just skip empty segments */
2054 crossing = gmt_crossing_dateline (GMT, S);
2055 n_split = (crossing) ? gmt_split_poly_at_dateline (GMT, S, &L) : 1;
2056 Dout->table[tbl]->n_segments += n_split;
2057 if (Dout->table[tbl]->n_segments > n_segs) { /* Must allocate more segment space */
2058 uint64_t old_n_segs = n_segs;
2059 n_segs = Dout->table[tbl]->n_segments;
2060 T->segment = gmt_M_memory (GMT, T->segment, n_segs, struct GMT_DATASEGMENT *); /* Allow more space for new segments */
2061 gmt_M_memset (&(T->segment[old_n_segs]), n_segs - old_n_segs, struct GMT_DATASEGMENT *); /* Set to NULL */
2062 }
2063 /* Here there are space for all segments if new ones are added via splitting */
2064 if (crossing) {
2065 n_split_tot++;
2066 for (kseg = 0; kseg < n_split; kseg++) {
2067 T->segment[seg_out++] = L[kseg]; /* Add the remaining segments to the end */
2068 }
2069 gmt_M_free (GMT, L);
2070 }
2071 else { /* Just duplicate */
2072 T->segment[seg_out++] = gmt_duplicate_segment (GMT, S);
2073 }
2074 }
2075 Dout->table[tbl]->n_segments = seg_out;
2076 Dout->n_segments += seg_out;
2077 }
2078 if (GMT_Write_Data (API, GMT_IS_DATASET, GMT_IS_FILE, GMT_IS_POLY, GMT_WRITE_SET, NULL, Ctrl->Out.file, Dout) != GMT_NOERROR) {
2079 Return (API->error);
2080 }
2081 GMT_Report (API, GMT_MSG_INFORMATION, "%" PRIu64 " segments split across the Dateline\n", n_split_tot);
2082 if (GMT_Destroy_Data (API, &Dout) != GMT_NOERROR) {
2083 Return (API->error);
2084 }
2085 }
2086 if (Ctrl->S.active && Ctrl->S.mode == POL_HOLE) { /* Flag polygons that are holes of others */
2087 uint64_t n_holes = 0, tbl1, seg1, tbl2, seg2, seg_out, k1, k2;
2088 uint64_t dim[GMT_DIM_SIZE] = {1, 0, 0, 0}; /* Only one output table */
2089 unsigned int side, *inside = NULL, *kase = NULL;
2090 int P_handedness, H_handedness;
2091 bool geo;
2092 double out[3];
2093 struct GMT_DATASET *Dout = NULL;
2094 struct GMT_DATATABLE *T1 = NULL, *T2 = NULL;
2095 struct GMT_DATASEGMENT *S1 = NULL, *S2 = NULL;
2096 struct GMT_DATASEGMENT_HIDDEN *SH = NULL;
2097 struct GMT_TBLSEG *K = NULL;
2098
2099 inside = gmt_M_memory (GMT, NULL, D->n_segments, unsigned int);
2100 kase = gmt_M_memory (GMT, NULL, D->n_segments, unsigned int);
2101 K = gmt_M_memory (GMT, NULL, D->n_segments, struct GMT_TBLSEG);
2102 geo = gmt_M_is_geographic (GMT, GMT_IN);
2103
2104 for (tbl1 = k1 = 0; tbl1 < D->n_tables; tbl1++) {
2105 T1 = D->table[tbl1];
2106 for (seg1 = 0; seg1 < T1->n_segments; seg1++, k1++) {
2107 K[k1].tbl = tbl1; K[k1].seg = seg1; /* Fill out the K lookup array */
2108 S1 = D->table[tbl1]->segment[seg1]; /* Current input segment */
2109 if (S1->n_rows == 0) continue; /* Just skip empty segments */
2110 if (S1->header && !strcmp (S1->header, "-Ph")) continue; /* Marked as a hole already */
2111 for (tbl2 = k2 = 0; tbl2 < D->n_tables; tbl2++) {
2112 T2 = D->table[tbl2];
2113 for (seg2 = 0; seg2 < T2->n_segments; seg2++, k2++) {
2114 if (tbl2 < tbl1) continue; /* Avoid duplication */
2115 if (tbl2 == tbl1 && seg2 <= seg1) continue; /* Avoid duplication */
2116 S2 = D->table[tbl2]->segment[seg2]; /* Current input segment */
2117 if (S2->n_rows == 0) continue; /* Just skip empty segments */
2118 if (S2->header && !strcmp (S2->header, "-Ph")) continue; /* Marked as a hole already */
2119 /* Here we determine if S1 is inside S2 or vice versa */
2120 side = gmt_inonout (GMT, S1->data[GMT_X][0], S1->data[GMT_Y][0], S2); /* Is S1 inside S2? */
2121 if (side == GMT_ONEDGE) {
2122 GMT_Report (API, GMT_MSG_ERROR, "Polygon A (tbl=%" PRIu64 ", seg=%" PRIu64 ") is tangent to Polygon B (tbl=%" PRIu64 ", seg=%" PRIu64 ") at (%g/%g). Skipping\n", tbl1, seg1, tbl2, seg2, S1->data[GMT_X][0], S1->data[GMT_Y][0]);
2123 continue;
2124 }
2125 else if (side == GMT_INSIDE) { /* S1 is inside S2 */
2126 n_holes++;
2127 kase[k1]++;
2128 inside[k1] = (unsigned int)k2 + 1; /* So 0 means not inside anything (a perimeter) */
2129 }
2130 side = gmt_inonout (GMT, S2->data[GMT_X][0], S2->data[GMT_Y][0], S1); /* Is S2 inside S1? */
2131 if (side == GMT_ONEDGE) {
2132 GMT_Report (API, GMT_MSG_ERROR, "Polygon B (tbl=%" PRIu64 ", seg=%" PRIu64 ") is tangent to Polygon A (tbl=%" PRIu64 ", seg=%" PRIu64 ") at (%g/%g). Skipping\n", tbl2, seg2, tbl1, seg1, S2->data[GMT_X][0], S2->data[GMT_Y][0]);
2133 continue;
2134 }
2135 else if (side == GMT_INSIDE) { /* S2 is inside S1 */
2136 n_holes++;
2137 kase[k2]++;
2138 inside[k2] = (unsigned int)k1 + 1; /* So 0 means not inside anything (a perimeter) */
2139 }
2140 }
2141 }
2142 }
2143 }
2144 for (k1 = 0; k1 < D->n_segments; k1++) { /* Make sure no polygon is inside more than one other polygon */
2145 if (kase[k1] > 1) {
2146 GMT_Report (API, GMT_MSG_ERROR, "Polygon # %d is inside more than one polygon or is both inside and contains other polygons\n", k1);
2147 gmt_M_free (GMT, kase);
2148 Return (API->error);
2149 }
2150 }
2151
2152 /* Create an output dataset with unallocated segments since rows = 0 */
2153
2154 dim[GMT_COL] = D->n_columns;
2155 if ((Dout = GMT_Create_Data (API, GMT_IS_DATASET, GMT_IS_POLY, 0, dim, NULL, NULL, 0, 0, NULL)) == NULL) Return (API->error);
2156 T1 = Dout->table[0]; /* Only one table used for output */
2157 T1->segment = gmt_M_memory (GMT, NULL, D->n_segments, struct GMT_DATASEGMENT *); /* Need this many segments */
2158 Dout->n_segments = T1->n_segments = D->n_segments;
2159 /* Shuffle segments so perimeters (inside[] == 0) are listed before their holes (inside[] > 0) */
2160 for (k1 = seg_out = 0; k1 < D->n_segments; k1++) { /* Loop over all polygons */
2161 if (inside[k1] == 0) { /* Perimeter polygon */
2162 tbl1 = K[k1].tbl; seg1 = K[k1].seg; /* Get the (tbl,seg) indices for the perimeter */
2163 S1 = D->table[tbl1]->segment[seg1]; /* Current input segment */
2164 /* Duplicate this polygon as next output polygon perimeter */
2165 T1->segment[seg_out++] = gmt_duplicate_segment (GMT, S1);
2166 /* Get perimeter handedness */
2167 P_handedness = gmtspatial_area_size (GMT, S1->data[GMT_X], S1->data[GMT_Y], S1->n_rows, out, geo);
2168 for (k2 = 0; k2 < D->n_segments; k2++) { /* Loop over all polygons */
2169 if (k2 == k1 || inside[k2] != (k1+1)) continue; /* Not a hole inside this perimeter */
2170 tbl2 = K[k2].tbl; seg2 = K[k2].seg; /* Get the (tbl,seg) indices for the hole */
2171 /* Duplicate this polygon as next output polygon hole */
2172 T1->segment[seg_out++] = S2 = gmt_duplicate_segment (GMT, D->table[tbl2]->segment[seg2]);
2173 /* Get hole handedness */
2174 H_handedness = gmtspatial_area_size (GMT, S2->data[GMT_X], S2->data[GMT_Y], S2->n_rows, out, geo);
2175 /* If same handedness then reverse order of polygon */
2176 if (H_handedness == P_handedness) {
2177 uint64_t row_f, row_l, col;
2178 for (row_f = 0, row_l = S2->n_rows - 1; row_f < S2->n_rows/2; row_f++, row_l--) {
2179 for (col = 0; col < S2->n_columns; col++) gmt_M_double_swap (S2->data[col][row_f], S2->data[col][row_l]);
2180 }
2181 }
2182 if (S2->header) { /* Must append -Ph to existing header - need to allocate more space */
2183 S2->header = realloc (S2->header, strlen (S2->header) + 5U);
2184 strncat (S2->header, " -Ph", 4U);
2185 }
2186 else
2187 S2->header = strdup ("-Ph");
2188 SH = gmt_get_DS_hidden (S2);
2189 SH->pol_mode = GMT_IS_HOLE;
2190 }
2191 }
2192 }
2193
2194 gmt_M_free (GMT, kase);
2195 gmt_M_free (GMT, inside);
2196 gmt_M_free (GMT, K);
2197 if (GMT_Write_Data (API, GMT_IS_DATASET, GMT_IS_FILE, GMT_IS_POLY, GMT_WRITE_SET, NULL, Ctrl->Out.file, Dout) != GMT_NOERROR) {
2198 Return (API->error);
2199 }
2200 if (GMT_Destroy_Data (API, &Dout) != GMT_NOERROR) {
2201 Return (API->error);
2202 }
2203 GMT_Report (API, GMT_MSG_INFORMATION, "%" PRIu64 " segments were holes in other polygons\n", n_holes);
2204 }
2205 #ifdef HAVE_GEOS
2206 if (Ctrl->S.active && Ctrl->S.mode == POL_BUFFER) { /* Compute buffer polygon */
2207 error = geos_methods(GMT, D, Ctrl->Out.file, Ctrl->S.width, "buffer");
2208 finishGEOS();
2209 if (error)
2210 Return (error);
2211 }
2212 if (Ctrl->S.active && Ctrl->S.mode == POL_CENTROID) { /* Compute centroid of polygons */
2213 error = geos_methods(GMT, D, Ctrl->Out.file, 0, "centroid");
2214 finishGEOS();
2215 if (error)
2216 Return (error);
2217 }
2218 #endif
2219
2220 if (Ctrl->F.active) { /* We read as polygons to force closure, now write out revised data */
2221 if (GMT_Write_Data (API, GMT_IS_DATASET, GMT_IS_FILE, GMT_IS_POLY, GMT_WRITE_SET, NULL, Ctrl->Out.file, D) != GMT_NOERROR) {
2222 Return (API->error);
2223 }
2224 Return (GMT_NOERROR);
2225 }
2226
2227 Return (GMT_NOERROR);
2228 }
2229
2230 #ifdef HAVE_GEOS
2231
geos_methods(struct GMT_CTRL * GMT,struct GMT_DATASET * D,char * fname,double buf_dist,char * method)2232 int geos_methods(struct GMT_CTRL *GMT, struct GMT_DATASET *D, char *fname, double buf_dist, char *method) {
2233 uint64_t dim[4] = {0,0,0,0};
2234 struct GMT_DATASET *Dout = NULL;
2235
2236 if (!strcmp(method, "buffer") && !strcmp(method, "centroid")) {
2237 GMT_Report (GMT->parent, GMT_MSG_NORMAL, "Unimplemented method -> %s\n", method);
2238 return -1;
2239 }
2240
2241 dim[GMT_TBL] = D->n_tables;
2242 dim[GMT_COL] = (D->n_columns == 2) ? 2 : 3;
2243 if ((Dout = GMT_Create_Data (GMT->parent, GMT_IS_DATASET, GMT_IS_PLP, 0, dim, NULL, NULL, 0, 0, NULL)) == NULL) {
2244 GMT_Report (GMT->parent, GMT_MSG_NORMAL, "Failed to create output dataset.\n");
2245 return GMT->parent->error;
2246 }
2247 Dout->n_segments = D->n_segments;
2248
2249 if (!strcmp(method, "centroid"))
2250 geos_method_polygon(GMT, D, Dout, "");
2251 else if (!strcmp(method, "buffer"))
2252 geos_method_linestring(GMT, D, Dout, buf_dist, "");
2253
2254 if (GMT_Write_Data (GMT->parent, GMT_IS_DATASET, GMT_IS_FILE, GMT_IS_PLP, GMT_WRITE_SET, NULL, fname, Dout) != GMT_NOERROR) {
2255 GMT_Report (GMT->parent, GMT_MSG_NORMAL, "Failed to write output dataset.\n");
2256 return (GMT->parent->error);
2257 }
2258 if (GMT_Destroy_Data (GMT->parent, &Dout) != GMT_NOERROR) {
2259 GMT_Report (GMT->parent, GMT_MSG_NORMAL, "Failed to destroy dataset container.\n");
2260 return (GMT->parent->error);
2261 }
2262
2263 return (GMT_NOERROR);
2264 }
2265
geos_method_polygon(struct GMT_CTRL * GMT,struct GMT_DATASET * Din,struct GMT_DATASET * Dout,char * method)2266 int geos_method_polygon(struct GMT_CTRL *GMT, struct GMT_DATASET *Din, struct GMT_DATASET *Dout, char *method) {
2267 /* This function calls GEOS functions that operate on POLYGON geometries */
2268 unsigned int nt, ns, nr, i, n_pts, n_col;
2269 bool is3D;
2270 GEOSCoordSequence *seq_in = NULL;
2271 const GEOSCoordSequence *seq_out = NULL;
2272 GEOSGeometry *geom = NULL, *geom_out = NULL, *shell = NULL;
2273 GEOSContextHandle_t handle = NULL;
2274
2275 is3D = (Din->n_columns >= 3);
2276 n_col = (Din->n_columns == 2) ? 2 : 3;
2277 handle = initGEOS_r(NULL, NULL);
2278
2279 for (nt = 0; nt < Din->n_tables; nt++) {
2280 Dout->table[nt]->segment = gmt_M_memory (GMT, NULL, 1, struct GMT_DATASEGMENT *);
2281 Dout->table[nt]->n_segments = 1;
2282
2283 Dout->table[nt]->segment[0] = GMT_Alloc_Segment (GMT->parent, GMT_NO_STRINGS, Din->table[nt]->n_segments, n_col, NULL, NULL);
2284 Dout->table[nt]->segment[0]->n_rows = Din->table[nt]->n_segments;
2285 Dout->table[nt]->n_records += Din->table[nt]->n_segments;
2286
2287 for (ns = 0; ns < Din->table[nt]->n_segments; ns++) {
2288 seq_in = GEOSCoordSeq_create_r(handle, (unsigned int)Din->table[nt]->segment[ns]->n_rows, n_col);
2289 if (!seq_in) {
2290 GMT_Report (GMT->parent, GMT_MSG_NORMAL, "Failed to create input GEOS sequence for table %d, segment %d\n", nt, ns);
2291 continue;
2292 }
2293 for (nr = 0; nr < Din->table[nt]->segment[ns]->n_rows; nr++) {
2294 GEOSCoordSeq_setX_r(handle, seq_in, nr, Din->table[nt]->segment[ns]->data[0][nr]);
2295 GEOSCoordSeq_setY_r(handle, seq_in, nr, Din->table[nt]->segment[ns]->data[1][nr]);
2296 if (is3D)
2297 GEOSCoordSeq_setY_r(handle, seq_in, nr, Din->table[nt]->segment[ns]->data[2][nr]);
2298 }
2299
2300 shell = GEOSGeom_createLinearRing_r(handle, seq_in);
2301 geom = GEOSGeom_createPolygon_r(handle, shell, NULL, 0);
2302 geom_out = GEOSGetCentroid_r(handle, geom);
2303
2304 if (!geom_out) {
2305 GMT_Report (GMT->parent, GMT_MSG_NORMAL, "Failed to compute %s for table %d, segment %d\n", method, nt, ns);
2306 continue;
2307 }
2308
2309 if ((n_pts = (unsigned int)GEOSGetNumCoordinates_r(handle, geom_out)) == 0) {
2310 GMT_Report (GMT->parent, GMT_MSG_NORMAL, "GEOS geometry is empty for table %d, segment %d\n", nt, ns);
2311 continue;
2312 }
2313
2314 if ((seq_out = GEOSGeom_getCoordSeq_r(handle, geom_out)) == NULL) {
2315 GMT_Report (GMT->parent, GMT_MSG_NORMAL, "Failed to create output GEOS sequence for table %d, segment %d.\n", nt, ns);
2316 continue;
2317 }
2318
2319 for (i = 0; i < n_pts; i++) {
2320 GEOSCoordSeq_getX_r(handle, seq_out, i, &Dout->table[nt]->segment[0]->data[0][ns]);
2321 GEOSCoordSeq_getY_r(handle, seq_out, i, &Dout->table[nt]->segment[0]->data[1][ns]);
2322 if (is3D)
2323 GEOSCoordSeq_getY_r(handle, seq_out, i, &Dout->table[nt]->segment[0]->data[2][ns]);
2324 }
2325
2326 GEOSGeom_destroy_r(handle, geom);
2327 GEOSGeom_destroy_r(handle, geom_out);
2328 }
2329 Dout->n_records += Dout->table[nt]->n_records;
2330 }
2331 return 0;
2332 }
2333
geos_method_linestring(struct GMT_CTRL * GMT,struct GMT_DATASET * Din,struct GMT_DATASET * Dout,double buf_dist,char * method)2334 int geos_method_linestring(struct GMT_CTRL *GMT, struct GMT_DATASET *Din, struct GMT_DATASET *Dout, double buf_dist, char *method) {
2335 /* This function calls GEOS functions that operate on LINESTRING geometries */
2336 unsigned int nt, ns, nr, i, n_pts, n_col;
2337 bool is3D;
2338 GEOSCoordSequence *seq_in = NULL;
2339 const GEOSCoordSequence *seq_out = NULL;
2340 GEOSGeometry *geom = NULL, *geom_out = NULL;
2341 GEOSContextHandle_t handle = NULL;
2342
2343 is3D = (Din->n_columns >= 3);
2344 n_col = (Din->n_columns == 2) ? 2 : 3;
2345 handle = initGEOS_r(NULL, NULL);
2346
2347 for (nt = 0; nt < Din->n_tables; nt++) {
2348 Dout->table[nt]->segment = gmt_M_memory (GMT, NULL, Din->table[nt]->n_segments, struct GMT_DATASEGMENT *);
2349 Dout->table[nt]->n_segments = Din->table[nt]->n_segments;
2350
2351 for (ns = 0; ns < Din->table[nt]->n_segments; ns++) {
2352 seq_in = GEOSCoordSeq_create_r(handle, (unsigned int)Din->table[nt]->segment[ns]->n_rows, n_col);
2353 if (!seq_in) {
2354 GMT_Report (GMT->parent, GMT_MSG_NORMAL, "Failed to create input GEOS sequence for table %d, segment %d\n", nt, ns);
2355 continue;
2356 }
2357 for (nr = 0; nr < Din->table[nt]->segment[ns]->n_rows; nr++) {
2358 GEOSCoordSeq_setX_r(handle, seq_in, nr, Din->table[nt]->segment[ns]->data[0][nr]);
2359 GEOSCoordSeq_setY_r(handle, seq_in, nr, Din->table[nt]->segment[ns]->data[1][nr]);
2360 if (is3D)
2361 GEOSCoordSeq_setY_r(handle, seq_in, nr, Din->table[nt]->segment[ns]->data[2][nr]);
2362 }
2363 geom = GEOSGeom_createLineString_r(handle, seq_in);
2364
2365 geom_out = GEOSBuffer_r(handle, geom, buf_dist, 30);
2366
2367 if (!geom_out) {
2368 GMT_Report (GMT->parent, GMT_MSG_NORMAL, "Failed to compute %s for table %d, segment %d\n", method, nt, ns);
2369 continue;
2370 }
2371
2372 if ((n_pts = (unsigned int)GEOSGetNumCoordinates_r(handle, geom_out)) == 0) {
2373 GMT_Report (GMT->parent, GMT_MSG_NORMAL, "GEOS geometry is empty for table %d, segment %d\n", nt, ns);
2374 continue;
2375 }
2376
2377 seq_out = GEOSGeom_getCoordSeq_r(handle, GEOSGetExteriorRing_r(handle, geom_out));
2378
2379 if (!seq_out) {
2380 GMT_Report (GMT->parent, GMT_MSG_NORMAL, "Failed to create output GEOS sequence for table %d, segment %d.\n", nt, ns);
2381 continue;
2382 }
2383
2384 Dout->table[nt]->segment[ns] = GMT_Alloc_Segment (GMT->parent, GMT_NO_STRINGS, (uint64_t)n_pts, n_col, NULL, NULL);
2385 Dout->table[nt]->segment[ns]->n_rows = (uint64_t)n_pts;
2386 Dout->table[nt]->n_records += (uint64_t)n_pts;
2387
2388 for (i = 0; i < n_pts; i++) {
2389 GEOSCoordSeq_getX_r(handle, seq_out, i, &Dout->table[nt]->segment[ns]->data[0][i]);
2390 GEOSCoordSeq_getY_r(handle, seq_out, i, &Dout->table[nt]->segment[ns]->data[1][i]);
2391 if (is3D)
2392 GEOSCoordSeq_getY_r(handle, seq_out, i, &Dout->table[nt]->segment[ns]->data[2][i]);
2393 }
2394
2395 GEOSGeom_destroy_r(handle, geom);
2396 GEOSGeom_destroy_r(handle, geom_out);
2397 }
2398 Dout->n_records += Dout->table[nt]->n_records;
2399 }
2400 return 0;
2401 }
2402
2403 #endif
2404