13

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?

2 Answers 2

15

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.

8
  • 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, 2017 at 17:03
  • 1
    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, 2019 at 18:25
  • 1
    @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, 2019 at 18:26
  • 2
    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, 2019 at 21:02
  • 2
    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, 2019 at 21:46
4

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")
}

enter image description here

2
  • 1
    nice solution, unfortunately with a large number of points 10000+ this method is quite slow. Oct 13, 2020 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, 2020 at 14:19

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.