3

Also asked here

Using R, I would like to overlay some spatial points and polygons in order to assign to the points some attributes of the geographic regions I have taken into consideration.

What I usually do is to use the command over of the sppackage. My problems is that I'm working with a large number of geo-referenced events that happened all over the globe and in some cases (especially in coastal areas), the longitude and latitude combination falls slightly outside the country/region border. Here a reproducible example based on in this very good question.

## example data
set.seed(1)
library(raster)
library(rgdal)
library(sp)
p <- shapefile(system.file("external/lux.shp", package="raster"))
p2 <- as(0.30*extent(p), "SpatialPolygons")
proj4string(p2) <- proj4string(p)
pts1 <- spsample(p2-p, n=3, type="random")
pts2<- spsample(p, n=10, type="random")
pts<-rbind(pts1, pts2)

## Plot to visualize
plot(pts, pch=16, cex=.5,col="red")
plot(p, col=colorRampPalette(blues9)(12), add=TRUE)

enter image description here

# overlay
pts_index<-over(pts, p)

# result
pts_index

   ID_1       NAME_1 ID_2           NAME_2 AREA
1    NA         <NA> <NA>             <NA>   NA
2    NA         <NA> <NA>             <NA>   NA
3    NA         <NA> <NA>             <NA>   NA
4     1     Diekirch    2         Diekirch  218
5     3   Luxembourg    8         Capellen  185
6     3   Luxembourg    9 Esch-sur-Alzette  251
7     3   Luxembourg    8         Capellen  185
8     1     Diekirch    1         Clervaux  312
9     1     Diekirch    2         Diekirch  218
10    3   Luxembourg    8         Capellen  185
11    1     Diekirch    2         Diekirch  218
12    1     Diekirch    3          Redange  259
13    2 Grevenmacher    7           Remich  129

Is there a way to give to the over function a sort of tolerance in order to capture also the points that are very close to the border?

NOTE:

Following this I could assign to the missing point the nearest polygon, but this is not exactly what I am after.

1

not sure how to do this in R, but this is the approach I use for this use case:-

  • buffer your polygons (I think rgeos has this option). This expands them by a certain distance outwards. This buffer amount should equal your tolerance.
  • test for point intersection with existing (unbuffered) polygons, just like you're doing now. if that works, use that polygon
  • if it doesn't, test for intersection on the buffered polygons. if that works, use that polygon.
  • failing that, you have an outlier.

It's important to do these tests in that order, to get around the problems where the buffering causes overlaps on the shared boundaries between polygons. The second test should catch those cases where the points fall slightly outside the outer boundary, but within your tolerance.

  • 1
    but do the test against the unbuffered polygons first so you don't waste time computing buffers you don't need if everything is in the source polygons! Also, it seems possible for a point to be in two (or more) buffers, so that case needs consideration and handling. – Spacedman Jul 18 '18 at 6:40
  • So, considering this, the nearest neighbor sounds like to most correct way to go. Am I wrong? – Nemesi Jul 18 '18 at 7:21
  • I agree with @Spacedman that using the buffers I would run the risk of having a point attributed to two different polygons. The example I reported above is only a small sample. In my real application, I have thousands of points, hundreds of them in coastal areas (often some fraction of a degree outside the polygon). – Nemesi Jul 18 '18 at 8:29

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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