# 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... – Skeboo 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"’ – Gregory Sep 6 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... – Spacedman Sep 6 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... – Spacedman Sep 6 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. – Spacedman Sep 6 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")