我可以通过多边形绑定 st_distance 调用吗？
Can I bound an st_distance call by a polygon?
我看过关于此主题的类似帖子（例如，参见 here and here），但不是专门针对 sf-tidyverse 生态系统的帖子。
我有一系列湖泊，每个湖泊内有一系列采样点，每个湖泊中都有一个“焦点”，表示船只下水的位置。
我想计算从船下水到每个采样点的“划船者的最短行程距离”。但是，我想以某种方式使用湖泊多边形“绑定”这些距离，这样就无法计算陆地上的距离。我可以想象这是通过让“直线”蛇沿着湖边一直延伸到它可以恢复为直线之前所需的时间来完成的。我还可以想象“直线”被分解成围绕任何中间土地弯曲的线段。我对各种实现持开放态度！
我在其他地方（例如 here）看到了将边界多边形转换为栅格的想法，使水成为一个值，使土地成为另一个更高的值，然后进行“最小-成本距离”，其中越过陆地的成本高得令人望而却步。但是，我不知道在 raster/sf 生态系统中实际如何做到这一点。
这是我用来制作这张地图的代码：
library(dplyr)
library(sf)
library(ggplot2)
Moose.ssw = sswMN.sf %>% filter(lake == "Moose")
Moose.lake = MN_lakes4 %>% filter(str_detect(map_label, "Moose")) %>% filter(cty_name == "Beltrami")
Moose.access = try3 %>% filter(LAKE_NAME == "Moose") %>% filter(COUNTYNAME == "Beltrami")
Moose.box = st_bbox(Moose.ssw)
ggplot() +
geom_sf(data = Moose.lake, color="lightblue") +
geom_sf(data = Moose.access, color = "red", size = 2) +
geom_sf(data = Moose.ssw, mapping = aes(color= Nitellopsis_obtusa_n)) +
coord_sf(xlim = c(Moose.box[1], Moose.box[3]), ylim = c(Moose.box[2], Moose.box[4]))
这里有一些代码可以很好地计算直线距离——也许它可以以某种方式包装？
Moose.pt.dists = st_distance(Moose.access, Moose.ssw, by_element = TRUE)
制作上述数据对象所需的文件可以从我的 Github 页面中提取（它们是通过 dput().Link to my Github.
生成的文件
我是一名扎实的 R 程序员，但我是地理空间工作的新手，所以如果我什至可以被指出一个富有成果的方向，我 可能 能够编写我自己的程序摆脱这个！
这是一个使用 sfnetworks 的解决方案，它非常适合 tidyverse。
下面的代码应该
- 定期对湖泊的边界框进行采样（创建均匀分布的点）
- 将他们与最近的邻居联系起来
- 删除跨陆地的连接
- 沿着剩下的线找到从船下水到样本位置的最短路径。
结果不准确，但非常接近。您可以通过增加采样点的数量来提高精度。下面使用1000点。
library(tidyverse)
library(sf)
library(sfnetworks)
library(nngeo)
# set seed to make the script reproducible,
# since there is random sampling used
set.seed(813)
# Getting your data:
x <- dget("https://raw.githubusercontent.com/BajczA475/random-data/main/Moose.lake")
# Subset to get just one lake
moose_lake <- x[5,]
boat_ramp <- dget("https://raw.githubusercontent.com/BajczA475/random-data/main/Moose.access")
sample_locations <- dget("https://raw.githubusercontent.com/BajczA475/random-data/main/Moose.ssw")
sample_bbox <- dget("https://raw.githubusercontent.com/BajczA475/random-data/main/Moose.box")
# sample the bounding box with regular square points, then connect each point to the closest 9 points
# 8 should've worked, but left some diagonals out.
sq_grid_sample <- st_sample(st_as_sfc(st_bbox(moose_lake)), size = 1000, type = 'regular') %>% st_as_sf() %>%
st_connect(.,.,k = 9)
# remove connections that are not within the lake polygon
sq_grid_cropped <- sq_grid_sample[st_within(sq_grid_sample, moose_lake, sparse = F),]
# make an sfnetwork of the cropped grid
lake_network <- sq_grid_cropped %>% as_sfnetwork()
# find the (approximate) distance from boat ramp to point 170 (far north)
pt170 <- st_network_paths(lake_network,
from = boat_ramp,
to = sample_locations[170,]) %>%
pull(edge_paths) %>%
unlist()
lake_network %>%
activate(edges) %>%
slice(pt170) %>%
st_as_sf() %>%
st_combine() %>%
st_length()
#> 2186.394 [m]
看起来从船下水点到湖最北端的样本位置 170 大约有 2186 米。
# Plotting all of the above, path from boat ramp to point 170:
ggplot() +
geom_sf(data = sq_grid_sample, alpha = .05) +
geom_sf(data = sq_grid_cropped, color = 'dodgerblue') +
geom_sf(data = moose_lake, color = 'blue', fill = NA) +
geom_sf(data = boat_ramp, color = 'springgreen', size = 4) +
geom_sf(data = sample_locations[170,], color = 'deeppink1', size = 4) +
geom_sf(data = lake_network %>%
activate(edges) %>%
slice(pt170) %>%
st_as_sf(),
color = 'turquoise',
size = 2) +
theme_void()
虽然上面只找到并绘制了从船下水点到一个样本点的距离，但网络可以找到所有距离。您可能需要使用 *apply 或 purrr，并且可能需要一个小的自定义函数来查找发射到所有样本点的 'one-to-many' 距离。
This page on sfnetworks 将有助于编写一对多解决方案。
编辑：
查找从船下水点到样本点的所有距离：
st_network_cost(lake_network,
from = boat_ramp,
to = sample_locations)
这对我来说比 for 循环或下面发布的 sp 解决方案快得多。采样中的一些代码可能需要调整，因为 st_network_cost 将删除任何重复的距离。 sample_locations（或@bajcz 答案中的 Moose.ssw）也需要清理以删除重复点。 322行有180个唯一点。
我本来打算今天回答我自己的问题，因为我刚刚设法使用链接到我的问题中的 gdistance 和 sp 方法获得了代码，但 mrhellmann 打败了我它！由于他们的答案也有效并且非常漂亮和优雅，我想我只是 post 代码在这里利用这两种方法并比较那些感兴趣的人的结果（以防其中一个或另一个不起作用在你的上下文中）。它们差不多快，虽然在 mrhellmann 的 sfnetworks 版本中有一个 for() 循环，所以如果那个位可以被矢量化可能会更快，我相信这是可能的。
#All the hooplah needed to get things started and talking the same language.
x = dget("https://raw.githubusercontent.com/BajczA475/random-data/main/Moose.lake")
# Subset to get just one lake
Moose.lake = x[5,]
Moose.access = dget("https://raw.githubusercontent.com/BajczA475/random-data/main/Moose.access")
Moose.ssw = dget("https://raw.githubusercontent.com/BajczA475/random-data/main/Moose.ssw")
Moose.box = dget("https://raw.githubusercontent.com/BajczA475/random-data/main/Moose.box")
library(sf)
library(sp)
library(raster)
library(gdistance)
library(tidyverse)
library(sfnetworks)
library(nngeo)
library(mosaic)
这是使用 sp 和 gdistance 工具的“我的”方式：
ptm <- proc.time()
#Make all the objects needed into spatial class objects.
test.pts = as(Moose.access, Class = "Spatial")
ssw.pts = as(Moose.ssw, Class = "Spatial")
test.bounds = as(Moose.lake, Class = "Spatial")
#Turn the lake into a raster of resolution 1000 x 1000 (which is arbitrary) and then make all points not in the lake = 0 so that no distances can "afford" to cross over land.
test.raster = raster(extent(test.bounds), nrow=1000, ncol=1000)
lake.raster = rasterize(test.bounds, test.raster)
lake.raster[is.na(lake.raster)] = 0
##For some lakes, although not this one, the public access was just off the lake, which resulted in distances of Inf. This code puts a circular buffer around the dock locations to prevent this. It makes a new raster with 1s at the dock location(s), finds all indices of cells within the buffer distance of each dock (here, 300, which is overkill), and makes the corresponding cells in the lake raster also 1 so that they are considered navigable. This makes the distances slightly inaccurate because it may allow some points to be more easily reachable than they should be, but it is better than the alternative!
access.raster = rasterize(test.pts, lake.raster, field = 1)
index.spots = raster::extract(access.raster, test.pts, buffer=300, cellnumbers = T)
index.spots2 = unlist(lapply(index.spots, "[", , 1))
lake.raster[index.spots2] = 1
#Make a transition matrix so that the cost of moving from one cell to the next can be understood. TBH, I don't understand this part well beyond that.
transition.mat1 = transition(lake.raster, transitionFunction=mean, directions=16) #This code does take a little while.
transition.mat1 = geoCorrection(transition.mat1,type="c")
#Calculates the actual cost-based distances.
dists.test = costDistance(transition.mat1, test.pts, ssw.pts)
proc.time() - ptm #About 55 seconds on my laptop.
然后为了比较，这里是mrhellmann提供的sfnetworks方式
ptm <- proc.time()
sq_grid_sample = st_sample(st_as_sfc(st_bbox(Moose.lake)), size = 1000, type = 'regular') %>% st_as_sf() %>%
st_connect(.,.,k = 9)
sq_grid_cropped = sq_grid_sample[st_within(sq_grid_sample, Moose.lake, sparse = F)] #I needed to delete a comma in this line to get it to work--I don't think we needed row,column indexing here.
lake_network = sq_grid_cropped %>% as_sfnetwork()
##Using their code to generate all the distances between the dock and the sample points. This might be vectorizable.
dists.test2 = rep(0, nrow(Moose.ssw))
for (i in 1:nrow(Moose.ssw)) {
dist.pt = st_network_paths(lake_network,
from = Moose.access,
to = Moose.ssw[i,]) %>%
pull(edge_paths) %>%
unlist() #This produces a vertices we must go through I think?
dists.test2[i] = lake_network %>%
activate(edges) %>%
slice(dist.pt) %>%
st_as_sf() %>%
st_combine() %>%
st_length()
}
proc.time() - ptm #About 58 seconds on my laptop.
我认为这里有两张图可以证明这两种方法都产生了相对相似的数字，没有很多偏差迹象，尽管对于 sfnetworks 来说，某些距离看起来确实有点长方法。第一个是使用这两种方法的距离估计的重叠密度图，第二个是相互绘制的距离的“一对一”散点图，您可以看到适合一对一的线会很不错
hist.df = data.frame(dists = c(dists.test, dists.test2), version = rep(c("sp", "sfnetworks"), each = 322))
gf_density(~dists, fill=~version, data=hist.df )
gf_point(dists.test2~dists.test)
我知道这个问题已经得到解答，但我想我还是会在这里提出另一种方法。而且，看起来（至少对我而言）这种方法比其他方法快一点。这使用 terra，它是 raster 包的“替代品”；许多函数具有相同的名称并完成相同的工作，但主要区别在于 terra 使用 SpatRaster 个对象，而 raster 使用 raster 个对象。
# import data ---
x = dget("https://raw.githubusercontent.com/BajczA475/random-data/main/Moose.lake")
# Subset to get just one lake
Moose.lake = x[5,]
Moose.access = dget("https://raw.githubusercontent.com/BajczA475/random-data/main/Moose.access")
Moose.ssw = dget("https://raw.githubusercontent.com/BajczA475/random-data/main/Moose.ssw")
Moose.box = dget("https://raw.githubusercontent.com/BajczA475/random-data/main/Moose.box")
我们做的第一件事是将湖泊制作成栅格，terra 将其称为 SpatRaster。我们使用 Moose.lake 的范围 (terra::ext())，1000 行和 1000 列，并将 crs 设置为与 Moose.lake 相同。最初，我们给每个单元格一个 1 的值，但是然后我们可以使用 terra::mask() 给 Moose.lake 之外的每个值一个 2 的值（这些将是“高成本”或“禁区”单元格）。
# step 1 - rasterize() the lake ---
# find the spatial extent of the lake
ext(Moose.lake) %>%
# make a raster with the same extent, and assign all cells to value of 1
rast(nrow = 1000, ncol = 1000, crs = st_crs(Moose.lake)[1], vals = 1) %>%
# set all cells outside the lake a value of 2
mask(vect(Moose.lake), updatevalue = 2) %>%
{. ->> moose_rast}
moose_rast
# class : SpatRaster
# dimensions : 1000, 1000, 1 (nrow, ncol, nlyr)
# resolution : 2.899889, 3.239947 (x, y)
# extent : 388474, 391373.9, 5265182, 5268422 (xmin, xmax, ymin, ymax)
# coord. ref. : +proj=utm +zone=15 +datum=NAD83 +units=m +no_defs
# source : memory
# name : lyr.1
# min value : 1
# max value : 2
plot(moose_rast)
然后，我们给 Moose.access 值为 3（或 1 或 2 以外的任何值）的单元格 - 这将是我们的起点。
# now, make the boat ramp a value of 3 (this is the starting point)
# find the cell number based on the coordinates of Moose.access
moose_rast %>%
cellFromXY(st_coordinates(Moose.access)) %>%
{. ->> access_cell}
access_cell
# [1] 561618
# assign that cell a value of 3 (or any number you want other than 1 or 2)
values(moose_rast)[access_cell] <- 3
# check it
moose_rast %>% plot
尽管由于单元格太小而无法在图上真正看到它，但我们可以看出单元格值已更改，因为我们的图例现在包含值 3.
接下来，我们使用 terra::gridDistance() 计算从起始单元格（3 的值）到所有其他单元格的距离。我们指定 origin = 3 将此单元格指定为起点，并告诉它不要使用 omit = 2 遍历任何值为 2 的单元格。此函数 returns a SpatRaster，但这次单元格值是到原点的距离。
# now, find distances from the access_cell to every other cell
moose_rast %>%
gridDistance(origin = 3, omit = 2) %>%
{. ->> moose_rast_distances}
moose_rast_distances
# class : SpatRaster
# dimensions : 1000, 1000, 1 (nrow, ncol, nlyr)
# resolution : 2.899889, 3.239947 (x, y)
# extent : 388474, 391373.9, 5265182, 5268422 (xmin, xmax, ymin, ymax)
# coord. ref. : +proj=utm +zone=15 +datum=NAD83 +units=m +no_defs
# source : memory
# name : lyr.1
# min value : 0
# max value : 2424.096
# check it
moose_rast_distances %>% plot
在此图中，距离接入小区最近的区域为白色，距离最远的区域为绿色。
所以现在我们要做的就是找到湖中样本点的单元格编号并提取它们的单元格值。我们使用 terra::cellFromXY() 根据一组 XY 坐标查找细胞编号。
# now, find the cell numbers of all the study sites
moose_rast %>%
cellFromXY(st_coordinates(Moose.ssw)) %>%
{. ->> site_cells}
# cell numbers
site_cells %>%
head(50)
# [1] 953207 953233 930156 930182 930207 930233 907156 907182 907207 884078 884130
# [12] 884156 884182 884207 861052 861078 861104 861130 861156 861182 861207 837026
# [23] 837052 837078 837104 837130 837156 837182 837207 814026 814052 814078 814104
# [34] 814130 814156 814182 814207 791026 791104 791130 791156 791182 768182 745182
# [45] 745207 722026 722233 699259 675052 675285
# and now their distance values
values(moose_rast_distances)[site_cells] %>% head(50)
# [1] 2212.998 2241.812 2144.020 2115.206 2138.479 2167.293 2069.501 2040.688 2063.960
# [10] 2081.424 2023.796 1994.983 1966.169 1989.441 2078.719 2006.905 1978.092 1949.278
# [19] 1920.464 1891.650 1914.923 2119.358 2043.960 1968.563 1900.333 1871.519 1842.705
# [28] 1813.891 1837.164 2086.047 2010.650 1935.253 1859.856 1797.000 1768.186 1739.372
# [37] 1762.645 2052.736 1826.545 1751.148 1693.667 1664.854 1590.335 1533.733 1488.110
# [46] 1952.805 1384.778 1281.445 1809.338 1174.872
为了使距离更加用户友好，我们可以将它们放入 sf collection 中的新列中。
# now, make a new column in the study sites which is the distance to the access_cell
Moose.ssw %>%
rowwise %>%
mutate(
distance_to_access = cellFromXY(moose_rast, st_coordinates(geometry)) %>% values(moose_rast_distances)[.]
) %>%
select(distance_to_access, everything())
# Simple feature collection with 322 features and 122 fields
# Geometry type: POINT
# Dimension: XY
# Bounding box: xmin: 388549 ymin: 5265332 xmax: 391324 ymax: 5268332
# CRS: +proj=utm +zone=15 +datum=NAD83 +unit=m
# # A tibble: 322 x 123
# # Rowwise:
# distance_to_access lake county dow lake_size_acres contact year_discovered first_year_trea~ survey_year
# <dbl> <chr> <chr> <int> <dbl> <chr> <int> <chr> <int>
# 1 2213. Moose Beltrami 4001100 601. Nicole Kovar 2016 NoTreat 2017
# 2 2242. Moose Beltrami 4001100 601. Nicole Kovar 2016 NoTreat 2017
# 3 2144. Moose Beltrami 4001100 601. Nicole Kovar 2016 NoTreat 2017
# 4 2115. Moose Beltrami 4001100 601. Nicole Kovar 2016 NoTreat 2017
# 5 2138. Moose Beltrami 4001100 601. Nicole Kovar 2016 NoTreat 2017
# 6 2167. Moose Beltrami 4001100 601. Nicole Kovar 2016 NoTreat 2017
# 7 2070. Moose Beltrami 4001100 601. Nicole Kovar 2016 NoTreat 2017
# 8 2041. Moose Beltrami 4001100 601. Nicole Kovar 2016 NoTreat 2017
# 9 2064. Moose Beltrami 4001100 601. Nicole Kovar 2016 NoTreat 2017
# 10 2081. Moose Beltrami 4001100 601. Nicole Kovar 2016 NoTreat 2017
# # ... with 312 more rows, and 114 more variables: survey_dates <chr>, surveyor <chr>, pointid <int>, depth_ft <dbl>,
# # Myriophyllum_spicatum_or_hybrid_n <int>, Potamogeton_crispus_n <int>, Bidens_beckii_n <int>,
# # Brasenia_schreberi_n <int>, Ceratophyllum_demersum_n <int>, Ceratophyllum_echinatum_n <int>, Chara_sp_n <int>,
# # Eleocharis_acicularis_n <int>, Eleocharis_palustris_n <int>, Elodea_canadensis_n <int>,
# # Elodea_nuttallii_n <int>, Heteranthera_dubia_n <int>, Lemna_minor_n <int>, Lemna_trisulca_n <int>,
# # Myriophyllum_heterophyllum_n <int>, Myriophyllum_sibiricum_n <int>, Myriophyllum_verticillatum_n <int>,
# # Najas_flexilis_n <int>, Najas_gracillima_n <int>, Najas_guadalupensis_n <int>, Najas_marina_n <int>, ...
距离以米为单位，因为 CRS 是 UTM。
此外，gridDistance() 也可用于查找从湖内每个点到湖边的最短距离。为此，我们说 origin = 2，而不是像我们之前的船坡道那样只是一个点，而是陆地上的每个单元格。
moose_rast %>%
gridDistance(origin = 2) %>%
plot
所以我运行你上面回答中的两种方法，并将结果和时间与这种方法进行比较。
Time-wise，terra 方法要快得多，无论如何在我的电脑上：
sp ~230 秒sfnetworks ~83 秒terra ~5 秒
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# sp
#All the hooplah needed to get things started and talking the same language.
x = dget("https://raw.githubusercontent.com/BajczA475/random-data/main/Moose.lake")
# Subset to get just one lake
Moose.lake = x[5,]
Moose.access = dget("https://raw.githubusercontent.com/BajczA475/random-data/main/Moose.access")
Moose.ssw = dget("https://raw.githubusercontent.com/BajczA475/random-data/main/Moose.ssw")
Moose.box = dget("https://raw.githubusercontent.com/BajczA475/random-data/main/Moose.box")
library(sf)
library(sp)
library(raster)
library(gdistance)
library(tidyverse)
library(sfnetworks)
library(nngeo)
library(mosaic)
ptm <- proc.time()
#Make all the objects needed into spatial class objects.
test.pts = as(Moose.access, Class = "Spatial")
ssw.pts = as(Moose.ssw, Class = "Spatial")
test.bounds = as(Moose.lake, Class = "Spatial")
#Turn the lake into a raster of resolution 1000 x 1000 (which is arbitrary) and then make all points not in the lake = 0 so that no distances can "afford" to cross over land.
test.raster = raster(extent(test.bounds), nrow=1000, ncol=1000)
lake.raster = rasterize(test.bounds, test.raster)
lake.raster[is.na(lake.raster)] = 0
##For some lakes, although not this one, the public access was just off the lake, which resulted in distances of Inf. This code puts a circular buffer around the dock locations to prevent this. It makes a new raster with 1s at the dock location(s), finds all indices of cells within the buffer distance of each dock (here, 300, which is overkill), and makes the corresponding cells in the lake raster also 1 so that they are considered navigable. This makes the distances slightly inaccurate because it may allow some points to be more easily reachable than they should be, but it is better than the alternative!
access.raster = rasterize(test.pts, lake.raster, field = 1)
index.spots = raster::extract(access.raster, test.pts, buffer=300, cellnumbers = T)
index.spots2 = unlist(lapply(index.spots, "[", , 1))
lake.raster[index.spots2] = 1
#Make a transition matrix so that the cost of moving from one cell to the next can be understood. TBH, I don't understand this part well beyond that.
transition.mat1 = transition(lake.raster, transitionFunction=mean, directions=16) #This code does take a little while.
transition.mat1 = geoCorrection(transition.mat1,type="c")
#Calculates the actual cost-based distances.
dists.test = costDistance(transition.mat1, test.pts, ssw.pts)
proc.time() - ptm
# 234 secs
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# sfnetworks
#All the hooplah needed to get things started and talking the same language.
x = dget("https://raw.githubusercontent.com/BajczA475/random-data/main/Moose.lake")
# Subset to get just one lake
Moose.lake = x[5,]
Moose.access = dget("https://raw.githubusercontent.com/BajczA475/random-data/main/Moose.access")
Moose.ssw = dget("https://raw.githubusercontent.com/BajczA475/random-data/main/Moose.ssw")
Moose.box = dget("https://raw.githubusercontent.com/BajczA475/random-data/main/Moose.box")
library(sf)
library(sp)
library(raster)
library(gdistance)
library(tidyverse)
library(sfnetworks)
library(nngeo)
library(mosaic)
ptm <- proc.time()
sq_grid_sample = st_sample(st_as_sfc(st_bbox(Moose.lake)), size = 1000, type = 'regular') %>% st_as_sf() %>%
st_connect(.,.,k = 9)
sq_grid_cropped = sq_grid_sample[st_within(sq_grid_sample, Moose.lake, sparse = F)] #I needed to delete a comma in this line to get it to work--I don't think we needed row,column indexing here.
lake_network = sq_grid_cropped %>% as_sfnetwork()
##Using their code to generate all the distances between the dock and the sample points. This might be vectorizable.
dists.test2 = rep(0, nrow(Moose.ssw))
for (i in 1:nrow(Moose.ssw)) {
dist.pt = st_network_paths(lake_network,
from = Moose.access,
to = Moose.ssw[i,]) %>%
pull(edge_paths) %>%
unlist() #This produces a vertices we must go through I think?
dists.test2[i] = lake_network %>%
activate(edges) %>%
slice(dist.pt) %>%
st_as_sf() %>%
st_combine() %>%
st_length()
}
proc.time() - ptm
# 83 secs
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# terra
#All the hooplah needed to get things started and talking the same language.
x = dget("https://raw.githubusercontent.com/BajczA475/random-data/main/Moose.lake")
# Subset to get just one lake
Moose.lake = x[5,]
Moose.access = dget("https://raw.githubusercontent.com/BajczA475/random-data/main/Moose.access")
Moose.ssw = dget("https://raw.githubusercontent.com/BajczA475/random-data/main/Moose.ssw")
Moose.box = dget("https://raw.githubusercontent.com/BajczA475/random-data/main/Moose.box")
library(tidyverse)
library(sf)
library(terra)
ptm <- proc.time()
# make rastre of moose lake
ext(Moose.lake) %>%
# make a raster with the same extent, and assign all cells to value of 1
rast(nrow = 1000, ncol = 1000, crs = st_crs(Moose.lake)[1], vals = 1) %>%
# set all cells outside the lake a value of 2
mask(vect(Moose.lake), updatevalue = 2) %>%
{. ->> moose_rast}
# find access cell
moose_rast %>%
cellFromXY(st_coordinates(Moose.access)) %>%
{. ->> access_cell}
# assign that cell a value of 3 (or any number you want other than 1 or 2)
values(moose_rast)[access_cell] <- 3
# find values to every other cell
moose_rast %>%
gridDistance(origin = 3, omit = 2) %>%
{. ->> moose_rast_distances}
# make column with distances to each sample site
distances <- Moose.ssw %>%
rowwise %>%
mutate(
distance_to_access = cellFromXY(moose_rast, st_coordinates(geometry)) %>% values(moose_rast_distances)[.]
) %>%
select(distance_to_access, everything())
proc.time() - ptm
# 5 secs
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
hist.df2 = data.frame(dists = c(dists.test, dists.test2, distances$distance_to_access), version = rep(c("sp", "sfnetworks", "terra"), each = 322))
gf_density(~dists, fill=~version, data=hist.df2)
并且 terra 方法的距离值与 sp 方法非常相似：
我看过关于此主题的类似帖子（例如，参见 here and here），但不是专门针对 sf-tidyverse 生态系统的帖子。
我有一系列湖泊，每个湖泊内有一系列采样点，每个湖泊中都有一个“焦点”，表示船只下水的位置。
我在其他地方（例如 here）看到了将边界多边形转换为栅格的想法，使水成为一个值，使土地成为另一个更高的值，然后进行“最小-成本距离”，其中越过陆地的成本高得令人望而却步。但是，我不知道在 raster/sf 生态系统中实际如何做到这一点。
这是我用来制作这张地图的代码：
library(dplyr)
library(sf)
library(ggplot2)
Moose.ssw = sswMN.sf %>% filter(lake == "Moose")
Moose.lake = MN_lakes4 %>% filter(str_detect(map_label, "Moose")) %>% filter(cty_name == "Beltrami")
Moose.access = try3 %>% filter(LAKE_NAME == "Moose") %>% filter(COUNTYNAME == "Beltrami")
Moose.box = st_bbox(Moose.ssw)
ggplot() +
geom_sf(data = Moose.lake, color="lightblue") +
geom_sf(data = Moose.access, color = "red", size = 2) +
geom_sf(data = Moose.ssw, mapping = aes(color= Nitellopsis_obtusa_n)) +
coord_sf(xlim = c(Moose.box[1], Moose.box[3]), ylim = c(Moose.box[2], Moose.box[4]))
这里有一些代码可以很好地计算直线距离——也许它可以以某种方式包装？
Moose.pt.dists = st_distance(Moose.access, Moose.ssw, by_element = TRUE)
制作上述数据对象所需的文件可以从我的 Github 页面中提取（它们是通过 dput().Link to my Github.
我是一名扎实的 R 程序员，但我是地理空间工作的新手，所以如果我什至可以被指出一个富有成果的方向，我 可能 能够编写我自己的程序摆脱这个！
这是一个使用 sfnetworks 的解决方案，它非常适合 tidyverse。
下面的代码应该
- 定期对湖泊的边界框进行采样（创建均匀分布的点）
- 将他们与最近的邻居联系起来
- 删除跨陆地的连接
- 沿着剩下的线找到从船下水到样本位置的最短路径。
结果不准确，但非常接近。您可以通过增加采样点的数量来提高精度。下面使用1000点。
library(tidyverse)
library(sf)
library(sfnetworks)
library(nngeo)
# set seed to make the script reproducible,
# since there is random sampling used
set.seed(813)
# Getting your data:
x <- dget("https://raw.githubusercontent.com/BajczA475/random-data/main/Moose.lake")
# Subset to get just one lake
moose_lake <- x[5,]
boat_ramp <- dget("https://raw.githubusercontent.com/BajczA475/random-data/main/Moose.access")
sample_locations <- dget("https://raw.githubusercontent.com/BajczA475/random-data/main/Moose.ssw")
sample_bbox <- dget("https://raw.githubusercontent.com/BajczA475/random-data/main/Moose.box")
# sample the bounding box with regular square points, then connect each point to the closest 9 points
# 8 should've worked, but left some diagonals out.
sq_grid_sample <- st_sample(st_as_sfc(st_bbox(moose_lake)), size = 1000, type = 'regular') %>% st_as_sf() %>%
st_connect(.,.,k = 9)
# remove connections that are not within the lake polygon
sq_grid_cropped <- sq_grid_sample[st_within(sq_grid_sample, moose_lake, sparse = F),]
# make an sfnetwork of the cropped grid
lake_network <- sq_grid_cropped %>% as_sfnetwork()
# find the (approximate) distance from boat ramp to point 170 (far north)
pt170 <- st_network_paths(lake_network,
from = boat_ramp,
to = sample_locations[170,]) %>%
pull(edge_paths) %>%
unlist()
lake_network %>%
activate(edges) %>%
slice(pt170) %>%
st_as_sf() %>%
st_combine() %>%
st_length()
#> 2186.394 [m]
看起来从船下水点到湖最北端的样本位置 170 大约有 2186 米。
# Plotting all of the above, path from boat ramp to point 170:
ggplot() +
geom_sf(data = sq_grid_sample, alpha = .05) +
geom_sf(data = sq_grid_cropped, color = 'dodgerblue') +
geom_sf(data = moose_lake, color = 'blue', fill = NA) +
geom_sf(data = boat_ramp, color = 'springgreen', size = 4) +
geom_sf(data = sample_locations[170,], color = 'deeppink1', size = 4) +
geom_sf(data = lake_network %>%
activate(edges) %>%
slice(pt170) %>%
st_as_sf(),
color = 'turquoise',
size = 2) +
theme_void()
虽然上面只找到并绘制了从船下水点到一个样本点的距离，但网络可以找到所有距离。您可能需要使用 *apply 或 purrr，并且可能需要一个小的自定义函数来查找发射到所有样本点的 'one-to-many' 距离。
This page on sfnetworks 将有助于编写一对多解决方案。
编辑：
查找从船下水点到样本点的所有距离：
st_network_cost(lake_network,
from = boat_ramp,
to = sample_locations)
这对我来说比 for 循环或下面发布的 sp 解决方案快得多。采样中的一些代码可能需要调整，因为 st_network_cost 将删除任何重复的距离。 sample_locations（或@bajcz 答案中的 Moose.ssw）也需要清理以删除重复点。 322行有180个唯一点。
我本来打算今天回答我自己的问题，因为我刚刚设法使用链接到我的问题中的 gdistance 和 sp 方法获得了代码，但 mrhellmann 打败了我它！由于他们的答案也有效并且非常漂亮和优雅，我想我只是 post 代码在这里利用这两种方法并比较那些感兴趣的人的结果（以防其中一个或另一个不起作用在你的上下文中）。它们差不多快，虽然在 mrhellmann 的 sfnetworks 版本中有一个 for() 循环，所以如果那个位可以被矢量化可能会更快，我相信这是可能的。
#All the hooplah needed to get things started and talking the same language.
x = dget("https://raw.githubusercontent.com/BajczA475/random-data/main/Moose.lake")
# Subset to get just one lake
Moose.lake = x[5,]
Moose.access = dget("https://raw.githubusercontent.com/BajczA475/random-data/main/Moose.access")
Moose.ssw = dget("https://raw.githubusercontent.com/BajczA475/random-data/main/Moose.ssw")
Moose.box = dget("https://raw.githubusercontent.com/BajczA475/random-data/main/Moose.box")
library(sf)
library(sp)
library(raster)
library(gdistance)
library(tidyverse)
library(sfnetworks)
library(nngeo)
library(mosaic)
这是使用 sp 和 gdistance 工具的“我的”方式：
ptm <- proc.time()
#Make all the objects needed into spatial class objects.
test.pts = as(Moose.access, Class = "Spatial")
ssw.pts = as(Moose.ssw, Class = "Spatial")
test.bounds = as(Moose.lake, Class = "Spatial")
#Turn the lake into a raster of resolution 1000 x 1000 (which is arbitrary) and then make all points not in the lake = 0 so that no distances can "afford" to cross over land.
test.raster = raster(extent(test.bounds), nrow=1000, ncol=1000)
lake.raster = rasterize(test.bounds, test.raster)
lake.raster[is.na(lake.raster)] = 0
##For some lakes, although not this one, the public access was just off the lake, which resulted in distances of Inf. This code puts a circular buffer around the dock locations to prevent this. It makes a new raster with 1s at the dock location(s), finds all indices of cells within the buffer distance of each dock (here, 300, which is overkill), and makes the corresponding cells in the lake raster also 1 so that they are considered navigable. This makes the distances slightly inaccurate because it may allow some points to be more easily reachable than they should be, but it is better than the alternative!
access.raster = rasterize(test.pts, lake.raster, field = 1)
index.spots = raster::extract(access.raster, test.pts, buffer=300, cellnumbers = T)
index.spots2 = unlist(lapply(index.spots, "[", , 1))
lake.raster[index.spots2] = 1
#Make a transition matrix so that the cost of moving from one cell to the next can be understood. TBH, I don't understand this part well beyond that.
transition.mat1 = transition(lake.raster, transitionFunction=mean, directions=16) #This code does take a little while.
transition.mat1 = geoCorrection(transition.mat1,type="c")
#Calculates the actual cost-based distances.
dists.test = costDistance(transition.mat1, test.pts, ssw.pts)
proc.time() - ptm #About 55 seconds on my laptop.
然后为了比较，这里是mrhellmann提供的sfnetworks方式
ptm <- proc.time()
sq_grid_sample = st_sample(st_as_sfc(st_bbox(Moose.lake)), size = 1000, type = 'regular') %>% st_as_sf() %>%
st_connect(.,.,k = 9)
sq_grid_cropped = sq_grid_sample[st_within(sq_grid_sample, Moose.lake, sparse = F)] #I needed to delete a comma in this line to get it to work--I don't think we needed row,column indexing here.
lake_network = sq_grid_cropped %>% as_sfnetwork()
##Using their code to generate all the distances between the dock and the sample points. This might be vectorizable.
dists.test2 = rep(0, nrow(Moose.ssw))
for (i in 1:nrow(Moose.ssw)) {
dist.pt = st_network_paths(lake_network,
from = Moose.access,
to = Moose.ssw[i,]) %>%
pull(edge_paths) %>%
unlist() #This produces a vertices we must go through I think?
dists.test2[i] = lake_network %>%
activate(edges) %>%
slice(dist.pt) %>%
st_as_sf() %>%
st_combine() %>%
st_length()
}
proc.time() - ptm #About 58 seconds on my laptop.
我认为这里有两张图可以证明这两种方法都产生了相对相似的数字，没有很多偏差迹象，尽管对于 sfnetworks 来说，某些距离看起来确实有点长方法。第一个是使用这两种方法的距离估计的重叠密度图，第二个是相互绘制的距离的“一对一”散点图，您可以看到适合一对一的线会很不错
hist.df = data.frame(dists = c(dists.test, dists.test2), version = rep(c("sp", "sfnetworks"), each = 322))
gf_density(~dists, fill=~version, data=hist.df )
gf_point(dists.test2~dists.test)
我知道这个问题已经得到解答，但我想我还是会在这里提出另一种方法。而且，看起来（至少对我而言）这种方法比其他方法快一点。这使用 terra，它是 raster 包的“替代品”；许多函数具有相同的名称并完成相同的工作，但主要区别在于 terra 使用 SpatRaster 个对象，而 raster 使用 raster 个对象。
# import data ---
x = dget("https://raw.githubusercontent.com/BajczA475/random-data/main/Moose.lake")
# Subset to get just one lake
Moose.lake = x[5,]
Moose.access = dget("https://raw.githubusercontent.com/BajczA475/random-data/main/Moose.access")
Moose.ssw = dget("https://raw.githubusercontent.com/BajczA475/random-data/main/Moose.ssw")
Moose.box = dget("https://raw.githubusercontent.com/BajczA475/random-data/main/Moose.box")
我们做的第一件事是将湖泊制作成栅格，terra 将其称为 SpatRaster。我们使用 Moose.lake 的范围 (terra::ext())，1000 行和 1000 列，并将 crs 设置为与 Moose.lake 相同。最初，我们给每个单元格一个 1 的值，但是然后我们可以使用 terra::mask() 给 Moose.lake 之外的每个值一个 2 的值（这些将是“高成本”或“禁区”单元格）。
# step 1 - rasterize() the lake ---
# find the spatial extent of the lake
ext(Moose.lake) %>%
# make a raster with the same extent, and assign all cells to value of 1
rast(nrow = 1000, ncol = 1000, crs = st_crs(Moose.lake)[1], vals = 1) %>%
# set all cells outside the lake a value of 2
mask(vect(Moose.lake), updatevalue = 2) %>%
{. ->> moose_rast}
moose_rast
# class : SpatRaster
# dimensions : 1000, 1000, 1 (nrow, ncol, nlyr)
# resolution : 2.899889, 3.239947 (x, y)
# extent : 388474, 391373.9, 5265182, 5268422 (xmin, xmax, ymin, ymax)
# coord. ref. : +proj=utm +zone=15 +datum=NAD83 +units=m +no_defs
# source : memory
# name : lyr.1
# min value : 1
# max value : 2
plot(moose_rast)
然后，我们给 Moose.access 值为 3（或 1 或 2 以外的任何值）的单元格 - 这将是我们的起点。
# now, make the boat ramp a value of 3 (this is the starting point)
# find the cell number based on the coordinates of Moose.access
moose_rast %>%
cellFromXY(st_coordinates(Moose.access)) %>%
{. ->> access_cell}
access_cell
# [1] 561618
# assign that cell a value of 3 (or any number you want other than 1 or 2)
values(moose_rast)[access_cell] <- 3
# check it
moose_rast %>% plot
尽管由于单元格太小而无法在图上真正看到它，但我们可以看出单元格值已更改，因为我们的图例现在包含值 3.
接下来，我们使用 terra::gridDistance() 计算从起始单元格（3 的值）到所有其他单元格的距离。我们指定 origin = 3 将此单元格指定为起点，并告诉它不要使用 omit = 2 遍历任何值为 2 的单元格。此函数 returns a SpatRaster，但这次单元格值是到原点的距离。
# now, find distances from the access_cell to every other cell
moose_rast %>%
gridDistance(origin = 3, omit = 2) %>%
{. ->> moose_rast_distances}
moose_rast_distances
# class : SpatRaster
# dimensions : 1000, 1000, 1 (nrow, ncol, nlyr)
# resolution : 2.899889, 3.239947 (x, y)
# extent : 388474, 391373.9, 5265182, 5268422 (xmin, xmax, ymin, ymax)
# coord. ref. : +proj=utm +zone=15 +datum=NAD83 +units=m +no_defs
# source : memory
# name : lyr.1
# min value : 0
# max value : 2424.096
# check it
moose_rast_distances %>% plot
在此图中，距离接入小区最近的区域为白色，距离最远的区域为绿色。
所以现在我们要做的就是找到湖中样本点的单元格编号并提取它们的单元格值。我们使用 terra::cellFromXY() 根据一组 XY 坐标查找细胞编号。
# now, find the cell numbers of all the study sites
moose_rast %>%
cellFromXY(st_coordinates(Moose.ssw)) %>%
{. ->> site_cells}
# cell numbers
site_cells %>%
head(50)
# [1] 953207 953233 930156 930182 930207 930233 907156 907182 907207 884078 884130
# [12] 884156 884182 884207 861052 861078 861104 861130 861156 861182 861207 837026
# [23] 837052 837078 837104 837130 837156 837182 837207 814026 814052 814078 814104
# [34] 814130 814156 814182 814207 791026 791104 791130 791156 791182 768182 745182
# [45] 745207 722026 722233 699259 675052 675285
# and now their distance values
values(moose_rast_distances)[site_cells] %>% head(50)
# [1] 2212.998 2241.812 2144.020 2115.206 2138.479 2167.293 2069.501 2040.688 2063.960
# [10] 2081.424 2023.796 1994.983 1966.169 1989.441 2078.719 2006.905 1978.092 1949.278
# [19] 1920.464 1891.650 1914.923 2119.358 2043.960 1968.563 1900.333 1871.519 1842.705
# [28] 1813.891 1837.164 2086.047 2010.650 1935.253 1859.856 1797.000 1768.186 1739.372
# [37] 1762.645 2052.736 1826.545 1751.148 1693.667 1664.854 1590.335 1533.733 1488.110
# [46] 1952.805 1384.778 1281.445 1809.338 1174.872
为了使距离更加用户友好，我们可以将它们放入 sf collection 中的新列中。
# now, make a new column in the study sites which is the distance to the access_cell
Moose.ssw %>%
rowwise %>%
mutate(
distance_to_access = cellFromXY(moose_rast, st_coordinates(geometry)) %>% values(moose_rast_distances)[.]
) %>%
select(distance_to_access, everything())
# Simple feature collection with 322 features and 122 fields
# Geometry type: POINT
# Dimension: XY
# Bounding box: xmin: 388549 ymin: 5265332 xmax: 391324 ymax: 5268332
# CRS: +proj=utm +zone=15 +datum=NAD83 +unit=m
# # A tibble: 322 x 123
# # Rowwise:
# distance_to_access lake county dow lake_size_acres contact year_discovered first_year_trea~ survey_year
# <dbl> <chr> <chr> <int> <dbl> <chr> <int> <chr> <int>
# 1 2213. Moose Beltrami 4001100 601. Nicole Kovar 2016 NoTreat 2017
# 2 2242. Moose Beltrami 4001100 601. Nicole Kovar 2016 NoTreat 2017
# 3 2144. Moose Beltrami 4001100 601. Nicole Kovar 2016 NoTreat 2017
# 4 2115. Moose Beltrami 4001100 601. Nicole Kovar 2016 NoTreat 2017
# 5 2138. Moose Beltrami 4001100 601. Nicole Kovar 2016 NoTreat 2017
# 6 2167. Moose Beltrami 4001100 601. Nicole Kovar 2016 NoTreat 2017
# 7 2070. Moose Beltrami 4001100 601. Nicole Kovar 2016 NoTreat 2017
# 8 2041. Moose Beltrami 4001100 601. Nicole Kovar 2016 NoTreat 2017
# 9 2064. Moose Beltrami 4001100 601. Nicole Kovar 2016 NoTreat 2017
# 10 2081. Moose Beltrami 4001100 601. Nicole Kovar 2016 NoTreat 2017
# # ... with 312 more rows, and 114 more variables: survey_dates <chr>, surveyor <chr>, pointid <int>, depth_ft <dbl>,
# # Myriophyllum_spicatum_or_hybrid_n <int>, Potamogeton_crispus_n <int>, Bidens_beckii_n <int>,
# # Brasenia_schreberi_n <int>, Ceratophyllum_demersum_n <int>, Ceratophyllum_echinatum_n <int>, Chara_sp_n <int>,
# # Eleocharis_acicularis_n <int>, Eleocharis_palustris_n <int>, Elodea_canadensis_n <int>,
# # Elodea_nuttallii_n <int>, Heteranthera_dubia_n <int>, Lemna_minor_n <int>, Lemna_trisulca_n <int>,
# # Myriophyllum_heterophyllum_n <int>, Myriophyllum_sibiricum_n <int>, Myriophyllum_verticillatum_n <int>,
# # Najas_flexilis_n <int>, Najas_gracillima_n <int>, Najas_guadalupensis_n <int>, Najas_marina_n <int>, ...
距离以米为单位，因为 CRS 是 UTM。
此外，gridDistance() 也可用于查找从湖内每个点到湖边的最短距离。为此，我们说 origin = 2，而不是像我们之前的船坡道那样只是一个点，而是陆地上的每个单元格。
moose_rast %>%
gridDistance(origin = 2) %>%
plot
所以我运行你上面回答中的两种方法，并将结果和时间与这种方法进行比较。
Time-wise，terra 方法要快得多，无论如何在我的电脑上：
sp~230 秒sfnetworks~83 秒terra~5 秒
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# sp
#All the hooplah needed to get things started and talking the same language.
x = dget("https://raw.githubusercontent.com/BajczA475/random-data/main/Moose.lake")
# Subset to get just one lake
Moose.lake = x[5,]
Moose.access = dget("https://raw.githubusercontent.com/BajczA475/random-data/main/Moose.access")
Moose.ssw = dget("https://raw.githubusercontent.com/BajczA475/random-data/main/Moose.ssw")
Moose.box = dget("https://raw.githubusercontent.com/BajczA475/random-data/main/Moose.box")
library(sf)
library(sp)
library(raster)
library(gdistance)
library(tidyverse)
library(sfnetworks)
library(nngeo)
library(mosaic)
ptm <- proc.time()
#Make all the objects needed into spatial class objects.
test.pts = as(Moose.access, Class = "Spatial")
ssw.pts = as(Moose.ssw, Class = "Spatial")
test.bounds = as(Moose.lake, Class = "Spatial")
#Turn the lake into a raster of resolution 1000 x 1000 (which is arbitrary) and then make all points not in the lake = 0 so that no distances can "afford" to cross over land.
test.raster = raster(extent(test.bounds), nrow=1000, ncol=1000)
lake.raster = rasterize(test.bounds, test.raster)
lake.raster[is.na(lake.raster)] = 0
##For some lakes, although not this one, the public access was just off the lake, which resulted in distances of Inf. This code puts a circular buffer around the dock locations to prevent this. It makes a new raster with 1s at the dock location(s), finds all indices of cells within the buffer distance of each dock (here, 300, which is overkill), and makes the corresponding cells in the lake raster also 1 so that they are considered navigable. This makes the distances slightly inaccurate because it may allow some points to be more easily reachable than they should be, but it is better than the alternative!
access.raster = rasterize(test.pts, lake.raster, field = 1)
index.spots = raster::extract(access.raster, test.pts, buffer=300, cellnumbers = T)
index.spots2 = unlist(lapply(index.spots, "[", , 1))
lake.raster[index.spots2] = 1
#Make a transition matrix so that the cost of moving from one cell to the next can be understood. TBH, I don't understand this part well beyond that.
transition.mat1 = transition(lake.raster, transitionFunction=mean, directions=16) #This code does take a little while.
transition.mat1 = geoCorrection(transition.mat1,type="c")
#Calculates the actual cost-based distances.
dists.test = costDistance(transition.mat1, test.pts, ssw.pts)
proc.time() - ptm
# 234 secs
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# sfnetworks
#All the hooplah needed to get things started and talking the same language.
x = dget("https://raw.githubusercontent.com/BajczA475/random-data/main/Moose.lake")
# Subset to get just one lake
Moose.lake = x[5,]
Moose.access = dget("https://raw.githubusercontent.com/BajczA475/random-data/main/Moose.access")
Moose.ssw = dget("https://raw.githubusercontent.com/BajczA475/random-data/main/Moose.ssw")
Moose.box = dget("https://raw.githubusercontent.com/BajczA475/random-data/main/Moose.box")
library(sf)
library(sp)
library(raster)
library(gdistance)
library(tidyverse)
library(sfnetworks)
library(nngeo)
library(mosaic)
ptm <- proc.time()
sq_grid_sample = st_sample(st_as_sfc(st_bbox(Moose.lake)), size = 1000, type = 'regular') %>% st_as_sf() %>%
st_connect(.,.,k = 9)
sq_grid_cropped = sq_grid_sample[st_within(sq_grid_sample, Moose.lake, sparse = F)] #I needed to delete a comma in this line to get it to work--I don't think we needed row,column indexing here.
lake_network = sq_grid_cropped %>% as_sfnetwork()
##Using their code to generate all the distances between the dock and the sample points. This might be vectorizable.
dists.test2 = rep(0, nrow(Moose.ssw))
for (i in 1:nrow(Moose.ssw)) {
dist.pt = st_network_paths(lake_network,
from = Moose.access,
to = Moose.ssw[i,]) %>%
pull(edge_paths) %>%
unlist() #This produces a vertices we must go through I think?
dists.test2[i] = lake_network %>%
activate(edges) %>%
slice(dist.pt) %>%
st_as_sf() %>%
st_combine() %>%
st_length()
}
proc.time() - ptm
# 83 secs
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# terra
#All the hooplah needed to get things started and talking the same language.
x = dget("https://raw.githubusercontent.com/BajczA475/random-data/main/Moose.lake")
# Subset to get just one lake
Moose.lake = x[5,]
Moose.access = dget("https://raw.githubusercontent.com/BajczA475/random-data/main/Moose.access")
Moose.ssw = dget("https://raw.githubusercontent.com/BajczA475/random-data/main/Moose.ssw")
Moose.box = dget("https://raw.githubusercontent.com/BajczA475/random-data/main/Moose.box")
library(tidyverse)
library(sf)
library(terra)
ptm <- proc.time()
# make rastre of moose lake
ext(Moose.lake) %>%
# make a raster with the same extent, and assign all cells to value of 1
rast(nrow = 1000, ncol = 1000, crs = st_crs(Moose.lake)[1], vals = 1) %>%
# set all cells outside the lake a value of 2
mask(vect(Moose.lake), updatevalue = 2) %>%
{. ->> moose_rast}
# find access cell
moose_rast %>%
cellFromXY(st_coordinates(Moose.access)) %>%
{. ->> access_cell}
# assign that cell a value of 3 (or any number you want other than 1 or 2)
values(moose_rast)[access_cell] <- 3
# find values to every other cell
moose_rast %>%
gridDistance(origin = 3, omit = 2) %>%
{. ->> moose_rast_distances}
# make column with distances to each sample site
distances <- Moose.ssw %>%
rowwise %>%
mutate(
distance_to_access = cellFromXY(moose_rast, st_coordinates(geometry)) %>% values(moose_rast_distances)[.]
) %>%
select(distance_to_access, everything())
proc.time() - ptm
# 5 secs
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
hist.df2 = data.frame(dists = c(dists.test, dists.test2, distances$distance_to_access), version = rep(c("sp", "sfnetworks", "terra"), each = 322))
gf_density(~dists, fill=~version, data=hist.df2)
并且 terra 方法的距离值与 sp 方法非常相似：