Search code examples
rlineintersection

Identify the intersected line between a path and multiple lines with R

I have two objects: one is a path and one contains multiple lines (used in this map).
The lines are spatial lines (GIS shapefile), the data looks like this and the path was simulated with R and then turned into a spatial line also .

I wanted to get the number of intersections between the path and the lines and the coordinates of the intersection points. I did that with this code:

intersection <- intersect(path, lines)
length(intersection)
icoords <- intersection@coords
icoords

My problem is that I also want to know which lines the path crossed. I have the coordinates of the intersection points, but not the lines to which they correspond.

Hence I want output like in the image above, but with the names of the crossed lines included on the right.

Solution

  • If you use the sf library, you can find the lines that cross with st_intersects. Here is a simplified example.

    library(sf)

# Create a line
ex1 <- st_as_sf(data.frame(id=c(1), lon=c(77,79), lat=c(37,40)), coords=c("lon","lat"), crs="epsg:4326")
l1 <- st_cast(st_combine(ex1$geometry), "LINESTRING")

# Another line. This will cross it.
ex2 <- st_as_sf(data.frame(id=c(1), lon=c(77,79), lat=c(40,37)), coords=c("lon","lat"), crs="epsg:4326")
l2 <- st_cast(st_combine(ex2$geometry), "LINESTRING")

# Another line. This will cross neither.
ex3 <- st_as_sf(data.frame(id=c(1), lon=c(74,76), lat=c(33,34)), coords=c("lon","lat"), crs="epsg:4326")
l3 <- st_cast(st_combine(ex3$geometry), "LINESTRING")

# Combine one of the crossing lines and the one that doesn't cross
ex4 <- st_cast(c(l1, l3), "MULTILINESTRING")

# Get results of the cross. List the multiline object first.
st_intersection(ex4, l2)
#Results of intersec*tion*
Geometry set for 1 feature 
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 78 ymin: 38.5355 xmax: 78 ymax: 38.5355
Geodetic CRS:  WGS 84
POINT (78 38.5355)

# See which cross of the combined line object. Again, list multiline object first.
st_intersects(ex4, l2)

    This is the single line object, l2.

    enter image description here

    This is the double line, ex4.

    enter image description here

    And these are the two objects plotted together.

    enter image description here

    So you can see that only one of the double lines intersects. With st_intersects(x4, l2), you get this:

    st_intersects(ex4, l2)
Sparse geometry binary predicate list of length 2, where the predicate
was `intersects'
 1: 1
 2: (empty)

    So you know it's the first line of the double line object that intersected the single line.

    EDIT: REPLY to @cfg comments.

    Note: normally, people don't edit their answers to reply to comments, but I worry this would just turn into a long thread of comments after the answer that I want to avoid. Instead, you should edit your question by adding more to the question (don't replace what's already there, add to it) after some kind of break like "Edit" if you have as much detail as you provided in your 3 comments. The comments are for asking for small clarifications or such, not a bunch of related code.

    You wrote in one comment:

    C<-st_as_sf(cables) 
c1 <- st_cast(C, "MULTILINESTRING") 
o<-st_as_sf(list.df$path1) 
o1 <- st_cast(st_combine(o$geometry), "LINESTRING") 
st_intersection(c1, o)

    First of all, you don't need st_cast() after the st_as_sf(). So you can just use "C" and "o".

    Secondly, st_intersects(), NOTE that's intersects not intersection, tells you which ones were intersected. It will tell you in the order the lines are listed in the object, so if the middle line is listed first for some reason, that's "1".

    st_intersection(c1, o) 
Simple feature collection with 2 features and 1 field 
Geometry type: POINT 
Dimension: XY 
Bounding box: xmin: -2.015635 ymin: 56.45792 xmax: -0.771279 ymax: 57.05269 Geodetic CRS: WGS 84 
fid geometry 
0 1 POINT (-0.771279 57.05269) 
2 3 POINT (-2.015635 56.45792)

    You have two intersections, but I don't know how they relate to anything in your objects, because I don't know what "C" or "o" actually are. You could guess by looking at the coordinates. Which one or two of you lines run near -2 Long. and 56.5 Lat.? Anyway, st_intersects() tells you, so you don't have to guess.