TransWikia.com

Creating multiple polygons from a list of coordinate matrices using sf in R

Geographic Information Systems Asked by jfalter on July 13, 2021

I’m trying to create an sf object representing the extents of geometry objects in an existing sf object. I would expect that st_polygon() could create a list of polygon objects directly from a list of coordinate matrices, but I can’t do it without introducing a clunky for loop.

Is there a faster and/or more elegant way of doing this?

st_bbox_to_coords <- function(bbox) {

    bbox_num <- as.numeric(bbox)
    coords_mat <- matrix(bbox_num[c(1, 1, 3, 3, 1, 2, 4, 4, 2, 2)], ncol = 2)
    return(coords_mat)
}

create_sf_from_extent <- function(sf_object) {

    bbox_list <- lapply(st_geometry(sf_object), st_bbox)
    bbox_coords_mat_list <- lapply(bbox_list, st_bbox_to_coords)

    # What should work
    # bbox_sfg_list <- lapply(bbox_coords_mat_list, st_polygon)

    # What actually works
    dummy_list <- bbox_coords_mat_list
    for (i in seq_along(bbox_coords_list)) { 
        dummy_list[[i]] = bbox_coords_mat_list[i]
    }
    bbox_sfg_list <- lapply(dummy_list, st_polygon)
    bbox_sfc <- st_sfc(bbox_sfg_list)
    bbox_sf <- sf_object
    st_geometry(bbox_sf) <- bbox_sfc

    return(bbox_sf)
}

One Answer

When generatings POLYGON from bboxes you should use st_as_sfc. I added an example, not able to say if it is more "elegant" or less clunky. Maybe the ugliest part is how to go from a list of POLYGON to an actual sfc object.

library(sf)
myshape <- st_read(system.file("shape/nc.shp", package = "sf"))
bbox_list <- lapply(st_geometry(myshape), st_bbox)
polys_list <- lapply(bbox_list, st_as_sfc)
endpol = polys_list[[1]]
for (i in 2:length(polys_list)) {
  endpol <- c(endpol, polys_list[[i]])
}
st_crs(endpol) <- st_crs(myshape)
endpol
#> Geometry set for 100 features 
#> geometry type:  POLYGON
#> dimension:      XY
#> bbox:           xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
#> epsg (SRID):    4267
#> proj4string:    +proj=longlat +datum=NAD27 +no_defs
#> First 5 geometries:
#> POLYGON ((-81.74107 36.23436, -81.23989 36.2343...
#> POLYGON ((-81.34754 36.36536, -80.90344 36.3653...
#> POLYGON ((-80.96577 36.23388, -80.43531 36.2338...
#> POLYGON ((-76.33025 36.07282, -75.77316 36.0728...
#> POLYGON ((-77.90121 36.16277, -77.07531 36.1627...

#Plots
plot(st_geometry(myshape))
plot(endpol, add=TRUE, border="red")

Created on 2020-02-22 by the reprex package (v0.3.0)

Answered by dieghernan on July 13, 2021

Add your own answers!

Ask a Question

Get help from others!

© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP