sf::st_as_sf() on a stars object crashes R

[mikoontz]

Hello,

I'm consistently crashing R when I try to use sf::st_as_sf() on a stars object when using the merge = TRUE argument (and not crashing with the default merge = FALSE). Is this something weird about what versions of GDAL, GEOS, and PROJ I have linked?

Here's a reproducible example:

library(stars)
library(sf)

df <- data.frame(x = c(295, 300, 305, 310, 295, 300, 305, 310),
                 y = c(310, 310, 310, 310, 315, 315, 315, 315),
                 d = c(2, 2, 2, 1, 2, 1, 1, 1))

outline <-
  df %>% 
  stars::st_as_stars() 

outline <-
  outline %>% 
  sf::st_as_sf(merge = TRUE)

And my sessionInfo:

> sessionInfo()
R version 3.6.3 (2020-02-29)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS Catalina 10.15.4

Matrix products: default
BLAS:   /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] stars_0.4-1 sf_0.9-1    abind_1.4-5

loaded via a namespace (and not attached):
 [1] compiler_3.6.3     magrittr_1.5       class_7.3-16       parallel_3.6.3    
 [5] DBI_1.1.0          tools_3.6.3        units_0.6-6        Rcpp_1.0.4.6      
 [9] KernSmooth_2.23-16 grid_3.6.3         e1071_1.7-3        lwgeom_0.2-1      
[13] classInt_0.4-3     rlang_0.4.5 

And my sf_extSoftVersion():

sf::sf_extSoftVersion()
          GEOS           GDAL         proj.4 GDAL_with_GEOS     USE_PROJ_H 
       "3.8.0"        "3.0.4"        "6.3.0"         "true"         "true" 

[mikoontz]

Note that using the merge = TRUE argument appears to work when converting from a raster as in this SO answer of mine (that I just updated to specify merge = TRUE in the call to sf::st_as_sf(): https://gis.stackexchange.com/a/313550/114075

Edit: But not for every raster. If I convert the object df above to a raster, R still crashes upon trying to use merge = TRUE:

## Convert your data.frame to a raster object
r <- raster::rasterFromXYZ(df)

outline <-
  r %>% 
  stars::st_as_stars() %>% 
  sf::st_as_sf(merge = TRUE)

The difference appears to be whether the raster has a CRS property associated with it. I can get R to crash if I take the working code from the SO answer above and strip away the CRS:

library(raster)
library(stars)
library(sf)
library(magrittr)

f <- system.file("external/test.grd", package="raster")
r <- raster(f)
r[r[] < 750] <- 0
r[r[] >= 750] <- 1

crs(r) <- NA

x <- st_as_stars(r) %>% 
  st_as_sf(merge = TRUE) %>% # this is the raster to polygons part
  st_cast("MULTILINESTRING") 

[edzer]

Still not fixed ... I'm pretty clueless!

[Athanasiamo]

I'm experiencing the same.

Perhaps not an ideal solution, but for my purposes I have no "real" crs, so I just applied an arbitrary one and then it works.

## Convert your data.frame to a raster object
r <- raster::rasterFromXYZ(df)
rsf  <– stars::st_as_stars(r)
sf::st_crs(rsf) <- 4326 

outline <- sf::st_as_sf(rsf, merge = TRUE)

[rsbivand]

Please NEVER use EPSG:4326 as an "arbitrary" CRS, it is firstly a GEOGCRS, and secondly swaps axes. @edzer - I think the wkt is a character NA, and is not checked for in sf/src/polygonize.cpp about line 80. It doesn't use the defensive behaviour of OGRsrs_from_crs() in sf/src/gdal.cpp at line 153. I haven't tried to fix it because my RCcpp is not active, but it looks like a plausible cause.

[edzer]

Very helpful - @rsbivand it is indeed NA, but also when it is st_crs(NA) it still broke. The NA CRS was then indeed not handled in sf/src/polygonize around line 80.