Finding euclidean distance in R{spatstat} between points, confined by an irregular polygon window

Question

I'm trying to find the euclidean distance between two points, confined by an irregular polygon. (ie. the distance would have to be calculated as a route through the window given)

Here is an reproducible example:

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

However we know just from plotting the points

plot(test.ppp)

that the distance when the route is confined to the polygon should be greater (in this case, 5.00).

Is there another function that I am not aware of in {spatstat} that would do this? Or does anybody have any other suggestions for another package that could do this?

I'm trying to find the distance between two points in a water body, so the irregular polygon in my actual data is more complex.

Any help is greatly appreciated!

Cheers

Mazonson

Interesting question. I might suggest converting bound to a RasterLayer object (perhaps first using maptools to supply as(bound, "SpatialPolygons") and as(test.ppp, "SpatialPoints")), and then using the gdistance package to compute a "least-cost distance" between the points, with the friction or cost of all grid points outside of bound set to infinity. gdistance ships with a nice vignette (run vignette("gdistance") to see it) that should give you a good start. – Ascot

Answer 1

OK, here's the gdistance-based approach I mentioned in comments yesterday. It's not perfect, since the segments of the paths it computes are all constrained to occur in one of 16 directions on a chessboard (king's moves plus knight's moves). That said, it gets within 2% of the correct values (always slightly overestimating) for each of the three pairwise distances in your example.

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))

Scaleboard

Great example using gdistance! I was trying to figure out a way to use a convex hull as well, but this is much better! Thanks! – Mazonson

Everything worked as needed, just the plotting to view the selected paths gave me the following error message: > SL12 <- shortestPath(tr1, Pts[1,], Pts[2,], output="SpatialLines") Error in validObject(.Object) : invalid class “CRS” object: invalid object for slot "projargs" in class "CRS": got class "logical", should be or extend class "character" – Mazonson

I won't be able to help you much without having your spatial data objects in front of me, and can only suggest that you carefully examine the projection metadata attached to both tr1 and Pts, using proj4string(), crs(), identicalCRS(), etc. If they don't match up, you may need to project your points to the raster's CRS or take some other fiddly steps to make them match up. Best of luck! – Ascot