# Calculate distance between points and nearest polygon in R

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)
W <- as.owin(shp_Poly_BAT)
W <- as.psp(W)
p <- ppp(lon, lat, W)
AVG <- nncross(p, W)
``````

How can i improve this?

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.

• That works really well on my end, just what i needed, thanks very much! I have about 25 polygons at the moment and that is not likely to change by much... Jan 19 '17 at 17:03
• When I try to reproduce your example, I'm getting and error at `apply(gDistance(spts, columbus,byid=TRUE),2,min)`: Error in (function (classes, fdef, mtable) : unable to find an inherited method for function ‘is.projected’ for signature ‘"sf"’ Sep 6 '19 at 18:25
• @Gregory you seem to be using `sf` spatial objects where my code is for `sp` spatial objects. Convert sp to sf with `as(spts, "Spatial")` or you'll need to rewrite using `sf_distance` functions somehow... Sep 6 '19 at 18:26
• Suggest you make a new question asking how to do this using `sf` objects. Make a reproducible example using sample data if possible, and include your code... Sep 6 '19 at 21:02
• Yup. The latest `spdep` example (I've just upgraded to 1.1.2) gives you `columbus` as an sf object. Convert that with `columbus = as(columbus,"Spatial")` and it should proceed. But a solution within `sf` using `st_distance` or `st_nearest_{feature/points}` might be better. Sep 6 '19 at 21:46

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 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")
}
`````` • nice solution, unfortunately with a large number of points 10000+ this method is quite slow. Oct 13 '20 at 18:11
• Hi @HermanToothrot, this doesn't surprises me. I think the current implementation computes all possible distances to all vertices of the polygons and then picks the minimum. If there are many vertices and/or many points from which you need to compute the distances, then that will also reduce a lot the computation speed. One work around would be reduce the number of the vertices in the polygons. Another, is simply brute force of many CPU-s, so parallelization. Sorry that I am not that helpful. Consider asking a separate question, link this post and feel free to ask for speed improvements. Oct 15 '20 at 14:19