Performance: import of many points is much faster in sp than in sf

[pbaylis]

I'd like to import many points as quickly as possible. It was taking longer than expected in sf, so I tried sp and was surprised to find that it was much quicker. I'm not sure if this is issue-appropriate, but perhaps this issue is unique to me? I'm running OS X, sf 0.6.0.

library(sf)
library(sp)

N <- 1000000

df <- data.frame(lng = runif(N, -120, -100), 
                 lat = runif(N, 30, 48))

tic <- Sys.time()
df_sf <- st_as_sf(df, coords = c("lng", "lat"), crs = "+proj=longlat +datum=WGS84")
print(Sys.time() - tic)
# Time difference of 22.33289 secs

tic <- Sys.time()
coordinates(df) <- ~lng+lat
proj4string(df) <- "+proj=longlat +datum=WGS84"
print(Sys.time() - tic)
# Time difference of 0.8412709 secs

[rsbivand]

In sp, SpatialPoints objects hold the coordinates in a matrix, so:

mat <- cbind(lng = runif(N, -120, -100), lat = runif(N, 30, 48))
system.time(sp_mat <- SpatialPoints(mat, proj4string=CRS("+proj=longlat +datum=WGS84")))
#   user  system elapsed 
#  0.071   0.000   0.071 

An sf MULTIPOINT geometry does the same:

 system.time(df_sf <- st_sfc(st_multipoint(mat), crs = "+proj=longlat +datum=WGS84"))
#   user  system elapsed 
#  0.029   0.000   0.029 

but casting to single POINT objects is costly:

system.time(df_sf1 <- st_cast(df_sf, "POINT"))
#   user  system elapsed 
# 22.808   0.000  22.881 

(Unchecked) sf objects will take less time in topological predicates than sp objects because they are already simple features.

Writing and reading are equivalent:

td <- tempdir()
system.time(st_write(df_sf1, file.path(td, "sf_write.gpkg"), driver="GPKG"))
#Writing layer `sf_write' to data source `/tmp/RtmploLAJz/sf_write.gpkg' using driver `GPKG'
#features:       1000000
#fields:         0
#geometry type:  Point
#   user  system elapsed 
# 15.345   3.132  18.512 
system.time(df_sf2 <- st_read(file.path(td, "sf_write.gpkg")))
#Reading layer `sf_write' from data source `/tmp/RtmploLAJz/sf_write.gpkg' using driver `GPKG'
#Simple feature collection with 1000000 features and 0 fields
#geometry type:  POINT
#dimension:      XY
#bbox:           xmin: -120 ymin: 30.00002 xmax: -100.0001 ymax: 47.99999
#epsg (SRID):    4326
#proj4string:    +proj=longlat +datum=WGS84 +no_defs
#   user  system elapsed 
# 11.500   0.018  11.543 

system.time(sp_read <- readOGR(file.path(td, "sf_write.gpkg")))
#OGR data source with driver: GPKG 
#Source: "/tmp/RtmploLAJz/sf_write.gpkg", layer: "sf_write"
#with 1000000 features
#It has 0 fields
#   user  system elapsed 
#  9.923   0.064  10.008 
system.time(writeOGR(sp_read, file.path(td, "sp_write.gpkg"), layer="1", driver="GPKG"))
#   user  system elapsed 
# 12.622   3.196  15.844 

[edzer]

This doesn't surprise me; I'm afraid the convenience of having points in a list-column comes with a cost. What we still could try is construct the object in C or C++, but I'm not sure that will help much.

[pbaylis]

Thanks for your responses. My actual task involves extracting raster values at these points using velox, so I was hoping that it would be able to use MULTIPOINT instead of a list of POINTs, but no dice there. Maybe something to ask the velox team. I'll close this now.

[mdsumner]

I did some work on this and got sidelined, it's a bit surprising still to me what is fast in R/Rcpp and vice versa, but I'm pretty limited in C++.

Part of the overhead is the structure() call, copying that object as a template and populating its (POINT) vector/matrix brought the time down from 20 to 6, which is good but I think the limit of what can be done at the R level.

I can't separate out the code properly right now but here's what I tried, it's quite a lot faster and I'm happy to follow up with a proper implementation iyi.

// [[Rcpp::export]]
List points_cpp(NumericMatrix pts, CharacterVector gdim = "XY") {
  int n = pts.nrow();
  List out(n);
  CharacterVector cls = CharacterVector::create(gdim[0], "POINT", "sfg");
  for (int i = 0; i < n; i++) {

    NumericVector lp = pts(i, _);
    lp.attr("class") = cls;
    out[i] = lp;
  }
  return out;
}

## x is a matrix, g a list of POINT sfg
  g <- sfcc::points_rcpp(x)

## add class and then everything else precision, n_empty, bbox etc. 
class(g) <- c("sfc_POINT", "sfc")
...

[Robinlovelace]

That will be awesome if st_as_sf() can be sped-up, used lots. Look forward to benchmarks!

[edzer]

yes, LGTM; I would also add somewhere a validator that the number of characters in gdim matches the number of columns in pts.

[rsbivand]

Does the point representation in PDAL relate to SF? It seems at least to have gdim as "XYZ", but may also have multiple attributes (Intensity, Time, Red, Green, and Blue for example)? I feel that SF and point clouds are probably orthogonal, aren't they?