problem with st_intersection in windows (between lines & polygon / multipolygon)

[LHMarshall]

Maybe I have done something funny with my lines or multi-polygon but I am getting a strange behaviour on my windows machine (this particular dataset finds no intersections in windows while for other datasets it finds some but not as expected - I am generating grids of line segments with a random start point). It all works fine on a mac.

Behaviour on mac (correctly identifies segments inside multipolygon):

Behaviour on windows (nothing):

Code to reproduce (apologies for clunky plotting!):

Links to the object files are rot_strata.robj and samps_list.robj

# load list of LINESTRING sfg objects
load(file = "samps_list.robj")
# load a MULTIPOLYGON sfg object
load(file = "rot_strata.robj")

# To allow plotting of list of LINESTRINGS via lapply
line.plot <- function(tmp.line, col, lwd){
  mat <- as.matrix(tmp.line)
  lines(mat, col = col, lwd = lwd)
}

# Create empty plot
plot(c(0,0), xlim = c(-2.5, 0), ylim = c(0, 3.8))
# Plot all the line segments
tmp <- lapply(samplers, line.plot, col = 2, lwd = 2)

# Find parts of line segments which intersect with rot.strata
to.keep <- lapply(samplers, sf::st_intersection, y = rot.strata)
# Plot intersecting line segments
tmp <- lapply(to.keep, line.plot, col = 3, lwd = 4)

# Plot rot.strata polygons over the top
tmp <- as.matrix(rot.strata[[1]])
lines(tmp, col = 1, lwd = 1)
tmp <- as.matrix(rot.strata[[2]])
lines(tmp, col = 1, lwd = 1)

[rsbivand]

Please report versions of sf on the two machines; for me succeeds (same as your Mac) with development versions of sf and the external software libraries it depends on :

> packageVersion("sf")
[1] '1.0.3'
> sf_extSoftVersion()
          GEOS           GDAL         proj.4 GDAL_with_GEOS     USE_PROJ_H 
   "3.10.0dev"        "3.3.2"        "8.1.1"         "true"         "true" 
          PROJ 
       "8.1.1" 

For sf >= 1.0.*, s2 is used for st_intersection() for geographical coordinates, but as no "crs" is declared here, that should be irrelevant.

[LHMarshall]

Thanks for your reply. I was using the latest CRAN version of sf, I see some of my dependencies could do with updating too though.

Windows machine:

> packageVersion("sf")
[1] ‘1.0.2’
> sf::sf_extSoftVersion()
          GEOS           GDAL         proj.4 GDAL_with_GEOS     USE_PROJ_H           PROJ 
       "3.9.0"        "3.2.1"        "7.2.1"         "true"         "true"        "7.2.1" 

Mac machine:

> packageVersion("sf")
[1] ‘1.0.2’
> sf::sf_extSoftVersion()
          GEOS           GDAL         proj.4 GDAL_with_GEOS     USE_PROJ_H           PROJ 
       "3.8.1"        "3.2.1"        "7.2.1"         "true"         "true"        "7.2.1" 

I have tried updating the windows machine to the development version of sf 1.0.3 and the issue persists. I don't remember manually installing the dependencies, is there an easy way to update them? When I look it's talking about installing conda to install more recent versions.

I develop my package on a mac so I only found this issue when I was setting up github workflows and the windows machine on there flagged it as a failure. Output from that can be found here. You'll see that the mac and ubuntu checks passed on github workflows too, it's just the windows setup.

[rsbivand]

The difference is not likely to be in sf. There is a difference in GEOS versions. Please try to avoid manually updating the external dependencies, things should work in the same way for the same version of sf; alternative environments and builds of external software just add noise. On Windows with the CRAN sf binary, I can reproduce your wrong output. On Mac (arm64, M1), with GEOS 3.9.1, the correct output is reproduced. @edzer - any suggestions? Is this a candidate for an additional test?

[edzer]

My computer has 3.9.0 installed, and I can reproduce your wrong figure. However,

plot(do.call(st_sfc, samplers)[rot.strata], add = TRUE, lwd = 4, col = 3)

creates the right plot. Odd!

[LHMarshall]

Hmm it is strange, thanks for your input. I will keep investigating. It is strange that this is not an issue with any of my other designs (continuous parallel lines or zigzag lines) just the segments and my code is very similar across the functions. The do.call statement doesn't quite do what I'm after as I do need any segments that are part in and part out to be clipped to the polygon boundary.

This is an example of what happens in a different study region, there are some samplers and you can see where they should be to complete the grid across the region but lots are dropped that should be retained.

The subsetting call retains all that should be retained but doesn't clip them.

meanwhile st_intersection retains some, infact all that have been clipped but misses some within the centre of the polygon. Note it does seem to do slightly better on non-rectangular regions.

This is an example of a rectangular region. I have run it multiple times and it only retains the segments which intersect the polygon boundary. I guess that explains why my other designs work though as all the transects in them always intersect the polygon boundaries.

library(dssd)  
region <- make.region()
design <- make.design(region = region, 
                      transect.type = "line",
                      design = "segmentedgrid",
                      spacing = 100, 
                      seg.length = 50,
                      seg.threshold = 0,
                      edge.protocol = "minus",
                      truncation = 50)
transects <- generate.transects(design)
plot(region, transects, lwd = 0.8, col = "blue")

As I say I will keep investigating. Just thought I'd post incase it is of use.

[LHMarshall]

It did occur to me that I ran simulations with this segmented grid design on this same windows machine that I am using now at the end of Feb / beginning of March. I have gone back and tried to re-run them and they now fail. After convincing myself that I haven't changed my code recently (Nov 2019 were my last edits to this design) I tried rolling back the sf package, the results are below:

Working as expected on windows with the following versions:

> packageVersion("sf")
[1] ‘0.9.6’
> sf::sf_extSoftVersion()
          GEOS           GDAL         proj.4 GDAL_with_GEOS     USE_PROJ_H 
       "3.8.0"        "3.0.4"        "6.3.1"         "true"         "true"

> packageVersion("sf")
[1] ‘0.9.7’
> sf::sf_extSoftVersion()
          GEOS           GDAL         proj.4 GDAL_with_GEOS     USE_PROJ_H 
       "3.8.0"        "3.0.4"        "6.3.1"         "true"         "true" 

The issue appears to first arise with the following installed:

> packageVersion("sf")
[1] ‘0.9.8’
> sf::sf_extSoftVersion()
          GEOS           GDAL         proj.4 GDAL_with_GEOS     USE_PROJ_H 
       "3.9.0"        "3.2.1"        "7.2.1"         "true"         "true"

[edzer]

Perhaps this is a regression in GEOS 3.9.0, corrected in 3.9.1? Quite a few things happened with the release of 3.9.0, IIRC.

[dr-jts]

There is a bug in GEOS 3.9.0 which causes overlay to essentially drop straight lines if they are axis-parallel. It very much sounds like this is what you are seeing. The bug was fixed in GEOS 3.9.1.

If it's not possible to upgrade, two workarounds are:

• rotate the lines a very small amount
• use intersects to detect if a line is wholly inside a polygon, and if so just include it in the result directly

[LHMarshall]

@dr-jts thanks! That sounds like the issue. I rotate the polygons and do all the transect creation parallel to the y-axis, so in many cases I can just move the step of finding the intersection back to after I have rotated back again.

@edzer I can implement a workaround but I was just wondering if there was a planned release of sf which would install GEOS 3.9.1? Thanks for all your help.

[dr-jts]

@LHMarshall good that you have a workflow that can avoid this issue.

[edzer]

@LHMarshall that depends not only on the release of sf but also on the update of https://github.com/rwinlib/gdal3 and potentially on the transition of R-windows to using the UCRT / MXE toolchain, and the geos version used there.

[rsbivand]

@dr-jts thanks very much for joining in, very helpful. @LHMarshall, you could adapt your workflow to rotate slightly if GEOS == 3.9.0 but not otherwise. The rwinlib GEOS version needs to be updated @jeroen could you please replace GEOS 3.9.0 with GEOS 3.9.1 in rwinlib/gdal3? Reason is locationtech/jts#685 fixing libgeos/geos#408, which was a regression from 3.8.* and earlier that occurred in moving to the new and better overlay engine.

[jeroen]

I have updated the rwinlib/gdal3, and also ported to the v3.2.1 tag. So new builds of sf, rgeos, etc will automatically use geos 3.9.1 now.

The CRAN winbuilder rebuilds packages approximately once per week, so it may tak a few days before the windows binaries on CRAN are fixed.

[edzer]

Thanks a lot, @jeroen !

[LHMarshall]

Thank you everyone for your help!
Looks like this issue is now resolved.
I just updated sf from CRAN on my windows machine and my GEOS version is now 3.9.1

sf::sf_extSoftVersion()
          GEOS           GDAL         proj.4 GDAL_with_GEOS     USE_PROJ_H           PROJ 
       "3.9.1"        "3.2.1"        "7.2.1"         "true"         "true"        "7.2.1" 

I ran one of the examples above (I hadn't had chance to implement a workaround yet) and all looks to be working as expected now.

library(dssd)  
region <- make.region()
design <- make.design(region = region, 
                      transect.type = "line",
                      design = "segmentedgrid",
                      spacing = 100, 
                      seg.length = 50,
                      seg.threshold = 0,
                      edge.protocol = "minus",
                      truncation = 50)
transects <- generate.transects(design)
plot(region, transects, lwd = 0.8, col = "blue")