CODEWITHSUNDEEP
rspatialr-sf
Alexa Fredston
Alexa Fredston

Reputation: 727

Measuring the distance of all points along a line in R (linestring in sf)

I have a smoothed line (a simplified abstraction of a coastline) that is a linestring and I want to measure the length of the line at frequent intervals along it. I can create the smoothed line, and measure its length: 

library(raster)
library(sf)
library(tidyverse)
library(rnaturalearth)
library(smoothr)

# create bounding box for the line 
xmin=-80
xmax=-66
ymin=24
ymax=45
bbox <- extent(xmin, xmax, ymin, ymax)

# get coarse outline of the US 
usamap <- rnaturalearth::ne_countries(scale = "small", country = "united states of america", returnclass = "sf")[1] %>% 
  st_cast("MULTILINESTRING")

# crop to extent of the bbox and get rid of a border line that isn't coastal 
bbox2 <- st_set_crs(st_as_sf(as(raster::extent(-80, -74, 42, 45.5), "SpatialPolygons")), st_crs(usamap))
cropmap <- usamap %>% 
  st_crop(bbox) %>% 
  st_difference(bbox2)

# smooth the line 
smoothmap <- cropmap %>% 
  smoothr::smooth(method="ksmooth", smoothness=8)

# measure the line length
st_length(smoothmap) # I get 1855956m

I'm going to be "snapping" sample sites to points along this line, and I need to know how far they are along the coastline. However, I can't figure out how to measure the length of the line at intervals along the coastline. The interval size isn't terribly important so long as it's at relatively fine resolution, perhaps every 1km or every 0.01 degree.

What I would like to produce is a dataframe with x,y columns containing the lat/lon of points along the line, and a "length" column containing the distance along the line (from the origin to that point). Here are some of the things I've tried:

Iterating over bounding boxes. I tried cropping the line with smaller and smaller bounding boxes (using regular intervals in lat/lon), but because the line bends "backward" around -76, 38, some of the cropped boxes don't encompass the complete line segment that I expected. This approach works for the top right half of the line, but not for the bottom left half--that just returns lengths of zero.

Cropping the extent before implementing the smoother, and then measuring the line. Since the smoother function does not produce the same shape if it is only measuring a segment of the original line, this doesn't actually measure distance along the same line.

Getting the coordinates of the linestring with st_coordinates, trimming off one row (one point on the line), and recasting the remaining coordinates as a linestring. This approach does not produce a single linestring but instead a chain of points (since st_cast doesn't know how to connect them again), so it can't be measured normally.

It would be ideal to "edit" the geometry of smoothmap to delete one row of coordinates at a time, repeatedly measure the line, and write out the end point coordinates and the line length to a dataframe. However, I'm not sure if it's possible to edit a sf object's coordinates without turning it into a dataframe.

Upvotes: 2

Views: 2396

Answers (2)

Kau&#234; Braga
Kau&#234; Braga

Reputation: 117

I believe you need sf::st_line_sample:

# Transform to metric
smoothmap_utm <- st_transform(smoothmap, 3857)

# Get samples at every kilometer
smoothmap_samples <- st_line_sample(smoothmap_utm, density = 1/1000)

# Transform back to a sf data.frame
smoothmaps_points <- map(smoothmap_samples, function(x) data.frame(geometry = st_geometry(x))) %>%
  map_df(as.data.frame) %>% st_sf() %>%
  st_cast("POINT") %>%
  st_set_crs(3857) %>%
  st_transform(4326)

mapview(smoothmaps_points) + mapview(smoothmap)

snap

Getting to your desired output:

# Function to transform sf to lon,lat

sfc_as_cols <- function(x, names = c("lon","lat")) {
  stopifnot(inherits(x,"sf") && inherits(sf::st_geometry(x),"sfc_POINT"))
  ret <- sf::st_coordinates(x)
  ret <- tibble::as_tibble(ret)
  stopifnot(length(names) == ncol(ret))
  x <- x[ , !names(x) %in% names]
  ret <- setNames(ret,names)
  ui <- dplyr::bind_cols(x,ret)
  st_set_geometry(ui, NULL)
}

smoothmaps_points_xy <- sfc_as_cols(smoothmaps_points) %>%
  mutate(dist = cumsum(c(0, rep(1000, times = n() - 1))))

smoothmaps_points_xy

        lon      lat dist
1 -67.06836 44.99794    0
2 -67.07521 44.99383 1000
3 -67.08206 44.98972 2000
4 -67.08891 44.98560 3000
5 -67.09575 44.98149 4000
6 -67.10260 44.97737 5000

Important

But if your ultimate goal is to get points distance in a path, I would recommend checking rgeos::gProject.

Upvotes: 0

Robert Hijmans
Robert Hijmans

Reputation: 47146

If I understand your question, I think you can do this:

x <- as(smoothmap, "Spatial")
g <- geom(x)
d <- pointDistance(g[-nrow(g), c("x", "y")], g[-1, c("x", "y")], lonlat=TRUE)
gg <- data.frame(g[, c('x','y')], seglength=c(d, 0))

gg$lengthfromhere <- rev(cumsum(rev(gg[,"seglength"])))

head(gg)

#          x        y seglength lengthfromhere
#1 -67.06494 45.00000 70850.765        1855956
#2 -67.74832 44.58805  2405.180        1785105
#3 -67.77221 44.57474  2490.175        1782700
#4 -67.79692 44.56095  2577.254        1780210
#5 -67.82248 44.54667  2666.340        1777633
#6 -67.84890 44.53189  2757.336        1774967

tail(gg)
#            x        y  seglength lengthfromhere
#539 -79.09383 33.34224   2580.481       111543.5
#540 -79.11531 33.32753   2542.557       108963.0
#541 -79.13648 33.31306   2512.564       106420.5
#542 -79.15739 33.29874   2479.949       103907.9
#543 -79.17802 33.28460 101427.939       101427.9
#544 -80.00000 32.68751      0.000            0.0

Upvotes: 2

Related Questions