10

I'm quite new to R and especially GIS with R. So recently I had to proceed to a quite simple analysis: counting health centers (points) per administrative unit (polygons). Quite common process in arcgis or qgis but I didn't found any similar tool in R, at least with the "sf" package that I'm trying to work with exclusively for now (for an easier learning of the way this package works).

So I googled it and found some people with the same problem but I didn't understand most of the methods they were suggesting, it looked oddly complicated for such a simple operation to me. So I created a function to solve my problem, which works now. But I would like to know in order to get better at this:

Is there an easier way that I missed?

Here's the code below for the function which works simply by: selecting polygons one by one, intersecting it with the point layer, counting the size of the extracted points, adding that count to a vector, repeating with the next polygon... And in the end, binding that "count" vector to my initial polygon dataframe.

CountPointsInPolygons <- function(pts, polygons){
  countPts = c()
  for (i in 1:nrow(polygons)) {
    polySelect <- polygons[i,]
    pts2 <- st_intersection(pts, polySelect)
    countPts[i] = nrow(pts2)

  }

  return(cbind(polygons,countPts))
}
3
  • It sounds like you are seeking a code review for which there is the Code Review Stack Exchange.
    – PolyGeo
    May 23 '19 at 23:03
  • @PolyGeo I found this question and answer useful (and helpful). It is aiming at techniques in using spatial tools in R, which involves writing code but specific to GIS. In my humble opinion, I think it could reopened. Oct 24 '19 at 15:33
  • @user3386170 to get the question re-opened simply edit it and then it will go through the review queue to see if the community agrees with you. Being closed does not prevent a Q&A from being useful and your upvote can be used to reflect your opinion.
    – PolyGeo
    Oct 24 '19 at 19:59
16

It's faster and easier to use st_intersects, no need to loop. The output is a bit obscure, essentially a classed list of feature IDs intersected, so we get the lengths() which is the number of points inside each feature.

First polygon has six points in it, in this example.

  library(sf)
#> Linking to GEOS 3.7.0, GDAL 2.4.0, PROJ 5.2.0
poly <- read_sf(system.file("shape/nc.shp", package = "sf"))
set.seed(77)
pts <- st_sample(poly, size = 750)

## add point count to each polygon
(poly$pt_count <- lengths(st_intersects(poly, pts)))
#>   [1]  6  8 12  3  9  2  3  4  9  9 10  5  6  5  4  9  4  7  4  1  4  4  4
#>  [24]  5  4 10 10 11  9  1  9  3  4  3  4  5 11  5  6  4  3  7 10  8  3  5
#>  [47]  7  9  6 14  9  9  6 11  6  9 12  7  5  7  9 11 10  6  6 12 12  6  4
#>  [70]  7  9  7  4  3  8  5  3  6 19  5  7 12  5 12  9 10  8 23  6  5  5  7
#>  [93] 12  9 12 16 18 17  2 14

plot(poly[1, 1], reset = FALSE, col = "grey")
plot(pts, add = TRUE)

Created on 2019-05-23 by the reprex package (v0.3.0)

1
  • Not the use of lengths not length (i.e. with an 's'). I'm putting this comment here because I always forget this and spend 30 minutes every time working out what I've done wrong
    – Phil
    Sep 8 at 21:02

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