GEOS-devel issues with OverlayNG

[rsbivand]

With GEOS devel version installed with --enable-overlayng, I see a failure in sf_intersection(). Here is the log and a zipfile of an rda file with objects sabre_x and sabre_y. Could you please check with released GEOS, then build sf and GDAL against GEOS built with --enable-overlayng, say from a recent nightly bz-ball, and report back? See also r-spatial/sf#794.

sabre_00check.log
GEOS_dev_issue.zip

Note that I'm also using GDAL 3.2.0 RC1, so the OGR: Corrupt data may be from GDAL rather than GEOS. I see a very similar error in the 'siland* package (not on github), logfile here:
siland_00check.log. If you get test containers up, debugging that error may help pin down the problem. If it is GDAL, feedback before 3.2.0 is released would help!

[rsbivand]

With:

> library(sf)
Linking to GEOS 3.9.0dev, GDAL 3.1.4, PROJ 7.2.0

siland and sabre pass, so it looks like a GDAL issue, perhaps? @edzer, will we be able to find out what happens before 3.2.0 is released?

[Nowosad]

@rsbivand Thank you for letting me know about the issue. I am now working on creating a testing docker image. Is is possible that the package data (sf and raster objects) need to be updated (due to PROJ6+ changes)?

[rsbivand]

Very curious, both now pass for sf 0.9-6 and GDAL 3.2.0 RC1:

> library(sf)
Linking to GEOS 3.9.0dev, GDAL 3.2.0, PROJ 7.2.0

As of now, the problem comes from running with the development version of sf. Installing sf from github seems to lead to failures in sabre and siland, but I'm re-checking.

[edzer]

siland and sabre pass, so it looks like a GDAL issue, perhaps? @edzer, will we be able to find out what happens before 3.2.0 is released?

Could you please elaborate the question?

[rsbivand]

Running revdep checks on 900+ packages with the nightly build of GEOS, GDAL 3.2.0 RC1, PROJ 7.2.0 RC1, and sf devel, I noticed failures in siland and sabre that were not there before, in sabre, example(mapcurves_calc) fails with:

> mc = mapcurves_calc(x = regions1, y = regions2, x_name = z, y_name = z)
OGR: Corrupt data
Error in CPL_geos_op2(op, x, y) : OGR error
Calls: mapcurves_calc -> mapcurves_calc.sf
Execution halted

In siland, example(landSiland) fails with:

> data(landSiland)
> plot(landSiland$geometry)
OGR: Corrupt data
Error in CPL_geos_is_empty(st_geometry(x)) : OGR error
Calls: plot ... plot.sfc_MULTIPOLYGON -> st_is_empty -> CPL_geos_is_empty
Execution halted

The failures are with remotes::install_github("r-spatial/sf"); with CRAN sf they pass, it turns out without any relationship to the version of GDAL. I've only checked GDAL 3.1.4 and 3.2.0RC1, and sf 0.9-6, 0.9-7 (current master), GEOS is 3.9.0dev --enable-overlayng (on one machine yesterday's nightly bz-ball, on the other git master of some days ago). So is there an interaction between sf and GEOS versions? I haven't varied GEOS, because I'd have to re-install GDAL.

[edzer]

Thanks; I will try to reproduce this, am now rebuilding the overlayNG docker image with Martin's / strk's recent updates.

[rsbivand]

Same errors with GEOS 3.8.1, GDAL 3.2.0 RC1, PROJ 7.2.0 RC1 and sf from github master, and with GEOS 3.9.0 devel without --enable-overlayng.

[edzer]

So GDAL RC is the smoking gun? Will build that against overlayNG.

[rsbivand]

No, the failures are only present with sf > 0.9-6, irrespective of GDAL 3.1.4 or 3.2.0 RC1, and GEOS 3.8.1, 3.9.0 with and without OverlayNG. I only tried the GEOS variants with GDAL 3.2.0 RC1. I don't see in sf/src/wkb.cpp where any GDAL code is involved, and I can't see how the GDAL error handler is being triggered.

[rsbivand]

I think @Nowosad found the problem - the extracts above are:

> st_crs(sabre_x)
Coordinate Reference System: NA
> st_crs(sabre_y)
Coordinate Reference System: NA

and r-spatial/sf@b8e4af4 is recent and changes how CRS are handled on export to GEOS, I think. And:

> st_crs(landSiland)
Coordinate Reference System: NA

Does the change take account of NA CRS? The sabre_[xy] objects also fail st_is_empty() with the same error.

[edzer]

the trouble in both cases were old-style crs that were NA.

[rsbivand]

attributes(st_geometry(sabre_x))

shows

$crs
Coordinate Reference System:
  No user input

So:

st_crs(sabre_x) <- NA
st_is_empty(sabre_x)
st_crs(sabre_y) <- NA
st_intersection(sabre_x, sabre_y)

work. As does:

> library(siland)
> data(landSiland)
> st_crs(landSiland) <- NA
> plot(st_geometry(landSiland))

@edzer could you tell the siland maintainer?

[edzer]

fix now in sf master.