How to measure length of a section sf linestring inserted between two features from sf multilinestri-CodePudding

I have a dataset in linestring format and an ocean bathymetry multilinestring (several lines with depth information). My aim is to measure the length section of these linestrings at each depth range (between two multilinestrings) and to have the sum lengths per interval.

EXPECTED OUTPUT EXAMPLE

       interval        length 
1  -2250 and -2500    5200.56 [m]
2  -2500 and -2750    xxxxxxx [m]
3  -2750 and -3000    xxxxxxx [m]
4  -3000 and -3250    xxxxxxx [m]

But, when I use st_intersection() I can get what depth this stretch of line is, but besides the information that comes about the depth being for example -1500 m instead of an interval, (example: -1500 m to -1750 m) I can't measure the length as st_intersection() only counts intersection points. Are there way to make this in r?

DATASET EXAMPLE

library(sf)
library(dplyr)


#creating dataset
id <- c("A","A", "B","B","C","C")
lat <- c(-25.31157, -25.42952, -25.4253, -25.19177, -25.18697, -25.12748)
long <- c(-41.39523, -39.99665, -41.00311, -41.29756, -41.30314, -39.37707)

df <- dplyr::tibble(id = as.factor(id), lat, long)

#convert sf
df.sf = sf::st_as_sf(df,
                      coords = c("long","lat"), 
                      crs    = 4326)

#creating linestrings
lines <- df.sf %>% 
  dplyr::group_by(id) %>%
  dplyr::summarise(do_union = FALSE) %>%
  sf::st_cast("LINESTRING")

######attempt at something similar:

#intersect
intersection <- sf::st_intersection (bathy, lines) %>% sf::st_as_sf() %>%
  dplyr::mutate(length=sf::st_length(.)) %>% 
  sf::st_drop_geometry()

#sum length per depth 
output <- intersection %>% 
  group_by(depth) %>% 
  summarise(length = sum(length)) %>%
  arrange(desc(length)) 

> output
   depth   length[m]
1  -2250     0 [m]
2  -2500     0 [m]
3  -2750     0 [m]
4  -3000     0 [m]

I unfortunately couldn't recreate a subset of the multilinestring object here, so I put in my github to download a subset in shapefile format if necessary, just click on "CODE" and then on "download ZIP" enter link description here. The object is like this:

> bathy
Simple feature collection with 15690 features and 1 field
Geometry type: MULTILINESTRING
Dimension:     XY
Bounding box:  xmin: -53.37185 ymin: -35.77762 xmax: -25 ymax: 7.081626
Geodetic CRS:  WGS 84
First 10 features:
   PROFUNDIDA                       geometry
1         -50 MULTILINESTRING ((-52.4267 ...
2         -50 MULTILINESTRING ((-52.77632...
3         -75 MULTILINESTRING ((-51.04274...
4         -75 MULTILINESTRING ((-52.38656...
5        -100 MULTILINESTRING ((-51.07005...
6        -100 MULTILINESTRING ((-52.18633...
7        -200 MULTILINESTRING ((-51.97665...
8        -300 MULTILINESTRING ((-51.95862...
9        -400 MULTILINESTRING ((-51.94465...
10       -500 MULTILINESTRING ((-51.93161...

CodePudding user response：

Here's a rough answer that might help. It would be easier to provide a reprex if your data was smaller. It helps to just pick the absolute minimum amount of data to reproduce your problem.

In this case, all you would need is one linestring for a path, and just a couple bathymetry linestrings. I suggest subsetting your data a lot, and then using the datapasta package or reprex to make an example that can be posted easily.

Assuming intersection is the intersection of just 1 LINESTRING with the bathy data, this might work:

intersection %>% 
  group_by(geomtype = st_geometry_type(geometry)) %>%
  group_modify(~ st_cast(.x, "POINT")) %>% 
  ungroup() %>% 
  select(-geomtype) %>% 
  st_as_sf() ->
  intersection_points

intersection_points %>% 
  mutate(geometry2 = lead(geometry)) %>%
  mutate(length = st_distance(geometry, geometry2, by_element = TRUE))