I've been experimenting with data.table + geos (see https://grantmcdermott.com/fast-geospatial-datatable-geos/ for the example I followed), but I've had issues converting the data.table+geos object back to sf (for mapping with tmap, for example).
In the reprex below, I take an sf object and convert it to a data.table with a geos column. I combine units on some criteria, and then map the result. This operation works fine with the nc data set (first section), but fails using the Digital Tax Map of NYC (original data, subset used in reprex). I am also unable to make the converted-back-to-sf object valid, or to recast it to MULTIPOLYGON. The error message is always Error: Not compatible with STRSXP: [type=NULL]..
As mentioned, I'm new to geos and mostly experimenting at this stage, so I might be missing something obvious, or doing something very unorthodox — please let me know if that's the case, how I could improve and how to make this work.
library(data.table)
library(geos)
library(sf)
#> Linking to GEOS 3.11.0, GDAL 3.5.3, PROJ 9.1.0; sf_use_s2() is TRUE
library(tmap)
#> The legacy packages maptools, rgdal, and rgeos, underpinning the sp package,
#> which was just loaded, will retire in October 2023.
#> Please refer to R-spatial evolution reports for details, especially
#> https://r-spatial.org/r/2023/05/15/evolution4.html.
#> It may be desirable to make the sf package available;
#> package maintainers should consider adding sf to Suggests:.
#> The sp package is now running under evolution status 2
#> (status 2 uses the sf package in place of rgdal)
# NC: maps fine ----
nc = st_read(system.file("shape/nc.shp", package = "sf"))
#> Reading layer `nc' from data source
#> `/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/sf/shape/nc.shp'
#> using driver `ESRI Shapefile'
#> Simple feature collection with 100 features and 14 fields
#> Geometry type: MULTIPOLYGON
#> Dimension: XY
#> Bounding box: xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
#> Geodetic CRS: NAD27
nc_dt = as.data.table(nc)[, geo := as_geos_geometry(geometry)][, geometry := NULL]
nc_div =
nc_dt[,
.(geo = geo |> geos_make_collection() |> geos_unary_union()),
by = .(region = ifelse(CNTY_ID<=1980, 'high', 'low'))]
tmap_mode("plot")
#> tmap mode set to plotting
tm_shape(st_as_sf(nc_div)) + tm_polygons()
tmap_mode("view")
#> tmap mode set to interactive viewing
tm_shape(st_as_sf(nc_div)) + tm_polygons()
# Sample of NYC Digital Tax Map: doesn't map ----
dtm_sf = readRDS("dtm_sf.Rds")
dtm_dt = as.data.table(dtm_sf)[, geo := as_geos_geometry(geometry)][, geometry := NULL]
dtm_consolidate =
dtm_dt[,
.(geo = geo |> geos_make_collection() |> geos_unary_union()),
by = .(boro, block, lot)]
dtm_consolidate_sf = st_as_sf(dtm_consolidate)
plot(dtm_consolidate[, geo])
tmap_mode("plot")
#> tmap mode set to plotting
tm_shape(dtm_consolidate_sf) + tm_polygons()
tmap_mode("view")
#> tmap mode set to interactive viewing
tm_shape(dtm_consolidate_sf) + tm_polygons()
#> Error: Not compatible with STRSXP: [type=NULL].
class(st_geometry(dtm_consolidate_sf))
#> [1] "sfc_GEOMETRY" "sfc"
table(st_geometry_type(dtm_consolidate_sf))
#>
#> GEOMETRY POINT LINESTRING POLYGON
#> 0 0 0 30
#> MULTIPOINT MULTILINESTRING MULTIPOLYGON GEOMETRYCOLLECTION
#> 0 0 36 0
#> CIRCULARSTRING COMPOUNDCURVE CURVEPOLYGON MULTICURVE
#> 0 0 0 0
#> MULTISURFACE CURVE SURFACE POLYHEDRALSURFACE
#> 0 0 0 0
#> TIN TRIANGLE
#> 0 0
dtm_consolidate_sf_mltp = dtm_consolidate_sf |> st_cast("MULTIPOLYGON")
#> Error in eval(expr, envir, enclos): Not compatible with STRSXP: [type=NULL].
Created on 2023-07-14 with reprex v2.0.2