Sf: Regression in revdep in st_sample when exact = TRUE

[Robinlovelace]

See #982

[Robinlovelace]

Reprex:

    remotes::install_github("robinlovelace/sf", "update-st_sample_exact2")
#> Skipping install of 'sf' from a github remote, the SHA1 (8c452d20) has not changed since last install.
#>   Use `force = TRUE` to force installation
    mtq = sf::st_read(system.file("gpkg/mtq.gpkg", package="cartography"))
#> Reading layer `mtq' from data source `/home/robin/R/x86_64-pc-linux-gnu-library/3.5/cartography/gpkg/mtq.gpkg' using driver `GPKG'
#> Simple feature collection with 34 features and 7 fields
#> geometry type:  MULTIPOLYGON
#> dimension:      XY
#> bbox:           xmin: 690574 ymin: 1592536 xmax: 735940.2 ymax: 1645660
#> epsg (SRID):    32620
#> proj4string:    +proj=utm +zone=20 +datum=WGS84 +units=m +no_defs
    mtq_pencil = cartography::getPencilLayer(x = mtq)
#> Error in st_coordinates.sfc(pt): not implemented for objects of class sfc_GEOMETRY

^{Created on 2019-02-23 by the reprex package (v0.2.1)}

[Robinlovelace]

Trying to reproduce but struggling to debug it. The reprex below shows preliminary attempts. Any ideas what is going on here @riatelab? Will plan to put the findings in a minimal test to ensure this does not happen again.

library(sf)
#> Linking to GEOS 3.7.0, GDAL 2.3.2, PROJ 4.9.3
mtq = sf::st_read(system.file("gpkg/mtq.gpkg", package="cartography"))
#> Reading layer `mtq' from data source `/home/robin/R/x86_64-pc-linux-gnu-library/3.5/cartography/gpkg/mtq.gpkg' using driver `GPKG'
#> Simple feature collection with 34 features and 7 fields
#> geometry type:  MULTIPOLYGON
#> dimension:      XY
#> bbox:           xmin: 690574 ymin: 1592536 xmax: 735940.2 ymax: 1645660
#> epsg (SRID):    32620
#> proj4string:    +proj=utm +zone=20 +datum=WGS84 +units=m +no_defs
mtq_pencil = cartography::getPencilLayer(x = mtq) # reproduces error
#> Error in st_coordinates.sfc(pt): not implemented for objects of class sfc_GEOMETRY
st_area(mtq)
#> Units: [m^2]
#>  [1] 699261 709840 697602 702229 698805 709840 722387 706811 706345 726831
#> [11] 697312 706811 716537 710977 697312 706811 734557 699261 697312 725227
#> [21] 723118 727316 722387 706345 692220 734605 725227 709622 702757 722913
#> [31] 715208 726610 702229 698805
x = mtq[1, ]
size = 3
a = 10000
buffer = 1
size <- round(sqrt(as.numeric(sf::st_area(x) * size/a)),
              0)
if (size <= 10) {
    size <- 10
}
pt <- sf::st_sample(sf::st_buffer(x, buffer), size = size)
class(pt)
#> [1] "sfc_POINT" "sfc"

^{Created on 2019-02-23 by the reprex package (v0.2.1)}

A separate issue is that I think there should be more tests on st_sample(). https://github.com/r-spatial/sf/commit/b7f7646fa2d3f9e392b83bf21aab33d05baa8c3c was just a starter.

[Robinlovelace]

Update: been looking at cartograpy source code. Trying to see under what conditions st_sample(..., exact = TRUE) is returning an object of sfc_GEOMETRY.

Great cartographic plotting package in any case. Heads-up @rcarto - any ideas/suggestions on this greatly appreciated. Thanks!

[rsbivand]

Might setting a seed help in seeing whether some outcomes fail?

[edzer]

Here is a clue:

> st_polygon(list(rbind(c(0,0), c(1,0), c(1,1), c(0,0))))
POLYGON ((0 0, 1 0, 1 1, 0 0))
> p = st_sfc(st_polygon(list(rbind(c(0,0), c(1,0), c(1,1), c(0,0)))))
> st_sample(p, 1)
Geometry set for 0 features 
bbox:           xmin: NA ymin: NA xmax: NA ymax: NA
epsg (SRID):    NA
proj4string:    NA
> x = st_sample(p, 1)
> x
Geometry set for 1 feature 
geometry type:  POINT
dimension:      XY
bbox:           xmin: 0.7386794 ymin: 0.5120839 xmax: 0.7386794 ymax: 0.5120839
epsg (SRID):    NA
proj4string:    NA
POINT (0.7386794 0.5120839)
> y = st_sample(p, 1)
> y
Geometry set for 0 features 
bbox:           xmin: NA ymin: NA xmax: NA ymax: NA
epsg (SRID):    NA
proj4string:    NA
> c(x, y)
Geometry set for 1 feature 
geometry type:  GEOMETRY
dimension:      XY
bbox:           xmin: 0.7386794 ymin: 0.5120839 xmax: 0.7386794 ymax: 0.5120839
epsg (SRID):    NA
proj4string:    NA
POINT (0.7386794 0.5120839)

[edzer]

and also

> c(y, x)
Geometry set for 1 feature 
geometry type:  GEOMETRY
dimension:      XY
bbox:           xmin: 0.7386794 ymin: 0.5120839 xmax: 0.7386794 ymax: 0.5120839
epsg (SRID):    NA
proj4string:    NA
POINT (0.7386794 0.5120839)

[rCarto]

I'll have a look at this issue on Monday.
About setting a seed, it would be great for reproducibility anyway. I remember reading somewhere it was not possible with 'sp'.

[edzer]

@rCarto - thanks, but the problem is clearly with sf, not cartography.

[Robinlovelace]

Aha, it happens when st_sample() returns nothing and concatenates the result with a valid point. I can imagine 2 solutions:

1. Preventing st_sample() from returning an empty result.

2. Make concatenating a point with an empty geometry return a point.

Plan to try strategy 1 on this branch: https://github.com/r-spatial/sf/compare/master...Robinlovelace:update-st_sample_exact2

Many thanks for the hints @edzer, interesting and hopefully usefully explorations of this thorny (or should I say pointy ; ) issue.

[edzer]

I vote for 2.

[Robinlovelace]

Which would involve changes to c.sfc() here I imagine: https://github.com/r-spatial/sf/blob/master/R/sfc.R#L139

Suspect it would involve changing this line but not sure:

https://github.com/r-spatial/sf/blob/a07168bdb316f4cc1bf5a1a959f204d473f8cdec/R/sfc.R#L146-L147

[rsbivand]

No, I think you have to leave c.sfc() as SF-compliant. IIRC we banned empty geometries in sp and later in rgeos as totally meaningless, but they are in the standard. I'd rather test before c() in st_sample_exact(), line 213 in sf/R/sample.R:

https://github.com/r-spatial/sf/blob/a07168bdb316f4cc1bf5a1a959f204d473f8cdec/R/sample.R#L213

Test that randon_pt_new isn't empty and catenate only if so. Here empty features are not needed, but they may be needed in c.

[rsbivand]

Looks as though:

if (length(random_pt_new) > 0L) random_pt = c(random_pt, random_pt_new)

should work, but this would need testing with much harder cases (use the off-shore sand bars in NC, for example), perhaps because I'm unsure how recursion will play out.

[edzer]

I'm entirely with you, and I have nothing against solving it at this level. Another problem however, which is the reason that this has to be solved at this level in the first place, is that

> c(st_sfc(), st_sfc(st_point(c(0,1))))
Geometry set for 1 feature 
geometry type:  GEOMETRY
dimension:      XY
bbox:           xmin: 0 ymin: 1 xmax: 0 ymax: 1
epsg (SRID):    NA
proj4string:    NA
POINT (0 1)

This has nothing to do with SFA compliancy, since SFA doesn't talk about the class of _sets of_ simple feature geometries, which class they should have when they're empty (= set of zero features, _not_ a feature without geometry = empty geometry).

I think it makes sense that an empty set has class GEOMETRY, but that this class once combined with something else should get the class of the something else, since the empty set essentially disappears right there.

Solving that would make solving it at this level obsolete.