diff --git a/src/gmt_shore.c b/src/gmt_shore.c index 30a7091538f..71390ae6613 100644 --- a/src/gmt_shore.c +++ b/src/gmt_shore.c @@ -1422,6 +1422,23 @@ int gmt_assemble_shore (struct GMT_CTRL *GMT, struct GMT_SHORE *c, int dir, bool if (c->ns > 0) p = gmt_M_memory (GMT, p, P, struct GMT_GSHHS_POL); /* Trim memory */ for (ku = 0; ku < P; ku++) gmtshore_path_shift2 (p[ku].lon, p[ku].n, west, east, c->leftmost_bin); /* Deal with possible longitude -/+360 issues */ + if (c->lon_sw >= east) { /* This bin sits at or beyond the east border, so its polygons were shifted west above. + * A point that lay exactly on that border is not shifted by the test above, which tears the polygon across + * the entire map: its bounding box then spans every column and the containment test paints unrelated nodes + * clear across the grid [issue #8235]. A polygon cannot be wider than its bin, so when one is we know the + * points left behind belong with the shifted part. */ + for (ku = 0; ku < P; ku++) { + double lon_min = DBL_MAX, lon_max = -DBL_MAX; + unsigned int kp; + for (kp = 0; kp < (unsigned int)p[ku].n; kp++) { + if (p[ku].lon[kp] < lon_min) lon_min = p[ku].lon[kp]; + if (p[ku].lon[kp] > lon_max) lon_max = p[ku].lon[kp]; + } + if ((lon_max - lon_min) <= c->bsize) continue; /* This one is intact */ + for (kp = 0; kp < (unsigned int)p[ku].n; kp++) + if (p[ku].lon[kp] >= east && (p[ku].lon[kp] - 360.0) >= west) p[ku].lon[kp] -= 360.0; + } + } *pol = p; return (P); /* Return list of polygons found */ diff --git a/src/grdlandmask.c b/src/grdlandmask.c index a215ac78009..03b7ead7621 100644 --- a/src/grdlandmask.c +++ b/src/grdlandmask.c @@ -439,12 +439,24 @@ EXTERN_MSC int GMT_grdlandmask (void *V_API, int mode, void *args) { if (p[k].lat[i] < ymin) ymin = p[k].lat[i]; if (p[k].lat[i] > ymax) ymax = p[k].lat[i]; } - col_min = MAX (0, irint (ceil (xmin * i_dx_inch - Grid->header->xy_off - GMT_CONV8_LIMIT))); - if (col_min > nx1) col_min = 0; - /* So col_min is in range [0,nx1] */ - col_max = MIN (nx1, irint (floor (xmax * i_dx_inch - Grid->header->xy_off + GMT_CONV8_LIMIT))); - if (col_max <= 0 || col_max < col_min) col_max = nx1; - /* So col_max is in range [1,nx1] */ + col_min = irint (ceil (xmin * i_dx_inch - Grid->header->xy_off - GMT_CONV8_LIMIT)); + col_max = irint (floor (xmax * i_dx_inch - Grid->header->xy_off + GMT_CONV8_LIMIT)); + if (wrap) { /* Keep the wrap-around treatment of a global grid */ + if (col_min < 0 || col_min > nx1) col_min = 0; + /* So col_min is in range [0,nx1] */ + if (col_max > nx1) col_max = nx1; + if (col_max <= 0 || col_max < col_min) col_max = nx1; + /* So col_max is in range [1,nx1] */ + } + /* For a non-global grid the same fallbacks would turn a round-off miss just outside the grid into + * the full grid width, so the winding test then ran over unrelated columns and made a mess of the + * mask [issue #8235]. Clip a genuine overlap instead, and skip a polygon that has none. */ + else if (col_min > nx1 || col_max < 0 || col_max < col_min) + continue; + else { + if (col_min < 0) col_min = 0; + if (col_max > nx1) col_max = nx1; + } row_min = MAX (0, irint (ceil ((GMT->current.proj.rect[YHI] - ymax) * i_dy_inch - Grid->header->xy_off - GMT_CONV8_LIMIT))); /* So row_min is in range [0,?] */ row_max = MIN (ny1, irint (floor ((GMT->current.proj.rect[YHI] - ymin) * i_dy_inch - Grid->header->xy_off + GMT_CONV8_LIMIT))); @@ -604,14 +616,15 @@ EXTERN_MSC int GMT_grdlandmask (void *V_API, int mode, void *args) { if (col_max < col_min) col_max += (int)Grid->header->n_columns; } else { /* Make sure we are inside our grid */ + /* lon_w and lon_e are the bin's edges measured east from the grid's west edge, so they must be + * wrapped against 360, not compared to the absolute wesn limits: doing the latter shifted a bin + * lying east of the west edge into negative offsets, so col_max came out below col_min and the + * bin painted nothing at all [issue #8235] */ double lon_w, lon_e; - lon_w = c.lon_sw - Grid->header->wesn[XLO]; lon_e = c.lon_sw + c.bsize - Grid->header->wesn[XLO]; - if (lon_w < Grid->header->wesn[XLO] && (lon_w+360.0) < Grid->header->wesn[XHI]) { - lon_w += 360.0; lon_e += 360.0; - } - else if (lon_e > Grid->header->wesn[XHI] && (lon_e-360.0) > Grid->header->wesn[XLO]) { - lon_w -= 360.0; lon_e -= 360.0; - } + lon_w = c.lon_sw - Grid->header->wesn[XLO]; lon_e = lon_w + c.bsize; + while (lon_w < 0.0) { lon_w += 360.0; lon_e += 360.0; } + while (lon_w >= 360.0) { lon_w -= 360.0; lon_e -= 360.0; } + if (lon_e > 360.0) { lon_w -= 360.0; lon_e -= 360.0; } /* Bin straddles the seam; use its western part */ col_min = irint (ceil (lon_w * HH->r_inc[GMT_X] - Grid->header->xy_off)); col_max = irint (floor (lon_e * HH->r_inc[GMT_X] - Grid->header->xy_off)); if (col_min < 0) col_min = 0;