裁剪SpatialPolygonsDataFrame

4 回答

• 5

  Streamlined method added 2014-10-9 ：

  raster::crop() 可用于裁剪 Spatial* （以及 Raster* ）个对象 .

  例如，以下是如何使用它来裁剪 SpatialPolygons* 对象：

  ## Load raster package and an example SpatialPolygonsDataFrame
  library(raster) 
  data("wrld_simpl", package="maptools")
  
  ## Crop to the desired extent, then plot
  out <- crop(wrld_simpl, extent(130, 180, 40, 70))
  plot(out, col="khaki", bg="azure2")
  

  ---

  Original (and still functional) answer:

  rgeos 函数 gIntersection() 使这非常简单 .

  使用mnel的漂亮示例作为跳跃点：

  library(maptools)
  library(raster)   ## To convert an "Extent" object to a "SpatialPolygons" object.
  library(rgeos)
  data(wrld_simpl)
  
  ## Create the clipping polygon
  CP <- as(extent(130, 180, 40, 70), "SpatialPolygons")
  proj4string(CP) <- CRS(proj4string(wrld_simpl))
  
  ## Clip the map
  out <- gIntersection(wrld_simpl, CP, byid=TRUE)
  
  ## Plot the output
  plot(out, col="khaki", bg="azure2")
  

  

  回复于 2024-05-05T05:04:07+08:00
• 1

  以下是使用世界 Map 作为示例如何使用 rgeos 执行此操作的示例

  这来自Roger Bivand，发表于R-sig-Geo mailing list . Roger是 sp 包的作者之一 .

  以世界 Map 为例

  library(maptools)
  
  data(wrld_simpl)
  
  # interested in the arealy bounded by the following rectangle
  # rect(130, 40, 180, 70)
  
  library(rgeos)
  # create  a polygon that defines the boundary
  bnds <- cbind(x=c(130, 130, 180, 180, 130), y=c(40, 70, 70, 40, 40))
  # convert to a spatial polygons object with the same CRS
  SP <- SpatialPolygons(list(Polygons(list(Polygon(bnds)), "1")),
  proj4string=CRS(proj4string(wrld_simpl)))
  # find the intersection with the original SPDF
  gI <- gIntersects(wrld_simpl, SP, byid=TRUE)
  # create the new spatial polygons object.
  out <- vector(mode="list", length=length(which(gI)))
  ii <- 1
  for (i in seq(along=gI)) if (gI[i]) {
    out[[ii]] <- gIntersection(wrld_simpl[i,], SP)
    row.names(out[[ii]]) <- row.names(wrld_simpl)[i]; ii <- ii+1
  }
  # use rbind.SpatialPolygons method to combine into a new object.
  out1 <- do.call("rbind", out)
  # look here is Eastern Russia and a bit of Japan and China.
  plot(out1, col = "khaki", bg = "azure2")
  

  

  回复于 2024-05-05T05:04:07+08:00
• 0

  您不能在sp多边形对象上使用裁剪 . 您将需要创建一个表示dat2的bbox坐标的多边形，然后可以在rgeos库中使用gIntersects .

  编辑：此评论与2012年提供的版本有关，现在已不再适用 .

  回复于 2024-05-05T05:04:07+08:00
• 54

  看？作物

  corp（x，y，filename =“”，snap ='near'，datatype = NULL，...）x Raster * object y Extent对象，或从中提取Extent对象的任何对象（请参阅Details

  您需要使用栅格包中的 rasterize 函数栅格化第一个SpatialPolygon

  我创建了一些数据来展示如何使用栅格化：

  n <- 1000
  x <- runif(n) * 360 - 180
  y <- runif(n) * 180 - 90
  xy <- cbind(x, y)
  vals <- 1:n
  p1 <- data.frame(xy, name=vals)
  p2 <- data.frame(xy, name=vals)
  coordinates(p1) <- ~x+y
  coordinates(p2) <- ~x+y
  

  如果我尝试：

  crop(p1,p2)
   unable to find an inherited method for function ‘crop’ for signature ‘"SpatialPointsDataFrame"’
  

  现在使用栅格化

  r <- rasterize(p1, r, 'name', fun=min)
  crop(r,p2)
  
  class       : RasterLayer 
  dimensions  : 18, 36, 648  (nrow, ncol, ncell)
  resolution  : 10, 10  (x, y)
  extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
  coord. ref. : +proj=longlat +datum=WGS84 
  data source : in memory
  names       : layer 
  values      : 1, 997  (min, max)
  

  回复于 2024-05-05T05:04:07+08:00