R sf Package – How to Count Points in Polygons

code-reviewcountrsf

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

Best Answer

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

Related Solutions

[GIS] Extract median value for polygons from multiresolution raster data, using R (velox, raster -package)

velox maintainer here.

The development version of velox (available on github) now implements a small option for the VeloxRaster_extract method. If small = TRUE, velox will return raster values for all polygons that intersect with the raster extent. Specifically, for all small or oddly shaped polygons that do not intersect with any cell center, velox will look for intersections with entire cell boxes.

Using the development version of velox and setting the small option to TRUE, your example now yields no NA values:

# Install latest version of velox
library(devtools)
install_github("hunzikp/velox")

# Reproducible example:   
library(velox)
library(raster)

## Make VeloxRaster
mat <- matrix(1:100, 10, 10)
extent <- c(0,1,0,1)
vx <- velox(mat, extent=extent, res=c(0.1,0.1), crs="+proj=longlat     +datum=WGS84 +no_defs")

## Make SpatialPolygonsDataFrame
library(sp)
library(rgeos)
set.seed(0)
coords <- cbind(runif(10, extent[1], extent[2]), runif(10, extent[3], extent[4]))
sp <- SpatialPoints(coords)

# Default example
# from https://cran.r-project.org/web/packages/velox/README.html
spol_norm <- gBuffer(sp, width=0.2, byid=TRUE)
spdf_norm <- SpatialPolygonsDataFrame(spol_norm,     data.frame(id=1:length(spol_norm)), FALSE)

# Smaller buffer
spol_small<- gBuffer(sp, width=0.05, byid=TRUE)
spdf_small <- SpatialPolygonsDataFrame(spol_small, data.frame(id=1:length(spol_small)), FALSE)

plot(spdf_norm); par(new=F)
plot(spdf_small)

## Extract values and calculate mean, see results
(ex.mat.norm <- vx$extract(spdf_norm, fun=median))
(ex.mat.small <- vx$extract(spdf_small, fun=median, small = TRUE)) # -> No NA values anymore.

R – Counting Points in Polygon by Time or Attribute Condition

If you want to go the R way I would suggest this workflow:

assign each polygon an unique id / such as the polyid in your example
perform sf::st_join(points, polygons) on your points and polygons objects, with - crucially - your points in the first place of the join; it will produce your points object, enriched by the data columns from the polygons object (such as the unique id / polyid).
at this step you can drop the geometry (= sf::st_drop_geometry()), and do simple aggregation over your time dimension and polygon id, as both are present in your points object; {dplyr} is your friend here as you are dealing with a regular data frame object

Since an example is worth 1000 of words :) consider this piece of code. It builds on the well known & much loved North Carolina shapefile that ships with {sf}.

It first produces five polygons of dissolved North Carolinas stacked on top of each other, with polygon id's from 1 to 5.

Then it produces three semi-random points, two with time character 1 and one with 3 (or whatever).

The third step is where the action happens: first the points are spatially joined to polygons (it will be a kind of a cross join, since all points are in NC) and then the result is filtered on the secondary criterion - I am using polygon id = time, which makes no sense but is there only for the sake of example. Then it is a question of group by and summarise / functionally equivalent to select polygon id, count(*) from yer data group by polygon id in SQL speak.

library(sf)
library(dplyr)

# a set of polygons in time
shape <- st_read(system.file("shape/nc.shp", package="sf")) %>%  # included with sf package
  summarise() %>% # single polygon over entire NC
  tidyr::expand_grid(poly_id = 1:5) %>% # sequence of ids from 1 to 5
  st_as_sf(crs = 4267) # this information was lost during the expand_grid (tidyr equivalent of cross join)

# a set of points in time
cities <- data.frame(name = c("Raleigh", "Greensboro", "Wilmington"),
                  x = c(-78.633333, -79.819444, -77.912222),
                  y = c(35.766667, 36.08, 34.223333),
                  time = c(1, 1, 3)) %>%  # time, to be matched with poly_id
  st_as_sf(coords = c("x", "y"), crs = 4326) %>% # CRS of the coordinates
  st_transform(st_crs(shape)) # transorm to CRS of the shape object

# this is the action!   
cities %>% 
  st_join(shape) %>% # spatial join only, without the second criterion (yet)
  st_drop_geometry() %>% # drop geometry (no longer necessary)
  filter(time == poly_id) %>%  # filter over the second condition
  group_by(poly_id) %>%  # group by
  tally(name = "point_totals") # summarise the result
# A tibble: 2 × 2
#  poly_id time_totals
#    <int>       <int>
#1       1           2
#2       3           1

Best Answer

Related Solutions

[GIS] Extract median value for polygons from multiresolution raster data, using R (velox, raster -package)

R – Counting Points in Polygon by Time or Attribute Condition

Related Question