From f82924741c403e38c869996685bf21dd6bd238ce Mon Sep 17 00:00:00 2001 From: Paul Wessel Date: Sun, 20 Aug 2023 17:58:12 +0200 Subject: [PATCH] Fix errors in gmtspatial related to -N (#7746) * Fix errors in gmtspatial related to -N See #7745 for background. * Do not write segments that are outside * Add a test test --- src/gmtspatial.c | 64 ++++++++++++++-------- test/gmtspatial/inout_features.sh | 89 +++++++++++++++++++++++++++++++ 2 files changed, 131 insertions(+), 22 deletions(-) create mode 100755 test/gmtspatial/inout_features.sh diff --git a/src/gmtspatial.c b/src/gmtspatial.c index b5a4b1b86a8..f1a15c20b57 100644 --- a/src/gmtspatial.c +++ b/src/gmtspatial.c @@ -38,7 +38,10 @@ #define THIS_MODULE_NEEDS "" #define THIS_MODULE_OPTIONS "-:RVabdefghijoqs" GMT_OPT("HMm") -#define GMT_W 3 +#define GMT_N_MODE_NOTSET 0 +#define GMT_N_MODE_REPORT 1 +#define GMT_N_MODE_ADD_ID 2 +#define GMT_W 3 #define POL_UNION 1 #define POL_INTERSECTION 2 @@ -49,10 +52,10 @@ #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 */ +#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 @@ -179,6 +182,7 @@ static void *New_Ctrl (struct GMT_CTRL *GMT) { /* Allocate and initialize a new 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); @@ -987,10 +991,10 @@ static int parse (struct GMT_CTRL *GMT, struct GMTSPATIAL_CTRL *Ctrl, struct GMT Ctrl->N.ID = (p[1]) ? atoi (&p[1]) : 1; break; case 'r': /* Just give a report */ - Ctrl->N.mode = 1; + Ctrl->N.mode = GMT_N_MODE_REPORT; break; - case 'z': /* Gave a new +s value */ - Ctrl->N.mode = 2; + case 'z': /* Add polygon ID as final column */ + Ctrl->N.mode = GMT_N_MODE_ADD_ID; break; } } @@ -1962,10 +1966,11 @@ EXTERN_MSC int GMT_gmtspatial (void *V_API, int mode, void *args) { } if (Ctrl->N.active) { /* Report the polygons that contain the given features */ + bool check_next; uint64_t tbl, row, first, last, n, p, np, seg, seg2, n_inside; unsigned int *count = NULL, nmode; int ID = -1; - char seg_label[GMT_LEN64] = {""}, record[GMT_BUFSIZ] = {""}, *kind[2] = {"Middle point", "All points"}; + char seg_label[GMT_LEN64] = {""}, record[GMT_BUFSIZ] = {""}, *kind[2] = {"%%d points", "All points"}; struct GMT_DATASET *C = NULL; struct GMT_DATATABLE *T = NULL; struct GMT_DATASEGMENT *S = NULL, *S2 = NULL; @@ -1980,7 +1985,7 @@ EXTERN_MSC int GMT_gmtspatial (void *V_API, int mode, void *args) { Return (GMT_DIM_TOO_SMALL); } gmt_reenable_bghio_opts (GMT); /* Recover settings provided by user (if -b -g -h -i were used at all) */ - nmode = (Ctrl->N.mode == 1) ? GMT_IS_NONE : GMT_IS_LINE; + nmode = (Ctrl->N.mode == GMT_N_MODE_REPORT) ? GMT_IS_NONE : GMT_IS_LINE; 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 */ Return (API->error); } @@ -1991,7 +1996,7 @@ EXTERN_MSC int GMT_gmtspatial (void *V_API, int mode, void *args) { Return (API->error); } - if (Ctrl->N.mode == 2) gmt_adjust_dataset (GMT, D, D->n_columns + 1); /* Add one more output column */ + if (Ctrl->N.mode == GMT_N_MODE_ADD_ID) gmt_adjust_dataset (GMT, D, D->n_columns + 1); /* Add one more output column */ T = C->table[0]; /* Only one input file so only one table */ count = gmt_M_memory (GMT, NULL, D->n_segments, unsigned int); @@ -2018,24 +2023,36 @@ EXTERN_MSC int GMT_gmtspatial (void *V_API, int mode, void *args) { for (seg = 0; seg < D->table[tbl]->n_segments; seg++, p++) { S = D->table[tbl]->segment[seg]; if (S->n_rows == 0) continue; - if (Ctrl->N.all) { first = 0; last = S->n_rows - 1; np = S->n_rows; } else { first = last = S->n_rows / 2; np = 1; } - for (row = first, n = 0; row <= last; row++) { /* Check one or all points if they are inside */ - n += (gmt_inonout (GMT, S->data[GMT_X][row], S->data[GMT_Y][row], S2) == GMT_INSIDE); + for (row = n = 0, check_next = true; check_next && row < S->n_rows; row++) { /* Check one or all points if they are inside */ + if (Ctrl->N.all && n < row) /* At least one point has been found outside so with +a we stop checking for more */ + check_next = false; + else /* Need to see if next point is inside (which always takes place for the first point row = 0) */ + n += (gmt_inonout (GMT, S->data[GMT_X][row], S->data[GMT_Y][row], S2) == GMT_INSIDE); + } + if (n == 0 || (Ctrl->N.all && n < S->n_rows)) { /* Nothing inside this polygon or not all of the points are inside (+a) */ + SH = gmt_get_DS_hidden (S); + SH->mode = GMT_WRITE_SKIP; /* Do not output this feature */ + continue; } - if (n < np) continue; /* Not inside this polygon */ if (count[p]) { GMT_Report (API, GMT_MSG_ERROR, "Segment %" PRIu64 "-%" PRIu64 " already inside another polygon; skipped\n", tbl, seg); continue; } count[p]++; - /* Here we are inside */ - if (Ctrl->N.mode == 1) { /* Just report on which polygon contains each feature */ - sprintf (record, "%s from table %" PRIu64 " segment %" PRIu64 " is inside polygon # %d", kind[Ctrl->N.all], tbl, seg, ID); + /* Here the feature was fully (+a) or partly inside the polygon */ + if (Ctrl->N.mode == GMT_N_MODE_REPORT) { /* Just report on which polygon contains each feature */ + if (Ctrl->N.all) + sprintf (record, "All points from table %" PRIu64 " segment %" PRIu64 " is inside polygon # %d", tbl, seg, ID); + else + sprintf (record, "%" PRIu64 " points from table %" PRIu64 " segment %" PRIu64 " is inside polygon # %d", n, tbl, seg, ID); GMT_Put_Record (API, GMT_WRITE_DATA, &Out); } - else if (Ctrl->N.mode == 2) { /* Add ID as last data column */ + else if (Ctrl->N.mode == GMT_N_MODE_ADD_ID) { /* Add ID as last data column */ for (row = 0, n = S->n_columns-1; row < S->n_rows; row++) S->data[n][row] = (double)ID; - GMT_Report (API, GMT_MSG_INFORMATION, "%s from table %" PRIu64 " segment %" PRIu64 " is inside polygon # %d\n", kind[Ctrl->N.all], tbl, seg, ID); + if (Ctrl->N.all) + GMT_Report (API, GMT_MSG_INFORMATION, "All points from table %" PRIu64 " segment %" PRIu64 " is inside polygon # %d\n", tbl, seg, ID); + else + GMT_Report (API, GMT_MSG_INFORMATION, "%" PRIu64 " points from table %" PRIu64 " segment %" PRIu64 " is inside polygon # %d\n", n, tbl, seg, ID); } else { /* Add ID via the segment header -Z */ if (gmt_parse_segment_item (GMT, S->header, "-Z", NULL)) @@ -2047,14 +2064,17 @@ EXTERN_MSC int GMT_gmtspatial (void *V_API, int mode, void *args) { sprintf (txt, " -Z%d", ID); strcat (buffer, txt); S->header = strdup (buffer); - GMT_Report (API, GMT_MSG_INFORMATION, "%s from table %" PRIu64 " segment %" PRIu64 " is inside polygon # %d\n", kind[Ctrl->N.all], tbl, seg, ID); + if (Ctrl->N.all) + GMT_Report (API, GMT_MSG_INFORMATION, "All points from table %" PRIu64 " segment %" PRIu64 " is inside polygon # %d\n", tbl, seg, ID); + else + GMT_Report (API, GMT_MSG_INFORMATION, "%" PRIu64 " points from table %" PRIu64 " segment %" PRIu64 " is inside polygon # %d\n", n, tbl, seg, ID); } } } } } for (p = n_inside = 0; p < D->n_segments; p++) if (count[p]) n_inside++; - if (Ctrl->N.mode != 1) { /* Write out results */ + if (n_inside && Ctrl->N.mode != GMT_N_MODE_REPORT) { /* Write out results */ if (GMT_Write_Data (API, GMT_IS_DATASET, GMT_IS_FILE, GMT_IS_POLY, GMT_WRITE_SET, NULL, Ctrl->Out.file, D) != GMT_NOERROR) { Return (API->error); } diff --git a/test/gmtspatial/inout_features.sh b/test/gmtspatial/inout_features.sh new file mode 100755 index 00000000000..482df301102 --- /dev/null +++ b/test/gmtspatial/inout_features.sh @@ -0,0 +1,89 @@ +#!/usr/bin/env bash +# +# Check that gmtspatial -N works for various modifiers +# -Npoly : Write out features who fully or partially are inside polygon +# -Npoly+z : Same, but place polygon ID as last column per row +# -Npoly+r : Same, just report if fully or partially are inside polygon +# -Npoly+a : Write out features who are fully inside polygon +# -Npoly+z+a : Same, but place polygon ID as last column per row +# -Npoly+r+a : Same, but just report if fully inside polygon +# The polygon +cat << EOF > test_poly.txt +> -Z999 +140 10 +140 20 +150 20 +150 10 +140 10 +EOF +# The features +cat << EOF > test_points.txt +> feature 1 crosses into the polygon from the outside +138 15 +139 15 +140 15 +141 15 +142 15 +143 15 +> feature 2 is completely outside +138 12 +139 12 +> feature 3 is completely inside +141 18 +142 18 +EOF + +echo "# Check -Ntest_poly.txt [Feature 1 and 3 expected]" > result +gmt spatial test_points.txt -Ntest_poly.txt >> result +echo "# Check -Ntest_poly.txt+z [Feature 1 and 3 expected]" >> result +gmt spatial test_points.txt -Ntest_poly.txt+z >> result +echo "# Check -Ntest_poly.txt+r [Feature 1 and 3 expected]" >> result +gmt spatial test_points.txt -Ntest_poly.txt+r >> result +echo "# Check -Ntest_poly.txt+a [Feature 3 expected]" >> result +gmt spatial test_points.txt -Ntest_poly.txt+a >> result +echo "# Check -Ntest_poly.txt+z+a [Feature 3 expected]" >> result +gmt spatial test_points.txt -Ntest_poly.txt+z+a >> result +echo "# Check -Ntest_poly.txt+r+a [Feature 3 expected]" >> result +gmt spatial test_points.txt -Ntest_poly.txt+r+a >> result + +# This is the verified expected outputs: +cat << EOF > answer +# Check -Ntest_poly.txt [Feature 1 and 3 expected] +> feature 1 crosses into the polygon from the outside -Z999 +138 15 +139 15 +140 15 +141 15 +142 15 +143 15 +138 15 +> feature 3 is completely inside -Z999 +141 18 +142 18 +# Check -Ntest_poly.txt+z [Feature 1 and 3 expected] +> feature 1 crosses into the polygon from the outside +138 15 999 +139 15 999 +140 15 999 +141 15 999 +142 15 999 +143 15 999 +138 15 999 +> feature 3 is completely inside +141 18 999 +142 18 999 +# Check -Ntest_poly.txt+r [Feature 1 and 3 expected] +3 points from table 0 segment 0 is inside polygon # 999 +2 points from table 0 segment 2 is inside polygon # 999 +# Check -Ntest_poly.txt+a [Feature 3 expected] +> feature 3 is completely inside -Z999 +141 18 +142 18 +# Check -Ntest_poly.txt+z+a [Feature 3 expected] +> feature 3 is completely inside +141 18 999 +142 18 999 +# Check -Ntest_poly.txt+r+a [Feature 3 expected] +All points from table 0 segment 2 is inside polygon # 999 +EOF +diff -q --strip-trailing-cr answer result > fail