I have two simulated linestrings, one that represents a satellite ground-track and the other representing a linear surface feature (e.g., a road or river). I am trying to calculate the minimum distance the two features.

From searching stack exchange, I have only found minimum distance methods for polygons to points, points to points, lines to points, etc. but cannot seem to find line to line examples.

My current code approach is as follows:


# Create 'coordinate'  dataframes -----------------------------------------------

#Simulate satellite transect coordinates
x0 <- runif(1, 1,10)
y0 <- runif(1, 1,10)
x1 <- runif(1, 1,10)
y1 <- runif(1, 1,10) 

#Simulate road coordinates
l_x0 <- runif(1, 1,10)
l_y0 <- runif(1, 1,10)
l_x1 <- runif(1, 1,10)
l_y1 <- runif(1, 1,10) 

# Set up a dataframe for the satellite transect
transect.df <- data.frame(
  'x' = c(x0, x1),
  'y' = c(y0, y1)

# Set up a dataframe for the simulated road info
road.df <- data.frame(
  'x' = c(l_x0, l_x1),
  'y' = c(l_y0, l_y1)

# Create 'sf'  objects -----------------------------------------------

road.sf <- st_sf(id = 'L01', st_sfc(st_linestring(as.matrix(road.df)))) #create a formal sf feature
st_crs(lead.sf) <- st_crs(4326) # assign crs
road.sf <- st_transform(road.sf, crs = 4326) # transform

transect.sf <- st_sf(id = 't01', st_sfc(st_linestring(as.matrix(transect.df)))) #create a formal sf feature
st_crs(transect.sf) <- st_crs(4326) # assign crs
transect.sf <- st_transform(transect.sf, crs = 4326) # transform

First, you do not need to "transform" your data if you are not projecting it to a different coordinate system. Your code is assigning the projection and then reprojecting to the same geographic coordinate space.

Now, is there a particular reason that you think that methods that apply to point or polygon to line distances do not apply here? You can derive a distance matrix for lines as well. For geographic projections, geod_inverse from PROJ, is used to calculate great circle distance. Here is a quick example, using your input data, where I add an extra line into road, return distances and subset the closest line from the road data. For more complex problems (ie., multiple line features in each object) you will have to work with a pairwise matrix but using which.min with apply will make it easy to return the feature index.


transect <- st_sf(id = 't01', st_sfc(st_linestring(
              cbind(c(runif(1, 1,10), runif(1, 1,10)),
              c(runif(1, 1,10), runif(1, 1,10)))))) 
  st_crs(transect) <- st_crs(4326) 
road <- rbind(
  st_sf(st_sfc(st_linestring(cbind(c(runif(1, 1,10), runif(1, 1,10)),
                             c(runif(1, 1,10), runif(1, 1,10)))))),
  st_sf(st_sfc(st_linestring(cbind(c(runif(1, 1,10), runif(1, 1,10)),
                             c(runif(1, 1,10), runif(1, 1,10)))))) )
    road$ID <- 1:2                       
      st_crs(road) <- st_crs(4326) 

Here is the distance between the lines.

( d <- st_distance(transect, road, by_element=TRUE) )

Here we can plot the roads (black dotted), transects (blue), and closest road (red).

plot(road[,1], lty=2)
  plot(transect[,1], col="blue", lwd=2, add=TRUE)
    plot(road[which.min(d),][,1], col="red", lwd=2, add=TRUE) 

Since the distances are based on nodes, you may need to densify nodes along the lines to increase accuracy.

  • Thanks so much @Jeffrey Evans. Your input was very helpful and I appreciate the time taken to assist. – Michael Joseph Feb 10 at 21:47

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.