I'm trying to create a multipolygon dataset where each row is identifiable.

Example below


# example raster data
r <- brick(nl = 2, nrows = 5, ncols = 10)
set.seed(54); r[] <- rpois(ncell(r)*nlayers(r), lambda = 3)
names(r) <- c("layer_1", "layer_2")

for (i in 1:nlayers(r)) {
  s <- r[[i]]
  # set values below x to NA, all else to 1
  s[s[] < quantile(r[[i]], 0.9)] <- NA
  s[!is.na(s[])] <- 1
  # convert from raster to multipolygon
  s <- st_union(as(rasterToPolygons(s, na.rm = TRUE), "sf"))
  # send to workspace
  ## ideally here I could bind to an existing empty dataset
  ## after naming to rows
  assign(x = paste0("poly", i), value = s)
# no way to tell which one is poly1 or poly2?
poly3 <- c(
  st_cast(poly1, "GEOMETRY"),
  st_cast(poly2, "GEOMETRY"))
names(poly3) <- names(r)
write_sf(poly3, "./test_poly_output.shp") ## works
rownames(poly3) <- names(r) ## error
poly3$name <- names(r)
write_sf(poly3, "./test_poly_output_2.shp") ## error

I want to be able to open the dataset in QGIS or Arc, and see which rows correspond to the appropriate layer. Opening the test_poly_output.shp shows the following which isn't particularly helpful, even though in this example FID_1 == 0 obviously corresponds with layer 1 enter image description here

Ideally I would be able to bind to an existing empty sfc class. The code below works for a MULTILINESTRING, but I cannot get it to work for MULTIPOLYGON or GEOMETRY types.

## Example where LINESTRINGS can be appended to blank sfc
rConts <- st_sf(st_sfc(crs = 4326))

for (i in 1:nlayers(r)) {
  contsR <- r[[i]]
  contsR <- rasterToContour(r[[i]], levels = quantile(r[[i]], 0.9))
  names(contsR) <- "Q90"
  contsR <- as(contsR, "sf")
  rownames(contsR) <- gsub("X_", "", names(r[[i]]))
  contsR$Year <- as.integer(gsub("layer_", "", names(r[[i]])))
  rConts <- rbind(rConts, contsR)
write_sf(rConts, "./test_contourLines_output.shp")

enter image description here

1 Answer 1


Your terminology is bit unclear but I think it boils down to you trying to save spatial features vectors and expecting each element to have some attributes. You can't do that with vectors, you have to make spatial features data frames.

If you have an object that is an sfc class:

> class(poly3)
[1] "sfc_GEOMETRY" "sfc"         

then saving it as a Shapefile will only save the geometry. To add extra information to each element of that geometry vector, create a spatial features data frame:

> sp3 = st_as_sf(data.frame(name=names(r), geometry=poly3))
> sp3
Simple feature collection with 2 features and 1 field
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: -180 ymin: -90 xmax: 180 ymax: 90
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs
     name                       geometry
1 layer_1 MULTIPOLYGON (((-180 18, -1...
2 layer_2 MULTIPOLYGON (((144 -54, 18...

Saving that to a Shapefile and loading into a GIS should give you a name column with the name in.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.