Calculating the distance between polygon and point in R

Question

I have a, not necessarily convex, polygon without intersections and a point outside this polygon. I'm wondering how calculate the Euclidian distance most efficiently in a 2-dimensional space. Is there a standard method in R?

My first idea was to calculate the minimum distance of all the lines of the polygon (extended infinitely so they are line, not line pieces) and then calculate the distance from the point to each individual line using the start of the line piece and Pythagoras.

Do you know about a package that implements an efficient algorithm?

Aphra

I don't know of anything prepackaged that does what you're asking for. Can you give us a little more information? Do you know anything about the the polygon before calculating distance, or can it be any shape? I assume since you mention Pythagoras that this is in a Cartesian space. – Nitride

Based on the proposed method, it sounds like the polygon will not allow internal intersections, which suggests that the vertices have a rotational order (going CW or CCW). Is this true? – Nitride

Have you looked at the spatstat package? I know they have, at least, some distance algorithms for line segment and point sets. – Fractious

Answer 1

You could use the rgeos package and the gDistance method. This will require you to prepare your geometries, creating spgeom objects from the data you have (I assume it is a data.frame or something similar). The rgeos documentation is very detailed (see the PDF manual of the package from the CRAN page), this is one relevant example from the gDistance documentation:

pt1 = readWKT("POINT(0.5 0.5)")
pt2 = readWKT("POINT(2 2)")
p1 = readWKT("POLYGON((0 0,1 0,1 1,0 1,0 0))")
p2 = readWKT("POLYGON((2 0,3 1,4 0,2 0))")
gDistance(pt1,pt2)
gDistance(p1,pt1)
gDistance(p1,pt2)
gDistance(p1,p2)

readWKT is included in rgeos as well.

Rgeos is based on the GEOS library, one of the de facto standards in geometric computing. If you don't feel like reinventing the wheel, this is a good way to go.

Bursa

Answer 2

I decided to return and write up a theoretical solution, just for posterity. This isn't the most concise example, but it is fully transparent for those who want to know how to go about solving a problem like this by hand.

The theoretical algorithm

First, our assumptions.

1. We assume the polygon's vertices specify the points of a polygon in a rotational order going clockwise or counter-clockwise and the lines of the polygon cannot intersect. This means we have a normal geometric polygon, and not some strangely-defined vector graphic shape.
2. We assume this is a set of Cartesian coordinates, using 'x' and 'y' values that represent location on a 2-dimensional plane.
3. We assume the point must be outside the internal area of the polygon.
4. Finally, we assume that the distance desired is the minimum distance between the point and all of the infinite number of points on the perimeter of the polygon.

Now before coding, we should write out in basic terms what we want to do. We can assume that the shortest distance between the polygon and the point outside the polygon will always be one of two things: a vertex of the polygon or a point on a line between two vertices. With this in mind, we do the following steps:

1. Calculate the distances between all vertices and the target point.
2. Find the two vertices closest to the target point.
3. If either: (a) the two closest vertices are not adjacent or (b) the inside angles of either of the two vertices is greater or equal to 90 degrees, then the closest vertex is the closest point. Calculate the distance between the closest point and the target point.
4. Otherwise, calculate the height of the triangle formed between the two points.

We're basically just looking to see if a vertex is closest to the point or if a point on a line is closest to the point. We have to use a few trig functions to make this work.

The code

To make this work properly, we want to avoid any 'for' loops and want to only use vectorized functions when looking at the entire list of polygon vertices. Luckily, this is pretty easy in R. We accept a data frame with 'x' and 'y' columns for our polygon's vertices, and we accept a vector with one 'x' and 'y' value for the point's location.

get_Point_Dist_from_Polygon <- function(.polygon, .point){

    # Calculate all vertex distances from the target point.
    vertex_Distance <- sqrt((.point[1] - .polygon$x)^2 + (.point[2] - .polygon$y)^2)

    # Select two closest vertices.
    min_1_Index <- which.min(vertex_Distance)
    min_2_Index <- which.min(vertex_Distance[-min_1_Index])

    # Calculate lengths of triangle sides made of
    # the target point and two closest points.
    a <- vertex_Distance[min_1_Index]
    b <- vertex_Distance[min_2_Index]
    c <- sqrt(diff(.polygon$x[c(min_1_Index, min_2_Index)])^2 + diff(.polygon$y[c(min_1_Index, min_2_Index)])^2)

    if(abs(min_1_Index - min_2_Index) != 1 |
        acos((b^2 + c^2 - a^2)/(2*b*c)) >= pi/2 | 
        acos((a^2 + c^2 - b^2)/(2*a*c)) >= pi/2
        ){
        # Step 3 of algorithm.
        return(vertex_Distance[min_1_Index])
    } else {
        # Step 4 of algorithm.
        # Here we are using the law of cosines.
        return(sqrt((a+b-c) * (a-b+c) * (-a+b+c) * (a+b+c)) / (2 * c))
    }
}

Demo

polygon <- read.table(text="
x,  y
0,  1
1,  0.8
2,  1.3
3,  1.4
2.5,0.3
1.5,0.5
0.5,0.1", header=TRUE, sep=",")

point <- c(3.2, 4.1)

get_Point_Dist_from_Polygon(polygon, point)
# 2.707397

Nitride

Perhaps I misunderstand this solution, but the facet of the polygon formed by the two closest vertices to the point need not be the closest facet of the polygon to the point. – Abolish

I completely agree with Richard. Dinre, your solution fails for a simple example where a query point is located below a quadrilateral with a very long edge at the bottom, where the two closest vertices to the query point are located right above the edge (not on the same side of the edge as the query point). – Decry

Answer 3

Otherwise:

p2poly <- function(pt, poly){
    # Closing the polygon
    if(!identical(poly[1,],poly[nrow(poly),])){poly<-rbind(poly,poly[1,])}
    # A simple distance function
    dis <- function(x0,x1,y0,y1){sqrt((x0-x1)^2 +(y0-y1)^2)}
    d <- c()   # Your distance vector
    for(i in 1:(nrow(poly)-1)){
        ba <- c((pt[1]-poly[i,1]),(pt[2]-poly[i,2])) #Vector BA
        bc <- c((poly[i+1,1]-poly[i,1]),(poly[i+1,2]-poly[i,2])) #Vector BC
        dbc <- dis(poly[i+1,1],poly[i,1],poly[i+1,2],poly[i,2]) #Distance BC
        dp <- (ba[1]*bc[1]+ba[2]*bc[2])/dbc          #Projection of A on BC
        if(dp<=0){ #If projection is outside of BC on B side
            d[i] <- dis(pt[1],poly[i,1],pt[2],poly[i,2])
            }else if(dp>=dbc){ #If projection is outside of BC on C side
                d[i] <- dis(poly[i+1,1],pt[1],poly[i+1,2],pt[2])
                }else{ #If projection is inside of BC
                    d[i] <- sqrt(abs((ba[1]^2 +ba[2]^2)-dp^2))
                    }
        }
    min(d)
    }

Example:

pt <- c(3,2)
triangle <- matrix(c(1,3,2,3,4,2),byrow=T, nrow=3)
p2poly(pt,triangle)
[1] 0.3162278

Sagesagebrush

There is an example at the bottom. Poly can be a matrix, an array or a data.frame, as long as each row is a vertex coordinates. This is the simplest algorithm I could think of, of the top of my head. – Sagesagebrush

I like your effort but the accepted answer was easier to implement. – Aphra

Answer 4

I used distm() function in geosphere package to calculate the distence when points and apexes are presented in coordinate system. Also, you can easily make some alternation by substitude dis <- function(x0,x1,y0,y1){sqrt((x0-x1)^2 +(y0-y1)^2)} for distm() .

algo.p2poly <- function(pt, poly){
  if(!identical(poly[1,],poly[nrow(poly),])){poly<-rbind(poly,poly[1,])}
  library(geosphere)
  n <- nrow(poly) - 1
  pa <- distm(pt, poly[1:n, ])
  pb <- distm(pt, poly[2:(n+1), ])
  ab <- diag(distm(poly[1:n, ], poly[2:(n+1), ]))
  p <- (pa + pb + ab) / 2
  d <- 2 * sqrt(p * (p - pa) * (p - pb) * (p - ab)) / ab
  cosa <- (pa^2 + ab^2 - pb^2) / (2 * pa * ab)
  cosb <- (pb^2 + ab^2 - pa^2) / (2 * pb * ab)
  d[which(cosa <= 0)] <- pa[which(cosa <= 0)]
  d[which(cosb <= 0)] <- pb[which(cosb <= 0)]
  return(min(d))
}

Example:

poly <- matrix(c(114.33508, 114.33616,
                 114.33551, 114.33824,
                 114.34629, 114.35053,
                 114.35592, 114.35951, 
                 114.36275, 114.35340, 
                 114.35391, 114.34715,
                 114.34385, 114.34349,
                 114.33896, 114.33917,
                 30.48271, 30.47791,
                 30.47567, 30.47356, 
                 30.46876, 30.46851,
                 30.46882, 30.46770, 
                 30.47219, 30.47356,
                 30.47499, 30.47673,
                 30.47405, 30.47723, 
                 30.47872, 30.48320),
               byrow = F, nrow = 16)
pt1 <- c(114.33508, 30.48271)
pt2 <- c(114.6351, 30.98271)
algo.p2poly(pt1, poly)
algo.p2poly(pt2, poly)

Outcome:

> algo.p2poly(pt1, poly)
[1] 0
> algo.p2poly(pt2, poly)
[1] 62399.81

Hydrography