1

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. enter image description here

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...
Caroline Portal
  • 195
  • 8
  • Not sure, but possibly `lwgeom::st_split` would help to split a line up by where it intersects the bathymetry line. – jmw Apr 27 '22 at 16:22
  • I think that like each multilinestring of bathy is a different feature, ```st_intersection()``` will not measure something between two features as it does not come from an intersection in single feature. – Caroline Portal Apr 27 '22 at 21:00

1 Answers1

2

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:

# st_intersects returns both POINT and MULTIPOINT geometries.
# It returns MULTIPOINT when the line x intersects the feature y 
# multiple times 
# We want just POINT geometries
# this is an unfortunately complex:
intersection %>% 
  group_by(geomtype = st_geometry_type(geometry)) %>%
  group_modify(~ st_cast(.x, "POINT")) %>% 
  ungroup() %>% 
  select(-geomtype) %>% 
  st_as_sf() ->
  intersection_points

# Once we have one POINT for each intesection, 
# along with the associated information 
# (Depth in this example) we can calculate the distance
# between adjacent points using dplyr::lead and sf::st_distance
intersection_points %>% 
  mutate(geometry2 = lead(geometry)) %>%
  mutate(length = st_distance(geometry, geometry2, by_element = TRUE)) ->
intersection_points

# Add Depth ranges:
intersection_points %>% 
   mutate(interval = paste(as.character(PROFUNDIDA),
                           "and",
                           as.character(lead(PROFUNDIDA))))
jmw
  • 443
  • 3
  • 9
  • 1
    Thank you very much, I am really grateful for your effort to answer me even with the difficulty to reproduce the data. The code worked perfectly, with the initial depth in output column I can know which range it refers to, this is no longer a problem. It just ignored the initial section of each linestring where there is obviously no intersection with the previous feature (of bathy). I didn't think of that when using ```st_intersection()``` I will measure the missed sections manually in QGIS. But, this already helps me a lot, otherwise I would have to manually measure each section per interval. – Caroline Portal Apr 27 '22 at 20:52
  • 1
    I tried many ways to reproduce the bathymetry object but i am little expertise in r. I confess that I did not know these packages, I will have to dedicate myself more to understand them for future questions here, thank you. – Caroline Portal Apr 27 '22 at 20:52
  • Glad you got it working! Looking at this problem led me to that problem where `st_cast` doesn't work normally for collections with both `MULTIPOINT` and `POINT`, so that's what the first part of my code is dealing with. I raised that issue here as well: https://github.com/r-spatial/sf/issues/1930 I think this issue should be re-opened. – jmw Apr 27 '22 at 21:21