我计划今天回答我自己的问题，因为我刚刚设法使用与我的问题相关的gdistance 和sp 方法获取代码，但mrhellmann 击败了我！由于他们的答案也有效并且非常好和优雅，我想我只是在这里发布代码，利用这两种方法并比较感兴趣的人的结果（如果其中一个或另一个在你的上下文中不起作用） .它们的速度差不多，尽管在 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 方法的一些距离最终会更长一些。第一个是使用这两种方法的距离估计的重叠密度图，第二个是相互绘制的距离的“1 对 1”散点图，您可以看到拟合到 1 对 1 线会很不错的。
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)