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

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
country = "CHE",
level = 3) %>%
subset(.@data\$NAME_1 == "Ticino")

# convert to sf format
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.

• I don't have time for a full answer now, but st_buffer to turn points into circles, then st_intersection to do the overlay, then st_area on the result? Feb 23 '17 at 5:31
• Thanks @andyteucher! I'll have a go at this and see where I get to. Feb 23 '17 at 5:58

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
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)) ``````# Summarize by NAME_3