Jim
Jim

Reputation: 578

Clipping surface polygon by a rather nasty large Spatial polygon data frame

I have come across a tricky issue.

I am trying to perfom a simple polygon clip using sp package, either using the fuction st_difference(st_union(x),st_union(y)) (or a variant there of) or st_intersectionfunction, whichever works best.

While this is easy with two surface polygons, I need to clip it to a horrible large downloaded Large SpatialPolygonsDataFrame, its just a shapefile for the UK, downloaded from: https://gadm.org/download_country_v3.html

The shapefiles are as as follows (plotted in leaflet): enter image description here

    > str(uk)
    Formal class 'SpatialPolygonsDataFrame' [package "sp"] with 5 slots
      ..@ data       :'data.frame': 1 obs. of  70 variables:
      .. ..$ ID_0      : Factor w/ 1 level "242": 1
      .. ..$ ISO       : Factor w/ 1 level "GBR": 1
      .. ..$ NAME_0    : Factor w/ 1 level "United Kingdom": 1
    # .....etc.
    #
    > str(box)
    sfc_POLYGON of length 1; first list element: List of 1
     $ : num [1:5, 1:2] -7.237 0.126 0.126 -7.237 -7.237 ...
     - attr(*, "class")= chr [1:3] "XY" "POLYGON" "sfg"

I want to clip the box (blue) to the uk, the reason for this is it takes too long to render the shapefile in of the UK (and france for that matter) in leaflet.

Upvotes: 0

Views: 350

Answers (1)

Wimpel
Wimpel

Reputation: 27732

This will probably do?

sf::st_intersection(UK, box)

Complete code

library(sf)

UK <- readRDS("./gadm36_GBR_0_sf.rds")

#create box since it was not provided in question
box <- c("POLYGON((-7.237 48,0 48,0 52,-7.237 52, -7.237 48))") %>% 
  st_as_sfc(crs = "+proj=longlat +datum=WGS84")

mapview::mapview(list(UK,box))

enter image description here

mapview::mapview( st_intersection(UK, box) )

enter image description here

update

If you want to cut the box with the shapefile of the UK, use st_difference()

mapview::mapview( st_difference (box, UK) )

enter image description here

Upvotes: 2

Related Questions