compatibility with sf library

I wonder if it would be possible to make data.table compatible with the new sf library. The library sf is promising to be a game changer for spatial analysis in R so it sounds like good idea to bring together the power of both libraries.

Currently, the class of an sf object is "sf" "data.frame" and it brings a column named geometry of class "sfc_MULTIPOLYGON" "sfc" .

Right now, I think the main incompatibility between dt and sf is that when we convert an sf data.frame into a data.table using setDT(), the geometry column is ruined and the object is not recognised anymore as an sf class.

I'm sure there are other points to take into account when making these two great libraries together so I just wanted to put the ball rolling if someone has not done this before.

Reproducible example please.

#1310 is related

Reproducible Example

library(data.table)
library(sf)

# load sf data
  nc <- st_read( system.file("shape/nc.shp", package="sf") )
  
  class(nc)
  > [1] "sf"         "data.frame"


head(nc)
  > Simple feature collection with 6 features and 14 fields
  > geometry type:  MULTIPOLYGON
  > dimension:      XY
  > bbox:           xmin: -81.74107 ymin: 36.07282 xmax: -75.77316 ymax: 36.58965
  > epsg (SRID):    4267
  > proj4string:    +proj=longlat +datum=NAD27 +no_defs
  >    AREA PERIMETER CNTY_ CNTY_ID        NAME  FIPS FIPSNO CRESS_ID BIR74 SID74 NWBIR74 BIR79 SID79 NWBIR79                        geometry
  > 1 0.114     1.442  1825    1825        Ashe 37009  37009        5  1091     1      10  1364     0      19 1 MULTIPOLYGON(((-81.47275543...
  > 2 0.061     1.231  1827    1827   Alleghany 37005  37005        3   487     0      10   542     3      12 2 MULTIPOLYGON(((-81.23989105...
  > 3 0.143     1.630  1828    1828       Surry 37171  37171       86  3188     5     208  3616     6     260 3 MULTIPOLYGON(((-80.45634460...
  > 4 0.070     2.968  1831    1831   Currituck 37053  37053       27   508     1     123   830     2     145 4 MULTIPOLYGON(((-76.00897216...
  > 5 0.153     2.206  1832    1832 Northampton 37131  37131       66  1421     9    1066  1606     3    1197 5 MULTIPOLYGON(((-77.21766662...
  > 6 0.097     1.670  1833    1833    Hertford 37091  37091       46  1452     7     954  1838     5    1237 6 MULTIPOLYGON(((-76.74506378...

Now when we try setting the sf object into data.table, note that the object looses both its "sf" class and the information in the geometry column

# set sf into data.table
  setDT(nc)
  
  class(nc)
  > [1] "data.table" "data.frame"
  
  head(nc)
  >     AREA PERIMETER CNTY_ CNTY_ID        NAME  FIPS FIPSNO CRESS_ID BIR74 SID74 NWBIR74 BIR79 SID79 NWBIR79 geometry
  > 1: 0.114     1.442  1825    1825        Ashe 37009  37009        5  1091     1      10  1364     0      19     <XY>
  > 2: 0.061     1.231  1827    1827   Alleghany 37005  37005        3   487     0      10   542     3      12     <XY>
  > 3: 0.143     1.630  1828    1828       Surry 37171  37171       86  3188     5     208  3616     6     260     <XY>
  > 4: 0.070     2.968  1831    1831   Currituck 37053  37053       27   508     1     123   830     2     145     <XY>
  > 5: 0.153     2.206  1832    1832 Northampton 37131  37131       66  1421     9    1066  1606     3    1197     <XY>
  > 6: 0.097     1.670  1833    1833    Hertford 37091  37091       46  1452     7     954  1838     5    1237     <XY>
  ```

For the record, I've also created an issue in the sf github page #428 because the compatibility between the two libraries might request a bit of collaboration from both sides.

If you lose the sf class you can still rely on the sfc class, but you have to take care of more stuff. Also, returning a geometry in a data.table seems tricky, but my DT skills are rusty, so there might be a better way. There is no reason it shouldn't work in theory, but in practice some things might have to be ironed. Looking at the tidy verbs implementation should give a good idea of what needs to be reconciled.

library(sf)
#> Linking to GEOS 3.5.1, GDAL 2.1.3, proj.4 4.9.2, lwgeom 2.2.1 r14555
library(data.table)
nc <- st_read(system.file("shape/nc.shp", package = "sf"))
#> Reading layer `nc' from data source `/usr/local/lib/R/site-library/sf/shape/nc.shp' (...)
#> Simple feature collection with 100 features and 14 fields
#> geometry type:  MULTIPOLYGON
#> dimension:      XY
#> bbox:           xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
#> epsg (SRID):    4267
#> proj4string:    +proj=longlat +datum=NAD27 +no_defs
nc <- setDT(nc)

nc[AREA > 0.1, st_area(geometry)]
#> Units: m^2
#>  [1] 1137388604 1423489919 1520740530 1179347051 1232769242 1136287383
#>  [...]
#> [61] 1264050755 2288682896 2181174167 2450241815 2165843695
nc[AREA > 0.1, sum(st_area(st_union(geometry))), by = SID74]
#>     SID74              V1
#>  1:     1  3598847644 m^2
#>  2:     5 11299175600 m^2
#> [...]
#> 20:    15  3850824500 m^2
#> 21:    29  1978619669 m^2
#> 22:    31  2439553215 m^2
#>     SID74              V1

But returning geometries is a problem

nc[, st_union(geometry), by = SID74]
#>     SID74
#>  1:     1
#> [...]
#> 83:    31
#>     SID74
#>                                                                                                                           V1
#>  1:                                                                                                                   <list>
#>  2:                                                                                                                   <list>
[...]
#> 83: -78.86451,-78.91947,-78.95074,-78.97536,-79.00224,-79.00642, 34.47720, 34.45364, 34.44938, 34.39917, 34.38804, 34.36627,
#>                                                                                                                           V1

Another path is to nc <- st_as_sf(setDT(nc)), but there are still problems.

A quick note that @SymbolixAU has been working on a project creating a spatial.data.table More info HERE. I'm not sure creating a new class is the way forward, but I'm sure there is a good overlap between all these projects.

@etiennebr

The problem @etiennebr is describing is because aggregation of geometry can lead to two different kinds of geometries (POLYGON and MULTIPOLYGON in this case). Try the following:

nc <- st_read(system.file("shape/nc.shp",package="sf"))
nc_DT <- as.data.table(nc)
nc %>% group_by(SID74) %>% summarise(geom = st_union(geometry))
nc_DT[,st_union(geometry),by=SID74][,table(SID74)] # indices where frequency > 1 are multipolygon geometries!

I still don't understand why data.table shows 83 rows for nc_DT[, st_union(geometry),by=SID74]. Regardless of the geometry aggregation issue shouldn't data.table still return only 23 rows here?

I have realized my error and I know that I have to do nc_DT[, .(st_union(geometry)), by=.(SID74)] but that still doesn't give me a valid object that I can plot!

nc %>% group_by(SID74) %>% summarise(geom=st_union(geometry)) %>% plot  # works!
plot(nc_DT[, .(st_union(geometry)), by=.(SID74)]) # does not work...some weird error!

Am I right in thinking it would still be nice if setDT() retained the sf in the class?
Or can this issue be closed now on the data.table side. The solutions above are workarounds?
I reran the reproducible example above with 1.12.0 but it's still the same. Note that geometry column is not really lost though. That's just a printing/display difference. It's still all there :

> nc$geometry
Geometry set for 100 features 
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
epsg (SRID):    4267
proj4string:    +proj=longlat +datum=NAD27 +no_defs
First 5 geometries:
MULTIPOLYGON (((-81.47276 36.23436, -81.54084 3...
MULTIPOLYGON (((-81.23989 36.36536, -81.24069 3...
MULTIPOLYGON (((-80.45634 36.24256, -80.47639 3...
MULTIPOLYGON (((-76.00897 36.3196, -76.01735 36...
MULTIPOLYGON (((-77.21767 36.24098, -77.23461 3...
>

Instead of putting sf last in the class I tried putting it first :
setattr(nc, "class", c("sf", "data.table", "data.frame"))
but this still resulted in problems. Perhaps you could try that and see if those problems can be overcome.
The root issue is perhaps that sf has its own [.sf method which either masks or is masked by [.data.table depending on whether sf is before or after data.table in class().
The class attribute being c("sf", "data.table", "data.frame") makes sense to me on first glance. data.table could change setDT() to retain sf like that. But then sf would likely need to become data.table-aware (import data.table) and call NextMethod() inside its [.sf if inherts(x,"data.table"). Something like that.

I think it may fit better to sf package, then it may be reused in other packages as well.

@jangorecki

@jangorecki Sounds good. I'll drop another note here once data.table v1.14.4 is out if my PR then gets pulled into the sf package.

Certainly my intention when writing the PR was that downstream packages would define methods for their own classes. sf is in a slightly unique case where it doesn't Depend/Import/Suggest data.table, so it's a bit more of a toss-up... still, I think hosting it in sf would be more appropriate.

If they are unwilling to merge we could do so here, since IINM we don't need to add any dependency, just to create an .sfg method.

I had a look at your PR, it looks quite elaborate for just adding a method! I am assuming because of the invisible-line dependency issue.

But the meat of your PR is very simple:

format_list_item.sfg = function(x, ...) {
    format.sfg(x)
}

Maybe format_list_item.default could check if there's a format method registered for the class of the object, and run that as a default, and only lacking that back up to the current <{class}[dim]> approach.

I'm not great with S3 registration stuff, but if this sounds easy to you, a PR would be great.

Hmmm. Looks I'm in the minority here, thinking that an internal data.table method is the simplest and cleanest thing to do. As Michael (and my earlier post) suggest, it would basically require a single (unexported) line of code.

format_list_item.sfg = function(x, ...) format(x)

Similarly, I don't think this method would change the behaviour for any other data frame(ish) object on the sf side. Both tibbles and vanilla data frames would keep their existing print behaviour. So I'm not sure there's much incentive on their end?

Final argument for keeping this on the data.table side: Converting an sf object to a data.table must not only be done explicitly, but also necessitates a slight change in workflow. (As documented in this thread, we need to refer to the geometry column explicitly when performing spatial operations, e.g. st_union(nc) vs nc_dt[, st_union(geometry)].) That's a perfectly acceptable trade-off from my perspective, for those of us that want to combine the power of the two packages. But the point is we're selecting into this modified workflow. IMHO the print method is just a small (but much appreciated!) convenience layer from the data.table side further enabling this integration.

@grantmcdermott

Now that I've thought about this a bit more, I think I'm with @grantmcdermott on this one.

It'd be really nice to be able to sidestep the elaborate hoops I go through in that sf-side PR just to be able to register the method without adding data.table to sf as an Imports or Depends. As both @grantmcdermott and @MichaelChirico note, adding the method on the data.table side would require just three lines of code and no additional dependencies.

FTR, I looked into @MichaelChirico suggestion that we could possibly modify the code of data.table:::format_list_item.default(). That solution (as implemented by the code below) works just fine for a simple feature geometry. Unfortunately, it's not a great general approach, since it'll cause all sorts of problems when list items have format methods that don't produce such nicely succinct output For example, it fails (throwing an error) when printing a data.table such as this: x <- data.table(list(mtcars, mtcars))

## Alternative implementation, that uses each list element's class to dispatch a format method
format_list_item.default <- function (x, ...)
{
    exists_format_method <- function(x) {
        !(is.null(getS3method("format", class = x, optional = TRUE)))
    }
    if (is.null(x))
        ""
    else if (is.atomic(x) || inherits(x, "formula"))
        paste(c(format(head(x, 6L), ...), if (length(x) > 6L) "..."),
               collapse = ",")
    else if (any(sapply(class(x), exists_format_method)))
        format(x)
    else
        paste0("<", class(x)[1L], paste_dims(x), ">")
}

Thanks for investigating Josh. point taken... I am only trying to avoid a maintenance headache if dozens of formatters are eventually requested in data.table.

the most scalable thing seems to me to be keeping a list of classes where the format method plays well.

but I think that can be kicked down the road a bit. I would be fine with a simple PR to data.table as described (with NEWS and mention it in the .Rd)

Add format_list_item() method for simple feature geometry columns

@MichaelChirico

@MichaelChirico That makes total sense, and I like your idea of, in the longer term, using a list to record classes with format methods that produce reasonable print.data.table-ready output.

In the PR I just pushed, I added this as a new item in NEWS, but could also see folding it into NEWS item 17, which announced the format_list_item() generic.

Also, I didn't add an example to the *.Rd file, only because I didn't want to induce any dependency on sf or trigger any CRAN-check warnings. Let me know if you think an example would be good, though, and I'll see how best to add one that'll satisfy CRAN's requirements.

Navigation Menu

Repository navigation

compatibility with sf library #2273

Reproducible Example

12 remaining items

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

compatibility with sf library #2273

Description

Activity

MichaelChirico commented on Jul 17, 2017

mattdowle commented on Jul 17, 2017

MichaelChirico commented on Jul 17, 2017

rafapereirabr commented on Jul 18, 2017

Reproducible Example

rafapereirabr commented on Jul 18, 2017

etiennebr commented on Jul 18, 2017

rafapereirabr commented on Jul 26, 2017

vlulla commented on Sep 7, 2017

vlulla commented on Jan 29, 2019

vlulla commented on Jan 29, 2019

mattdowle commented on Jan 29, 2019

mattdowle commented on Jan 30, 2019

12 remaining items

jangorecki commented on Oct 14, 2021

JoshOBrien commented on Oct 14, 2021

MichaelChirico commented on Oct 15, 2021

grantmcdermott commented on Oct 15, 2021

JoshOBrien commented on Oct 15, 2021

MichaelChirico commented on Oct 15, 2021

JoshOBrien commented on Oct 16, 2021

Metadata

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Issue actions