7

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.

  • 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? – andyteucher Feb 23 '17 at 5:31
  • Thanks @andyteucher! I'll have a go at this and see where I get to. – Nick Tierney Feb 23 '17 at 5:58
7

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(.))
  • This got me to where I wanted to get to, thank you very much! Amazing that simple features plays so well with dplyr, this is magical. I've got a couple of extra steps that I've done afterwards to calculate the proportion of area covered, am I best to post these as another comment/answer below? – Nick Tierney Feb 24 '17 at 2:01
  • Sure, I think another answer showing how you put it all together would great. – andyteucher Feb 24 '17 at 4:30

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.