Create a circle of defined radius around a point and then find the overlapping area on a shapefile

Question

My use case is related to calculating the area of a polygon in a shapefile that is covered by another polygon.

Specifically, I want to do the following:

1. Create circles that have a 100m radius, from a data.frame of long/lat points
2. Overlay these circles onto a shapefile
3. Calculate the area of each of these polygons that are covered by these circles.
4. Ideally return back a dataframe where each row contains the area covered by all circles for each polygon.

Here I create some simulated data points.

mean_lat <- 46.07998
sd_lat <- 0.1609196
mean_long <- 8.931849
sd_long <-  0.1024659

dat_sim <- data.frame(lat = rnorm(500, mean = mean_lat, sd = sd_lat),
                      long = rnorm(500, mean = mean_long, sd = sd_long))

Then we get the shapefiles

library(raster)
library(tidyverse)
library(sf)

# regular sp format
dat_ticino_sp <- getData(name = "GADM", 
                         country = "CHE", 
                         level = 3) %>% 
  subset(.@data$NAME_1 == "Ticino")


# convert to sf format
dat_ticino_sf <- getData(name = "GADM", 
                         country = "CHE", 
                         level = 3) %>%
  # convert to simple features
  sf::st_as_sf() %>%
  # Filter down to Ticino
  dplyr::filter(NAME_1 == "Ticino")

So what I want to do is :

• create circles with a 100m radius in the dat_sim, to overlay onto the shapefile dat_ticino_sf or dat_ticino_sp files.

• For each polygon in the shapefile, calculate how much of that polygon is covered by circles.

• Return a data frame where each row (polygon) contains the area covered by the circles.

• Plot this.

Bonus points if this can be done in Simple Features.

Answer 1

I think this should do what you want.

Essentially what I mentioned in my comment. One thing I didn't mention is that you will want to transform to an appropriate equal-area projection that uses metres so that you can buffer by 100m and be confident that the areas calculated across your study are equivalent. I don't know the correct one for Switzerland, but I found this SO answer and used that.

The beauty of sf is that it integrates (almost magically!) so well with dplyr - see the last step.

library(raster)
library(tidyverse)
library(sf)
#> Linking to GEOS 3.5.0, GDAL 2.1.1, proj.4 4.9.3

# convert to sf format
dat_ticino_sf <- getData(name = "GADM", 
                         country = "CHE", 
                         level = 3) %>%
  # convert to simple features
  sf::st_as_sf() %>%
  # Filter down to Ticino
  dplyr::filter(NAME_1 == "Ticino") %>% 
  st_transform(3035) # Reproject to EPSG:3035 as mentioned above

mean_lat <- 46.07998
sd_lat <- 0.1609196
mean_long <- 8.931849
sd_long <-  0.1024659

set.seed(42)
dat_sim <- data.frame(lat = rnorm(500, mean = mean_lat, sd = sd_lat),
                      long = rnorm(500, mean = mean_long, sd = sd_long))

# Convert to sf, set the crs to EPSG:4326 (lat/long), 
# and transform to EPSG:3035
dat_sf <- st_as_sf(dat_sim, coords = c("long", "lat"), crs = 4326) %>% 
  st_transform(3035)

# Buffer circles by 100m
dat_circles <- st_buffer(dat_sf, dist = 100)

# Intersect the circles with the polygons
ticino_int_circles <- st_intersection(dat_ticino_sf, dat_circles)
#> Warning: attribute variables are assumed to be spatially constant
#> throughout all geometries

## Plot a zoomed in map of polygons and overlay intersected circles to double check:
## (I assumed NAME_3 is the best unique identifier for the polygons)
bb <- st_bbox(ticino_int_circles)

plot(dat_ticino_sf[, "NAME_3"], 
     xlim = c(mean(c(bb["xmin"], bb["xmax"])) - 1000, 
              mean(c(bb["xmin"], bb["xmax"])) + 1000), 
     ylim = c(mean(c(bb["ymin"], bb["ymax"])) - 1000, 
              mean(c(bb["ymin"], bb["ymax"])) + 1000))

plot(ticino_int_circles[, "NAME_3"], add = TRUE)

# Summarize by NAME_3
# which aggregates all of the circles by NAME_3, 
# then calculate the area
ticino_int_circles_summary <- group_by(ticino_int_circles, NAME_3) %>% 
  summarise() %>% 
  mutate(area = st_area(.))