3

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?

Is my function efficiently written or could I improve it somehow?

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

closed as too broad by PolyGeo May 23 at 23:04

Please edit the question to limit it to a specific problem with enough detail to identify an adequate answer. Avoid asking multiple distinct questions at once. See the How to Ask page for help clarifying this question. If this question can be reworded to fit the rules in the help center, please edit the question.

  • It sounds like you are seeking a code review for which there is the Code Review Stack Exchange. – PolyGeo May 23 at 23:03
4

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)

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