-
Notifications
You must be signed in to change notification settings - Fork 375
/
Copy pathgmtspatial.c
2603 lines (2423 loc) · 122 KB
/
gmtspatial.c
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
/*--------------------------------------------------------------------
*
* Copyright (c) 1991-2025 by the GMT Team (https://www.generic-mapping-tools.org/team.html)
* See LICENSE.TXT file for copying and redistribution conditions.
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU Lesser General Public License as published by
* the Free Software Foundation; version 3 or any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU Lesser General Public License for more details.
*
* Contact info: www.generic-mapping-tools.org
*--------------------------------------------------------------------*/
/*
* gmtspatial performs miscellaneous geospatial operations on polygons, such
* as truncating them against a clipping polygon, calculate areas, find
* crossings with other polygons, etc.
*
* Author: Paul Wessel
* Date: 10-Jun-2009
* Version: 6 API
*
* Note on KEYS: DD(=f mean -D takes an optional input Dataset as argument via the +f modifier.
* ND(= means -N takes a input Dataset as argument which may be followed by optional modifiers.
*/
#include "gmt_dev.h"
#include "longopt/gmtspatial_inc.h"
#define THIS_MODULE_CLASSIC_NAME "gmtspatial"
#define THIS_MODULE_MODERN_NAME "gmtspatial"
#define THIS_MODULE_LIB "core"
#define THIS_MODULE_PURPOSE "Geospatial operations on points, lines and polygons"
#define THIS_MODULE_KEYS "<D{,DD(=f,ND(=,TD(,>D}"
#define THIS_MODULE_NEEDS ""
#define THIS_MODULE_OPTIONS "-:RVabdefghijoqs" GMT_OPT("HMm")
#define GMT_N_MODE_NOTSET 0
#define GMT_N_MODE_REPORT 1
#define GMT_N_MODE_ADD_ID 2
#define GMT_N_MODE_CLOUD 3
#define GMT_W 3
#define POL_UNION 1
#define POL_INTERSECTION 2
#define POL_CLIP 3
#define POL_SPLIT 4
#define POL_JOIN 5
#define POL_HOLE 6
#define POL_BUFFER 7
#define POL_CENTROID 8
#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 */
#define MIN_SEPARATION 0 /* If the two closest points for two features are > 0 units apart then they are not the same feature */
#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 */
#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 */
#ifdef HAVE_GEOS
#include <geos_c.h>
int geos_methods(struct GMT_CTRL *GMT, struct GMT_DATASET *D, char *fname, double buf_dist, char *method);
int geos_method_polygon(struct GMT_CTRL *GMT, struct GMT_DATASET *Din, struct GMT_DATASET *Dout, char *method);
int geos_method_linestring(struct GMT_CTRL *GMT, struct GMT_DATASET *Din, struct GMT_DATASET *Dout, double buf_dist, char *method);
#endif
struct GMTSPATIAL_DUP { /* Holds information on which single segment is closest to the current test segment */
uint64_t point;
uint64_t segment;
unsigned int table;
int mode;
bool inside;
double distance;
double mean_distance;
double closeness;
double setratio;
double a_threshold;
double d_threshold;
double c_threshold;
double s_threshold;
};
struct DUP_INFO {
uint64_t point;
int mode;
double distance;
double closeness;
double setratio;
};
struct GMTSPATIAL_CTRL {
struct GMTSPATIAL_Out { /* -> */
bool active;
char *file;
} Out;
struct GMTSPATIAL_A { /* -Aa<min_dist>, -A */
bool active;
unsigned int mode;
int smode;
double min_dist;
char unit;
} A;
struct GMTSPATIAL_C { /* -C */
bool active;
} C;
struct GMTSPATIAL_D { /* -D[pol] */
bool active;
int mode;
char unit;
char *file;
struct GMTSPATIAL_DUP I;
} D;
struct GMTSPATIAL_E { /* -E+n|p */
bool active;
unsigned int mode;
} E;
struct GMTSPATIAL_F { /* -F */
bool active;
unsigned int geometry;
} F;
struct GMTSPATIAL_I { /* -I[i|e] */
bool active;
unsigned int mode;
} I;
struct GMTSPATIAL_L { /* -L */
bool active;
char unit;
double s_cutoff, path_noise, box_offset;
} L;
struct GMTSPATIAL_N { /* -N<file>[+a][+p>ID>][+r][+z] */
bool active;
bool all; /* All points in lines and polygons must be inside a polygon for us to report ID */
bool p_set; /* True if +p was used */
unsigned int mode; /* 0 for reporting ID in -Z<ID> header, 1 via data column, 2 just as a report */
unsigned int ID; /* If 1 we use running numbers */
char *file;
} N;
struct GMTSPATIAL_Q { /* -Q[<unit>][+c<min>/<max>][+h][+l][+p][+s[a|d]] */
bool active;
bool header; /* Place dimension and centroid in segment headers */
bool area; /* Apply range test on dimension */
bool sort; /* Sort segments based on dimension */
unsigned int mode; /* 0 use input as is, 1 force to line, 2 force to polygon */
unsigned int dmode; /* for geo data: 1 = flat earth, 2 = great circle, 3 = geodesic (for distances) */
int dir; /* For segment sorting: -1 is descending, +1 is ascending [Default] */
double limit[2]; /* Min and max area or length for output segments */
char unit;
} Q;
struct GMTSPATIAL_S { /* -S[u|i|c|j|h] */
bool active;
unsigned int mode;
double width;
} S;
struct GMTSPATIAL_T { /* -T[pol] */
bool active;
char *file;
} T;
struct GMTSPATIAL_W { /* -W<length>[unit][+f|l] */
bool active;
bool extend[2];
unsigned int smode[2];
double length[2];
char unit[2];
} W;
};
struct GMTSPATIAL_PAIR {
double node;
uint64_t pos;
};
static void *New_Ctrl (struct GMT_CTRL *GMT) { /* Allocate and initialize a new control structure */
struct GMTSPATIAL_CTRL *C;
C = gmt_M_memory (GMT, NULL, 1, struct GMTSPATIAL_CTRL);
/* Initialize values whose defaults are not 0/false/NULL */
C->A.unit = 'X'; /* Cartesian units as default */
C->L.box_offset = 1.0e-10; /* Minimum significant amplitude */
C->L.path_noise = 1.0e-10; /* Minimum significant amplitude */
C->D.unit = 'X'; /* Cartesian units as default */
C->D.I.a_threshold = MIN_AREA_DIFF;
C->D.I.d_threshold = MIN_SEPARATION;
C->D.I.c_threshold = MIN_CLOSENESS;
C->D.I.s_threshold = MIN_SUBSET;
C->N.mode = GMT_N_MODE_NOTSET;
C->Q.mode = GMT_IS_POINT; /* Undecided on line vs poly */
C->Q.dmode = GMT_GREATCIRCLE; /* Great-circle distance if not specified */
return (C);
}
static void Free_Ctrl (struct GMT_CTRL *GMT, struct GMTSPATIAL_CTRL *C) { /* Deallocate control structure */
if (!C) return;
gmt_M_str_free (C->Out.file);
gmt_M_str_free (C->D.file);
gmt_M_str_free (C->N.file);
gmt_M_str_free (C->T.file);
gmt_M_free (GMT, C);
}
GMT_LOCAL unsigned int gmtspatial_area_size (struct GMT_CTRL *GMT, double x[], double y[], uint64_t n, double *out, int geo) {
double size = gmt_centroid_area (GMT, x, y, n, geo, out);
out[GMT_Z] = fabs (size);
return ((size < 0.0) ? GMT_POL_IS_CCW : GMT_POL_IS_CW);
}
#if 0
GMT_LOCAL unsigned int gmtspatial_gmtspatial_area_size_old (struct GMT_CTRL *GMT, double x[], double y[], uint64_t n, double *out, int geo) {
uint64_t i;
double wesn[4], xx, yy, size, ix, iy;
double *xp = NULL, *yp = NULL;
wesn[XLO] = wesn[YLO] = DBL_MAX;
wesn[XHI] = wesn[YHI] = -DBL_MAX;
gmt_centroid (GMT, x, y, n, out, geo); /* Get mean location */
for (i = 0; i < n; i++) {
wesn[XLO] = MIN (wesn[XLO], x[i]);
wesn[XHI] = MAX (wesn[XHI], x[i]);
wesn[YLO] = MIN (wesn[YLO], y[i]);
wesn[YHI] = MAX (wesn[YHI], y[i]);
}
xp = gmt_M_memory (GMT, NULL, n, double); yp = gmt_M_memory (GMT, NULL, n, double);
if (geo == 1) { /* Initializes GMT projection parameters to the -JA settings */
GMT->current.proj.projection_GMT = GMT->current.proj.projection = GMT_LAMB_AZ_EQ;
GMT->current.proj.unit = 1.0;
GMT->current.proj.pars[3] = 39.3700787401574814;
GMT->common.R.oblique = false;
GMT->common.J.active = true;
GMT->current.setting.map_line_step = 1.0e7; /* To avoid nlon/nlat being huge */
GMT->current.proj.pars[0] = out[GMT_X];
GMT->current.proj.pars[1] = out[GMT_Y];
if (gmt_M_err_pass (GMT, gmt_proj_setup (GMT, wesn), "")) {
gmt_M_free (GMT, xp); gmt_M_free (GMT, yp);
return (0);
}
ix = 1.0 / GMT->current.proj.scale[GMT_X];
iy = 1.0 / GMT->current.proj.scale[GMT_Y];
for (i = 0; i < n; i++) {
gmt_geo_to_xy (GMT, x[i], y[i], &xx, &yy);
xp[i] = (xx - GMT->current.proj.origin[GMT_X]) * ix;
yp[i] = (yy - GMT->current.proj.origin[GMT_Y]) * iy;
}
}
else { /* Just take out mean coordinates */
for (i = 0; i < n; i++) {
xp[i] = x[i] - out[GMT_X];
yp[i] = y[i] - out[GMT_Y];
}
}
size = gmt_pol_area (xp, yp, n);
gmt_M_free (GMT, xp);
gmt_M_free (GMT, yp);
if (geo) size *= (GMT->current.map.dist[GMT_MAP_DIST].scale * GMT->current.map.dist[GMT_MAP_DIST].scale);
out[GMT_Z] = fabs (size);
return ((size < 0.0) ? GMT_POL_IS_CCW : GMT_POL_IS_CW);
}
#endif
GMT_LOCAL void gmtspatial_length_size (struct GMT_CTRL *GMT, double x[], double y[], uint64_t n, double *out) {
uint64_t i;
double length = 0.0, mid, f, *s = NULL;
assert (n > 0);
/* Estimate 'average' position */
s = gmt_M_memory (GMT, NULL, n, double);
for (i = 1; i < n; i++) {
length += gmt_distance (GMT, x[i-1], y[i-1], x[i], y[i]);
s[i] = length;
}
mid = 0.5 * length;
i = 0;
while (s[i] <= mid) i++; /* Find mid point length-wise */
f = (mid - s[i-1]) / (s[i] - s[i-1]);
out[GMT_X] = x[i-1] + f * (x[i] - x[i-1]);
out[GMT_Y] = y[i-1] + f * (y[i] - y[i-1]);
out[GMT_Z] = length;
gmt_M_free (GMT, s);
}
GMT_LOCAL int gmtspatial_comp_pairs (const void *a, const void *b) {
const struct GMTSPATIAL_PAIR *xa = a, *xb = b;
/* Sort on node value */
if (xa->node < xb->node) return (-1);
if (xa->node > xb->node) return (+1);
return (0);
}
GMT_LOCAL void gmtspatial_write_record (struct GMT_CTRL *GMT, double **R, uint64_t n, uint64_t p) {
uint64_t c;
double out[GMT_MAX_COLUMNS];
struct GMT_RECORD Out;
Out.data = out; Out.text = NULL;
for (c = 0; c < n; c++) out[c] = R[c][p];
GMT_Put_Record (GMT->parent, GMT_WRITE_DATA, &Out);
}
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) {
/* Given single line segment S and a dataset of many line segments in D, determine the closest neighbor
* to S in D (call it S'), and if "really close" it might be a duplicate or slight revision to S.
* There might be several features S' in D close to S so we return how many near or exact matches we
* found and pass the details via L. We return 0 if no duplicates are found. If some are found then
* these are the possible types of result per match:
* +1 : S has identical duplicate in D (S == S')
* -1 : S has identical duplicate in D (S == S'), but is reversed
* +2 : S is very close to S', probably a slight revised version
* -2 : S is very close to S', probably a slight revised version, but is reversed
* +3 : Same as +2 but S is likely a subset of S'
* -3 : Same as -2 but S is likely a reversed subset of S'
* +4 : Same as +2 but S is likely a superset of S'
* -4 : Same as -2 but S is likely a reversed superset of S'
* 5 : Two lines split exactly at the Dateline
* The algorithm first finds the smallest separation from any point on S to the line Sp.
* We then reverse the situation and check the smallest separation from any point on any
* of the Sp candidates to the line S.
* If the smallest distance found exceeds I->d_threshold then we move on (not close enough).
* We next compute the mean (or median, with +Ccmax) closeness, defined as the ratio between
* the mean (or median) separation between points on S and the line Sp (or vice versa) and
* the length of S (or Sp). We compute it both ways since the segments may differ in lengths
* and degree of overlap (i.e., subsets or supersets).
* NOTE: Lines that intersect obviously has a zero min distance but since we only compute
* distances from the nodes on one line to the nearest point along the other line we will
* not detect such crossings. For most situations this should not matter much (?).
*/
bool status;
unsigned int k, tbl, n_close = 0, n_dup = 0, mode1, mode3;
uint64_t row, seg, pt, np, sno, *n_sum = NULL;
int k_signed;
double dist, f_seg, f_pt, d1, d2, closest, length[2], separation[2], close[2], *d_mean = NULL;
double med_separation[2], med_close[2], high = 0, low = 0, use_length, use_sep, use_close, *sep = NULL;
struct GMT_DATASEGMENT *Sp = NULL;
struct GMT_DATASEGMENT_HIDDEN *SH = NULL;
GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Determine the segments in D closest to our segment\n");
I->distance = I->mean_distance = DBL_MAX;
mode3 = 3 + 10 * I->inside; /* Set gmt_near_lines modes */
mode1 = 1 + 10 * I->inside; /* Set gmt_near_lines modes */
d_mean = gmt_M_memory (GMT, NULL, D->n_segments, double); /* Mean distances from points along S to other lines */
n_sum = gmt_M_memory (GMT, NULL, D->n_segments, uint64_t); /* Number of distances from points along S to other lines */
/* We first want to obtain the mean distance from one segment to all others. We do this by computing the
* nearest distance from each point along our segment S to the other segments and compute the average of
* all those distances. Then we reverse the search and continue to accumulate averages */
/* Process each point along the trace in S and find sum of nearest distances for each segment in table */
for (row = 0; row < S->n_rows; row++) {
for (sno = tbl = 0; tbl < D->n_tables; tbl++) { /* For each table to compare it to */
for (seg = 0; seg < D->table[tbl]->n_segments; seg++, sno++) { /* For each segment in current table */
if (D->table[tbl]->segment[seg]->n_rows == 0) continue; /* Skip segments with no records (may be itself) */
dist = DBL_MAX; /* Reset for each line to find distance to that line */
(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);
d_mean[sno] += dist; /* Sum of distances to this segment */
n_sum[sno]++; /* Number of such distances */
}
}
}
/* Must also do the reverse test: for each point for each line in the table, compute distance to S; it might be shorter */
for (sno = tbl = 0; tbl < D->n_tables; tbl++) { /* For each table to compare it to */
for (seg = 0; seg < D->table[tbl]->n_segments; seg++, sno++) { /* For each segment in current table */
Sp = D->table[tbl]->segment[seg]; /* This is S', one of the segments that is close to S */
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 */
dist = DBL_MAX; /* Reset for each line to find distance to that line */
(void) gmt_near_a_line (GMT, Sp->data[GMT_X][row], Sp->data[GMT_Y][row], seg, S, mode3, &dist, &f_seg, &f_pt);
d_mean[sno] += dist; /* Sum of distances to this segment */
n_sum[sno]++; /* Number of such distances */
}
}
}
/* Compute the average distances */
for (sno = 0; sno < D->n_segments; sno++) {
d_mean[sno] = (n_sum[sno] > 0) ? d_mean[sno] / n_sum[sno] : DBL_MAX;
}
gmt_M_free (GMT, n_sum);
/* Process each point along the trace in S and find nearest distance for each segment in table */
for (row = 0; row < S->n_rows; row++) {
for (sno = tbl = 0; tbl < D->n_tables; tbl++) { /* For each table to compare it to */
for (seg = 0; seg < D->table[tbl]->n_segments; seg++, sno++) { /* For each segment in current table */
SH = gmt_get_DS_hidden (D->table[tbl]->segment[seg]);
dist = DBL_MAX; /* Reset for each line to find distance to that line */
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);
if (!status && I->inside) continue; /* Only consider points that project perpendicularly within the line segment */
pt = lrint (f_pt); /* We know f_seg == seg so no point assigning that */
if (dist < I->distance && d_mean[sno] < I->mean_distance) { /* Keep track of the single closest feature */
I->point = pt;
I->segment = seg;
I->distance = dist;
I->mean_distance = d_mean[sno];
I->table = tbl;
}
if (dist > I->d_threshold) continue; /* Not close enough for duplicate consideration */
if (SH->mode == 0) {
n_close++;
L[tbl][seg].distance = DBL_MAX;
}
SH->mode = 1; /* Use mode to flag segments that are close enough */
if (dist < L[tbl][seg].distance) { /* Keep track of the closest feature */
L[tbl][seg].point = pt;
L[tbl][seg].distance = dist;
}
}
}
}
/* Must also do the reverse test: for each point for each line in the table, compute distance to S; it might be shorter */
for (sno = tbl = 0; tbl < D->n_tables; tbl++) { /* For each table to compare it to */
for (seg = 0; seg < D->table[tbl]->n_segments; seg++, sno++) { /* For each segment in current table */
Sp = D->table[tbl]->segment[seg]; /* This is S', one of the segments that is close to S */
SH = gmt_get_DS_hidden (Sp);
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 */
dist = DBL_MAX; /* Reset for each line to find distance to that line */
status = gmt_near_a_line (GMT, Sp->data[GMT_X][row], Sp->data[GMT_Y][row], seg, S, mode3, &dist, &f_seg, &f_pt);
if (!status && I->inside) continue; /* Only consider points that project perpendicularly within the line segment */
pt = lrint (f_pt);
if (dist < I->distance && d_mean[sno] < I->mean_distance) { /* Keep track of the single closest feature */
I->point = pt;
I->segment = seg;
I->distance = dist;
I->mean_distance = d_mean[sno];
I->table = tbl;
}
if (dist > I->d_threshold) continue; /* Not close enough for duplicate consideration */
if (SH->mode == 0) {
n_close++;
L[tbl][seg].distance = DBL_MAX;
}
SH->mode = 1; /* Use mode to flag segments that are close enough */
if (dist < L[tbl][seg].distance) { /* Keep track of the closest feature */
L[tbl][seg].point = pt;
L[tbl][seg].distance = dist;
}
}
}
}
gmt_M_free (GMT, d_mean);
if (n_close == 0)
GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "No other segment found within dmax [probably due to +p requirement]\n");
else
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);
/* Here we have found the shortest distance from a point on S to another segment; S' is segment number seg */
if (I->distance > I->d_threshold) return (0); /* We found no duplicates or slightly different versions of S */
/* Here we need to compare one or more S' candidates and S a bit more. */
for (row = 1, length[0] = 0.0; row < S->n_rows; row++) { /* Compute length of S once */
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]);
}
for (tbl = 0; tbl < D->n_tables; tbl++) { /* For each table to compare it to */
for (seg = 0; seg < D->table[tbl]->n_segments; seg++) { /* For each segment in current table */
SH = gmt_get_DS_hidden (D->table[tbl]->segment[seg]);
if (SH->mode != 1) continue; /* Not one of the close segments we flagged earlier */
SH->mode = 0; /* Remove this temporary flag */
Sp = D->table[tbl]->segment[seg]; /* This is S', one of the segments that is close to S */
if (S->n_rows == Sp->n_rows) { /* Exactly the same number of data points; check for identical duplicate (possibly reversed) */
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 */
d1 += gmt_distance (GMT, S->data[GMT_X][row], S->data[GMT_Y][row], Sp->data[GMT_X][row], Sp->data[GMT_Y][row]);
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]);
}
GMT_Report (GMT->parent, GMT_MSG_DEBUG, "S to Sp of same length and gave d1 = %g and d2 = %g\n", d1, d2);
if (gmt_M_is_zero (d1)) { /* Exact duplicate */
L[tbl][seg].mode = +1;
n_dup++;
continue;
}
else if (gmt_M_is_zero (d2)) { /* Exact duplicate but reverse order */
L[tbl][seg].mode = -1;
n_dup++;
continue;
}
}
/* We get here when S' is not an exact duplicate of S, but approximate.
* We compute the mean [and possibly median] separation between S' and S
* and the closeness ratio (separation/length). */
separation[0] = 0.0;
if (I->mode) { /* Various items needed to compute median separation */
sep = gmt_M_memory (GMT, NULL, MAX (S->n_rows, Sp->n_rows), double);
low = DBL_MAX;
high = -DBL_MAX;
}
for (row = np = 0; row < S->n_rows; row++) { /* Process each point along the trace in S */
dist = DBL_MAX;
status = gmt_near_a_line (GMT, S->data[GMT_X][row], S->data[GMT_Y][row], seg, Sp, mode1, &dist, NULL, NULL);
if (!status && I->inside) continue; /* Only consider points that project perpendicularly within the Sp segment */
separation[0] += dist;
if (I->mode) { /* Special processing for median calculation */
sep[np] = dist;
if (dist < low) low = dist;
if (dist > high) high = dist;
}
np++; /* Number of points within the overlap zone */
}
separation[0] = (np > 1) ? separation[0] / np : DBL_MAX; /* Mean distance between S and S' */
use_length = (np) ? length[0] * np / S->n_rows : length[0]; /* ~reduce length to overlap section assuming equal point spacing */
close[0] = (np > 1) ? separation[0] / use_length : DBL_MAX; /* Closeness as viewed from S */
use_sep = (separation[0] == DBL_MAX) ? GMT->session.d_NaN : separation[0];
use_close = (close[0] == DBL_MAX) ? GMT->session.d_NaN : close[0];
GMT_Report (GMT->parent, GMT_MSG_INFORMATION,
"S has length %.3f km, has mean separation to Sp of %.3f km, and a closeness ratio of %g [n = %" PRIu64 "/%" PRIu64 "]\n",
length[0], use_sep, use_close, np, S->n_rows);
if (I->mode) {
if (np > 1) {
gmt_median (GMT, sep, np, low, high, separation[0], &med_separation[0]);
med_close[0] = med_separation[0] / use_length;
}
else med_close[0] = DBL_MAX;
use_sep = (np == 0 || med_separation[0] == DBL_MAX) ? GMT->session.d_NaN : med_separation[0];
use_close = (np == 0 || med_close[0] == DBL_MAX) ? GMT->session.d_NaN : med_close[0];
GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "S has median separation to Sp of %.3f km, and a robust closeness ratio of %g\n",
use_sep, use_close);
}
/* Must now compare the other way */
separation[1] = length[1] = 0.0;
if (I->mode) { /* Reset for median calculation */
low = +DBL_MAX;
high = -DBL_MAX;
}
for (row = np = 0; row < Sp->n_rows; row++) { /* Process each point along the trace in S' */
dist = DBL_MAX; /* Reset for each line to find distance to that line */
status = gmt_near_a_line (GMT, Sp->data[GMT_X][row], Sp->data[GMT_Y][row], seg, S, mode1, &dist, NULL, NULL);
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]);
if (!status && I->inside) continue; /* Only consider points that project perpendicularly within the Sp segment */
separation[1] += dist;
if (I->mode) { /* Special processing for median calculation */
sep[np] = dist;
if (dist < low) low = dist;
if (dist > high) high = dist;
}
np++; /* Number of points within the overlap zone */
}
separation[1] = (np > 1) ? separation[1] / np : DBL_MAX; /* Mean distance between S' and S */
use_length = (np) ? length[1] * np / Sp->n_rows : length[1]; /* ~reduce length to overlap section assuming equal point spacing */
close[1] = (np > 1) ? separation[1] / use_length : DBL_MAX; /* Closeness as viewed from S' */
use_sep = (separation[1] == DBL_MAX) ? GMT->session.d_NaN : separation[1];
use_close = (close[1] == DBL_MAX) ? GMT->session.d_NaN : close[1];
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",
length[1], use_sep, use_close, np, Sp->n_rows);
if (I->mode) {
if (np > 1) {
gmt_median (GMT, sep, np, low, high, separation[1], &med_separation[1]);
med_close[1] = med_separation[1] / use_length;
}
else med_close[1] = DBL_MAX;
gmt_M_free (GMT, sep);
use_sep = (med_separation[1] == DBL_MAX) ? GMT->session.d_NaN : med_separation[1];
use_close = (med_close[1] == DBL_MAX) ? GMT->session.d_NaN : med_close[1];
GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Sp has median separation to S of %.3f km, and a robust closeness ratio of %g\n",
use_sep, use_close);
k = (med_close[0] <= med_close[1]) ? 0 : 1; /* Pick the setup with the smallest robust closeness */
closest = med_close[k]; /* The longer the segment and the closer they are, the smaller the closeness */
}
else {
k = (close[0] <= close[1]) ? 0 : 1; /* Pick the setup with the smallest robust closeness */
closest = close[k]; /* The longer the segment and the closer they are, the smaller the closeness */
}
if (closest > I->c_threshold) continue; /* Not a duplicate or slightly different version of S */
L[tbl][seg].closeness = closest; /* The longer the segment and the closer they are the smaller the I->closeness */
/* 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 */
d1 = gmt_distance (GMT, S->data[GMT_X][0], S->data[GMT_Y][0], Sp->data[GMT_X][0], Sp->data[GMT_Y][0]);
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]);
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) */
if (length[0] < (length[1] / I->s_threshold))
k_signed = 3; /* S is a subset of Sp */
else if (length[1] < (length[0] / I->s_threshold))
k_signed = 4; /* S is a superset of Sp */
else
k_signed = 2; /* S and S' are practically the same feature */
L[tbl][seg].mode = (d1 < d2) ? k_signed : -k_signed; /* Negative means Sp is reversed w.r.t. S */
/* Last check: If two lines are split at 180 the above will find them to be approximate duplicates even
* though they just share one point (at 180 longitude) */
if (L[tbl][seg].closeness < GMT_CONV8_LIMIT && L[tbl][seg].distance < GMT_CONV8_LIMIT) {
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))
L[tbl][seg].mode = 5;
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))
L[tbl][seg].mode = 5;
}
n_dup++;
}
}
return (n_dup);
}
struct NN_DIST {
double data[4]; /* Up to x,y,z,weight */
double distance; /* Distance to nearest neighbor */
int64_t ID; /* Input ID # of this point (original input record number 0,1,...)*/
int64_t neighbor; /* Input ID # of this point's neighbor */
};
struct NN_INFO {
int64_t sort_rec; /* Input ID # of this point's neighbor */
int64_t orig_rec; /* Rec # of this point */
};
GMT_LOCAL int gmtspatial_compare_nn_points (const void *point_1v, const void *point_2v) {
/* Routine for qsort to sort NN data structure on distance.
*/
const struct NN_DIST *point_1 = point_1v, *point_2 = point_2v;
if (gmt_M_is_dnan (point_1->distance)) return (+1);
if (gmt_M_is_dnan (point_2->distance)) return (-1);
if (point_1->distance < point_2->distance) return (-1);
if (point_1->distance > point_2->distance) return (+1);
return (0);
}
GMT_LOCAL struct NN_DIST *gmtspatial_NNA_update_dist (struct GMT_CTRL *GMT, struct NN_DIST *P, uint64_t *n_points) {
/* Return array of NN results sorted on smallest distances */
int64_t k, k2, np;
double *distance = gmt_M_memory (GMT, NULL, *n_points, double);
#ifdef DEBUG
static int iteration = 0;
#endif
np = *n_points;
for (k = 0; k < (np-1); k++) {
if (gmt_M_is_dnan (P[k].distance)) continue; /* Skip deleted point */
P[k].distance = DBL_MAX;
/* We split the loop over calculation of distance from the loop over assignments since
* if OpenMP is used then we cannot interchange k and k2 as there may be overprinting.
*/
#ifdef _OPENMP
#pragma omp parallel for private(k2) shared(k,np,GMT,P,distance)
#endif
for (k2 = k + 1; k2 < np; k2++) {
if (gmt_M_is_dnan (P[k2].distance)) continue; /* Skip deleted point */
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]);
}
for (k2 = k + 1; k2 < np; k2++) {
if (gmt_M_is_dnan (P[k2].distance)) continue; /* Skip deleted point */
if (distance[k2] < P[k].distance) {
P[k].distance = distance[k2];
P[k].neighbor = P[k2].ID;
}
if (distance[k2] < P[k2].distance) {
P[k2].distance = distance[k2];
P[k2].neighbor = P[k].ID;
}
}
}
gmt_M_free (GMT, distance);
/* Prefer mergesort since qsort is not stable for equalities */
mergesort (P, np, sizeof (struct NN_DIST), gmtspatial_compare_nn_points); /* Sort on small to large distances */
for (k = np; k > 0 && gmt_M_is_dnan (P[k-1].distance); k--); /* Skip the NaN distances that were placed at end */
*n_points = k; /* Update point count */
#ifdef DEBUG
if (gmt_M_is_verbose (GMT, GMT_MSG_DEBUG)) {
GMT_Report (GMT->parent, GMT_MSG_DEBUG, "===> Iteration = %d\n", iteration);
for (k = 0; k < (int64_t)(*n_points); k++)
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);
}
iteration++;
#endif
return (P);
}
GMT_LOCAL struct NN_DIST *gmtspatial_NNA_init_dist (struct GMT_CTRL *GMT, struct GMT_DATASET *D, uint64_t *n_points) {
/* Return array of NN results sorted on smallest distances */
uint64_t tbl, seg, row, col, n_cols;
int64_t k, np = 0; /* Must be signed due to Win OpenMP retardedness */
double *distance = NULL;
struct GMT_DATASEGMENT *S = NULL;
struct NN_DIST *P = gmt_M_memory (GMT, NULL, D->n_records, struct NN_DIST);
n_cols = MIN (D->n_columns, 4); /* Expects lon,lat and makes room for at least z, w and other columns */
np = (int64_t)(D->n_records * (D->n_records - 1)) / 2;
distance = gmt_M_memory (GMT, NULL, np, double);
for (tbl = 0, np = 0; tbl < D->n_tables; tbl++) {
for (seg = 0; seg < D->table[tbl]->n_segments; seg++) {
S = D->table[tbl]->segment[seg];
for (row = 0; row < S->n_rows; row++) {
for (col = 0; col < n_cols; col++) P[np].data[col] = S->data[col][row]; /* Duplicate coordinates */
if (n_cols < 4) P[np].data[GMT_W] = 1.0; /* No weight provided, set to unity */
P[np].ID = (uint64_t)np; /* Assign ID based on input record # from 0 */
P[np].distance = DBL_MAX;
/* We split the loop over calculation of distance from the loop over assignments since
* if OpenMP is used then we cannot interchange k and np as there may be overprinting.
*/
#ifdef _OPENMP
#pragma omp parallel for private(k) shared(np,GMT,P,distance)
#endif
for (k = 0; k < np; k++) { /* Get distances to other points */
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]);
}
for (k = 0; k < np; k++) { /* Compare this point to all previous points */
if (distance[k] < P[k].distance) { /* Update shortest distance so far, and with which neighbor */
P[k].distance = distance[k];
P[k].neighbor = np;
}
if (distance[k] < P[np].distance) { /* Update shortest distance so far, and with which neighbor */
P[np].distance = distance[k];
P[np].neighbor = k;
}
}
np++;
}
}
}
gmt_M_free (GMT, distance);
/* Prefer mergesort since qsort is not stable for equalities */
mergesort (P, np, sizeof (struct NN_DIST), gmtspatial_compare_nn_points);
*n_points = (uint64_t)np;
#ifdef DEBUG
if (gmt_M_is_verbose (GMT, GMT_MSG_DEBUG)) {
GMT_Report (GMT->parent, GMT_MSG_DEBUG, "===> Initialization\n");
for (k = 0; k < (int64_t)(*n_points); k++)
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);
}
#endif
return (P);
}
GMT_LOCAL int gmtspatial_compare_nn_info (const void *point_1v, const void *point_2v) {
/* Routine for qsort to sort NN rec numbers structure on original record order.
*/
const struct NN_INFO *point_1 = point_1v, *point_2 = point_2v;
if (point_1->orig_rec < point_2->orig_rec) return (-1);
if (point_1->orig_rec > point_2->orig_rec) return (+1);
return (0);
}
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) {
/* Return revised array of NN ID lookups via sorting on neighbor IDs */
uint64_t k;
struct NN_INFO *info = (I) ? I : gmt_M_memory (GMT, NULL, n_points, struct NN_INFO);
for (k = 0; k < n_points; k++) {
info[k].sort_rec = k;
info[k].orig_rec = int64_abs (NN_dist[k].ID);
}
/* Prefer mergesort since qsort is not stable for equalities */
mergesort (info, n_points, sizeof (struct NN_INFO), gmtspatial_compare_nn_info);
/* Now, I[k].sort_rec will take the original record # k and return the corresponding record in the sorted array */
return (info);
}
static int usage (struct GMTAPI_CTRL *API, int level) {
const char *name = gmt_show_name_and_purpose (API, THIS_MODULE_LIB, THIS_MODULE_CLASSIC_NAME, THIS_MODULE_PURPOSE);
if (level == GMT_MODULE_PURPOSE) return (GMT_NOERROR);
GMT_Usage (API, 0, "usage: %s [<table>] [-A[a<min_dist>]] [-C] [-D[+a<amax>][+c|C<cmax>][+d<dmax>][+f<file>][+p][+s<factor>]] [-E+n|p] "
"[-F[l]] [-I[i|e]] [-L%s/<noise>/<offset>] [-N<pfile>[+a][+i][+p[<ID>]][+r][+z]] [-Q[<unit>][+c<min>[/<max>]][+h][+l][+p][+s[a|d]]] [%s] "
"[-Sb<width>|h|s] [-T[<cpol>]] [-W<dist>[<unit>][+f|l]] [%s] [%s] [%s] [%s] [%s] [%s] [%s] [%s] [%s] [%s] [%s] [%s] [%s] [%s] [%s]\n", name, GMT_DIST_OPT, GMT_Rgeo_OPT,
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);
if (level == GMT_SYNOPSIS) return (GMT_MODULE_SYNOPSIS);
GMT_Message (API, GMT_TIME_NONE, " OPTIONAL ARGUMENTS:\n");
GMT_Option (API, "<");
GMT_Usage (API, 1, "\n-A[a<min_dist>]");
GMT_Usage (API, -2, "Nearest Neighbor (NN) Analysis. Compute minimum distances between NN point pairs. "
"Append unit used for NN distance calculation. Returns minimum distances and point IDs for pairs. "
"Use -Aa to replace close neighbor pairs with their weighted average location until "
"no point pair has a NN distance less than the specified <min_dist> distance [0]. "
"Considers 3rd column as z (if present) and 4th as w, if present [weight = 1].");
GMT_Usage (API, 1, "\n-C Clip polygons to the given region box (requires -R), possibly yielding new closed polygons. "
"For truncation instead (possibly yielding open polygons, i.e., lines), see -T.");
GMT_Usage (API, 1, "\n-D[+a<amax>][+c|C<cmax>][+d<dmax>][+f<file>][+p][+s<factor>]");
GMT_Usage (API, -2, "Look for (near-)duplicate lines or polygons. Duplicate lines have a minimum point separation less than <dmax> and a closeness "
"ratio (mean separation/length) less than <cmax>. "
"If near-duplicates have lengths that differ by <sfact> or more then they are subsets or supersets. "
"To flag duplicate polygons, the fractional difference in areas must be less than <amax>. "
"By default, we consider all points when comparing two lines. Change these settings via modifiers:");
GMT_Usage (API, 3, "+a Set limit of fractional difference of polygon areas [0.01].");
GMT_Usage (API, 3, "+c Set closeness ratio (mean separation/length) [0.01]. "
"Use +C to use median separation instead.");
GMT_Usage (API, 3, "+d Set minimum mean point separation [0].");
GMT_Usage (API, 3, "+f Compare <table> against <file> instead of itself.");
GMT_Usage (API, 3, "+p limit comparison to points that project perpendicularly on to the other line [all points].");
GMT_Usage (API, 3, "+s Set length ratio difference threshold factor [2].");
GMT_Usage (API, 1, "\n-E+n|p");
GMT_Usage (API, -2, "Orient all polygons to have the same handedness:");
GMT_Usage (API, 3, "+n Impose clockwise (negative) handedness.");
GMT_Usage (API, 3, "+p Impose counter-clockwise (positive) handedness.");
GMT_Usage (API, 1, "\n-F[l]");
GMT_Usage (API, -2, "Force all input segments to become closed polygons on output by adding repeated point if needed. "
"Use -Fl instead to ensure input lines are not treated as polygons.");
GMT_Usage (API, 1, "\n-I[i|e]");
GMT_Usage (API, -2, "Compute Intersection locations between input polygon(s). Optionally append a directive:");
GMT_Usage (API, 3, "e: Compute external crossings only [Default is both].");
GMT_Usage (API, 3, "i: Compute internal crossings only [Default is both].");
GMT_Usage (API, -2, "Use uppercase E or I to consider all segments within the same table as one entity [separate].");
GMT_Usage (API, 1, "\n-L%s/<noise>/<offset>", GMT_DIST_OPT);
GMT_Usage (API, -2, "Remove tile Lines. These are superfluous lines along the -R border. "
"Append <dist> (in m) [0], coordinate noise [1e-10], and max offset from gridline [1e-10].");
GMT_Usage (API, 1, "\n-N<pfile>[+a][+i][+p[<ID>]][+r][+z]");
GMT_Usage (API, -2, "Determine ID of polygon (in <pfile>) enclosing each input feature (or point if +i). By default, "
"the ID starts at 0 and auto-increments. Alternatively, the ID is set as follows:");
GMT_Usage (API, 3, "%s If shapefile or OGR/GMT polygons, get polygon ID via -a for Z column.", GMT_LINE_BULLET);
GMT_Usage (API, 3, "%s Interpret segment value (-Z<value>) as polygon ID.", GMT_LINE_BULLET);
GMT_Usage (API, 3, "%s Interpret segment label (-L<label>) as polygon ID, else we increment IDs from 0.", GMT_LINE_BULLET);
GMT_Usage (API, -2, "Additional modifiers are available:");
GMT_Usage (API, 3, "+a All points of a feature (line, polygon) must be inside the ID polygon [any point].");
GMT_Usage (API, 3, "+i Determine polygon IDs of all individual input points, append ID to each output row as a new data column.");
GMT_Usage (API, 3, "+p Change origin for auto-incrementing polygon IDs [0].");
GMT_Usage (API, 3, "+r No table output; just report which polygon a feature is inside.");
GMT_Usage (API, 3, "+z Append the ID as a new output data column [Default adds -Z<ID> to segment header].");
GMT_Usage (API, 1, "\n-Q[<unit>][+c<min>[/<max>]][+h][+l][+p][+s[a|d]]");
GMT_Usage (API, -2, "Measure area and handedness of polygon(s) or length of line segments. If -fg is used "
"you may append a <unit> from %s [k]; otherwise it will be based on the input Cartesian data unit. "
"We also compute polygon centroid or line mid-point. See documentation for more information. Optional modifiers:", GMT_LEN_UNITS_DISPLAY);
GMT_Usage (API, 3, "+c Limit output segments to those with area or length within specified range [output all segments]. "
"If <max> is not given then it defaults to infinity.");
GMT_Usage (API, 3, "+h Place the (area, handedness) or length result in the segment header on output.");
GMT_Usage (API, 3, "+p Consider all input as polygons and close them if necessary "
"[only closed polygons are considered polygons].");
GMT_Usage (API, 3, "+l Consider all input as lines even if closed [closed polygons are considered polygons].");
GMT_Usage (API, 3, "+s Sort segments based on area or length; append a for ascending or d for descending [ascending].");
GMT_Usage (API, -2, "[Default only reports results to standard output].\n");
GMT_Option (API, "R");
GMT_Usage (API, 1, "\n-Sb<width>|h|i|j|s|u");
GMT_Usage (API, -2, "Spatial manipulation of polygons; choose a directive:");
#ifdef HAVE_GEOS
GMT_Usage (API, 3, "b: Compute buffer polygon around line/polygon. Append <width> of buffer zone. "
"Note: this is a purely Cartesian operation so <width> must be in data units.");
#endif
GMT_Usage (API, 3, "h: Detect holes and reverse them relative to perimeters.");
//GMT_Usage (API, 3, "i: Find intersection [Not implemented yet].");
//GMT_Usage (API, 3, "j: Join polygons that were split by the Dateline [Not implemented yet].");
GMT_Usage (API, 3, "s: Split polygons that straddle the Dateline.");
//GMT_Usage (API, 3, "u: Find union [Not implemented yet].");
GMT_Usage (API, 1, "\n-T[<cpol>]");
GMT_Usage (API, -2, "Truncate polygons against the clip polygon <cpol>; if <cpol> is not given we require -R "
"and clip against a polygon derived from the region border.");
GMT_Option (API, "V");
GMT_Usage (API, 1, "\n-W<dist>[<unit>][+f|l]");
GMT_Usage (API, -2, "Extend all segments with extra first and last points that are <dist> units away from the "
"original end points in the directions implied by the line ends. For geographic data append a unit from %s [k]. "
" To extend the two ends by different amounts, give <distf>[<unit>]/<distl>[<unit>] instead. "
"By default we extend both ends of the segments. Choose one modifier to change this:", GMT_LEN_UNITS_DISPLAY);
GMT_Usage (API, 3, "+f Only extend the start of the segments.");
GMT_Usage (API, 3, "+l Only extend the end of the segments.");
GMT_Option (API, "a,bi2,bo,d,e,f,g,h,i,j,o,q,s,:,.");
return (GMT_MODULE_USAGE);
}
static int parse (struct GMT_CTRL *GMT, struct GMTSPATIAL_CTRL *Ctrl, struct GMT_OPTION *options) {
/* This parses the options provided to grdsample and sets parameters in CTRL.
* Any GMT common options will override values set previously by other commands.
* It also replaces any file names specified as input or output with the data ID
* returned when registering these sources/destinations with the API.
*/
bool got_i = false;
unsigned int pos, n_errors = 0, got_n = 0;
int n;
char txt_a[GMT_LEN64] = {""}, txt_b[GMT_LEN64] = {""}, txt_c[GMT_LEN64] = {""}, p[GMT_LEN256] = {""};
char *s = NULL, *c = NULL, *q = NULL;
struct GMT_OPTION *opt = NULL;
struct GMTAPI_CTRL *API = GMT->parent;
if (gmt_M_is_geographic (GMT, GMT_IN)) Ctrl->Q.unit = 'k'; /* Default geographic distance unit is km */
for (opt = options; opt; opt = opt->next) {
switch (opt->option) {
case '<': /* Skip input files after checking the exist */
if (GMT_Get_FilePath (API, GMT_IS_DATASET, GMT_IN, GMT_FILE_REMOTE, &(opt->arg))) n_errors++;
break;
case '>': /* Got named output file */
n_errors += gmt_M_repeated_module_option (API, Ctrl->Out.active);
n_errors += gmt_get_required_file (GMT, opt->arg, opt->option, 0, GMT_IS_DATASET, GMT_OUT, GMT_FILE_LOCAL, &(Ctrl->Out.file));
break;
/* Processes program-specific parameters */
case 'A': /* Do nearest neighbor analysis */
n_errors += gmt_M_repeated_module_option (API, Ctrl->A.active);
if (opt->arg[0] == 'a' || opt->arg[0] == 'A') { /* Spatially average points until minimum NN distance is less than given distance */
Ctrl->A.mode = (opt->arg[0] == 'A') ? 2 : 1; /* Slow mode is an undocumented test mode */
Ctrl->A.smode = gmt_get_distance (GMT, &opt->arg[1], &(Ctrl->A.min_dist), &(Ctrl->A.unit));
}
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]))) {
/* Just compute NN distances and return the unique ones */
if (opt->arg[0] == '-') /* Flat earth calculation */
Ctrl->A.smode = 1;
else if (opt->arg[0] == '+') /* Geodetic */
Ctrl->A.smode = 3;
else
Ctrl->A.smode = 2; /* Great circle */
Ctrl->A.unit = (Ctrl->A.smode == 2) ? opt->arg[0] : opt->arg[1];
if (gmt_M_is_cartesian (GMT, GMT_IN)) { /* Data was not geographic, revert to Cartesian settings */
Ctrl->A.smode = 0;
Ctrl->A.unit = 'X';
}
}
else if (opt->arg[0]) {
GMT_Report (API, GMT_MSG_ERROR, "Bad modifier in %c in -A option\n", opt->arg[0]);
n_errors++;
}
break;
case 'C': /* Clip to given region */
n_errors += gmt_M_repeated_module_option (API, Ctrl->C.active);
n_errors += gmt_get_no_argument (GMT, opt->arg, opt->option, 0);
break;
case 'D': /* Look for duplications */
n_errors += gmt_M_repeated_module_option (API, Ctrl->D.active);
pos = 0;
while (gmt_strtok (opt->arg, "+", &pos, p)) {
switch (p[0]) {
case 'a': /* Gave a new +a<dmax> value */
GMT_Report (API, GMT_MSG_ERROR, "+a not implemented yet\n");
Ctrl->D.I.a_threshold = atof (&p[1]);
break;
case 'd': /* Gave a new +d<dmax> value */
Ctrl->D.mode = gmt_get_distance (GMT, &p[1], &(Ctrl->D.I.d_threshold), &(Ctrl->D.unit));
break;
case 'C': /* Gave a new +C<cmax> value */
Ctrl->D.I.mode = 1; /* Median instead of mean */
/* Intentionally fall through */
case 'c': /* Gave a new +c<cmax> value */
if (p[1]) Ctrl->D.I.c_threshold = atof (&p[1]); /* This allows +C by itself just to change to median */
break;
case 's': /* Gave a new +s<fact> value */
Ctrl->D.I.s_threshold = atof (&p[1]);
break;
case 'p': /* Consider only inside projections */
Ctrl->D.I.inside = 1;
break;
case 'f': /* Gave a file name */
Ctrl->D.file = strdup (&p[1]);
if (GMT_Get_FilePath (API, GMT_IS_DATASET, GMT_IN, GMT_FILE_REMOTE, &(Ctrl->D.file))) n_errors++;
break;
}
}
break;
case 'E': /* Orient polygons -E+n|p (old -E-|+) */
n_errors += gmt_M_repeated_module_option (API, Ctrl->E.active);
if (opt->arg[0] == '-' || strstr (opt->arg, "+n"))
Ctrl->E.mode = GMT_POL_IS_CW;
else if (opt->arg[0] == '+' || strstr (opt->arg, "+p"))
Ctrl->E.mode = GMT_POL_IS_CCW;
else
n_errors++;
break;
case 'F': /* Force polygon or line mode */
n_errors += gmt_M_repeated_module_option (API, Ctrl->F.active);
Ctrl->F.geometry = (opt->arg[0] == 'l') ? GMT_IS_LINE : GMT_IS_POLY;
break;
case 'I': /* Compute intersections between polygons */
n_errors += gmt_M_repeated_module_option (API, Ctrl->I.active);
if (opt->arg[0] == 'i') Ctrl->I.mode = 1;
if (opt->arg[0] == 'I') Ctrl->I.mode = 2;
if (opt->arg[0] == 'e') Ctrl->I.mode = 4;
if (opt->arg[0] == 'E') Ctrl->I.mode = 8;
break;
case 'L': /* Remove tile lines */
n_errors += gmt_M_repeated_module_option (API, Ctrl->L.active);
n = sscanf (opt->arg, "%[^/]/%[^/]/%s", txt_a, txt_b, txt_c);
if (n >= 1) Ctrl->L.s_cutoff = atof (txt_a);
if (n >= 2) Ctrl->L.path_noise = atof (txt_b);
if (n == 3) Ctrl->L.box_offset = atof (txt_c);
break;
case 'N': /* Determine containing polygons for features */
n_errors += gmt_M_repeated_module_option (API, Ctrl->N.active);
if (gmt_validate_modifiers (GMT, opt->arg, 'N', "aiprz", GMT_MSG_ERROR)) n_errors++;
if ((s = strchr (opt->arg, '+')) == NULL) { /* No modifiers */
Ctrl->N.file = strdup (opt->arg);
continue;
}
else { /* Hide modifiers until we duplicate the polygon name */
s[0] = '\0';
Ctrl->N.file = strdup (opt->arg);
s[0] = '+';
}
q = gmt_first_modifier (GMT, opt->arg, "aiprz");
pos = 0; txt_a[0] = 0;
while (gmt_getmodopt (GMT, 'N', q, "aiprz", &pos, txt_a, &n_errors) && n_errors == 0) {
/* Note: +i can only be combined with +p */
switch (txt_a[0]) {