5

我试图找到由不规则多边形限制的两点之间的欧几里德距离。（即，必须将距离计算为通过给定窗口的路线）

这是一个可重现的示例：

library(spatstat)

#Simple example of a polygon and points.
ex.poly <- data.frame(x=c(0,5,5,2.5,0), y=c(0,0,5,2.5,5))
points <- data.frame(x=c(0.5, 2.5, 4.5), y=c(4,1,4))

bound <- owin(poly=data.frame(x=ex.poly$x, y=ex.poly$y))

test.ppp <- ppp(x=points$x, y=points$y, window=bound)

pairdist.ppp(test.ppp)#distance between every point
#The distance result from this function between point 1 and point 3, is given as 4.0

但是，我们仅通过绘制点就知道

plot(test.ppp)

路线被限制在多边形内时的距离应该更大（在本例中为 5.00）。

在 {spatstat} 中是否有另一个我不知道的功能可以做到这一点？或者有人对另一个可以做到这一点的包有任何其他建议吗？

我试图找出水体中两点之间的距离，所以我的实际数据中的不规则多边形更复杂。

任何帮助是极大的赞赏！

干杯

4

1 回答 1

7

好的，这是我昨天在评论中提到的基于gdistance的方法。它并不完美，因为它计算的路径段都被限制在棋盘上的 16 个方向之一（国王的移动加上骑士的移动）。也就是说，对于您的示例中的三个成对距离中的每一个，它都在正确值的 2% 以内（总是略微高估）。

library(maptools)  ## To convert spatstat objects to sp objects
library(gdistance) ## Loads raster and provides cost-surface functions

## Convert *.ppp points to SpatialPoints object
Pts <- as(test.ppp, "SpatialPoints")

## Convert the lake's boundary to a raster, with values of 1 for
## cells within the lake and values of 0 for cells on land
Poly <- as(bound, "SpatialPolygons")           ## 1st to SpatialPolygons-object
R <- raster(extent(Poly), nrow=100,  ncol=100) ## 2nd to RasterLayer ...
RR <- rasterize(Poly, R)                       ## ...
RR[is.na(RR)]<-0                               ## Set cells on land to "0"

## gdistance requires that you 1st prepare a sparse "transition matrix"
## whose values give the "conductance" of movement between pairs of
## adjacent and next-to-adjacent cells (when using directions=16)
tr1 <- transition(RR, transitionFunction=mean, directions=16)
tr1 <- geoCorrection(tr1,type="c")

## Compute a matrix of pairwise distances between points
## (These should be 5.00 and 3.605; all are within 2% of actual value).  
costDistance(tr1, Pts)
##          1        2
## 2 3.650282         
## 3 5.005259 3.650282

## View the selected paths
plot(RR)
plot(Pts, pch=16, col="gold", cex=1.5, add=TRUE)
SL12 <- shortestPath(tr1, Pts[1,], Pts[2,], output="SpatialLines")
SL13 <- shortestPath(tr1, Pts[1,], Pts[3,], output="SpatialLines")
SL23 <- shortestPath(tr1, Pts[2,], Pts[3,], output="SpatialLines")
lapply(list(SL12, SL13, SL23), function(X) plot(X, col="red", add=TRUE, lwd=2))

在此处输入图像描述

于 2014-04-09T22:58:31.150 回答