R - Compute density around a point

Question

I have a dataset with cities in a French département (more or less a county) and the population in each of the cities. I have another dataset with geographic coordinates of nursing homes in this département.

I would like to compute the population around a radius of 20km of each of my nursing homes.

But the thing is, I do not know how to do that...

Here is the list of the problems I should resolve:

1. Which cities are inside the circle?
2. Which percentage of a city is there in the circle?
3. Under the assumption that the population is homogeneously evenly distributed, how many people are there in the part of city inside the circle?
4. Compute the sum of the population inside the circle

Here is a map of my département, the nursing homes are the points in black, and the buffers in yellow.

Here is my code (well, it does not work...):

library(sf)
# To know which cities are in the circle:
cities_in_circle <- st_intersects(departement, buffer_nh)

It gives a list and I do not know how to use it...

To know the area overlapped I should use st_intersection and st_area but it does not work with my list above.

Any help would be gladly welcomed!

---

Here is a sample of my data : For the cities :

              NOM                 POP               geometry
1          Cuire-Le-Haut   230.2751 MULTIPOLYGON (((841575.8 65...
2           Cuire-Le-Bas   219.7269 MULTIPOLYGON (((842816.9 65...
3         Pierre-Brunier   164.8633 MULTIPOLYGON (((842340.5 65...
4     Coste-Croix-Rousse   214.6620 MULTIPOLYGON (((842691.5 65...
5             Coste-Nord   105.7979 MULTIPOLYGON (((842469 6522...
6           Centre-Bourg   220.3814 MULTIPOLYGON (((843779 6524...
7                 Vernay   109.1020 MULTIPOLYGON (((843308.8 65...
8  Jean-Moulin de-Gaulle   169.4209 MULTIPOLYGON (((843332.1 65...
9             Montchoisi   121.2074 MULTIPOLYGON (((844903.3 65...
10     Magnolles Pasteur   495.6844 MULTIPOLYGON (((842804.5 65...

Unfortunately, I cannot manage to give the whole geometry...

And for my nursing homes :

Simple feature collection with 7 features and 1 field
geometry type:  POINT
dimension:      XY
bbox:           xmin: 841863.4 ymin: 6522364 xmax: 843572.9 ymax: 6524230
projected CRS:  RGF93 / Lambert-93
                              name                 geometry
1        EHPAD Cercle de la Carette POINT (843446.3 6522377)
2                 EHPAD La Rochette POINT (841890.1 6522749)
3                   EHPAD Le Manoir POINT (841863.4 6522792)
4        EHPAD Residence des Canuts   POINT (842821 6522520)
5            RESIDENCE LE VAL FORON POINT (843022.7 6523452)
6              Residence Marie Lyan POINT (843572.9 6524230)
7 USLD du docteur Frederic Dugoujon   POINT (842758 6522364)

Answer 1

I think a method that will work well is to

1. Assume your population is evenly distributed within each department
2. Calculate area and density for each department
3. Rasterize department polygons using the density attribute
4. Extract values within your already computed buffer_nh polygons.

It looks like you are using EPSG:2154, which is projected in meters. I'm assuming your departments are already in the same projection. It looks like your polygons are called departments and buffers are already computed and called buffers-nh.

Without sample data I can't test this, but I think it should work:

library(sf)
library(raster)

res = 100  # Resolution to use for raster, set to appropriate value

## calculate population density per raster pixel
departement$density = departement$pop/ (st_area(departement) / res^2)

## Making raster from polygons. 
d_geo = st_bbox(departement)  # extract bounding box
num_x = (d_geo$xmax - d_geo$xmin) / res 
num_y = (d_geo$ymax - d_geo$ymin) / res

r = raster(ncol= num_x , 
           nrow= num_y, 
           xmn = d_geo$xmin ,
           xmx= d_geo$xmax, 
           ymn= d_geo$ymin , 
           ymx= d_geo$ymax, 
           crs = crs(departement)
           )
r = rasterize(departement, r, field = "density")

plot(r) # should look like your departement, with color to population density

buffer_nh$pop = raster::extract(r, buffer_nh, fun = sum, na.rm = TRUE)

Another option is to use the velox package, which was recently taken off CRAN, but previous versions remain available. The velox package is built on top of the raster package and handles raster extraction much faster, which can be helpful if you are using rasters with a higher resolution. To use this method, you would replace the last line of the method above with:

library(velox)
vs = velox(r)
buffer_nh$pop2 <- vs$extract(buffer_nh, fun = function(x){sum(x, na.rm= TRUE)})