I'm trying to calculate density within a shapefile, but I'm fairly confident I'm doing it wrong. The idea is to figure out which geographical regions there have been the most sales by density.

Here is a link to the file that I use (testdata.shp)


sample <- st_read("testdata.shp")

sample$area <- st_area(sample$geometry)

density_calc <-sample %>% st_buffer(0) %>% group_by(areas) %>% summarise(`Sales (density)` = sum(sales)/sum(area))

Here are the details of the shapefile:

Geometry set for 2106 features 
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: -120.0065 ymin: 35.00184 xmax: -114.0396 ymax: 42.00221
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs

I guess my issue is, I don't really know what is right and wrong, so I have no clue if I did it correctly.

Sorry if it's not the most extensive question, I just don't remember my high school geometry that well!


If you want to calculate the number of sales per area unit, your code works fine. Although you could also just run

sample$density <- sample$sales / sample$area

instead. Word of caution: there are NA values in your data, be careful you interpret them correctly.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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