# R/sf - Transform a geometry column made of POINTs and MULTIPOINTs into xy-coordinates columns (x1, y1, x2, y2, ...)

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 ...
``````

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
[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...

• 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. Apr 6, 2021 at 12:40
• Its the kind of problem I like and importantly you gave sample data and desired solution! A lesson to everyone! Apr 6, 2021 at 12:58