0

This question is about how while using sf::st_intersection on a LINESTRING and a set of POLYGONs to get several LINESTRINGs in the CORRECT order. sf::st_intersection automatically orders the newly created LINESTRINGs based on the names of the POLYGONs you have. To demonstrate this behavior is the following example:

# sample data for linestring
df <- data.frame(
  id = c(2,2),
  lon = c(4,9),
  lat = c(9,5)
)

df <- df %>% st_as_sf(coords = c("lon", "lat"), crs = 4326) %>% 
  group_by(id) %>% 
  summarise(geometry = st_union(geometry), do_union = FALSE) %>% 
  st_cast("LINESTRING")

# function to create set of polygons
create_grid <- function(lonmin, lonmax, latmin, latmax, no_squares) {
  # a matrix defines the bounding box as required by st_make_grid()
  globe_bb <- matrix(c(lonmin,  latmax,
                   lonmax,  latmax,
                   lonmax, latmin,
                   lonmin, latmin,
                   lonmin,  latmax), byrow = TRUE, ncol = 2) %>%
  list() %>% 
  st_polygon() %>% 
  st_sfc(., crs = 4236)

grid <- st_make_grid(globe_bb, n = c(no_squares, no_squares),
                   crs = 4326, what = 'polygons')  %>%
  st_sf('geometry' = ., data.frame('ID' = 1:length(.)))
}
grid <- create_grid(0, 10, 1, 10, 4) %>% st_transform(4326)

# applying `sf::st_intersection` to linestring (df) and set of polygons (grid)
df_in <- df %>% 
  group_by(id) %>% st_intersection(grid)

the current and expected output are explained in the following image:

enter image description here

Based on where the original LINESTRING starts and ends the order of the newly created LINESTRINGs is not correct since it is ordered by the names of the POLYGONs.

Is there a way to set an argument in sf::st_intersection to turn off ordering by names of polygons by default, or is there maybe another approach to this issue ?

  • I don't think so, these structures get completely decomposed in the underlying libs and there's no control over it - it will come out in the order that the algorithm rebuilds the features. Simple features pointedly does not provide identity for sub-feature elements, and most tools that use the standard work around this in their own ways. – mdsumner Sep 13 '18 at 10:05
  • @mdsumner maybe there is a way to get the right sequence and reorder the rows afterwards ? – adl Sep 13 '18 at 10:17
  • If you st_intersection(df_in) you'll see four POINT geometries which are the intersections of the line segments. Using the origins column you can pull out the line segment rows and then with a bit of graph theory work out the order... I'm not sure this is the best way to do it so I've not worked it through yet... – Spacedman Sep 13 '18 at 12:13
  • Yep, it's possible - I think @Spacedman is the world expert on it actually! I find it hard and unworkable usually but am still keen to find ways – mdsumner Sep 13 '18 at 15:26
  • I'm not sure its even guaranteed that the line segments will be in the same direction as the original source line... – Spacedman Sep 13 '18 at 20:46
2

This function takes a data frame of line segments that, when connected, form a continuous chain. This is not tested for and it will fail if this isn't true.

library(igraph) # for connecting the points

reorder_segments <- function(src, parts){

    joints = st_intersection(parts)
    joints = joints[st_is(joints,"POINT"),]
    jgraph = igraph::graph_from_edgelist(do.call(rbind, joints$origins), directed=FALSE)

    ends = which(igraph::degree(jgraph) == 1)

    sps = igraph::shortest_paths(jgraph, ends[1], ends[2])

    path = sps$vpath[[1]]
    return(parts[path,])
}

Here's a helpful function for plotting the order of line segments in a data frame:

show_order <- function(lines,...){
    plot(st_geometry(lines))
    pts = st_coordinates(st_centroid(lines))
    text(pts[,1],pts[,2], 1:nrow(lines),...)
}

If you plot your data like this you'll see the order as in your image:

show_order(df_in,cex=3)

but then...

rs = reorder_segments(df, df_in)
show_order(rs,cex=3)

will show it in the new order - 1,2,3,4,5. There is a problem in that it can be in reverse order - 5,4,3,2,1. This could be tested for by overlaying the segments with the first point of the original data to see which is the first one. Left as an exercise....

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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