Pangkas objek fitur sederhana di R

20

Apakah ada fungsi untuk memotong objek peta sf, mirip dengan yang maptools::pruneMap(lines, xlim= c(4, 10), ylim= c(10, 15))digunakan untuk SpatialPolygon atau SpatialLine?

Saya sedang mempertimbangkan st_intersection()tetapi mungkin ada cara yang tepat.

Kazuhito
sumber

Jawaban:

17

st_intersectionmungkin cara terbaik. Temukan cara apa pun yang terbaik untuk mendapatkan sfobjek yang bersinggungan dengan input Anda. Inilah cara menggunakan kenyamanan raster::extentdan campuran lama dan baru. ncdibuat oleh example(st_read):

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

Saya tidak berpikir Anda bisa membujuk st_intersectionuntuk tidak membutuhkan CRS yang cocok, jadi atur keduanya menjadi NA atau pastikan keduanya sama. Tidak ada alat yang mudah untuk bbox / afaik tingkat, jadi menggunakan raster adalah cara yang baik untuk membuat semuanya menjadi mudah.

mdsumner
sumber
Terima kasih banyak @mdsumner, ini bekerja seperti pesona. Saya menghabiskan waktu berjam-jam st_intersectiontetapi tidak bisa menyelesaikannya sendiri.
Kazuhito
Anda sekarang dapat menggunakan spex::spexuntuk mengganti st_as_sf(as(...))panggilan. Juga, tmaptools::crop_shape()bisa melakukan ini.
AF7
1
sfsekarang termasuk st_crop, lihat jawaban saya untuk detail.
AF7
23

Sejak hari ini , ada st_cropfungsi dalam versi github sf( devtools::install_github("r-spatial/sf"), mungkin pada CRAN dalam waktu dekat juga).

Hanya masalah:

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

Vektor harus dinamai dengan xmin xmax ymin ymax(dalam urutan apa pun).

Anda juga dapat menggunakan objek apa pun yang dapat dibaca st_bboxsebagai batas pemangkasan, yang sangat berguna.

AF7
sumber
5

Solusi lain, bagi saya lebih cepat untuk shapefile yang lebih besar:

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)
pbaylis
sumber
Terima kasih. Ini alur kerja yang menarik, kombinasi raster :: crop () dan st_as_sf () ... +1 dari saya. Saya berharap kita dapat memiliki jenis fungsi yang mudah diakses seperti crop () di versi sf mendatang . Mengenai kecepatan, menjalankan sistem.waktu dengan fungsi Anda melaporkan pengguna: 5,42, sistem: 0,09, berlalu 5,52 , sementara st_intersection()pendekatan adalah pengguna: 1,18, sistem: 0,05, berlalu 1,23 pada dataset Anda. (Mungkin lingkungan saya berbeda dengan lingkungan Anda ... tidak yakin.)
Kazuhito
Itu menarik - pendekatan st_intersection membutuhkan sekitar 80-an untuk saya.
pbaylis
Perlu diingat bahwa fungsi raster :: crop, ketika diterapkan pada objek geometri sp, bertindak sebagai pembungkus untuk fungsi rgeos. Meskipun, pembungkus yang sangat nyaman. API GEOS beroperasi pada objek WKT sehingga, akan selalu menjadi standar untuk operasi overlay sf.
Jeffrey Evans
1
Dan itu berubah dari waktu ke waktu, sekarang sf memiliki "pengindeksan spasial" built-in di 0,5-1 cran.r-project.org/web/packages/sf/news.html jadi mungkin sekarang lebih cepat daripada sp / rgeos.
mdsumner
1
sfsekarang termasuk st_crop, lihat jawaban saya untuk detail.
AF7
1

Solusi @ mdsumner sebagai fungsi. Bekerja jika rastaRasterBrick, luas, bbox, dll.

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

Itu membuang informasi crs dari raster karena saya tidak tahu bagaimana mengubah raster crs () menjadi st_crs ()

Di mesin saya dan untuk sampel data saya ini memiliki kinerja setara raster::cropdengan versi data SpatialLinesDataFrame.

Solusi @ pbaylis sekitar 2,5 kali lebih lambat:

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

Sunting: Komentar Somebodies menyarankan spex , yang menghasilkan SpatialPolygons dengan cr dari rasta, jika memiliki crs.

Kode ini menggunakan metode yang sama dengan 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))
}
cmc
sumber
sf sekarang memiliki st_cropfungsi yang mungkin perlu dicoba.
cmc