Cómo convertir de forma eficiente origen-destino de coordenadas en líneas en R

Question

Reproducible ejemplo para mostrar lo que quiero decir:

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

El reto es convertir esto en líneas rectas conexión de los orígenes (fx y fy) a los destinos (tx y ty). Aquí es la mejor manera que se me ocurrió, pero parece torpe e ineficiente - ¿hay una manera más rápida de hacerlo, sin que los bucles?

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

Dos lapply llamadas - es más rápido? No sé.

> 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())
}))

Para funsies (yo esperaba que el siguiente resultado:

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)

Una cosa a tener en cuenta (si tienes un montón de estas líneas) es hacer un CSV y dejar que el ogr cmdline utils hacer el trabajo por usted:

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)

(rasgado de http://www.gdal.org/drv_csv.html)

Que se puede hacer a través de:

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)

}

Desde que se estipula que este era para el mayor punto de listas, tenemos que igualar el campo de juego:

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

Esto muchas pts toma un tiempo no importa qué, así que terminó la limitación de la times a 10 y tiene esto para la 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

No sé por qué la vapply es más lento (que no me metiera en eso, tho). Pero mi hipótesis de speedup método es, de hecho, no lo suficientemente rápido como para justificar el uso de (OMI) - especialmente conmigo no hacer la data.frame de transformación como parte del proceso de prueba :-) Permitir que se ejecute 100x podría mostrar una profunda diferencia entre los 3 restantes, sino que es un ejercicio izquierda para el lector :-)