for #2063; reproducing rgeos workflow for assigning interior to exterior rings

[rsbivand]

r-spatial/evolution#9, rstudio/leaflet#833

This is not optimized at all, but since it is unlikely to be entered unless someone has canned, stale sp polygon objects in their workflows, speed doesn't matter. Here the workhorse is: st_contains_properly(raw0[!holes], raw0[holes]), where holes are got from the input object. If no holes were assigned initially, the flow of control reverts effectively to the previous code. st_contains_properly() may fail if a quasi-hole shares a border segment with its containing exterior ring. Holes are assigned as interior rings to exterior rings in order of exterior ring area from smallest to largest (the nested island/holes case), and assigned holes the dropped from the larger exterior rings.

[edzer]

According to the simple feature standard, holes do not share a border segment with the exterior ring. Does this PR bring in non-SF logic into the package? What speaks against making such objects (SF-)valid at this point?

[rsbivand]

I'll check and report back - st_contains_properly() checks that all of B is in the interior of A, so perhaps using st_contains() would catch the case. Now I think that such a hole would warn that no exterior ring could be found, but if st_contains() was used, it would assign the interior ring to the correct exterior ring, and st_make_valid() used later should restore SF conformance.

With current a9f1b98 and adding:

  p9a <- Polygon(cbind(x=c(1, 1, 2, 2, 1), y=c(6, 7, 7, 6, 6)), hole=TRUE) # H
  p7a <- Polygon(cbind(x=c(14, 14, 15, 15, 14), y=c(13, 14, 14, 13, 13)), hole=TRUE) # H
  lp <- list(p1, p2, p13, p7, p7a, p6, p5, p4, p3, p8, p11, p12, p9, p9a, p10)
# 
plot(spls, col="white", bg="grey")

to tests/testthat/test_sp.R:

> spls_sfc <- sf::st_as_sfc(spls)
Warning message:
In FUN(X[[i]], ...) : orphaned hole, cannot find containing polygon
# 
plot(spls_rt, col="white", bg="grey")

as expected (holes sharing borders are merged, hole on edge of central exterior ring not found).

[rsbivand]

With st_contains() instead of st_contains_properly(), no warning of orphaned hole and the hole on the edge of the central island removed and the edge notched:

Will update the PR; in neither case would coercion to "sf" create an invalid object, but some geometry details would have been lost with st_contains_properly().

[edzer]

Great, so this is ready to go, right?

[rsbivand]

Yes, as far as I can see. Should resolve r-spatial/evolution#9, and rstudio/leaflet#833, where rgeos::createPolygonsComment() isn't exactly patched by as(sf::st_as_sfc(x), "Spatial") because x is a "Spatial*", so this round trip corresponds to rgeos::createSPComment(). I don't think it is necessary to create coercions from "Polygons" to "sfg" and back, but https://github.com/rstudio/leaflet/blob/51daa695419fe2bd9bab3d9fadf69d59bbf795ee/R/normalize-sp.R#L85-L109 might well be re-configured to coerce input "Spatial*" objects into "sf" objects then use normalize-sf.R, probably for all the vector objects.

[edzer]

Strange enough this passes locally, but fails on windows, see here. Is it OK if I outcomment the test for the 1.0-10 release?

[rsbivand]

Which sp version were you using? Winbuilder is on Rtools43, but I'm sure with released sp. My local testing was on Fedora with sp fixed for the main() problem - do we also need a release https://cran.r-project.org/web/checks/check_results_roger.bivand_at_nhh.no.html to remove the rgdal and rgeos warnings. I'll submit rgdal shortly with updated src/Makevars.ucrt for the Rtools43 error.

[rsbivand]

Did you also update to: https://svn.r-project.org/R-dev-web/trunk/WindowsBuilds/winutf8/ucrt3/r_packages/patches/CRAN/sf.diff for GDAL 3.6.2 ? I think Win-builder has an intermediate version with older GEOS - maybe a cause? Trying now ...

[rsbivand]

With R Under development (unstable) (2023-01-03 r83550 ucrt), matching Rtools with GDAL 3.6.2 and GEOS 3.11.1:

00check.log
00install.zip
testthat.zip

[rsbivand]

With released CRAN sf and GEOS 3.9.3, I see (in debug and after pasting in the new process_pl_comment()):

>   spls_rt <- as(spls_sfc, "Spatial")
> comment(slot(spls_rt, "polygons")[[1]])
[1] "0 1 1 1"
> sf_extSoftVersion()
          GEOS           GDAL         proj.4 GDAL_with_GEOS     USE_PROJ_H 
       "3.9.3"        "3.5.2"        "8.2.1"         "true"         "true" 
          PROJ 
       "8.2.1" 

so maybe a dependency on GEOS versions later than 3.9.3? Shall I try to establish which? Should we condition this part of test_sp.R on say GEOS >= 3.11? For example:

skip_if_not(package_version(sf_extSoftVersion()["GEOS"]) >= "3.11.0")

[rsbivand]

Confirmed that with GEOS 3.9.4 on Fedora, st_make_valid() returns two GEOMETRYCOLLECTION features and four POLYGON features, with 3-.11.1, six POLYGON features.