/*-------------------------------------------------------------------- * * Copyright (c) 1991-2020 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 */ #include "gmt_dev.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}" #define THIS_MODULE_NEEDS "" #define THIS_MODULE_OPTIONS "-:RVabdefghijoqs" GMT_OPT("HMm") #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 PW_TESTING #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 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, -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[+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 */ unsigned int mode; /* 0 for reporting ID in -Z 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[+c/][+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_PAIR { double node; uint64_t pos; }; #ifdef __APPLE__ /* macOX has it built in, so ensure we define this flag */ #define HAVE_MERGESORT #endif #ifndef HAVE_MERGESORT #include "mergesort.c" #endif 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->Q.mode = GMT_IS_POINT; /* Undecided on line vs poly */ C->Q.dmode = 2; /* 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); #ifdef PW_TESTING GMT_Message (API, GMT_TIME_NONE, "usage: %s [] [-A[a]] [-C]\n\t[-D[+f][+a][+d%s][+c|C][+l][+s][+p]]\n\t[-E+p|n] [-F[l]] [-I[i|e]] [-L%s//] [-N[+a][+p][+r][+z]]\n\t[-Q[+c[/]][+h][+l][+p][+s[a|d]]] [%s] [-Sb|h|i|j|s|u]\n", name, GMT_DIST_OPT, GMT_DIST_OPT, GMT_Rgeo_OPT); #else GMT_Message (API, GMT_TIME_NONE, "usage: %s [
] [-A[a]] [-C]\n\t[-D[+f][+a][+d%s][+c|C][+l][+s][+p]]\n\t[-E+p|n] [-F[l]] [-I[i|e]] [-N[+a][+p][+r][+z]]\n\t[-Q[+c[/]][+h][+l][+p][+s[a|d]]] [%s] [-Sb|h|i|j|s|u]\n", name, GMT_DIST_OPT, GMT_Rgeo_OPT); #endif GMT_Message (API, GMT_TIME_NONE, "\t[-T[]] [%s] [%s] [%s] [%s] [%s]\n\t[%s] [%s]\n\t[%s] [%s] [%s]\n\t[%s] [%s] [%s] [%s]\n\n", GMT_V_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, "\tOPTIONS:\n"); GMT_Option (API, "<"); GMT_Message (API, GMT_TIME_NONE, "\t-A Nearest Neighbor (NN) Analysis. Compute minimum distances between NN point pairs.\n"); GMT_Message (API, GMT_TIME_NONE, "\t Append unit used for NN distance calculation. Returns minimum distances and point IDs for pairs.\n"); GMT_Message (API, GMT_TIME_NONE, "\t Use -Aa to replace close neighbor pairs with their weighted average location until\n"); GMT_Message (API, GMT_TIME_NONE, "\t no point pair has a NN distance less than the specified distance [0].\n"); GMT_Message (API, GMT_TIME_NONE, "\t Considers 3rd column as z (if present) and 4th as w, if present [weight = 1].\n"); GMT_Message (API, GMT_TIME_NONE, "\t-C Clip polygons to the given region box (requires -R), possibly yielding new closed polygons.\n"); GMT_Message (API, GMT_TIME_NONE, "\t For truncation instead (possibly yielding open polygons, i.e., lines), see -T.\n"); GMT_Message (API, GMT_TIME_NONE, "\t-D Look for (near-)duplicates in
, or append +f to compare
against .\n"); GMT_Message (API, GMT_TIME_NONE, "\t Near-duplicates have a minimum point separation less than [0] and a closeness\n"); GMT_Message (API, GMT_TIME_NONE, "\t ratio (mean separation/length) less than [0.01]. Use +d and +c to change these.\n"); GMT_Message (API, GMT_TIME_NONE, "\t Use +C to use median separation instead [+c uses mean separation].\n"); GMT_Message (API, GMT_TIME_NONE, "\t If near-duplicates have lengths that differ by or more then they are subsets or supersets [2].\n"); GMT_Message (API, GMT_TIME_NONE, "\t To flag duplicate polygons, the fractional difference in areas must be less than [0.01].\n"); GMT_Message (API, GMT_TIME_NONE, "\t By default we consider all points when comparing two lines. Use +p to limit\n"); GMT_Message (API, GMT_TIME_NONE, "\t the comparison to points that project perpendicularly on to the other line.\n"); GMT_Message (API, GMT_TIME_NONE, "\t-E Orient all polygons to have the same handedness.\n"); GMT_Message (API, GMT_TIME_NONE, "\t Append +p for counter-clockwise (positive) or +n for clockwise (negative) handedness.\n"); GMT_Message (API, GMT_TIME_NONE, "\t-F Force all input segments to become closed polygons on output by adding repeated point if needed.\n"); GMT_Message (API, GMT_TIME_NONE, "\t Use -Fl instead to ensure input lines are not treated as polygons.\n"); GMT_Message (API, GMT_TIME_NONE, "\t-I Compute Intersection locations between input polygon(s).\n"); GMT_Message (API, GMT_TIME_NONE, "\t Append e or i for external or internal crossings only [Default is both].\n"); GMT_Message (API, GMT_TIME_NONE, "\t Use uppercase E or I to consider all segments within the same table as one entity [separate].\n"); #ifdef PW_TESTING GMT_Message (API, GMT_TIME_NONE, "\t-L Remove tile Lines. These are superfluous lines along the -R border.\n"); GMT_Message (API, GMT_TIME_NONE, "\t Append (in m) [0], coordinate noise [1e-10], and max offset from gridline [1e-10].\n"); #endif GMT_Message (API, GMT_TIME_NONE, "\t-N Determine ID of polygon (in ) enclosing each input feature. The ID is set as follows:\n"); GMT_Message (API, GMT_TIME_NONE, "\t a) If OGR/GMT polygons, get polygon ID via -a for Z column, else\n"); GMT_Message (API, GMT_TIME_NONE, "\t b) Interpret segment labels (-Z) as polygon IDs, else\n"); GMT_Message (API, GMT_TIME_NONE, "\t c) Interpret segment labels (-L