是否存在用于裁剪sf映射对象的函数，类似于用于SpatialPolygon或SpatialLine的maptools::pruneMap(lines, xlim= c(4, 10), ylim= c(10, 15))？

我正在考虑st_intersection()，但可以采取适当的方法。

#1 楼

st_intersection可能是最好的方法。寻找使sf对象与您的输入相交的最佳方法。这是使用raster::extent的便利以及新旧混合的一种方式。 nc是由example(st_read)创建的：

st_intersection(nc, st_set_crs(st_as_sf(as(raster::extent(-82, -80, 35, 36), "SpatialPolygons")), st_crs(nc)))


我不认为您可以哄骗st_intersection不需要完全匹配的CRS，因此请将它们都设置为NA或确保它们是相同的。没有用于bbox / extent afaik的简便工具，因此使用栅格是使事情变得简单的好方法。

评论

非常感谢@mdsumner，它就像一个魅力。我在st_intersection上花了几个小时，但我自己无法解决。

– Kazuhito
17 Mar 6 '17 at 11:48
现在，您可以使用spex :: spex替换st_as_sf（as（...））调用。另外，tmaptools :: crop_shape（）可以做到这一点。

–AF7
18 Mar 5 '18 at 10:36
旧金山现在包括st_crop，有关详细信息，请参见我的答案。

–AF7
18年4月23日在7:34

#2 楼

从今天开始，在st_crop的github版本中有一个sf函数（devtools::install_github("r-spatial/sf")，不久的将来也可能在CRAN上）。

只是问题：

st_crop(nc, c(xmin=-82, xmax=-80, ymin=35, ymax=36))


向量必须以xmin xmax ymin ymax命名（以任意顺序）。

您也可以使用st_bbox可以读取的任何对象作为裁剪限制，这非常方便。

#3 楼

另一个解决方法，对我来说，较大的shapefile更快：

library(sf)
library(raster)
library(rgeos)
library(ggplot2)

# Load National Forest shapefile
# https://data.fs.usda.gov/geodata/edw/edw_resources/shp/S_USA.AdministrativeForest.zip
nf.poly <- st_read("S_USA.AdministrativeForest"), "S_USA.AdministrativeForest")

crop_custom <- function(poly.sf) {
  poly.sp <- as(poly.sf, "Spatial")
  poly.sp.crop <- crop(poly.sp, extent(c(-82, -80, 35, 36)))
  st_as_sf(poly.sp.crop)
}

cropped <- crop_custom(nf.poly)


评论

谢谢。这是有趣的工作流程，是我的raster :: crop（）和st_as_sf（）... + 1的组合。我希望我们在将来的sf版本中可以拥有这类易于访问的函数，例如crop（）。至于速度，在您的数据集上，以您的函数报告的用户：5.42，系统：0.09，运行时间5.52快速运行system.time，而用户使用st_intersection（）方法：1.18，系统：0.05，运行时间1.23。 （可能我的环境与您的环境有所不同...不确定。）

– Kazuhito
17年4月1日在7:33


这很有趣-st_intersection方法对我来说大约需要80秒钟。

– pbaylis
17年4月3日，下午3:36
请记住，将raster :: crop函数应用于sp几何对象时，它充当rgeos函数的包装器。虽然，这是一个非常方便的包装器。 GEOS API在WKT对象上运行，因此，它将始终是sf覆盖操作的标准。

–杰弗里·埃文斯（Jeffrey Evans）
17年4月14日在15:46
而且它会随着时间而变化，SF现在在0.5-1 cran.r-project.org/web/packages/sf/news.html中内置了“空间索引”，因此现在可能比sp / rgeos更快。

–mdsumner
17年7月3日在10:16
旧金山现在包括st_crop，有关详细信息，请参见我的答案。

–AF7
18年4月23日在7:34

#4 楼

@mdsumner的解决方案作为函数。如果rasta是RasterBrick，范围，bbox等，则可以使用。

# Crop a Simple Features Data Frame to the extent of a raster
crop.sf = function(sfdf, rasta) {
  st_intersection(sfdf, st_set_crs(st_as_sf(as(extent(rasta), "SpatialPolygons")), st_crs(sfdf)))
}


它丢弃了栅格的crs信息，因为我不知道如何转换栅格将crs（）转换为st_crs（）

在我的计算机上以及对于我的数据示例，具有与SpatialLinesDataFrame版本的数据一样的raster::crop性能。

@ pbaylis的解决方案是慢约2.5倍：

# Slower option.
crop.sf2 = function(sfdf, rasta) {
  st_as_sf(crop(as(sfdf, "Spatial"), rasta))
}


编辑：有人评论建议spex，如果有crs，它会从rasta中产生crs，产生SpatialPolygons。

此代码使用与spex相同的方法：

# Crop a Simple Features Data Frame to the extent of a raster
crop.sf3 <- function(sfdf, rasta) {
  # Get extent and crs
  ext.sp <- as(extent(rasta), "SpatialPolygons")
  crs(ext.sp) <- crs(rasta)

  # crop
  st_intersection(sfdf, st_as_sf(ext.sp))
}


评论

sf现在有一个st_crop函数，可能值得一试。

– cmc
19年7月23日在16:03