CODEWITHSUNDEEP
rr-sf
CzechInk
CzechInk

Reputation: 473

Make sf object out of polygons with holes and set crs

With contourLines() I have extracted the 95% contour for my data. I would like to make an sf object with the correct crs. While I can't share my actual data set, I've adapted an example from this SO post to illustrate where I am stuck.

The main issue is that there is a hole in one of my polygons, but I can't figure out how to make an sf object where this hole is recognized and where I can specify the correct crs.

I am still fairly new to the sf-package and haven't been able to figure this out. Below is what I've tried so far. I've updated the output from contourLines() so that it looks like pts, where each list element is a matrix of point coordinates for a polygon. Using st_polygon() will remove the hole... but then I can't specify a 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")

image of pl1

Alternatively, I could try make each list element a polygon and specify the correct crs... but then I can't figure out how to remove the hole:

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")

image of pl2

Upvotes: 3

Views: 1612

Answers (1)

nniloc
nniloc

Reputation: 4243

In this example they use st_polygon then st_sfc. This seems to work

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')

enter image description here

Upvotes: 4

Related Questions