Create line_sample function about stplanr HOT 18 CLOSED

This works:

# Aim: create points proportional to some value

#' @examples
l = flowlines[2:5,]
l$All
npoints = 100
p = line_samplepoints(l, 50)
plot(p)
p = line_samplepoints(l, 50, prob = 1:length(l))
plot(p)
line_samplepoints <- function(l, npoints, prob = 1, type = "regular") {
  # divide points between l elements
  if(length(prob) < length(l))
    prob <- rep(prob, length.out = length(l))
  npoints_lines = round(npoints) * prob / mean(prob) / nrow(l)
  uvals = unique(npoints_lines)
  for(i in 1:length(uvals)) {
    line_sel = npoints_lines == uvals[i]
    if(i == 1)
      p = spsample(l[line_sel,], n = sum(npoints_lines[line_sel]), type = type) else
        p = tmap::sbind(p, spsample(l[line_sel,], n = sum(npoints_lines[line_sel]), type = type) )
  }
  p
}



# test spsample
set.seed(99)
p = spsample(l, 100, "regular")
plot(p)
p2 = spsample(l, 100, "regular", prob = c(1, 2, 100, 1000))
plot(p2)

from stplanr.

@Robinlovelace , sf now has a st_line_sample, which randomly samples a set of LINESTRINGs (sfc), with density possibly set per element. It returns a MULTIPOINT for each LINESTRING. Let me know whether it is useful (I might change the name and interface of this function later on.)

from stplanr.

Nice - sounds useful and faster than my hacked code above.

We may have to import sf sooner than expected!

from stplanr.

Quick question: can you do regular sampling, as in the example above?

Please try to reproduce if you get time.

from stplanr.

Without variable prob field:

With it as prob = 1:length(l)

p = line_samplepoints(l, 50, prob = 1:length(l))
plot(p)

from stplanr.

library(sf)
ls = st_sfc(st_linestring(rbind(c(0,0),c(0,1))),
    st_linestring(rbind(c(0,0),c(.2,0))), crs = 4326)
ls = st_transform(ls, 3857)

png("line.png", 800, 800)
par(mfrow = c(2,2))
x = st_line_sample(ls, density = 1/3000)
plot(ls)
plot(x, add = TRUE, col = 'red')

x = st_line_sample(ls, density = c(1/3000, 1/1000))
plot(ls)
plot(x, add = TRUE, col = 'red')

x = st_line_sample(ls, density = 1/3000, "random")
plot(ls)
plot(x, add = TRUE, col = 'red')

x = st_line_sample(ls, density = c(1/10000, 1/1000), "random")
plot(ls)
plot(x, add = TRUE, col = 'red')

from stplanr.

Great work!

from stplanr.

Makes my code redundant though :( hehe

from stplanr.

It's still unresolved how to do this

• for MULTILINESTRING geometries, and
• for lines on the sphere (longlat coordinates; the GDAL function I use does not provide this, geosphere also seems not to, but has gcIntermediate)

from stplanr.

Good to know - it's a pretty good start though. Solid - would you like me to do some benchmarks on this, e.g. in sfr/tests/benchmarks/st_sample.R?

from stplanr.

Your benchmark will not affect what I do in this particular case. I was more curious whether you can get around with the MULTIPOINT geometries that come out.

from stplanr.

As long as they can go back into a SpatialPoints object that should be fine - this is the kind of application we're aiming for btw, potentially to present to the UK's Department for Transport and then hopefully go international: http://rpubs.com/RobinLovelace/232396

from stplanr.

On a slightly different note, @mem48 is looking at the direct line to raster option but it's being slow on the giant files we have. Any ideas what the timeline is on C implementation of that in raster?

from stplanr.

Are you looking at doing this with the raster package?

from stplanr.

@edzer Yes, but the performance is very poor.

We have lines which represent people traveling to work along the road network, each line has a value for the number of people.

We would like to rasterize those lines with the value of number of people. So we would have a kind of heat map showing busy and quite roads. (See example image, n.b. this also shows a base map, the colours were done in QGIS)

The problems we are having are:

1. If two lines overlap we would like to sum the number of people, rasterize only takes one value per cell.
   We solve this by creating groups of non-overlapping lines, and rasterizing each group separately, then adding the rasters together
2. Due to the large nature of the dataset (about 4 million lines) we end up with many groups (around 5000) each covering the whole of England. To create the raster for one group takes around 45 min, so the total process would take around 150 days!

I had looked at using gdalUtils::gdal_rasterize which was faster but I couldn't get it to work reliably.
I also looked at a new package velox which claims to be 35x faster than raster but currently only rasterizes polygons and not lines.

from stplanr.

Note that gdal_rasterize is now part of GDAL's C++ api, since gdal 2.1.

from stplanr.

@mem48 @Robinlovelace do you have an example for the rasterize slowness? I'd explore a reasonable sized example that was relevant to you. I've been doing exactly this for fisheries impacts, we use raster's cell abstraction to build a "cell index" pool and usually spatstats::pixellate for the local rasterization itself, it's adaptable to direct sums like "time along line" or line intersection counts, and specifically for us repeated impacts within a local cell. Also, cellFromLine is probably fine if density and tricky weighting isn't required. I never actually instantiate the entire high resolution raster, which is suggested by your 45min creation time, so was thinking that might help here.

from stplanr.

Very interesting - imagine Malcolm will be keen - hope is to eventually implement what we did (on 4 million routes using the HPC cluster) in a function and it sounds like your approach could be useful.

from stplanr.