From 8b76b268e519677d434ccbd2b1b0691ba2dd0a2e Mon Sep 17 00:00:00 2001 From: Joaquim Date: Thu, 14 May 2026 13:43:00 +0100 Subject: [PATCH 1/4] coast: DCW painting inside/outside issue Fix #5151 Assisted-by: Claude Opus 4.7 --- src/gmt_map.c | 54 ++++++++++++++++++++++++++++++++++++++++----------- 1 file changed, 43 insertions(+), 11 deletions(-) diff --git a/src/gmt_map.c b/src/gmt_map.c index 86355bcf89e..ddbec8d6b8e 100644 --- a/src/gmt_map.c +++ b/src/gmt_map.c @@ -1976,18 +1976,50 @@ uint64_t gmt_map_wesn_clip (struct GMT_CTRL *GMT, double *lon, double *lat, uint } /*! . */ -GMT_LOCAL uint64_t gmtmap_radial_boundary_arc (struct GMT_CTRL *GMT, int this_way, double end_x[], double end_y[], double **xarc, double **yarc) { +GMT_LOCAL bool gmtmap_arc_long_way (struct GMT_CTRL *GMT, double end_x[], double end_y[], double *lon, double *lat, uint64_t np) { + /* Decide which horizon arc connects a pair of polygon crossings (issue #5151). + * end_x[]/end_y[] are the two crossing points expressed as offsets from the map center. + * We take the midpoint of the SHORT arc, invert-project it back to lon/lat, and test + * whether that midpoint lies inside the original spherical polygon. If it does NOT, + * the short arc sits outside the polygon -- so the polygon's interior must be reached + * via the LONG (complementary) arc, e.g. when the polygon encloses the projection + * center. */ + double az1, az2, d_az, mid_az, xr, yr, mx, my, mlon, mlat, lon_ref; + if (np < 3) return false; + az1 = d_atan2d(end_y[0], end_x[0]); + az2 = d_atan2d(end_y[1], end_x[1]); + gmt_M_set_delta_lon(az1, az2, d_az); + mid_az = az1 + 0.5 * d_az; + sincosd(mid_az, &yr, &xr); + mx = GMT->current.proj.r * (1.0 + xr); + my = GMT->current.proj.r * (1.0 + yr); + gmt_xy_to_geo(GMT, &mlon, &mlat, mx, my); + /* Bring mlon into the polygon's longitude range to avoid 360 wrap mismatches */ + lon_ref = lon[0]; + while (mlon - lon_ref > 180.0) mlon -= 360.0; + while (mlon - lon_ref < -180.0) mlon += 360.0; + return (gmt_non_zero_winding(GMT, mlon, mlat, lon, lat, np) == GMT_OUTSIDE); +} + +/*! . */ +GMT_LOCAL uint64_t gmtmap_radial_boundary_arc (struct GMT_CTRL *GMT, int this_way, double end_x[], double end_y[], double **xarc, double **yarc, bool long_way) { uint64_t n_arc, k, pt; double az1, az2, d_az, da, xr, yr, da_try, *xx = NULL, *yy = NULL; /* When a polygon crosses out then in again into the circle we need to add a boundary arc * to the polygon where it is clipped. We simply sample the circle as finely as the arc - * length and the current line_step demands */ + * length and the current line_step demands. If long_way is true we instead trace the + * complementary (long) arc around the horizon -- needed when the polygon encloses the + * projection center so that the polygon's interior wraps around the visible disk. */ da_try = (GMT->current.setting.map_line_step * 360.0) / (TWO_PI * GMT->current.proj.r); /* Angular step in degrees */ az1 = d_atan2d (end_y[0], end_x[0]); /* azimuth from map center to 1st crossing */ az2 = d_atan2d (end_y[1], end_x[1]); /* azimuth from map center to 2nd crossing */ - gmt_M_set_delta_lon (az1, az2, d_az); /* Insist we take the short arc for now */ + gmt_M_set_delta_lon(az1, az2, d_az); /* Short arc delta in [-180, 180] */ + if (long_way) { /* Take the complementary arc instead */ + if (d_az > 0.0) d_az -= 360.0; + else d_az += 360.0; + } n_arc = lrint (ceil (fabs (d_az) / da_try)); /* Get number of integer increments of da_try degree... */ if (n_arc < 2) n_arc = 2; /* ...but minimum 2 */ da = d_az / (n_arc - 1); /* Reset da to get exact steps */ @@ -2055,13 +2087,12 @@ GMT_LOCAL uint64_t gmtmap_radial_clip (struct GMT_CTRL *GMT, double *lon, double add_boundary = !this_side; /* We only add boundary arcs if we first exited and now entered the circle again */ } if (add_boundary) { /* Crossed twice. Now add arc between the two crossing points */ - /* PW: Currently, we make the assumption that the shortest arc is the one we want. However, - * extremely large polygons could cut the boundary so that it is the longest arc we want. - * The way to improve this algorithm in the future is to find the two opposite points on - * the circle boundary that lies on the bisector of az1,az2, and see which point lies - * inside the polygon. This would require that gmt_inonout_sphpol be called. - */ - if ((n_arc = gmtmap_radial_boundary_arc (GMT, this_side, &end_x[nx-2], &end_y[nx-2], &xarc, &yarc)) > 0) { + /* Choose between short and long arc by testing the midpoint of the short arc: + * project it back to lon/lat and check if it lies inside the original polygon. + * If outside, the polygon's interior is on the far side of the disk and we must + * take the long arc instead (issue #5151). */ + bool long_way = gmtmap_arc_long_way(GMT, &end_x[nx-2], &end_y[nx-2], lon, lat, np); + if ((n_arc = gmtmap_radial_boundary_arc(GMT, this_side, &end_x[nx-2], &end_y[nx-2], &xarc, &yarc, long_way)) > 0) { if ((n + n_arc) >= n_alloc) gmt_M_malloc2 (GMT, xx, yy, n + n_arc, &n_alloc, double); gmt_M_memcpy (&xx[n], xarc, n_arc, double); /* Copy longitudes of arc */ gmt_M_memcpy (&yy[n], yarc, n_arc, double); /* Copy latitudes of arc */ @@ -2084,7 +2115,8 @@ GMT_LOCAL uint64_t gmtmap_radial_clip (struct GMT_CTRL *GMT, double *lon, double } if (nx == 2) { /* Must close polygon by adding boundary arc */ - if ((n_arc = gmtmap_radial_boundary_arc (GMT, this_side, end_x, end_y, &xarc, &yarc)) > 0) { + bool long_way = gmtmap_arc_long_way(GMT, end_x, end_y, lon, lat, np); + if ((n_arc = gmtmap_radial_boundary_arc(GMT, this_side, end_x, end_y, &xarc, &yarc, long_way)) > 0) { if ((n + n_arc) >= n_alloc) gmt_M_malloc2 (GMT, xx, yy, n + n_arc, &n_alloc, double); gmt_M_memcpy (&xx[n], xarc, n_arc, double); /* Copy longitudes of arc */ gmt_M_memcpy (&yy[n], yarc, n_arc, double); /* Copy latitudes of arc */ From bb783c30dd3517fa5edc208447f49026d612f5fc Mon Sep 17 00:00:00 2001 From: Joaquim Date: Sat, 16 May 2026 12:30:29 +0100 Subject: [PATCH 2/4] New version that still fixes on Windows and doesn't break many other tests. --- src/gmt_map.c | 47 +++++++++++++++++++++++++++++++++++++++-------- 1 file changed, 39 insertions(+), 8 deletions(-) diff --git a/src/gmt_map.c b/src/gmt_map.c index ddbec8d6b8e..185c8ab4a21 100644 --- a/src/gmt_map.c +++ b/src/gmt_map.c @@ -1976,7 +1976,35 @@ uint64_t gmt_map_wesn_clip (struct GMT_CTRL *GMT, double *lon, double *lat, uint } /*! . */ -GMT_LOCAL bool gmtmap_arc_long_way (struct GMT_CTRL *GMT, double end_x[], double end_y[], double *lon, double *lat, uint64_t np) { +GMT_LOCAL bool gmtmap_polygon_encloses_center(struct GMT_CTRL *GMT, double *lon, double *lat, uint64_t np) { + /* Returns true iff (central_meridian, proj.pole) lies inside the spherical polygon + * defined by (lon,lat,np). Used as a gate before deciding the short/long arc closure + * for radial clipping (issue #5151); we only need to consider flipping the arc if the + * polygon actually encloses the projection center -- otherwise the historical short-arc + * default is correct. */ + struct GMT_DATASEGMENT *S = NULL; + bool inside, saved_sph_inside; + double clon, clat; + if (np < 3) return false; + clon = GMT->current.proj.central_meridian; + clat = GMT->current.proj.pole; + if (gmt_M_is_dnan(clon) || gmt_M_is_dnan(clat)) return false; + S = gmt_get_segment(GMT, 2); + if (S == NULL) return false; + gmt_alloc_segment(GMT, S, np, 2, 0, true); + gmt_M_memcpy(S->data[GMT_X], lon, np, double); + gmt_M_memcpy(S->data[GMT_Y], lat, np, double); + gmt_set_seg_minmax(GMT, GMT_IS_POLY, 2, S); + saved_sph_inside = GMT->current.proj.sph_inside; + GMT->current.proj.sph_inside = true; /* Force spherical in-polygon test */ + inside = (gmt_inonout(GMT, clon, clat, S) == GMT_INSIDE); + GMT->current.proj.sph_inside = saved_sph_inside; + gmt_free_segment(GMT, &S); + return inside; +} + +/*! . */ +GMT_LOCAL bool gmtmap_arc_long_way(struct GMT_CTRL *GMT, double end_x[], double end_y[], double *lon, double *lat, uint64_t np) { /* Decide which horizon arc connects a pair of polygon crossings (issue #5151). * end_x[]/end_y[] are the two crossing points expressed as offsets from the map center. * We take the midpoint of the SHORT arc, invert-project it back to lon/lat, and test @@ -2059,7 +2087,7 @@ GMT_LOCAL uint64_t gmtmap_radial_clip (struct GMT_CTRL *GMT, double *lon, double uint64_t n = 0, n_arc; unsigned int i, nx; unsigned int sides[4]; - bool this_side = false, add_boundary = false; + bool this_side = false, add_boundary = false, may_need_long = false; double xlon[4], xlat[4], xc[4], yc[4], end_x[3], end_y[3], xr, yr; double *xx = NULL, *yy = NULL, *xarc = NULL, *yarc = NULL; @@ -2067,6 +2095,11 @@ GMT_LOCAL uint64_t gmtmap_radial_clip (struct GMT_CTRL *GMT, double *lon, double if (np == 0) return (0); + /* Only consider flipping to the long horizon arc if the polygon spherically encloses the + * projection center. Polygons that do not enclose the center keep the historical + * short-arc behavior, which avoids regressing existing tests (issue #5151). */ + may_need_long = gmtmap_polygon_encloses_center(GMT, lon, lat, np); + if (!gmt_map_outside (GMT, lon[0], lat[0])) { gmt_M_malloc2 (GMT, xx, yy, n, &n_alloc, double); gmt_geo_to_xy (GMT, lon[0], lat[0], &xx[0], &yy[0]); @@ -2087,11 +2120,9 @@ GMT_LOCAL uint64_t gmtmap_radial_clip (struct GMT_CTRL *GMT, double *lon, double add_boundary = !this_side; /* We only add boundary arcs if we first exited and now entered the circle again */ } if (add_boundary) { /* Crossed twice. Now add arc between the two crossing points */ - /* Choose between short and long arc by testing the midpoint of the short arc: - * project it back to lon/lat and check if it lies inside the original polygon. - * If outside, the polygon's interior is on the far side of the disk and we must - * take the long arc instead (issue #5151). */ - bool long_way = gmtmap_arc_long_way(GMT, &end_x[nx-2], &end_y[nx-2], lon, lat, np); + /* Only run the short/long arc midpoint test for polygons that enclose the + * projection center; otherwise stay with the historical short arc default. */ + bool long_way = may_need_long && gmtmap_arc_long_way(GMT, &end_x[nx-2], &end_y[nx-2], lon, lat, np); if ((n_arc = gmtmap_radial_boundary_arc(GMT, this_side, &end_x[nx-2], &end_y[nx-2], &xarc, &yarc, long_way)) > 0) { if ((n + n_arc) >= n_alloc) gmt_M_malloc2 (GMT, xx, yy, n + n_arc, &n_alloc, double); gmt_M_memcpy (&xx[n], xarc, n_arc, double); /* Copy longitudes of arc */ @@ -2115,7 +2146,7 @@ GMT_LOCAL uint64_t gmtmap_radial_clip (struct GMT_CTRL *GMT, double *lon, double } if (nx == 2) { /* Must close polygon by adding boundary arc */ - bool long_way = gmtmap_arc_long_way(GMT, end_x, end_y, lon, lat, np); + bool long_way = may_need_long && gmtmap_arc_long_way(GMT, end_x, end_y, lon, lat, np); if ((n_arc = gmtmap_radial_boundary_arc(GMT, this_side, end_x, end_y, &xarc, &yarc, long_way)) > 0) { if ((n + n_arc) >= n_alloc) gmt_M_malloc2 (GMT, xx, yy, n + n_arc, &n_alloc, double); gmt_M_memcpy (&xx[n], xarc, n_arc, double); /* Copy longitudes of arc */ From 8a1177c2eea8f09c363ec5c97e49b96a5d47e331 Mon Sep 17 00:00:00 2001 From: Joaquim Date: Sat, 16 May 2026 12:58:38 +0100 Subject: [PATCH 3/4] A fix to avoid numerical instability and make this fix work on Linux (and likely Mac too) --- src/gmt_map.c | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/src/gmt_map.c b/src/gmt_map.c index de83d637576..22f5e58785d 100644 --- a/src/gmt_map.c +++ b/src/gmt_map.c @@ -2019,8 +2019,13 @@ GMT_LOCAL bool gmtmap_arc_long_way(struct GMT_CTRL *GMT, double end_x[], double gmt_M_set_delta_lon(az1, az2, d_az); mid_az = az1 + 0.5 * d_az; sincosd(mid_az, &yr, &xr); - mx = GMT->current.proj.r * (1.0 + xr); - my = GMT->current.proj.r * (1.0 + yr); + /* Pull the test point slightly inside the horizon (0.99 * r) instead of placing it on + * the horizon edge. At exactly r the inverse projection is numerically singular + * (involves asin(1)/sqrt(0)) and ULP-level libm differences between Windows (MSVC) and + * Linux (glibc) produce different (mlon,mlat) for the same (mid_az), which then flips + * the gmt_non_zero_winding result and the short/long arc decision (issue #5151). */ + mx = GMT->current.proj.r * (1.0 + 0.99 * xr); + my = GMT->current.proj.r * (1.0 + 0.99 * yr); gmt_xy_to_geo(GMT, &mlon, &mlat, mx, my); /* Bring mlon into the polygon's longitude range to avoid 360 wrap mismatches */ lon_ref = lon[0]; From 60b9f316b8ac263beeff88186b5cac60b97b04ad Mon Sep 17 00:00:00 2001 From: Joaquim Date: Sat, 16 May 2026 14:58:15 +0100 Subject: [PATCH 4/4] Another fix to fix a pssolar breaking. --- src/gmt_map.c | 91 ++++++++++++++++++++++++++------------------------- 1 file changed, 46 insertions(+), 45 deletions(-) diff --git a/src/gmt_map.c b/src/gmt_map.c index 22f5e58785d..291032de729 100644 --- a/src/gmt_map.c +++ b/src/gmt_map.c @@ -1976,62 +1976,53 @@ uint64_t gmt_map_wesn_clip (struct GMT_CTRL *GMT, double *lon, double *lat, uint } /*! . */ -GMT_LOCAL bool gmtmap_polygon_encloses_center(struct GMT_CTRL *GMT, double *lon, double *lat, uint64_t np) { - /* Returns true iff (central_meridian, proj.pole) lies inside the spherical polygon - * defined by (lon,lat,np). Used as a gate before deciding the short/long arc closure - * for radial clipping (issue #5151); we only need to consider flipping the arc if the - * polygon actually encloses the projection center -- otherwise the historical short-arc - * default is correct. */ - struct GMT_DATASEGMENT *S = NULL; - bool inside, saved_sph_inside; - double clon, clat; - if (np < 3) return false; - clon = GMT->current.proj.central_meridian; - clat = GMT->current.proj.pole; - if (gmt_M_is_dnan(clon) || gmt_M_is_dnan(clat)) return false; +GMT_LOCAL struct GMT_DATASEGMENT *gmtmap_build_test_segment(struct GMT_CTRL *GMT, double *lon, double *lat, uint64_t np) { + /* Wrap (lon,lat,np) in a temporary GMT_DATASEGMENT so gmt_inonout (spherical mode) can + * be used for in-polygon tests during radial clipping (issue #5151). Caller frees. */ + struct GMT_DATASEGMENT *S; + if (np < 3) return NULL; S = gmt_get_segment(GMT, 2); - if (S == NULL) return false; + if (S == NULL) return NULL; gmt_alloc_segment(GMT, S, np, 2, 0, true); gmt_M_memcpy(S->data[GMT_X], lon, np, double); gmt_M_memcpy(S->data[GMT_Y], lat, np, double); gmt_set_seg_minmax(GMT, GMT_IS_POLY, 2, S); - saved_sph_inside = GMT->current.proj.sph_inside; - GMT->current.proj.sph_inside = true; /* Force spherical in-polygon test */ - inside = (gmt_inonout(GMT, clon, clat, S) == GMT_INSIDE); - GMT->current.proj.sph_inside = saved_sph_inside; - gmt_free_segment(GMT, &S); - return inside; + return S; +} + +/*! . */ +GMT_LOCAL bool gmtmap_polygon_encloses_center(struct GMT_CTRL *GMT, struct GMT_DATASEGMENT *S) { + /* Returns true iff (central_meridian, proj.pole) lies inside the spherical polygon S. + * Used as a gate before deciding the short/long arc closure for radial clipping. */ + double clon, clat; + if (S == NULL) return false; + clon = GMT->current.proj.central_meridian; + clat = GMT->current.proj.pole; + if (gmt_M_is_dnan(clon) || gmt_M_is_dnan(clat)) return false; + return (gmt_inonout(GMT, clon, clat, S) == GMT_INSIDE); } /*! . */ -GMT_LOCAL bool gmtmap_arc_long_way(struct GMT_CTRL *GMT, double end_x[], double end_y[], double *lon, double *lat, uint64_t np) { +GMT_LOCAL bool gmtmap_arc_long_way(struct GMT_CTRL *GMT, double end_x[], double end_y[], struct GMT_DATASEGMENT *S) { /* Decide which horizon arc connects a pair of polygon crossings (issue #5151). * end_x[]/end_y[] are the two crossing points expressed as offsets from the map center. - * We take the midpoint of the SHORT arc, invert-project it back to lon/lat, and test - * whether that midpoint lies inside the original spherical polygon. If it does NOT, - * the short arc sits outside the polygon -- so the polygon's interior must be reached - * via the LONG (complementary) arc, e.g. when the polygon encloses the projection - * center. */ - double az1, az2, d_az, mid_az, xr, yr, mx, my, mlon, mlat, lon_ref; - if (np < 3) return false; + * Take the midpoint of the SHORT arc, pull it inward to 0.99 * r so the inverse map is + * well-conditioned (avoids the asin(1) singularity that causes libm-level Win/Linux + * divergence at the horizon edge), invert-project to (lon,lat) and ask gmt_inonout in + * SPHERICAL mode (handles dateline-spanning / hemispheric polygons like pssolar's + * terminator) whether it lies inside the polygon. If it does NOT, the short arc sits + * outside the polygon -- so we close along the LONG (complementary) arc instead. */ + double az1, az2, d_az, mid_az, xr, yr, mx, my, mlon, mlat; + if (S == NULL) return false; az1 = d_atan2d(end_y[0], end_x[0]); az2 = d_atan2d(end_y[1], end_x[1]); gmt_M_set_delta_lon(az1, az2, d_az); mid_az = az1 + 0.5 * d_az; sincosd(mid_az, &yr, &xr); - /* Pull the test point slightly inside the horizon (0.99 * r) instead of placing it on - * the horizon edge. At exactly r the inverse projection is numerically singular - * (involves asin(1)/sqrt(0)) and ULP-level libm differences between Windows (MSVC) and - * Linux (glibc) produce different (mlon,mlat) for the same (mid_az), which then flips - * the gmt_non_zero_winding result and the short/long arc decision (issue #5151). */ mx = GMT->current.proj.r * (1.0 + 0.99 * xr); my = GMT->current.proj.r * (1.0 + 0.99 * yr); gmt_xy_to_geo(GMT, &mlon, &mlat, mx, my); - /* Bring mlon into the polygon's longitude range to avoid 360 wrap mismatches */ - lon_ref = lon[0]; - while (mlon - lon_ref > 180.0) mlon -= 360.0; - while (mlon - lon_ref < -180.0) mlon += 360.0; - return (gmt_non_zero_winding(GMT, mlon, mlat, lon, lat, np) == GMT_OUTSIDE); + return (gmt_inonout(GMT, mlon, mlat, S) == GMT_OUTSIDE); } /*! . */ @@ -2092,18 +2083,24 @@ GMT_LOCAL uint64_t gmtmap_radial_clip (struct GMT_CTRL *GMT, double *lon, double uint64_t n = 0, n_arc; unsigned int i, nx; unsigned int sides[4]; - bool this_side = false, add_boundary = false, may_need_long = false; + bool this_side = false, add_boundary = false, may_need_long = false, saved_sph_inside; double xlon[4], xlat[4], xc[4], yc[4], end_x[3], end_y[3], xr, yr; double *xx = NULL, *yy = NULL, *xarc = NULL, *yarc = NULL; + struct GMT_DATASEGMENT *S_test = NULL; *total_nx = 0; /* Keep track of total of crossings */ if (np == 0) return (0); - /* Only consider flipping to the long horizon arc if the polygon spherically encloses the - * projection center. Polygons that do not enclose the center keep the historical - * short-arc behavior, which avoids regressing existing tests (issue #5151). */ - may_need_long = gmtmap_polygon_encloses_center(GMT, lon, lat, np); + /* Build a temporary segment wrapping the polygon and force spherical in-polygon tests + * so dateline-spanning / hemispheric polygons (e.g., pssolar terminator) are handled + * correctly. Only consider flipping to the long horizon arc if the polygon spherically + * encloses the projection center -- polygons that do not enclose the center keep the + * historical short-arc behavior (issue #5151, no regressions on other tests). */ + saved_sph_inside = GMT->current.proj.sph_inside; + S_test = gmtmap_build_test_segment(GMT, lon, lat, np); + if (S_test) GMT->current.proj.sph_inside = true; + may_need_long = gmtmap_polygon_encloses_center(GMT, S_test); if (!gmt_map_outside (GMT, lon[0], lat[0])) { gmt_M_malloc2 (GMT, xx, yy, n, &n_alloc, double); @@ -2127,7 +2124,7 @@ GMT_LOCAL uint64_t gmtmap_radial_clip (struct GMT_CTRL *GMT, double *lon, double if (add_boundary) { /* Crossed twice. Now add arc between the two crossing points */ /* Only run the short/long arc midpoint test for polygons that enclose the * projection center; otherwise stay with the historical short arc default. */ - bool long_way = may_need_long && gmtmap_arc_long_way(GMT, &end_x[nx-2], &end_y[nx-2], lon, lat, np); + bool long_way = may_need_long && gmtmap_arc_long_way(GMT, &end_x[nx-2], &end_y[nx-2], S_test); if ((n_arc = gmtmap_radial_boundary_arc(GMT, this_side, &end_x[nx-2], &end_y[nx-2], &xarc, &yarc, long_way)) > 0) { if ((n + n_arc) >= n_alloc) gmt_M_malloc2 (GMT, xx, yy, n + n_arc, &n_alloc, double); gmt_M_memcpy (&xx[n], xarc, n_arc, double); /* Copy longitudes of arc */ @@ -2151,7 +2148,7 @@ GMT_LOCAL uint64_t gmtmap_radial_clip (struct GMT_CTRL *GMT, double *lon, double } if (nx == 2) { /* Must close polygon by adding boundary arc */ - bool long_way = may_need_long && gmtmap_arc_long_way(GMT, end_x, end_y, lon, lat, np); + bool long_way = may_need_long && gmtmap_arc_long_way(GMT, end_x, end_y, S_test); if ((n_arc = gmtmap_radial_boundary_arc(GMT, this_side, end_x, end_y, &xarc, &yarc, long_way)) > 0) { if ((n + n_arc) >= n_alloc) gmt_M_malloc2 (GMT, xx, yy, n + n_arc, &n_alloc, double); gmt_M_memcpy (&xx[n], xarc, n_arc, double); /* Copy longitudes of arc */ @@ -2162,6 +2159,10 @@ GMT_LOCAL uint64_t gmtmap_radial_clip (struct GMT_CTRL *GMT, double *lon, double if (n == n_alloc) gmt_M_malloc2 (GMT, xx, yy, n, &n_alloc, double); xx[n] = xx[0]; yy[n] = yy[0]; n++; /* Close the polygon */ } + if (S_test) { + gmt_free_segment(GMT, &S_test); + GMT->current.proj.sph_inside = saved_sph_inside; + } n_alloc = n; gmt_M_malloc2 (GMT, xx, yy, 0U, &n_alloc, double); *x = xx;