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:
- Create circles that have a 100m radius, from a data.frame of long/lat points
- Overlay these circles onto a shapefile
- Calculate the area of each of these polygons that are covered by these circles.
- 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 shapefiledat_ticino_sf
ordat_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
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(.))
No comments:
Post a Comment