5

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)
3
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)

enter image description here

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

enter image description here

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 :-)

| improve this answer | |
  • Any idea of how to get that back into R again, using the latter method? The mean of all 3 (4 including Barry Rowlingson's) seem to be within 10% of each other so there must be one to split the pack! – RobinLovelace Jun 4 '16 at 6:22
  • Although both answers are not incorrect I'm ticking yours as 'right' due to the amazing amount of information transmitted in your answer. Your solution using ogr is awesome. Thanks so much! – RobinLovelace Jun 9 '16 at 15:54
4

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)
| improve this answer | |

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.