[R-sig-Geo] 'Customizing' spsample

[Paul Hiemstra]

Jean-Paul Kibambe Lubamba wrote:
> Hello everybody,
>
> I am using 'spsample' to derive points along a line and it works perfectly
> by doing the following :
>
> xx <- readShapeLines("line2",proj4string=CRS("+proj=utm +zone=34
> +datum=WGS84"))
> sppts <- spsample(xx, n=50, type="random")
> plot(sppts)
>
> However, I was wondering if the user can somehow 'customize' the sampling
> process, as for example, by separating sampling points by a given
> distance. And this distance being in the same units as the line. The 'n'
> argument in spsample will therefore vary according with the length of each
> line.
>
> Any help is welcome!
>
>
> Jean-Paul
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at stat.math.ethz.ch
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>   
Hi Jean-Paul,

My suggestion would be to take the function sample.Line (from sp) as a 
starting point and modify it to your specifications. You can find the 
source in the /R directory (spsample.R) of the source package, which you 
can download from CRAN.

cheers,
Paul

sample.Line
function (x, n, type, offset = runif(1), proj4string = 
CRS(as.character(NA)),
    ...)
{
    offset = offset[1]
    if (missing(n))
        n <- as.integer(NA)
    if (!is.finite(n) || n < 1)
        return(NULL)
    cc = coordinates(x)
    lengths = LineLength(cc, longlat = FALSE, sum = FALSE)
    if (any(abs(lengths) < .Machine$double.eps)) {
        wl <- which(abs(lengths) < .Machine$double.eps)
        cc <- cc[-(wl), ]
        lengths <- lengths[-(wl)]
    }
    csl = c(0, cumsum(lengths))
    maxl = csl[length(csl)]
    if (type == "random")
        pts = runif(n) * maxl
    else if (type == "stratified")
        pts = ((1:n) - runif(n))/n * maxl
    else if (type == "regular")
        pts = ((1:n) - (1 - offset))/n * maxl
    else stop(paste("type", type, "not available for Line"))
    int = findInterval(pts, csl, all.inside = TRUE)
    where = (pts - csl[int])/diff(csl)[int]
    xy = cc[int, , drop = FALSE] + where * (cc[int + 1, , drop = FALSE] -
        cc[int, , drop = FALSE])
    if (nrow(xy) < 1)
        return(NULL)
    SpatialPoints(xy, proj4string)
}
<environment: namespace:sp>



-- 
Drs. Paul Hiemstra
Department of Physical Geography
Faculty of Geosciences
University of Utrecht
Heidelberglaan 2
P.O. Box 80.115
3508 TC Utrecht
Phone:  +3130 274 3113 Mon-Tue
Phone:  +3130 253 5773 Wed-Fri
http://intamap.geo.uu.nl/~paul