R中到道路（空间线）的距离答案

【问题标题】：Distance to roads (spatial lines) in RR中到道路（空间线）的距离
【发布时间】：2021-09-25 14:02:37
【问题描述】：

几天来，我一直在努力生成一个表示到道路距离的栅格。我的研究区域中有一个道路空间线形状文件，我想以与其他协变量相同的分辨率（30 米像素）生成到道路要素的距离。在过去的几天里，我尝试了几件事，但最近我将模板栅格转换为空间点并使用了gDistance，但这似乎一直持续下去。
我的景观有点大，但不是很大（~25000 平方公里），其他空间分析不需要这么长时间，即使有 100 多个栅格。在尝试了我在网上找到的大多数建议后，我现在有点迷失了。

【问题讨论】：

你能提供行数据或类似格式的东西吗？是否可以投影到平面坐标，然后使用 spatstat 包？
当然，我可以提供线路数据。我没有使用 spatstat 包，但这似乎应该没问题。但是，我不确定在堆栈溢出时共享数据的最佳方式是什么。
在下面的答案中，我只是使用了 R 中 sf 包中的示例数据。如果您有很多行，读取行数据并转换为 spatstat 格式有时会非常缓慢。或许通过maptools中的as.linnet.SpatialLines可以加快速度，如果路网很详细，也许可以简化。

标签： r raster spatial

【解决方案1】：

用sf读入shapefile数据（注意自定义坐标参考用于投影到平面空间的系统 CRS 被发现使用 https://projectionwizard.org):

library(sf)
#> Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1
dat_path <- system.file("shape/olinda1.shp", package="sf")
dat <- st_read(dat_path)
#> Reading layer `olinda1' from data source `/usr/lib/R/site-library/sf/shape/olinda1.shp' using driver `ESRI Shapefile'
#> Simple feature collection with 470 features and 6 fields
#> Geometry type: POLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -34.91692 ymin: -8.044467 xmax: -34.82779 ymax: -7.954672
#> Geodetic CRS:  GRS 1980(IUGG, 1980)
poly <- st_geometry(dat)
new_crs <- st_crs("+proj=tcea +lon_0=-34.87 +datum=WGS84 +units=m +no_defs")
flat_poly <- st_transform(poly, new_crs)
flat_lines <- st_boundary(flat_poly)
plot(flat_lines)

转换为spatatstat 格式（注意spatstat 不做任何自动单位转换或命名，所以部分代码只是为了美观目的。坐标也是移位的，所以原点在左下角以获得更好的打印输出，但如果您曾经转换回其他格式，则应省略以避免在地图上放错区域。）

library(spatstat)
#> Loading required package: spatstat.data
#> Loading required package: spatstat.geom
#> spatstat.geom 2.2-2.002
#> Loading required package: spatstat.core
#> Loading required package: nlme
#> Loading required package: rpart
#> spatstat.core 2.3-0.003
#> Loading required package: spatstat.linnet
#> spatstat.linnet 2.3-0
#> 
#> spatstat 2.2-0       (nickname: 'That's not important right now') 
#> For an introduction to spatstat, type 'beginner'
flat_lines_psp <- as.psp(flat_lines)
flat_lines_psp
#> planar line segment pattern: 12235 line segments
#> window: rectangle = [-5172.766, 4653.619] x [-895506, -885510.1] units
L <- shift(flat_lines_psp, origin = "bottomleft")
unitname(L) <- "m"
L
#> planar line segment pattern: 12235 line segments
#> window: rectangle = [0, 9826.385] x [0, 9995.897] m

128x128 像素的默认 distmap：

system.time(D128pix <- distmap(L))
#>    user  system elapsed 
#>   0.527   0.000   0.527
D128pix
#> real-valued pixel image
#> 128 x 128 pixel array (ny, nx)
#> enclosing rectangle: [0, 9826.4] x [0, 9995.9] m
plot(D128pix, log = TRUE)

您可以使用dimyx 控制距离图的分辨率（产生的像素数）或eps（像素分辨率/大小）：

system.time(D512pix <- distmap(L, dimyx = 512))
#>    user  system elapsed 
#>   8.083   0.020   8.104
plot(D512pix, log = TRUE)

system.time(D20m <- distmap(L, eps = 20))
#>    user  system elapsed 
#>   7.704   0.000   7.704
plot(D20m, log = TRUE)

我认为代码在像素数上或多或少呈线性增长，并且线段的数量，所以如果你有很多小线段可能需要很长时间。至少在这个例子中减少了大约 2 倍的行数将计算时间减少了 2 倍：

Lsmall <- L[owin(c(2500,7500),c(2500,7500))]
Lsmall
#> planar line segment pattern: 5988 line segments
#> window: rectangle = [2500, 7500] x [2500, 7500] m
system.time(Dsmall512pix <- distmap(Lsmall, dimyx = 512))
#>    user  system elapsed 
#>   4.179   0.004   4.184
plot(Dsmall512pix, log = TRUE)

【讨论】：