# R - Compute density around a point

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)
``````
• What format is your population data in? The process that occurs to me is to buffer the points to get polygons, and rasterize the population data, then use raster::extract or velox::velox\$extract to pull the population within the buffer polygons. Can you provide some sample data, particularly showing how your population data is structured? Apr 7, 2021 at 17:04
• You can get a grid of 100m population counts for France from Worldpop: worldpop.org - create 20km polygon buffers with st_buffer around your nursing homes and extract total population count from the WorldPop raster. This way you don't need to use city boundary or population data at all. Apr 7, 2021 at 17:34
• I thank you @Spacedman, but I would prefer to keep my cities since I have precise data on the number of people aged 80years old and over, people aged 75 years old and over, the housing rents... Apr 8, 2021 at 8:21

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

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)})
``````
• Thank you for your answer ! However, the package `velox` has been removed from CRAN and cannot be installed anymore... Apr 20, 2021 at 14:07
• @Emeline I edited the answer to separate the `velox` approach. The `raster::extract()` route should work fine, it might take a little while if you have a dense raster layer. When I tested the code on a ~2500 x 4000 cell raster (about 2.5mb on disk as tiff) and 9 polygons the extract step took about 3 minutes, while the velox approach was around 15 seconds. Apr 21, 2021 at 21:17
• This is EXCELLENT ! I really tried to find how `raster::extract()` worked but I never managed to make it actually work ! Thank you !!! Apr 22, 2021 at 7:29