2

After intersecting a set of routes (LINESTRING) with a perimeter (say the border of a country) using st_intersection from sf package, I get a list of points and multipoints: there is one point when the route crossed the border once, and there are many points (multipoints) when the route crossed the border several times. My goal would be then to transform this output into columns of xy-coordinates, with NA when necessary.

Example of input:

library(sf)
g = st_sfc(st_multipoint(matrix(1:6,,2)),
           st_point(c(1, 2)),
           st_multipoint(matrix(11:18,,2)))
df = data.frame(a = 1:3, b = 5:7)
st_geometry(df) = g

With df being:

  a b                       geometry
1 1 5 MULTIPOINT ((1 4), (2 5), (...
2 2 6                    POINT (1 2)
3 3 7 MULTIPOINT ((11 15), (12 16...

The desired output would look like:

  a b x1 y1 x2 y2 ...
1 1 5  1  4  2  5 ...
2 2 6  1  2 NA NA ...
3 3 7 11 15 12 16 ...

1 Answer 1

3

Convert the geometry into a data frame of X,Y,L1. Have to cast to MULTIPOINT for consistency in the geometry (it just means a POINT gets cast to a single point MULTIPOINT then all the geoms are MULTIPOINT):

xyL = data.frame(st_coordinates(st_cast(df$geometry,"MULTIPOINT")))

That looks like this, which we can then do some kind of split/process/combine on:

> xyL
   X  Y L1
1  1  4  1
2  2  5  1
3  3  6  1
4  1  2  2
5 11 15  3
6 12 16  3
7 13 17  3
8 14 18  3

As we are going to make a matrix with NA's padding out the rows, we need to know how many points we have in the biggest multipoint so we can pad it out:

maxPoints = max(table(xyL$L1))

and here is a padding function:

padNA = function(v,nmax){c(v, rep(NA,nmax-length(v)))}

which works like this:

> padNA(c(3,2,1),3)
[1] 3 2 1
> padNA(c(3,2,1),5)
[1]  3  2  1 NA NA

Now split the data frame by L1 (the original line in the geometry data frame), convert the X and Y columns in each part to a vector in the right order, pad it out to the right length, and then do an rbind over the lot to get:

do.call(rbind,lapply(split(xyL, xyL$L1), function(m){padNA(c(rbind(m$X,m$Y)), maxPoints*2)}))

  [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
1    1    4    2    5    3    6   NA   NA
2    1    2   NA   NA   NA   NA   NA   NA
3   11   15   12   16   13   17   14   18

That's a matrix you can save and splat onto the original data frame:

mxy = do.call(rbind,lapply(split(xyL, xyL$L1), function(m){padNA(c(rbind(m$X,m$Y)), maxPoints*2)}))

> dfx = cbind(df, mxy)
> dfx
Simple feature collection with 3 features and 10 fields
geometry type:  GEOMETRY
dimension:      XY
bbox:           xmin: 1 ymin: 2 xmax: 14 ymax: 18
CRS:            NA
  a b X1 X2 X3 X4 X5 X6 X7 X8                       geometry
1 1 5  1  4  2  5  3  6 NA NA MULTIPOINT ((1 4), (2 5), (...
2 2 6  1  2 NA NA NA NA NA NA                    POINT (1 2)
3 3 7 11 15 12 16 13 17 14 18 MULTIPOINT ((11 15), (12 16...

If you are fussy about the exact column names you can set them easily enough...

2
  • Thank you very much! You've answered my question and very well explained how it works. I applied it to my real case and it works very well.
    – Antonin
    Apr 6, 2021 at 12:40
  • 1
    Its the kind of problem I like and importantly you gave sample data and desired solution! A lesson to everyone!
    – Spacedman
    Apr 6, 2021 at 12:58

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.