在 R 中裁剪加拿大北部多边形时出现问题 - 几何体可能不包含其他几何体

Question

我正在用 R 制作地图。我在加拿大北部发现了一些严重的多边形撕裂。 问题似乎出在 gClip 函数和加拿大北部地区的几何结构上。如何剪切几何图形以避免撕裂并在地图中显示北部边界？

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)

我收到此错误。

> # 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 

任何解决问题的帮助将不胜感激。

Answer 1

您正在使用旧的和已弃用的软件包。您可以使用“sf”或“terra”来代替。例如

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))