Search code examples
rrasterspatialr-sf

How to split network in equal piece line segments in R


I create a SpatialLines from x and y coordinates. I was wondering is there a way to pick 5, 10 or m number points on the network which are equal distance to each other at least (2:m-1) ones have equal distance.

I was thinking maybe I can compute the distance for each point to the previous one and get the cumulative length and use seq(., ., length.out=m) would give me equal distance but I cannot get x and y coordinates from there.

xy <- cbind(xco=c(172868.6, 172891, 172926.8, 172949.2, 172985, 173007.3, 173029.7, 173065.5, 173087.9, 173123.7, 173146, 173181.9, 173204.2, 173226.6, 173262.4, 173320.7, 173343, 173365.4, 173401.3, 173423.7, 173459.6, 173482, 173504.3, 173541.2, 173563.5, 173601.2, 173623.4, 173644.9, 173684.7, 173706.7, 173747.3, 173769.3, 173811, 173832.8, 173875.2, 173896.9, 173918.5, 173962, 173983.6, 174027.1, 174048.6, 174092.1, 174113.7, 174157.2, 174178.8, 174222.2, 174243.7, 174287.1, 174308.7, 174330.3, 174373.6, 174395.2, 174438.9, 174460.4, 174504.2, 174525.7, 174569.4, 174590.9, 174634.5, 174656, 174700.4, 174722, 174742.4, 174790.4, 174860.9, 174914.5, 174932.6, 174988.4, 175005.4, 175063.6, 175079.2, 175138.5, 175200.4, 175213.2, 175276, 175340.7, 175405.8, 175415.9, 175481, 175546.1, 175611.2, 175621.3, 175686.4, 175751.5, 175761.6, 175826.3, 175893.1, 175960.8, 176029.4, 176098.8, 176168.1, 176237.6, 176307.1, 176375.5, 176442.9, 176509.2, 176517.7, 176584.2, 176648.6, 176712.4),
yco=c(3376152, 3376130, 3376096, 3376074, 3376040, 3376019, 3375997, 3375963, 3375941, 3375907, 3375886, 3375851, 3375830, 3375808, 3375774, 3375718, 3375697, 3375676, 3375641, 3375620, 3375586, 3375564, 3375543, 3375509, 3375489, 3375454, 3375434, 3375415, 3375381, 3375362, 3375328, 3375310, 3375276, 3375259, 3375226, 3375209, 3375192, 3375159, 3375142, 3375109, 3375093, 3375060, 3375043, 3375010, 3374994, 3374960, 3374944, 3374911, 3374894, 3374878, 3374844, 3374828, 3374795, 3374778, 3374745, 3374729, 3374696, 3374680, 3374646, 3374630, 3374597, 3374581, 3374567, 3374536, 3374494, 3374464, 3374455, 3374428, 3374420, 3374395, 3374389, 3374366, 3374346, 3374341, 3374324, 3374307, 3374292, 3374289, 3374274, 3374258, 3374243, 3374240, 3374225, 3374209, 3374207, 3374192, 3374181, 3374172, 3374166, 3374162, 3374161, 3374163, 3374167, 3374175, 3374184, 3374197, 3374198, 3374213, 3374230, 3374248))

library(sp)
myLine190 <- Line(xy)
myL190 <- Lines(list(myLine190), ID = 1)
Spt1 <- SpatialLines(list(myL190))
Y1 <- coordinates (Spt1)  # the same as xy above

Cumulative distance

DD1 = sqrt(rowSums(do.call("rbind", lapply(1:(dim(Y1)[1]-1), function(i){(Y1[i,])-(Y1[i+1,])}))^2))
cd1 = cumsum(c(0,DD1))
grid = seq(min(cd1), max(cd1), length.out = 20)

Solution

  • Example data

    pts <- cbind(xco=c(172868.6, 172891, 172926.8, 172949.2, 172985, 173007.3, 173029.7, 173065.5, 173087.9, 173123.7, 173146, 173181.9, 173204.2, 173226.6, 173262.4, 173320.7, 173343, 173365.4, 173401.3, 173423.7, 173459.6, 173482, 173504.3, 173541.2, 173563.5, 173601.2, 173623.4, 173644.9, 173684.7, 173706.7, 173747.3, 173769.3, 173811, 173832.8, 173875.2, 173896.9, 173918.5, 173962, 173983.6, 174027.1, 174048.6, 174092.1, 174113.7, 174157.2, 174178.8, 174222.2, 174243.7, 174287.1, 174308.7, 174330.3, 174373.6, 174395.2, 174438.9, 174460.4, 174504.2, 174525.7, 174569.4, 174590.9, 174634.5, 174656, 174700.4, 174722, 174742.4, 174790.4, 174860.9, 174914.5, 174932.6, 174988.4, 175005.4, 175063.6, 175079.2, 175138.5, 175200.4, 175213.2, 175276, 175340.7, 175405.8, 175415.9, 175481, 175546.1, 175611.2, 175621.3, 175686.4, 175751.5, 175761.6, 175826.3, 175893.1, 175960.8, 176029.4, 176098.8, 176168.1, 176237.6, 176307.1, 176375.5, 176442.9, 176509.2, 176517.7, 176584.2, 176648.6, 176712.4),
    yco=c(3376152, 3376130, 3376096, 3376074, 3376040, 3376019, 3375997, 3375963, 3375941, 3375907, 3375886, 3375851, 3375830, 3375808, 3375774, 3375718, 3375697, 3375676, 3375641, 3375620, 3375586, 3375564, 3375543, 3375509, 3375489, 3375454, 3375434, 3375415, 3375381, 3375362, 3375328, 3375310, 3375276, 3375259, 3375226, 3375209, 3375192, 3375159, 3375142, 3375109, 3375093, 3375060, 3375043, 3375010, 3374994, 3374960, 3374944, 3374911, 3374894, 3374878, 3374844, 3374828, 3374795, 3374778, 3374745, 3374729, 3374696, 3374680, 3374646, 3374630, 3374597, 3374581, 3374567, 3374536, 3374494, 3374464, 3374455, 3374428, 3374420, 3374395, 3374389, 3374366, 3374346, 3374341, 3374324, 3374307, 3374292, 3374289, 3374274, 3374258, 3374243, 3374240, 3374225, 3374209, 3374207, 3374192, 3374181, 3374172, 3374166, 3374162, 3374161, 3374163, 3374167, 3374175, 3374184, 3374197, 3374198, 3374213, 3374230, 3374248))
    

    Here is a solution that can probably be improved upon

    # interval
    x <- 500
    
    d <- raster::pointDistance(pts[-nrow(pts),], pts[-1,], lonlat=FALSE)
    cd <- cumsum(d)
    n <- trunc(sum(d) / x)
    s <- (1:n) * x
    
    xy <- matrix(ncol=2, nrow=length(s))  
    for (i in 1:length(s)) {
        j <- which(cd >= s[i])[1]
        pd <- raster::pointDistance(pts[j,], pts[j+1,], lonlat=FALSE)
        r <- (s[i] - cd[j-1]) / pd
        xy[i,] <- c(((1 - r) * pts[j,1] + r * pts[j+1,1]), 
                    ((1 - r) * pts[j,2] + r * pts[j+1,2]))
    }
    

    Have a look

    plot(pts, type="l")
    points(xy)
    

    This could be adjusted for lon/lat data by using geosphere::destPoint in the last row of the loop.

    spsample("regular") and findInterval can get you close.

    sp <- raster::spLines(xy)
    s <- sp::spsample(sp, 5, "regular")
    plot(sp)
    points(s)
    

    But that selects nodes, not points inbetween nodes.