Calculating the distance between a polygon and a point in R

Question

I have an optionally convex polygon without intersections and a point outside this polygon. I am wondering how to most efficiently calculate the Euclidean distance in two-dimensional space. Is there a standard method in R ?

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

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

Answer 1

You can use the rgeos package and the gDistance method. This will require you to prepare your geometries by creating spgeom objects from the data you have (I assume this is data.frame or something similar). The rgeos documentation is very detailed (see the PDF Package Guide on the CRAN page), this is one example given in 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 also included in rgeos.

Rgeos is based on the GEOS library, one of the actual standards in geometric computing. If you don't want to reinvent the wheel, this is a good way to go.

Answer 2

I decided to go back and write a theoretical solution, only for posterity. This is not the shortest example, but it is completely transparent to those who want to know how to solve this manually.

Theoretical algorithm

First, our assumptions.

• Assume that the vertices of the polygon define the points of the polygon in rotational order, going clockwise or counterclockwise, and the lines of the polygon cannot intersect. This means that we have a normal geometric polygon, and not some weirdly defined vector graphic shape.
• Suppose this is a set of Cartesian coordinates using the values ​​"x" and "y", which represent a location on a two-dimensional plane.
• Suppose that the point must be outside the inner region of the polygon.
• Finally, we assume that the desired distance is the minimum distance between a point and the entire infinite number of points on the perimeter of the polygon.

Now, before coding, we must write down the basic words that we want to make. We can assume that the shortest distance between a polygon and a point outside the polygon will always be one of two things: the vertices of the polygon or the points on the line between two vertices. With this in mind, we take the following steps:

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

We basically just see if the vertex is closer to the point or if the point on the line is closest to the point. To do this, we need to use several trigger functions.

Code

To make this work properly, we want to avoid any "for" loops and want to use only vectorized functions when viewing the entire list of polygon vertices. Fortunately, this is pretty easy in R. We take a data frame with columns “x” and “y” for our vertices of the polygons, and we accept a vector with one value “x” and “y” for the location of the point.

 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+bc) * (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 

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 

Answer 4

I used the distm() function in the geosphere package to calculate the distance when points and vertices are represented in the coordinate system. In addition, you can easily do some interlacing by substituting 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) 

Result:

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