r-spatial / sf

Simple Features for R
https://r-spatial.github.io/sf/
Other
1.35k stars 299 forks source link

st_collection_extract() removes everything if an empty geometry collection is present #2367

Open dschlaep opened 8 months ago

dschlaep commented 8 months ago

Describe the bug The presence of an empty geometrycollection causes st_collection_extract()to return an empty object even if geometries of the requested type are present.

For geometrycollections, st_cast_sfc_default() checks for all(lengths(x)) https://github.com/r-spatial/sf/blame/7bf3b66cec65ef21ee580926cb69444af7f3a931/R/cast_sfc.R#L125 which is false if there is an empty geometrycollection. Thus, the object is not unlisted and everything gets removed because nothing appears to be correctly identified.

This appears to be related to issue https://github.com/r-spatial/sf/issues/1767 and commit https://github.com/r-spatial/sf/commit/92376396eb5f3f2f82f43f14d1579856d8eb5db0

Thanks!

To Reproduce

library("sf")
#> Linking to GEOS 3.12.1, GDAL 3.8.4, PROJ 9.3.1; sf_use_s2() is TRUE
packageVersion("sf")
#> [1] '1.0.15'

# Examples based on ?st_geometrycollection
pt <- st_point(c(1, 0))
ls <- st_linestring(matrix(c(4, 3, 0, 0), ncol = 2))
poly1 <- st_polygon(list(matrix(c(5.5, 7, 7, 6, 5.5, 0, 0, -0.5, -0.5, 0), ncol = 2)))
poly2 <- st_polygon(list(matrix(c(6.6, 8, 8, 7, 6.6, 1, 1, 1.5, 1.5, 1), ncol = 2)))

i <- st_geometrycollection(list(pt, ls, poly1, poly2))
(f <- st_as_sf(st_sfc(i, st_geometrycollection(list(pt)))))
#> Simple feature collection with 2 features and 0 fields
#> Geometry type: GEOMETRYCOLLECTION
#> Dimension:     XY
#> Bounding box:  xmin: 1 ymin: -0.5 xmax: 8 ymax: 1.5
#> CRS:           NA
#>                                x
#> 1 GEOMETRYCOLLECTION (POINT (...
#> 2 GEOMETRYCOLLECTION (POINT (...

# st_collection_extract() works as expected, i.e., returns polygons
st_collection_extract(f, "POLYGON")
#> Simple feature collection with 2 features and 0 fields
#> Geometry type: POLYGON
#> Dimension:     XY
#> Bounding box:  xmin: 5.5 ymin: -0.5 xmax: 8 ymax: 1.5
#> CRS:           NA
#>                                x
#> 1 POLYGON ((5.5 0, 7 0, 7 -0....
#> 2 POLYGON ((6.6 1, 8 1, 8 1.5...

# Issue: object with empty geometrycollection
(f2 <- st_as_sf(st_sfc(i, st_geometrycollection())))
#> Simple feature collection with 2 features and 0 fields (with 1 geometry empty)
#> Geometry type: GEOMETRYCOLLECTION
#> Dimension:     XY
#> Bounding box:  xmin: 1 ymin: -0.5 xmax: 8 ymax: 1.5
#> CRS:           NA
#>                                x
#> 1 GEOMETRYCOLLECTION (POINT (...
#> 2       GEOMETRYCOLLECTION EMPTY

# st_collection_extract() returns an empty object
# -- expected behavior: returns polygons (as with extraction of `f`)
st_collection_extract(f2, "POLYGON")
#> Warning in st_collection_extract.sf(f2, "POLYGON"): x contains no geometries of
#> specified type
#> Simple feature collection with 0 features and 0 fields
#> Bounding box:  xmin: NA ymin: NA xmax: NA ymax: NA
#> CRS:           NA
#> [1] x
#> <0 rows> (or 0-length row.names)

Created on 2024-03-22 with reprex v2.1.0

Additional context

sf::sf_extSoftVersion()
#>           GEOS           GDAL         proj.4 GDAL_with_GEOS     USE_PROJ_H 
#>       "3.12.1"        "3.8.4"        "9.3.1"         "true"         "true" 
#>           PROJ 
#>        "9.3.1"

Created on 2024-03-22 with reprex v2.1.0