Calculate distance between points and nearest polygon in R

Question

I have a data frame with among others longitude and latitude as variables for about 17,000 points. I also have a shapefile consisting of multiple polygons. The points fall inside or outside the polygons. I would like to add to the data frame the distance to nearest polygon for each point but i have been struggling quite a lot.

This is what i have so far after reading multiple posts on the topic but this only works for points inside the polygons apparently.

library(spatstat)
library(maptools)
shp_Poly <- readShapeSpatial("POLYGON.shp")
W <- as.owin(shp_Poly_BAT)
W <- as.psp(W)
p <- ppp(lon, lat, W)
AVG <- nncross(p, W)

How can i improve this?

Answer 1

Let's set up some sample data:

library(sp)
library(spdep)
example(columbus)
plot(columbus)

We'll now get some points. Click 20 times on the map, some inside and some outside polygons:

pts = locator(20,type="p")

Convert to Spatial data type:

spts = SpatialPoints(pts)

Now use rgeos to compute point-polygon distances, and take the minimum for each point:

library(rgeos)
> apply(gDistance(spts, columbus,byid=TRUE),2,min)
         1          2          3          4          5          6          7 
0.25947439 0.03898401 1.27156515 1.50490316 0.00000000 0.00000000 0.00000000 
         8          9         10         11         12         13         14 
0.00000000 0.00000000 0.31312329 0.26742466 0.53934325 0.07764322 0.03909773 
        15         16         17         18         19         20 
0.11156343 0.29243322 0.08872334 0.00000000 0.00000000 0.00000000 

There you go. That's the minimum distance from each of the 20 points to any of the polygons. The zeroes are for points inside one of the polygons.

Note you should use data in a projected spatial coordinate system and not lat-long.

If you have 17,000 points then you might want to do this in smaller subsets to avoid creating a 17,000 x 5,000 matrix if you have 5,000 polygons. You didn't say.

Answer 2

As mentioned here, one can also use the geosphere::dist2Line for unprojected coordinates (lat-long). Below is an example:

library(sp)
library(geosphere)

# some country polygons to try on
data(wrld_simpl, package = "maptools")
wrld_subset <- wrld_simpl[wrld_simpl@data$ISO2 %in% c("RO","HU","AT","DE","FR"),]

# Generate random points (in and out)
set.seed(2017)
pts <- sp::makegrid(wrld_subset, n = 5)

# compute the shortest distance between points and polygons
# (from ?dist2Line): "returns matrix with distance and lon/lat of the nearest point" & 
# "the ID (index) of (one of) the nearest objects"; distance is in meters (default)
dist.mat <- geosphere::dist2Line(p = pts, line = wrld_subset)

# bind results with original points
pts.wit.dist <- cbind(pts, dist.mat)
pts.wit.dist[1:3,]
##      x1   x2 distance        lon      lat ID
## 1 -10.0 40.2 767133.6 -1.7808770 43.35992  2
## 2  -0.2 40.2 282022.2  0.1894124 42.71846  2
## 3   9.6 40.2 134383.0  9.1808320 41.36472  2

Plot some results to get an idea

pts.sp <- sp::SpatialPoints(coords      = pts[,c("x1","x2")], # order matters
                            proj4string = wrld_subset@proj4string)
plot(pts.sp, col="red")
plot(wrld_subset, add=TRUE)
# plot arrows to indicate the direction of the great-circle-distance
for (i in 1:nrow(pts.wit.dist)) {
    arrows(x0 = pts.wit.dist[i,1], 
           y0 = pts.wit.dist[i,2], 
           x1 = pts.wit.dist[i,4], 
           y1 = pts.wit.dist[i,5],
           length = 0.1,
           col = "green")
}