[SOLVED] Problem with clipping northern Canada polygons in R- geometry may not contain other geometry

Problem with clipping northern Canada polygons in R- geometry may not contain other geometry

I'm making a map in R. I get some intense tearing of polygons in northern Canada. It appears the issue is with this function gClip and the structure of the geometries in northern areas of Canada. How can I clip the geometry to avoid tearing and also show northern Cborders in my map?

library(geodata)
library(raster)
library(rgeos)
library(terra)
library(sp)
library(sf)

provinces <- c("British Columbia", "Alberta", "Saskatchewan", "Manitoba", "Ontario", "Yukon", "Northwest Territories","Nunavut")
    
# download from GADM
canada <- getData("GADM",country="CAN",level=1, path = ".")
    
# select named places
ca.provinces <- canada[canada$NAME_1 %in% provinces,]
    
# change everything to sf
 can<-st_as_sf(ca.provinces)
    
 can<-st_transform(can,4326)
 can<-as_Spatial(can)
    
 # canada polygons need to be clipped. if the polygons are outside the extent of the background raster it will create a distorted map
 # second argument in function is bounding box
    
 gClip <- function(shp, bb){
      if (class(bb) == "matrix") b_poly <- as(extent(as.vector(t(bb))), "SpatialPolygons")
      else b_poly <- as(extent(bb), "SpatialPolygons")
      proj4string(b_poly) <- proj4string(shp)
      gIntersection(shp, b_poly, byid = T)
 }
    
 # clip the provinces to the extent  
 can <- gClip(can, matrix(c(-150, -95,35, 65)))
    
 # polygon df
 can <- fortify(can)

I get this error.

> # clip the provinces to the extent  
> can <- gClip(can, matrix(c(-150, -95,35, 55)))
output subgeometry 3, row.name: 13 1
subsubgeometry 0: Polygon
subsubgeometry 1: Point
Error in RGEOSBinTopoFunc(spgeom1, spgeom2, byid, id, drop_lower_td, unaryUnion_if_byid_false,  : 
  Geometry collections may not contain other geometry collections
In addition: Warning messages:
1: In if (class(bb) == "matrix") b_poly <- as(extent(as.vector(t(bb))),  :
  the condition has length > 1 and only the first element will be used
2: In proj4string(shp) : CRS object has comment, which is lost in output
3: In RGEOSBinTopoFunc(spgeom1, spgeom2, byid, id, drop_lower_td, unaryUnion_if_byid_false,  :
 
 Error in RGEOSBinTopoFunc(spgeom1, spgeom2, byid, id, drop_lower_td, unaryUnion_if_byid_false,  : 
  Geometry collections may not contain other geometry collections

Any help with a work around would be greatly apreciated.

Solution

You are using old and deprecated packages. You can do this with "sf" or "terra" instead. For example

library(terra)
library(geodata)
provinces <- c("British Columbia", "Alberta", "Saskatchewan", "Manitoba", "Ontario", "Yukon", "Northwest Territories","Nunavut")
canada <- geodata::gadm("CAN", level=1, path = ".")
ca.provinces <- canada[canada$NAME_1 %in% provinces,]

can <- crop(canada, c(-150, -95,35, 65))