Cortar objeto de recursos simples em R

20

Existe uma função para cortar o objeto sf map, semelhante ao maptools::pruneMap(lines, xlim= c(4, 10), ylim= c(10, 15))usado para SpatialPolygon ou SpatialLine?

Estou considerando, st_intersection()mas pode haver maneira adequada.

Kazuhito
fonte

Respostas:

17

st_intersectioné provavelmente o melhor caminho. Encontre a melhor maneira de fazer com que um sfobjeto se cruze com sua entrada. Aqui está uma maneira de usar a conveniência raster::extente uma mistura de antigo e novo. ncé criado por example(st_read):

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

Eu não acho que você pode convencer st_intersectiona não precisar de um CRS correspondente exato; portanto, defina-os como NA ou verifique se são iguais. Não há ferramentas fáceis para o bbox / extension afaik, portanto, usar o raster é uma boa maneira de facilitar as coisas.

mdsumner
fonte
Muito obrigado @mdsumner, funcionou como um encanto. Passei horas, st_intersectionmas não consegui resolver sozinho.
Kazuhito
Agora você pode usar spex::spexpara substituir a st_as_sf(as(...))chamada. Além disso, tmaptools::crop_shape()pode fazer isso.
AF7 5/0318
1
sfagora inclui st_crop, veja minha resposta para detalhes.
AF7 23/04/19
23

Desde hoje , existe uma st_cropfunção na versão github de sf( devtools::install_github("r-spatial/sf"), provavelmente também no CRAN no futuro próximo).

Basta emitir:

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

O vetor deve ser nomeado com xmin xmax ymin ymax(em qualquer ordem).

Você também pode usar qualquer objeto que possa ser lido st_bboxcomo limites de corte, o que é muito útil.

AF7
fonte
5

Outra solução alternativa, para mim, foi mais rápida para arquivos de forma maiores:

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
fonte
Obrigado. É um fluxo de trabalho interessante, combinação de raster :: crop () e st_as_sf () ... + 1 de mim. Eu gostaria que pudéssemos ter esse tipo de função facilmente acessível como crop () em versões futuras do sf . Quanto à velocidade, uma execução rápida de system.time com sua função relatou usuário: 5,42, sistema: 0,09, decorrido 5,52 , enquanto a st_intersection()abordagem foi user: 1,18, sistema: 0,05, decorrido 1,23 no seu conjunto de dados. (Provavelmente o meu ambiente é diferente com o seu ... não tenho certeza.)
Kazuhito
Isso é interessante - a abordagem st_intersection leva cerca de 80 anos para mim.
Pbaylis
Lembre-se de que a função raster :: crop, quando aplicada a objetos de geometria sp, atua como um invólucro para funções rgeos. Embora, um invólucro muito conveniente. A API do GEOS opera em objetos WKT, portanto, será invariavelmente um padrão para operações de sobreposição de sf.
Jeffrey Evans
1
E isso muda com o tempo. Agora, o sf possui "indexação espacial" embutida em 0.5-1 cran.r-project.org/web/packages/sf/news.html, portanto, provavelmente agora é mais rápido que o sp / rgeos.
Mdsumner 03/07
1
sfagora inclui st_crop, veja minha resposta para detalhes.
AF7 23/04/19
1

A solução do @ mdsumner como uma função. Funciona se rastafor um RasterBrick, extensão, bbox, etc.

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

Ele joga fora as informações crs da varredura, porque não sei como converter uma varredura crs () em st_crs ()

Na minha máquina e para minha amostra de dados, isso tem desempenho equivalente ao de raster::cropuma versão SpatialLinesDataFrame dos dados.

A solução do @ pbaylis é cerca de 2,5 vezes mais lenta:

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

Edit: Somebodies comment sugere spex , que produz SpatialPolygons com os crs do rasta, se ele tiver um crs.

Este código usa o mesmo método que 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
fonte
sf agora tem uma st_cropfunção que provavelmente vale a pena conferir.
cmc