3

I have a set of points (latitude and longitude) and a set of polygons with the same CRS. Since some of the polygons are interescting, it may be that some of the points are in more than one polygon.

How can I count how many times a point lies in all the set of polygons?

inside <- data.frame(as.numeric(!is.na(over(points, as(get(polygons), "SpatialPolygons")))))

This line of code works only to generate a dummy variable equal to 1 if the point is contained in at least a polygon and 0 otherwise.

  • @mdsumner The code in that post counts how many points are inside a certain polygon (the polygons - Italy's regions - do not intercept). Here, what I am trying to find is how many times a point lies inside a different polygon (since the polygons can intersect). – Matteo Ruzzante Feb 19 '17 at 20:54
  • 1
    Ok, whoops. Use the code there to build up an example so we don't have to do it for you. Feel free to ask about making a reproducible example if you need. – mdsumner Feb 19 '17 at 21:22
  • 1
    For example, if a point intersects 4 overlapping polygons, that point should be assigned a value of 4? For each point, count the number of intersecting polygons? – Aaron Feb 19 '17 at 21:56
  • 1
    over has an argument, returnList = TRUE, that gives you all intersections instead of the first. Alternatively, swap x and y in over and see what that gives you. – Edzer Pebesma Feb 20 '17 at 9:15
1

You can count the times one point is in a polygon or in different polygons with simple R functions using base and sp packages. With the reproducible example below you can:

  1. Identify and count the times one point is inside a set of polygons
  2. Identify which point/polygon intersects

Load libraries and make some example data:

# Load libraries
library("sp")
library("raster")

# Example coords
coordsPoly1 <- data.frame('x' = c(-63.8, -40.9, -40.7, -63.5, -63.8), 'y' = c(-30.1, -30.1, -37.4, -37.1, -30.1))
coordsPoly2 <- data.frame('x' = c(-47.1, -38.4, -37.3, -46.6, -47.1), 'y' = c(-28.7, -28.6, -47.1, -46.6, -28.7))
coordsPoly3 <- data.frame('x' = c(-48.4, -28.0, -27.4, -48.1, -48.4), 'y' = c(-40.9, -40.9, -49.6, -49.4, -40.9))
coordsPoly4 <- data.frame('x' = c(-49.7, -29.4, -28.6, -48.9, -49.7), 'y' = c(-25.6, -25.8, -39.1, -38.8, -25.6)) 

# Coords to SpatialPolygons objects    
polygons <- SpatialPolygons(Srl = list(Polygons(srl = list(Polygon(coords = coordsPoly1)), ID = "1"), 
                                       Polygons(srl = list(Polygon(coords = coordsPoly2)), ID = "2"), 
                                       Polygons(srl = list(Polygon(coords = coordsPoly3)), ID = "3"),
                                       Polygons(srl = list(Polygon(coords = coordsPoly4)), ID = "4")))

# Make sample points 
points <- spsample(polygons, n = 10, type = "random")

Plot polygons and points:

plot(polygons, border = c("#EF2929", "#729FCF", "#8AE234", "#AD7FA8"), lwd = 3, 
     xlab = "Longitude", ylab = "Latitude", main = "Plot polygons")
text(polygons, labels = seq_along(polygons), col = c("#EF2929", "#729FCF", "#8AE234", "#AD7FA8"), cex = 4)
plot(points, pch = 19, add = TRUE)
text(points, labels = seq_along(points), pos = 3)
box()

examplePlot

Spatial Query: points in polygon

pointsInPolygons <- sp::over(x = points, y = polygons, returnList = TRUE)

# Count times in polygon
counting <- lapply(pointsInPolygons, FUN = function(x) length(x))

df <- data.frame("point" = rownames(t(do.call("cbind", counting))), 
                 "count" = t(do.call("cbind", counting)), 
                 "polygon" = paste(pointsInPolygons))

> df
   point count polygon
1      1     1       2
2      2     1       4
3      3     1       3
4      4     2 c(2, 4)
5      5     2 c(1, 4)
6      6     2 c(1, 4)
7      7     1       2
8      8     1       4
9      9     1       1
10    10     1       4

Note: your plot and example points will be different because the spsample function generates random sampling each time.

  • Thank you very much @Guzmán for the detailed reply! In the meanwhile, I came out with a possible solution to the problem I had using the function "returnList=T". list <- over(points, polygons, "SpatialPolygons"), returnList = T) count <- data.frame(sapply(list,length)) – Matteo Ruzzante Mar 3 '17 at 14:16

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.