How to efficiently convert origin-destination coordinates into lines in R

Question

Reproducible example to show what I mean:

# Create origins and destinations
odf = data.frame(fx = 1:5, fy = 1:5, tx = 0:4, ty = rep(1, 5))

The challenge is to convert this into straight lines connecting the origins (fx and fy) to the destinations (tx and ty). Here's the best way I could come up with, but it seems clunky and inefficient - is there a faster way to do it, without for loops?

l <- vector("list", nrow(odf))
for(i in 1:nrow(odf)){
  o = c(odf$fx[i], odf$fy[i])
  d = c(odf$tx[i], odf$ty[i])
  l[[i]] <- sp::Lines(list(sp::Line(rbind(o, d))), as.character(i))
}
l <- sp::SpatialLines(l)
plot(l)

Answer 1

Two lapply calls - is it faster? Dunno.

> V = lapply(1:nrow(odf), 
   function(i){return(list(L=Line(matrix(unlist(odf[i,]),ncol=2,byrow=TRUE)),i=i))})
> VV = SpatialLines(lapply(V, function(E){Lines(list(E$L),as.character(E$i))}))
> plot(VV)

Answer 2

library(uuid)
library(sp)

sp::SpatialLines(apply(odf, 1, function(r) {
  sp::Lines(list(sp::Line(cbind(r[c(1,3)], r[c(2,4)]))), uuid::UUIDgenerate())
}))

For funsies (I expected the following result):

library(sp)
library(uuid)
library(microbenchmark)
library(ggplot2)

odf <- data.frame(fx = 1:100,
                  fy = 1:100, 
                  tx = 0:99, 
                  ty = rep(1, 100))

f_apply <- function() {
  sp::SpatialLines(apply(odf, 1, function(r) {
    sp::Lines(list(sp::Line(cbind(r[c(1,3)], r[c(2,4)]))), uuid::UUIDgenerate())
  })) -> l
}

f_for <- function() {
  l <- vector("list", nrow(odf))
  for(i in 1:nrow(odf)){
    o = c(odf$fx[i], odf$fy[i])
    d = c(odf$tx[i], odf$ty[i])
    l[[i]] <- sp::Lines(list(sp::Line(rbind(o, d))), as.character(i))
  }
  l <- sp::SpatialLines(l)
}

f_vapply <- function() {
  V <- lapply(1:nrow(odf), 
             function(i){return(list(L=Line(matrix(unlist(odf[i,]),ncol=2,byrow=TRUE)),i=i))})
  VV <- SpatialLines(lapply(V, function(E){Lines(list(E$L),as.character(E$i))}))
}

mb <- microbenchmark(f_apply(), f_for(), f_vapply())

autoplot(mb)

One thing to consider (if you have alot of these lines) is to make a CSV and let the ogr cmdline utils do the work for you:

Consider the following CSV file (test.csv):

way_id,pt_id,x,y
1,1,2,49
1,2,3,50
2,1,-2,49
2,2,-3,50

With a GDAL build with Spatialite enabled, `ogrinfo test.csv -dialect SQLite -sql "SELECT way_id, MakeLine(MakePoint(CAST(x AS float),CAST(y AS float))) FROM test GROUP BY way_id"` will return :

OGRFeature(SELECT):0
  way_id (String) = 1
  LINESTRING (2 49,3 50)

OGRFeature(SELECT):1
  way_id (String) = 2
  LINESTRING (-2 49,-3 50)

(ripped from http://www.gdal.org/drv_csv.html)

That can be done via:

library(rgdal)
library(readr)
library(tools)

f_ogr <- function() {

  csv <- tempfile(fileext=".csv")
  shp <- tempfile(fileext=".shp")

  readr::write_csv(odf2, csv)

  # there's a way to get the following to work with gdalUtils::ogr2ogr, but 
  # it ends up calling the same system command, so this is just more explicit

  system(sprintf('ogr2ogr -f "ESRI Shapefile" %s -dialect SQLite -sql "SELECT way_id,MakeLine(MakePoint(CAST(x AS float),CAST(y AS float))) FROM %s GROUP BY way_id" %s',
                 shp, tools::file_path_sans_ext(basename(csv)), csv))

  lines <- rgdal::readOGR(shp, tools::file_path_sans_ext(basename(shp)), verbose=FALSE)

  unlink(csv)
  unlink(shp)

}

Since I stipulated that this was for larger point lists, we need to even the playing field:

odf <- data.frame(fx = 1:100000,
                  fy = 1:100000,
                  tx = 0:99999,
                  ty = rep(1, 100000))

odf2 <- data.frame(x = c(1:100000, 0:99999),
                   y = c(1:100000, rep(1, 100000)),
                   way_id = c(1:100000, 1:100000),
                   pt_id = c(rep(1, 100000), rep(2, 100000)))

This many pts takes a while no matter what so I ended up limiting the times to 10 and got this for the microbenchmark:

Unit: seconds
       expr      min       lq     mean   median       uq      max neval cld
    f_ogr() 49.02599 50.13780 52.45603 50.99081 53.40040 61.10121    10  a 
  f_apply() 49.56835 51.02472 53.79773 52.73900 57.15918 58.54167    10  a 
    f_for() 52.26334 53.20995 55.51910 55.61787 57.80937 58.76571    10  a 
 f_vapply() 78.09049 80.62011 84.16231 83.60744 84.87795 93.89976    10   b

Not sure why the vapply is slower (I didn't poke at it, tho). But my hypothesized speedup method is, in fact, not fast enough to warrant using it (IMO) — especially with me not doing the data.frame transform as part of the test process :-) Letting it run 100x might show a more profound difference between the remaining 3, but that is an exercise left to the reader :-)