网格 r 的每个单元内的河流长度

Question

我正在尝试测量网格的每个单元包含多少公里的河流。

我的河流 shapefile 是单行 shapefile（不是每条河流都有一行）。

我曾尝试修改this 和this 代码但失败了。

下面是两个对象的描述

> > grid5km Simple feature collection with 4075 features and 5 fields geometry type: MULTIPOLYGON dimension: XY bbox: xmin: 225320.3 ymin: 481277.6 xmax: 681165.7 ymax: 945385.3 epsg (SRID): 32629 proj4string: +proj=utm +zone=29 +datum=WGS84 +units=m +no_defs First 10 features: id xmin xmax ymin ymax geometry 0 28 -10.21953 -10.17431 8.50657 8.55179 MULTIPOLYGON (((369972.9 94... 1 29 -10.17431 -10.12909 8.50657 8.55179 MULTIPOLYGON (((370749.6 94...

> river Simple feature collection with 1 feature and 1 field geometry type: MULTILINESTRING dimension: XY bbox: xmin: 231870.6 ymin: 483505.6 xmax: 680763.9 ymax: 945087.4 epsg (SRID): 32629 proj4string: +proj=utm +zone=29 +datum=WGS84 +units=m +no_defs Strahler geometry 0 3 MULTILINESTRING ((232019.2 ...

Answer 1

您的代码不能完全重现，所以请允许我继续我的示例；请注意，它是建立在一组多边形（{sf} 包附带的标准北卡罗来纳 shapefile）和一个线串对象之上的。

你发现你的 multilinestring 对象不太好用；在这种情况下，请考虑 sf::st_line_merge() 将多线串转换为线串（我无法检查）。

关键方面是使用sf::st_intersection()，后跟sf::st_length()；结果将以您的 CRS 为单位 - 在这个例子中，我使用的是一个古朴的本地 CRS（让弗隆再次变得伟大......）

library(sf)
library(dplyr)

# included with sf package
shape <- st_read(system.file("shape/nc.shp", package="sf")) %>% 
  select(NAME) %>%  # just geometry, no data
  st_transform(crs = 6543)  # note: the result will be in US survey feet

# a slice of the 36th parallel
parallel36 <- st_linestring(matrix(c(-84, 36, -75, 36), 
                                   nrow = 2, byrow = T), 
                            dim = "XY") %>% 
  st_sfc(crs = 4326) %>% 
  st_transform(crs = 6543)  # again, a quaint local projection

# overview of result
plot(st_geometry(shape), col = "grey50")
plot(parallel36, add = T, col = "red")

# this is where the action happens ...
xsection <- st_intersection(shape, parallel36)  %>% # cross section
  mutate(length = st_length(.)) %>%  # calculate length
  st_drop_geometry() %>%  # geometry is no longer required
  units::drop_units() # drop units denomination

# to add the data back to the shape object
shape <- shape %>% 
  left_join(xsection, by = c("NAME"))

plot(shape["length"])