Skip to content

Commit

Permalink
Fix errors in gmtspatial related to -N (#7746)
Browse files Browse the repository at this point in the history
* Fix errors in gmtspatial related to -N

See #7745 for background.

* Do not write segments that are outside

* Add a test test
  • Loading branch information
PaulWessel authored Aug 20, 2023
1 parent cbb8051 commit f829247
Show file tree
Hide file tree
Showing 2 changed files with 131 additions and 22 deletions.
64 changes: 42 additions & 22 deletions src/gmtspatial.c
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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 <geos_c.h>
Expand Down Expand Up @@ -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);
Expand Down Expand Up @@ -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<fact> value */
Ctrl->N.mode = 2;
case 'z': /* Add polygon ID as final column */
Ctrl->N.mode = GMT_N_MODE_ADD_ID;
break;
}
}
Expand Down Expand Up @@ -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;
Expand All @@ -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);
}
Expand All @@ -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);
Expand All @@ -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))
Expand All @@ -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);
}
Expand Down
89 changes: 89 additions & 0 deletions test/gmtspatial/inout_features.sh
Original file line number Diff line number Diff line change
@@ -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

0 comments on commit f829247

Please sign in to comment.