Sf: Feature request: st_dissolve(x, by = 'field')

[wgrundlingh]

There currently doesn't seem to be a dissolve feature as part of sf.

Here's a view of the menu option in QGIS:

Here's a visual of the operation:

Related references:

• Dissolve polygons in R
• Polygon Dissolve Operations using R
• ArcMap help on dissolve operation

[tim-salabim]

I wanted to raise something similar so I am adding to this.
With the latest mapview develop commit we can do:

library(mapview)
library(sf)
library(dplyr)

franconia %>% 
  group_by(district)

which creates a grouped_df that should be st_union-able to what is shown above with

franconia %>% 
  group_by(district) %>%
  st_union() %>%
  mapview(zcol = "district")

but what we gt is a union of all polygons instead of a grouped union.

[edzer]

franconia %>% 
  group_by(district) %>% 
  summarise(m = mean(SHAPE_LEN)) %>% 
  st_cast() %>% 
  mapview(zcol = "m")

[tim-salabim]

Cool! We obviously cannot infer how to summarise attributes automatically but shouldn't something like the above (group_by %>% st_union) behave more like

franconia %>%
  split(.$district) %>%
  lapply(st_union) %>%
  mapview()

but preserve the district attribute and obviously be a sf / data.frame rather than a list (which I couldn't figure out how to achieve not having much time before work)?

[edzer]

In your example, st_union doesn't know what to do with non-geometry columns, so they get dropped, and nothing binds the list together. This works:

franconia %>%  
  split(.$district) %>% 
  lapply(st_union) %>% 
  do.call(c, .) %>% # bind the list element to a single sfc
  st_cast() %>% # mapview doesn't like GEOMETRY -> cast to MULTIPOLYGON
  mapview()

[edzer]

@wgrundlingh : please remove the copyrighted figure from your question, but retain a link to it.

[tim-salabim]

Thanks for the solution @edzer. I was trying with do.call(rbind, .) rather than c. And good catch here with the GEOMETRY. I will look into this.

Sorry for side-tracking this issue a little.

[edzer]

I think this exactly addresses the question.

I was also trying rbind, but that combines rows, not elements (geometries).

I'm still somewhat torn whether to put st_cast() for instance inside st_union, or assume downstream pkgs to take care of it.

[tim-salabim]

I will investigate how I can handle this best in mapview and report back. Though in st_read things are cast by default, aren't they?

[edzer]

In the sense of {single,multi} become {multi,multi}, by default yes; not in unfolding GEOMETRYCOLLECTIONs (which rarely come from st_read) - sth st_cast also does.

[tim-salabim]

Ok, the {single,multi} is exactly the problem above (mixed POLYGON and MULTIPOLYGON). I think I have mixed types (mix of *POINT, *LINE, *POLYGON) under control by splitting them by st_dimension and then plotting one layer for each type. Though this is not what the above example should do, I think, as the resulting geometries from st_union all originate from the same layer (all MULTIPOLYGONS).

[wgrundlingh]

@edzer While I don't see anything wrong in using that image here with a link to the original source, I've updated my original post with my own creation showing the dissolve operation.

[edzer]

The page you copied from has, at the end, this message: "Copyright © Environmental Systems Research Institute, Inc.". You _copied_ content from that page without confirming that you had obtained rights to do so, so I assumed that was a copyright infringement, even if you attribute the source. _Linking_ to such a page is something different.

The figure you replaced it with, did you create that yourself?

[wgrundlingh]

@edzer The use of the image falls under _fair use_ and its intent is for educational / non-profit / explanatory / research use, without taking anything away from the original copyright owner:

"...the fair use of a copyrighted work, including such use by reproduction in copies or phonorecords or by any other means specified by that section, for purposes such as criticism, comment, news reporting, teaching (including multiple copies for classroom use), scholarship, or research, is not an infringement of copyright."

Source: 17 USC Section 107

The image currently displayed in my post is one I created.

[edzer]

I think that's up to the copyright holder, but thank you in any case.

[edzer]

@tim-salabim looks like the need to st_cast() before mapview() was caused by another bug in sf...

[wgrundlingh]

@edzer Excuse my confusion here, but I'm still not clear on how this could be implemented using only sf. From the closure I assume this is not something you'll consider implementing as part of sf?

Here's an example of what I'm trying to do (which also replicates the updated image I posted in my original request):

• Download Canada Census Subdivisions from Statistics Canada's website.
• Dissolve under the field ERUID

Using

library(sf)
csd <- st_read('gcsd000b11a_e.shp') # Read Census Subdivision (CSD) shape file
temp <- csd %>% group_by(ERUID) %>% st_union()

leaves me with this error:

gcsd000b11a_e.shp is a MULTIPOLYGON. Is that because of the geography files or my use of st_union or something else altogether?

[tiernanmartin]

@wgrundlingh - is there a specific need for an sf verb here?

What would it provide that the group_by() %>% summarize() pattern doesn't?

library(sf) 
library(tidyverse) 
library(purrr)
library(magrittr) 

nc <- st_read(system.file("shape/nc.shp", package="sf")) 

nc_grp <-
    nc %>%
    mutate(MAX_X = map_dbl(geometry, ~ st_bbox(.x) %>% extract2("xmax")), # extract the bounding box's max x coordinate
           GRP = factor(ntile(MAX_X, 3)),
           STATUS = 'distinct')

nc_grp_distinct <-
    nc_grp %>%
    group_by(NAME) %>%
    summarise(GRP = first(GRP),
              STATUS = 'distinct')  %>% # make a facetting variable for ggplot
    ungroup() %>%
    select(-NAME)

nc_grp_dissolved <-
    nc_grp %>%
    group_by(GRP) %>%
    summarise(STATUS = 'dissolved') # make a facetting variable for ggplot

nc_combined <- rbind(nc_grp_distinct, nc_grp_dissolved)

# Have a look
ggplot(data = nc_combined) +
    ggplot2::geom_sf(aes(fill = GRP)) +
    facet_wrap(~STATUS, ncol = 1)

[edzer]

@wgrundlingh you still need to look into how split-apply-combine works; either by hand, as I explained above (using split, lapply, c) or by dplyr, as @tiernanmartin explains. This worked for me:

temp <- csd %>% group_by(ERUID) %>% summarize()

giving

of which I of course have no idea whether it makes sense.

The error you see refers to an invalid geometry; it's quite suspect that I don't see it while working on the same file. You may want to try to update sf and retry, run st_make_valid (#67) on the object and retry, or if st_make_valid is not available to you try to run st_buffer(csd, 0) on it instead.