new segfault with GEOS 3.10.0beta2, OK in 3.9.1

[rsbivand]

> library(sf)
Linking to GEOS 3.10.0beta2, GDAL 3.3.2, PROJ 8.1.1
>   p1 <- st_as_sfc("POLYGON((0 0, 0 10, 10 10, 10 0, 0 0))")
>   p2 <- st_as_sfc("POLYGON((0 0, 0 10, 10 0))")
>   set1 <- c(p1, p2, p1, p2, p1, p2)
> st_is_valid(set1)

 *** caught segfault ***
address (nil), cause 'memory not mapped'

Traceback:
 1: CPL_geos_is_valid(x, as.logical(NA_on_exception))
 2: st_is_valid.sfc(set1)
 3: st_is_valid(set1)

Possible actions:
1: abort (with core dump, if enabled)
2: normal R exit
3: exit R without saving workspace
4: exit R saving workspace
Selection: 2
Save workspace image? [y/n/c]: n

and

> library(sf)
Linking to GEOS 3.10.0beta2, GDAL 3.3.2, PROJ 8.1.1
> library(rgeos)
Loading required package: sp
rgeos version: 0.5-8, (SVN revision 680)
 GEOS runtime version: 3.10.0beta2-CAPI-1.16.0 
 Please note that rgeos will be retired by the end of 2023,
plan transition to sf functions using GEOS at your earliest convenience.
 GEOS using OverlayNG
 Linking to sp version: 1.4-6 
 Polygon checking: TRUE 

> p1 <- st_as_sfc("POLYGON((0 0, 0 10, 10 10, 10 0, 0 0))")
> p2 <- st_as_sfc("POLYGON((0 0, 0 10, 10 0))")
> set1 <- c(p1, p2, p1, p2, p1, p2)
> sp_set1 <- as(set1, "Spatial")
Error in CPL_geos_is_empty(st_geometry(x)) : 
  Evaluation error: IllegalArgumentException: Points of LinearRing do not form a closed linestring.
# as in 3.9.1, unchanged
> p1_sp <- readWKT("POLYGON((0 0, 0 10, 10 10, 10 0, 0 0))")
> p2_sp <- readWKT("POLYGON((0 0, 0 10, 10 0))")
Error: Unable to parse: POLYGON((0 0, 0 10, 10 0))
GEOS reported: "rgeos_readWKT: unable to read wkt"
# as in 3.9.1, unchanged

Simplified:

> p2 <- st_as_sfc("POLYGON((0 0, 0 10, 10 0))")
> st_is_valid(p2)

 *** caught segfault ***
address (nil), cause 'memory not mapped'

Traceback:
 1: CPL_geos_is_valid(x, as.logical(NA_on_exception))
 2: st_is_valid.sfc(p2)
 3: st_is_valid(p2)

Possible actions:
1: abort (with core dump, if enabled)
2: normal R exit
3: exit R without saving workspace
4: exit R saving workspace

[rsbivand]

Looks like NEWS -> - Fix IsValidOp to correctly report certain kinds of invalid LinearRings (Martin Davis) @dr-jts ?
NEWS.txt

[rsbivand]

> library(sf)
Linking to GEOS 3.10.0beta2, GDAL 3.3.2, PROJ 8.1.1
> p2 <- st_as_sfc("POLYGON((0 0, 0 10, 10 0))")
> st_is_empty(p2)
Error in CPL_geos_is_empty(st_geometry(x)) : 
  Evaluation error: IllegalArgumentException: Points of LinearRing do not form a closed linestring.

[rsbivand]

p2 maybe should have failed to be created by CPL_sfc_from_wkt(); by GDALs OGRGeometryFactory::createFromWkt():

// [[Rcpp::export]]
Rcpp::List CPL_sfc_from_wkt(Rcpp::CharacterVector wkt) {
	std::vector<OGRGeometry *> g(wkt.size());
	OGRGeometryFactory f;
	for (int i = 0; i < wkt.size(); i++) {
		char *wkt_str = wkt(i);
#if GDAL_VERSION_MAJOR <= 2 && GDAL_VERSION_MINOR <= 2
		handle_error(f.createFromWkt(&wkt_str, NULL, &(g[i])));
#else
		handle_error(f.createFromWkt( (const char*) wkt_str, NULL, &(g[i])));
#endif
	}
	return sfc_from_ogr(g, true);
}

[rsbivand]

GDB:

#0  0x00007fffe40cd156 in geos::operation::valid::IsValidOp::isValidGeometry(geos::geom::Geometry const*) () from /usr/local/lib64/libgeos.so.3.10.0
#1  0x00007fffe40cd22c in geos::operation::valid::IsValidOp::getValidationError() () from /usr/local/lib64/libgeos.so.3.10.0
#2  0x00007fffe4ebc59d in GEOSisValid_r () from /usr/local/lib64/libgeos_c.so.1
#3  0x00007fffe6a1d58a in CPL_geos_is_valid (sfc=..., 
    NA_on_exception=NA_on_exception@entry=true) at geos.cpp:566
#4  0x00007fffe69ed2b6 in _sf_CPL_geos_is_valid (sfcSEXP=0x4444970, 
    NA_on_exceptionSEXP=<optimized out>) at RcppExports.cpp:592
#5  0x00007ffff7c09dec in R_doDotCall (ofun=<optimized out>, 
    nargs=nargs@entry=2, cargs=cargs@entry=0x7fffffff3d50, 
    call=call@entry=0x3f11c98) at ../../../src/main/dotcode.c:601
#6  0x00007ffff7c0a36d in do_dotcall (call=0x3f11c98, op=<optimized out>, 
    args=<optimized out>, env=<optimized out>)
    at ../../../src/main/dotcode.c:1281
#7  0x00007ffff7c451e5 in bcEval (body=<optimized out>, rho=<optimized out>, 
    useCache=<optimized out>) at ../../../src/main/eval.c:7117
#8  0x00007ffff7c5ed60 in Rf_eval (e=0x3f11d40, rho=rho@entry=0x3f119f8)
    at ../../../src/main/eval.c:729
...

[rsbivand]

In base sf, you need an extra non-default argument to create a non-closed polygon (sf:::MtrxSet() is called by polygon() with needClosed=TRUE):

sf:::MtrxSet(list(cbind(x=c(0, 0, 10), y=c(0, 10, 0))), dim="XY", type="POLYGON", needClosed=TRUE)
Error in sf:::MtrxSet(list(cbind(x = c(0, 0, 10), y = c(0, 10, 0))), dim = "XY",  : 
  polygons not (all) closed
p2g <- sf:::MtrxSet(list(cbind(x=c(0, 0, 10), y=c(0, 10, 0))), dim="XY", type="POLYGON")
POLYGON ((0 0, 0 10, 10 0))

but in the WKT case, unlike rgeos::readWKT(), p2g leads to:

> st_is_valid(p2g)

 *** caught segfault ***
address (nil), cause 'memory not mapped'

Traceback:
 1: CPL_geos_is_valid(x, as.logical(NA_on_exception))
 2: st_is_valid.sfc(st_geometry(x), ...)
 3: st_is_valid(st_geometry(x), ...)
 4: st_is_valid.sfg(p2g)
 5: st_is_valid(p2g)

I don't think GEOS 3.10.0 should cause a segfault when 3.9.1 didn't, but something has changed in trapping the exception, and arguably CPL_sfc_from_wkt() has a different default on permitting unclosed ring polygons than those just created in R.

[rsbivand]

Interestingly:

> library(sf)
Linking to GEOS 3.10.0beta2, GDAL 3.3.2, PROJ 8.1.1
> p2 <- st_as_sfc("POLYGON((0 0, 0 10, 10 0))")
> st_is_valid(p2, NA_on_exception = FALSE)
Error in CPL_geos_is_valid(x, as.logical(NA_on_exception)) : 
  Evaluation error: IllegalArgumentException: Points of LinearRing do not form a closed linestring.

so is the error handler not working and the error not being trapped? Is HAVE350 being detected? Adding an Rprintf() in src/geos.cpp l. 542 in #ifdef HAVE350 #endif gives:

> library(sf)
Linking to GEOS 3.10.0beta2, GDAL 3.3.2, PROJ 8.1.1
> p2 <- st_as_sfc("POLYGON((0 0, 0 10, 10 0))")
> st_is_valid(p2, NA_on_exception = FALSE)
HAVE350 detected
Error in CPL_geos_is_valid(x, as.logical(NA_on_exception)) : 
  Evaluation error: IllegalArgumentException: Points of LinearRing do not form a closed linestring.

so yes, detected.

[rsbivand]

Commenting out the __countErrorHandler and __emptyNoticeHandler (in HAVE350) in sf/src/geos.cpp in CPL_geos_is_valid(), we get no segfault:

> library(sf)
Linking to GEOS 3.10.0beta2, GDAL 3.3.2, PROJ 8.1.1
> p2 <- st_as_sfc("POLYGON((0 0, 0 10, 10 0))")
> st_is_valid(p2, NA_on_exception = FALSE)
Error in CPL_geos_is_valid(x, as.logical(NA_on_exception)) : 
  Evaluation error: IllegalArgumentException: Points of LinearRing do not form a closed linestring.
> st_is_valid(p2, NA_on_exception = TRUE)
Error in CPL_geos_is_valid(x, as.logical(NA_on_exception)) : 
  Evaluation error: IllegalArgumentException: Points of LinearRing do not form a closed linestring.

but also no NA set. The GEOS error is thrown immediately GEOSisValid_r() is called.

[paleolimbot]

I wonder if the error handler is doing a longjmp and/or throwing an exception. I don't think the error handler is supposed to do either and I wonder if that's the problem (or part of it). FWIW, in 'geos' I just record the last error message and call error() after whatever was called returns an error code:

https://github.com/paleolimbot/geos/blob/master/src/geos-common.c#L5-L18

https://github.com/paleolimbot/geos/blob/master/src/geos-io.c#L25

[rsbivand]

@paleolimbot thanks for joining in. I cannot see where else to try things here, __countErrorHandler() is used elsewhere too, and there is nothing in NEWS to indicate changes to its arguments that would change the behaviour of st_is_valid() with an input polygon that should cause an exception. Will you try to see if Edzer is online or should I? If this is a GEOS regression (I'm unsure here), now is the time to raise this before the beta becomes an RC. Really it is unfortunate that the GDAL WKT reader permits an unclosed polygon to be created.

[edzer]

Thanks; building a docker image now, not sure if I get to this before tomorrow evening.

[paleolimbot]

It looks like the message handler gets called in isValid on exception:

https://github.com/libgeos/geos/blob/main/capi/geos_ts_c.cpp#L707

...and that the userdata that's passed is stack-allocated memory:

https://github.com/r-spatial/sf/blob/master/src/geos.cpp#L1142-L1144

Could you pass NULL as userdata since it's empty?

[rsbivand]

Right, but NOTICE_MESSAGE() is being called, not ERROR_MESSAGE(), so why is it segfaulting? I see a lot of recent commits to geos_ts_c.cpp. @pramsey, does it seem possible that previous behaviour in

sf/src/geos.cpp

Lines 546 to 576 in 30cccf5

	#ifdef HAVE350
	GEOSContext_setNoticeMessageHandler_r(hGEOSCtxt,
	(GEOSMessageHandler_r) __emptyNoticeHandler, (void *) &notice);
	GEOSContext_setErrorMessageHandler_r(hGEOSCtxt,
	(GEOSMessageHandler_r) __countErrorHandler, (void *) &notice);
	#endif
	}
	Rcpp::LogicalVector out(sfc.size());
	for (int i = 0; i < out.length(); i++) {
	// get geometry i:
	Rcpp::List geom_i = Rcpp::List::create(sfc[i]);
	geom_i.attr("precision") = sfc.attr("precision");
	geom_i.attr("class") = sfc.attr("class");
	geom_i.attr("crs") = sfc.attr("crs");
	SEXP classes = sfc.attr("classes");
	if (classes != R_NilValue) {
	Rcpp::CharacterVector cl = sfc.attr("classes");
	geom_i.attr("classes") = cl[i];
	}
	std::vector<GeomPtr> gmv = geometries_from_sfc(hGEOSCtxt, geom_i, NULL); // where notice might be set!
	int ret = GEOSisValid_r(hGEOSCtxt, gmv[0].get());
	if (NA_on_exception && (ret == 2 || notice != 0))
	out[i] = NA_LOGICAL; // no need to set notice back here, as we only consider 1 geometry #nocov
	else
	out[i] = chk_(ret);
	notice = 0; // reset notice.
	}
	#ifdef HAVE350
	GEOSContext_setNoticeHandler_r(hGEOSCtxt, __warningHandler);
	GEOSContext_setErrorHandler_r(hGEOSCtxt, __errorHandler);
	#endif

was wrong, and cleanups have revealed its problem, or have recent cleanups led to a possible regression with a hard exception being thrown where the intention was simply to signal invalidity? The case here is an unclosed polygon. I know that:

If your geometries are invalid, the problem is with your geometry, not with GEOS.

from https://trac.osgeo.org/geos/wiki/TopologyExceptions, but checking validity is not the same as a predicate or operation.

[pramsey]

I'm finding it hard to understand how you are calling things. Is it possible to make a standalone that dies? The was a big re-write in IsValidOp, so that's where things changed most likely.

[rsbivand]

We see that where before an unclosed polygon threw an exception (correctly) which we could trap, now the same code segfaults with memory not mapped.

We'll try to create a standalone, but two of us are CEST, so it will take time.

Might the code fail even before getting to validity checking, has the polygon creation code also changed?

@edzer the same difference is found with released sf as with master latest, which was what I started by testing.

[pramsey]

So, I made up a little unit test

    GEOSCoordSequence* shell_seq = GEOSCoordSeq_create(4, 2);
    double shell_coords[] = {0,0, 0,10, 10,10, 10,0};
    for (unsigned int i = 0; i < 4; i++) {
        GEOSCoordSeq_setXY(shell_seq, i, shell_coords[2*i], shell_coords[2*i+1]);
    }

    GEOSGeometry* shell = GEOSGeom_createLinearRing(shell_seq);
    ensure(shell != nullptr);
    GEOSGeometry* polygon = GEOSGeom_createPolygon(shell, nullptr, 0);
    ensure(polygon != nullptr);
    char isvalid = GEOSisValid(polygon);
    ensure_equals(0, isvalid);
    GEOSGeom_destroy(polygon);

I get IllegalArgumentException: Points of LinearRing do not form a closed linestring when GEOSGeom_createLinearRing() is called. If I fail to catch that null, then I end up passing a null into IsValid which segfaults, so I guess I'm wondering if you've got a null further down the line?

[pramsey]

Yeah, seems odd, you should be getting a notice and a null return when you attempt to create an unclosed LinearRing.

[paleolimbot]

It looks like there's no NULL check when creating the pointer from WKB:

https://github.com/r-spatial/sf/blob/master/src/geos.cpp#L168

There are a lot of other places in geos.cpp without NULL checks, too, although there is a function to check for a null geometry:

https://github.com/r-spatial/sf/blob/master/src/geos.cpp#L701-L706

I haven't tested this for memory leaks, but I think that sf is leaking the context object whenever this throws! My strategy in the 'geos' package is a global context object so that I don't have to worry about leaking that particular object.

[edzer]

It looks like @paleolimbot found the right spot; checking here for a NULL pointer returned gives me

> library(sf)
Linking to GEOS 3.10.0beta2, GDAL 3.3.2, PROJ 8.0.1
> p1 <- st_as_sfc("POLYGON((0 0, 0 10, 10 10, 10 0, 0 0))")
> p2 <- st_as_sfc("POLYGON((0 0, 0 10, 10 0))")
> set1 <- c(p1, p2, p1, p2, p1, p2)
> st_is_valid(set1)
Error in CPL_geos_is_valid(x, as.logical(NA_on_exception)) : 
  GEOS returned NULL pointer on invalid geometry
Calls: st_is_valid -> st_is_valid.sfc -> CPL_geos_is_valid

It seems that this is new behavior of GEOS, but makes complete sense IMO.

[rsbivand]

I have in the GEOSIsValid_r() call:

		int ret;
                if (gmv[0].get() == NULL) 
                    ret = 2;
                else
                    ret = GEOSisValid_r(hGEOSCtxt, gmv[0].get());

where both nullptr and NULL seem to work. Curiously, GEOSisValidReason_r() never had a problem, but it doesn't swap out the handlers, so perhaps some interaction?

[rsbivand]

I also see that test_valid.R was modified 24 August, so after 1.0-2 was released; 1.0-2 st_is_valid(p2) still segfaults with GEOS 3.10.0beta2, so that change isn't relevant.

[edzer]

Thanks everyone, very helpful - should work now!

[rsbivand]

Looks good here too, no segfault with GEOS 3.10.0beta2.

[pramsey]

So, an FYI. The fact that you could pass nullptr into GEOSIsValid() in version 3.9 was interesting so I had a look at that, and found that (a) you could pass null and (b) the return value would be 1 (aka valid) which was a little shocking. I think I see in your code that you ignore the return value if you have extant notices, which I guess was a workaround for something like this. In general the issue of unclosed rings is a little hard to deal with. PostGIS actually closes rings before converting into GEOS as a way of working around the GEOS limitation on not accepting unclosed rings in the LinearRing constructor. I'm not sure if you would find that more reliable.