How to split network in equal piece line segments in R

Question

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)

It would be much easier to help you if you provided `Y1` as a matrix. E.g. `Y1 <- cbind(x=c( ), y = c( ))` — Robert Hijmans, Nov 13 '20 at 04:12
You can also use this code `Y1<- read.table(text = 'copy-paste data from the post', header = TRUE)` on first `Y1` it will work — iHermes, Nov 13 '20 at 04:25
Thanks, but that still does not work. I should be able to copy and paste your code into R and run it. Of course I could figure out a way to do this if I really wanted to, but I am not the one asking for help. — Robert Hijmans, Nov 13 '20 at 06:05
`Y134<- read.table(text = ' xco yco 1 172868.6 3376152 2 172891.0 3376130 3 172926.8 3376096 4 172949.2 3376074 5 172985.0 3376040 6 173007.3 3376019 7 173029.7 3375997 8 173065.5 3375963 9 173087.9 3375941 10 173123.7 3375907 11 173146.0 3375886', header = TRUE)` when I just copied this way, it worked on me. Can I email you the code? or problem. Yes, I do need help. And I appreciate for your time. — iHermes, Nov 13 '20 at 06:22
Yes, and if you edit your question so that others can copy that with the data.... Sorry to insist but I am trying to be helpful. — Robert Hijmans, Nov 13 '20 at 07:23
I do not understand you question. Especially what oyu mean with ```but I cannot get x and y coordinates from there```. What should the x and y coordinates exactly represent? What are their purpose? Should they split the two dimensional surface such that each grid cell contains the same number of points? Or should the grid cells just be equal in their x and y length? — MacOS, Nov 13 '20 at 09:09
Also, what you are doing here ```cumsum(c(0,DD1))``` is including the distance from 0 to the first point. Why is this necessary? And here ```seq(min(cd1),max(cd1),length.out = 20)``` you are creating a sequence from 0 to the biggest entry that cumsum returns, which is naturally its last entry. So what you are effectively doing is creating a grid on a two dimensional surface (a table, for example) and splitting this table into equally long grid cells, but the length of the grid cells depends on the distances of the points? So it can happen that some grid cells are not on the table? Correct? — MacOS, Nov 13 '20 at 09:16
@MacOS After picking some points which are equally distant to other points, I will use those x and y coordinates in future analysis. `cumsum()` is the way I thought I could solve the problem. My initial idea was if I can get cumsum for the points the I can figure where the last point located and create seq points between zero to last distance. I couldn't convert back to distance to X and Y coordinates. — iHermes, Nov 13 '20 at 15:32

Robert Hijmans · Accepted Answer · 2020-11-14T03:39:46.490

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.

`sp <- spLines(xy)` this code does not work on my end. Which library it is coming from? — iHermes, Nov 14 '20 at 02:48

score 0 · Answer 2 · answered Nov 13 '20 at 10:04

Here is a posssible solution.

df coordinates is your Y1 df.

dist.matrix <- as.matrix(dist(coordinates, method = "euclidean"))

Create a distance matrix and split it into its diagonals. See here for details.

d <- row(dist.matrix) - col(dist.matrix)
diagonals <- split(dist.matrix, d)

This is the first diagonal below the diagonal

dist.points <- diagonals$`1`

Our cum sum vector.

dist.cum <- cumsum(0, dist.points)

Our distance segments vector.

dist.segments <- seq(min(dist.cum), max(dist.cum), length.out = 20)

This list of boolean vectors tell you up to which position in the dist.points diagonal you can go such that it is within a segment. Consequently, you also know at which point from the original data frame you are. Because entry (2, 1) tells you the distance from point two to one. In other words, this is the key.

segments <- lapply(dist.segments,
                   FUN = function(dist.segment)
                   {
                     return(dist.cum <= dist.segment)
                   })

From there you just have to write a function that returns the point and the next point (by using the segments list of boolean vectors) and calculate the extra distance you need to get to the exact cumsum. And because you know the next point, you can easily calculate the (x, y) coordinates.

HTH!

This is pretty much what I did. Well, I could not get from distance to x and y coordinates, that's the first place why I posted the question. — iHermes, Nov 13 '20 at 15:47

How to split network in equal piece line segments in R

Cumulative distance

2 Answers2