Faça o objeto sf de polígonos com buracos e defina crs

Sep 25 2020

Com contourLines()eu extraí o contorno de 95% para meus dados. Eu gostaria de fazer um objeto sf com o crs correto. Embora eu não possa compartilhar meu conjunto de dados real, adaptei um exemplo deste post do SO para ilustrar onde estou preso.

O principal problema é que há um buraco em um dos meus polígonos, mas não consigo descobrir como fazer um objeto sf onde esse buraco é reconhecido e onde posso especificar o crs correto.

Eu ainda sou bastante novo no pacote sf e não fui capaz de descobrir isso. Abaixo está o que tentei até agora. Atualizei a saída de contourLines()para que pareça pts, onde cada elemento da lista é uma matriz de coordenadas de ponto para um polígono. Usando st_polygon()irá remover o buraco ... mas então eu não posso especificar um crs:

library(sf)
library(dplyr)
library(purrr)

# data example
outer1 <- matrix(c(0,0,4,0,4,4,0,4,0,0), ncol = 2, byrow = TRUE)
hole   <- matrix(c(1,1,1,2,2,2,2,1,1,1), ncol = 2, byrow = TRUE)
outer2 <- matrix(c(5,5,5,6,6,6,6,5,5,5), ncol = 2, byrow = TRUE)
pts <- list(outer1, hole, outer2)

# removes hole, but can't add crs 
# - nothing to transform with st_transform() 
# - st_set_crs() throws error:
#  > Error in UseMethod("st_crs<-") : 
#  >  no applicable method for 'st_crs<-' applied to an object of class "c('XY', 'POLYGON', 'sfg')"
pl1 <- st_polygon(pts)
plot(pl1, col = "red")

Como alternativa, eu poderia tentar tornar cada elemento da lista um polígono e especificar o crs correto ... mas não consigo descobrir como remover o buraco:

pl2 <- pts %>% map(function(pts.x) {
  pts.x %>%
    as.data.frame() %>%
    set_colnames(c("x", "y")) %>%
    st_as_sf(coords = c("x", "y"), crs = 32611) %>%
    summarise(geometry = st_combine(geometry)) %>%
    st_cast("POLYGON")
}) %>%
  bind_rows
plot(pl2, col = "red")

Respostas

2 nniloc Sep 25 2020 at 07:00

No presente exemplo, eles usam st_polygonentão st_sfc. Isso parece funcionar

pl1 <- st_polygon(list(outer1, hole)) %>%
  st_sfc(crs = 32611)

#----------
Geometry set for 1 feature 
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: 0 ymin: 0 xmax: 4 ymax: 4
projected CRS:  WGS 84 / UTM zone 11N
POLYGON ((0 0, 4 0, 4 4, 0 4, 0 0), (1 1, 1 2, ...

plot(pl1, col = 'red')