Dissolving large shapefiles of protected areas in R

Question

My overall aim is to measure the number of GPS points that overlap with protected areas. I am using the protected planet database to source my shapefiles for Ireland and the UK.

I don't need any detail of the shapefiles aside from the outline and so I try to dissolve them. Indeed, if I don't do this then I encounter issues with multiparts.

My attempt seems to work for Ireland but not the UK, presumably due to the larger size of the latter. To clarify, R hangs for hours when I run it for the UK. For this code to run you'll need to download the zip from protected planet and extract the three individual files for either UK or Ireland.

library(sf)
library(tidyverse)

spath <- "data"
(shps <- list.files(spath, "polygons.shp$"))
mydata <- do.call(rbind, lapply(shps, read_sf))
marine_only <- mydata %>% filter(MARINE == "2") 

marine_only$GIS_M_AREA <- st_area(marine_only)
marine_only

marine_dissolved <-
  marine_only %>%
  summarise(REP_AREA = sum(REP_AREA))

Answer 1

I found that st_within() is a bit slow. Something like the procedure below might work (watch out for typos)

library (sf)

# extract the geometry of the conservation areas to a new variable
# to get rid of all other attributes

polygon_outlines <- st_geometry (conservation_area_polygons) %>% st_as_sf()

# add a new field to the polygon_outlines to identify points within or without
# so points that have within==1 are ... within boundaries

polygon_outlines$within <- 1   

# create a join between the points and the polygons. this adds the field 'within' 
# to the points - points within a polygon will have a value of within==1 otherwise NA

dt <- st_join (point_data, polygon_outlines)

# select the points where within!=NA

dt <- dt[!(is.na(dt$within), ]

print (nrow(dt))

The number of records in dt should reflect the number of points within conservation areas.