TransWikia.com

Converting points to polygons by group

Geographic Information Systems Asked by cschwab98 on April 3, 2021

I have a dataframe with around 40 points. They are grouped by an ID variable (their name), which four points in each group. I’m trying to turn these points into polygons, but I’ve had a lot of trouble with different methods in both sf and sp — currently I have them as a large sf dataframe, but I can’t seem to find any function that would use a grouping variable to create polygons.

I’ve tried using the st_cast method, the general idea being

polypoints <- polypoints %>% group_by("ID") %>% st_cast("MULTIPOINT") %>% st_cast("MULTILINESTRING") %>% st_cast("MULTIPOLYGON")

but the step from multipoint to multiline string leaves the geography empty — what are the steps for a process like this?

2 Answers

Make a reproducible example:

> set.seed(999)
> xy = data.frame(x=runif(24), y=runif(24))
> xy$ID = rep(c(1:6), rep(4,6))
> head(xy)
           x         y ID
1 0.38907138 0.8260703  1
2 0.58306072 0.8195141  1
3 0.09466569 0.5684927  1
4 0.85263123 0.6196068  1
5 0.78674676 0.8308805  2
6 0.11934226 0.4588336  2

Make the data frame into an sf data frame:

> xys = st_as_sf(xy, coords=c("x","y"))

Then aggregate by ID, combine the points and cast to POLYGON, turn the whole thing into an sf data frame. In a one-liner:

> polys = st_sf(
  aggregate(
    xys$geometry,
    list(xys$ID),
    function(g){
       st_cast(st_combine(g),"POLYGON")
    }
   ))

> polys
Simple feature collection with 6 features and 1 field
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: 0.03014558 ymin: 0.01902308 xmax: 0.9875201 ymax: 0.8690149
epsg (SRID):    NA
proj4string:    NA
  Group.1                       geometry
1       1 POLYGON ((0.3890714 0.82607...
2       2 POLYGON ((0.7867468 0.83088...
3       3 POLYGON ((0.3907724 0.52724...
4       4 POLYGON ((0.03014558 0.8162...
5       5 POLYGON ((0.1665847 0.10961...
6       6 POLYGON ((0.9074913 0.35951...

Looks pretty awful but that's random data for you...

enter image description here

Correct answer by Spacedman on April 3, 2021

You need to summarise after your group_by statement then your approach works perfectly fine. Directly from POINTS --> POLYGON, and keeping the crs (if there is one).

set.seed(999)
xy = data.frame(x=runif(24), y=runif(24))
xy$ID = rep(1:6, each = 4)
head(xy)

xys = st_as_sf(xy, coords=c("x","y"))

polys = xys %>% 
  dplyr::group_by(ID) %>% 
  dplyr::summarise() %>%
  st_cast("POLYGON")

plot(polys)

enter image description here

If you want the outer polygon you can add st_convex_hull()

polys = polys %>% 
  st_convex_hull()

plot(polys)

enter image description here

Answered by Latrunculia on April 3, 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