[R-sig-Geo] Shortest distance between points within a polygon

[Robert J. Hijmans]

Ruari, You might be able to use gCrosses together with optim or
another optimization function to find the best route. Below is a very
simple example. For more complex polygons you may need to allow for
several intermediate points (perhaps following the polygon boundary at
times). Unless you want to spend a lot of time on this, gDistance is
probably the safer way to go.  Robert


library(raster)
library(rgeos)

opfun <- function(p) {
if (p[1] < xmin(e) | p[1] > xmax(e) | p[2] < ymin(e) | p[2] > ymax(e)) {
return(Inf)
}
via <- cbind(p[1], p[2])
line <- SpatialLines(list(Lines(list(Line(rbind(from, via, to))), "1")))
if (gCrosses(line, pol)) {
return(Inf)
}
d1 <- pointDistance(from, via, lonlat=TRUE)
d2 <- pointDistance(via, to, lonlat=TRUE)
d1 + d2
}


pol <- matrix(c(-79.49102, 11.152694, -60.77844, 50.074852, -33.08383,
5.913173, 43.26348, 2.919161, -27.84431, -29.266470, -48.80240,
-68.188628, -73.50300, -32.260482, -108.68264, -42.739525, -84.73054,
-14.296409, -110.17965, 26.122755, -79.49102, 11.152694), ncol=2,
byrow=TRUE)
pol <- Polygons(list(Polygon(pol)), 1)
pol <- SpatialPolygons( list( pol) )
e <- extent(pol)
from <- c(-95, 11)
to <- c(-59, 40)

op <- optim(as.vector(coordinates(gCentroid(pol))), fun)

via <- rbind(op$par)
opline <- SpatialLines(list(Lines(list(Line(rbind(from, via, to))), "1")))

plot(pol, border='red')
x <- rbind(from, via, to)
points(x, pch=20)
lines(x, col='blue', lwd=2)

distance <- op$value
distance


On Mon, Oct 13, 2014 at 11:13 AM, Roman Luštrik <roman.lustrik at gmail.com> wrote:
> I was thinking along the lines of this example from `gCrosses`. If
> line crosses the shape, it returns TRUE. If not, you've got your
> distance. I can't speak for the performance part, you better do some
> field testing.
>
> require(rgeos)
>
> l1 = readWKT("LINESTRING(0 3,1 1,2 2,3 0)")
> l2 = readWKT("LINESTRING(0 0.5,1 1,2 2,3 2.5)")
> l3 = readWKT("LINESTRING(1 3,1.5 1,2.5 2)")
>
> pt1 = readWKT("MULTIPOINT(1 1,3 0)")
> pt2 = readWKT("MULTIPOINT(1 1,3 0,1 2)")
>
> p1 = readWKT("POLYGON((0 0,0 2,1 3.5,3 3,4 1,3 0,0 0))")
> p2 = readWKT("POLYGON((2 2,3 4,4 1,4 0,2 2))")
>
>
> plot(p1,border='blue',col='blue');plot(l1,add=TRUE)
> title(paste("Crosses:",gCrosses(p1,l1),
>             "\nOverlaps:",gOverlaps(p1,l1)))
>
> Cheers,
> Roman
>
>
> On Mon, Oct 13, 2014 at 7:12 PM, Barry Rowlingson
> <b.rowlingson at lancaster.ac.uk> wrote:
>> On Mon, Oct 13, 2014 at 12:12 PM, Roman Luštrik <roman.lustrik at gmail.com> wrote:
>>> One way would be to calculate the distance and see if the line connecting
>>> the dots intersects with the polygon edge. `rgeos` package would be the
>>> tool of trade for this.
>>>
>>
>>  I can't see how that works...
>>
>>  A solution would be to rasterise your polygon and use `gdistance` to
>> compute a cost-weighted minimum distance on the raster - you'd set the
>> cost on the transition surface outside your polygon to +Inf and inside
>> the polygon to +1, or maybe its zero cost inside the polygon and
>> anything finite outside.
>>
>> The gdistance vignette has an example of travelling between places in
>> Europe where travel is restricted to the land mass only.
>>
>> The precision of your raster will determine the precision of your
>> result, but finer rasters will take longer.
>>
>> Barry
>
>
>
> --
> In God we trust, all others bring data.
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo