Function to aggregate polygons

[mpadge]

This will likely be a very useful post-processing function for cases when the package is used to identify a whole lot of polygons purportedly representing the same thing/area, and this mess of polygons need to be converted to some kind of heat map. A raster environment will provide likely the most straightforward method of aggregation, and could be subsequently converted to vector contours. Alternatively, it should also be possible to do sequential polygon intersections, and aggregate increasing numbers of tiny intersected polygons.

@mdsumner Any thoughts on this: Simply put, consider a Venn diagram of 3 circles, each represented as a spatial polygon. I want a function that gives a value of 3 in the middle, 2 for each overlap of two polygons, 1 in the outer polygon reaches, and 0 for background. Any insights?

[mdsumner]

Spacebucket!

remotes::install_github("mdsumner/spacebucket")

  library(spacebucket)
library(sf)
#> Linking to GEOS 3.7.0, GDAL 2.4.0, PROJ 5.2.0
x <- st_sf(a = 1:3, geometry = st_sfc(list(st_point(cbind(0, 0)), st_point(cbind(0, 1)), st_point(cbind(1, 0)))))
bx <- st_buffer(x, 0.75)
bucket <- do.call(spacebucket, split(bx, seq_len(nrow(bx))))

plot(bx)

## get various n-intersections
sf_geom <- function(x) st_geometry(st_union(x))
plot(st_geometry(bx), reset = FALSE)
bucket %>% n_intersections(3) %>% sf_geom() %>% plot(add = TRUE, lwd = 5)

plot(st_geometry(bx), reset = FALSE)
bucket %>% n_intersections(2) %>% sf_geom() %>% plot(add = TRUE, lwd = 5)

plot(st_geometry(bx), reset = FALSE)
bucket %>% n_intersections(1) %>% sf_geom() %>% plot(add = TRUE, lwd = 5)

^{Created on 2019-07-04 by the reprex package (v0.3.0)}

(gad, sorry 50 reproducible-edits later ...)

I have to refamiliarize with how the inputs are recorded in the list-column, but n_intersections is a helper to get out the underlying pieces.

[mdsumner]

## how many n-pieces? 
n_types <- bucket$index %>% group_by(path_) %>% n_groups()
sf_df <- function(x, n) st_sf(n = n, geometry = st_union(x))

## note that these are now overlapping polygons, with a record of n-pieces, 
## so we order in increasing n for the plot
a <- do.call(rbind, purrr::map(1:n_types, ~{bucket %>% n_intersections(.x) %>% sf_df(n = .x)}))
plot(a)

ggplot(a) + geom_sf(aes(colour = as.factor(n))) + facet_wrap(~n)

Along the right track? We literally have an index in the bucket for every single fragment, but it's stored as densely as I can understand - pulling this summary out could be much more efficient but it'll take some thought.

[mpadge]

Ah, you never let me down! That's awesome. Only problem at the outset would be spacebucket's absence from cran...

[mdsumner]

Heh, it's these two lines - neither too hard to arrange

  path <- silicate::PATH(sfall)
  RTri <- pfft::edge_RTriangle(path)

hmm, got an idea ...

[mpadge]

It's the middle of the night, dude!

[mdsumner]

It is late but I've been watching Happy Valley bits of S01 I missed, so good.

What I've lost sight of is why sb gives >= n and this seems to give exclusive == n ... but at any rate this is pretty swish and CRAN-friendly ;)

library(sf)
#> Linking to GEOS 3.7.0, GDAL 2.4.0, PROJ 5.2.0
library(sfdct)
library(dplyr)
#> 
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag
#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union
  library(ggplot2)
x <- st_sf(a = 1:3, geometry = st_sfc(list(st_point(cbind(0, 0)), st_point(cbind(0, 1)), st_point(cbind(1, 0)))))
bx <- st_buffer(x, 0.75)

mesh <- st_cast(ct_triangulate(st_union(st_cast(bx, "MULTILINESTRING"))))
mesh_centroids <- st_centroid(st_cast(mesh))


intscts <- tibble::tibble(triangle = unlist(st_intersects(bx, mesh_centroids))) %>% 
  mutate(feature = purrr::map(st_intersects(mesh_centroids[triangle], bx), 
                              ~tibble::tibble(feature = .x)))

a <- do.call(rbind, intscts %>% tidyr::unnest() %>% group_by(triangle, feature) %>% mutate(n = n()) %>% 
          ungroup() %>% group_by(n) %>% split(.$n) %>% 
          purrr::map(~st_sf(geometry = st_union(mesh[.x$triangle]), n = .x$n[1]))) 
#> Warning: `cols` is now required.
#> Please use `cols = c(feature)`

plot(a)

ggplot(a) + geom_sf() + facet_wrap(~n)

^{Created on 2019-07-04 by the reprex package (v0.3.0)}

[mpadge]

That's nearly perfect! May I then make a request of you? Do you wanna do a direct push or PR (you should have invite already, and I'll admin you) for basically this code in order to to:

1. convert a collection of polygons to your a object above?
2. Extract contours as polygons -- not multipolygons -- for defined levels of aggregation?

The latter requires some slight change, because your multipolygon objects exclude the higher-level polygons contained within, while the desired polygon objects (in your example, for n = 2) should include the higher levels (n = 3). That just means it'll need some kind of sequential st_union operation along the n values. The only real change to the former would be that there is no purrr dep at present, and no real need for that one function, so slight de-purrring would also be good there.

That'd be awesome, and then we'd be equal co-authors of the whole shebang, which I'm hoping to have ready for an rOS sub within a coupla weeks. This could then be a proper joint effort between us, for which I'd be mighty appreciative!

[mpadge]

Thanks loads for the code Mike, but ... but ... but ... it's not scalable:

Standard usage according to the peeps for whom I'm doing this little project will be n~O(several dozen) or even hundreds, so that would mean processing time of up to half an hour or so 😭 Although sfdct absolutely smokes, the killer is the st_intersects call. So I've got another idea: Just bundle clipper to access it's point-in-polygon routine and do the following, noting that clipper for a polygon layer returns a vector with index for each polygon in layer:

# initialisation:
for (p in polygons)
  for (i in npts (p))
    p [i].count = 1
# count is an extra variable appended to each vertex, starting at 1 for self-count

for (p in polygons) {
  p_rest = clipper::polygon_layer (polygons [-p]) # all polygons as one layer
  for (i in npts (p)) {
    p[i].count += sum (point_in_polygon (p_rest)) # +extra code to not count clipper's "-1" for *on*
  }
}

Each vertex then just counts the number of polygons it appears within. The the vertices just need to be grouped by count, and polygons re-traced around them. Simple! No triangulation necessary, all can be done by simple point-in-polygon calls.

The only trick: The polygon tracing. That amounts to then requiring ... convex polygon tracing at the limit of maximum convexity. But I've messed around with that enough to realise how simple this whole game really is, and I'm sure I can just code up my own, custom built and highly efficient version. Then we'd have it.

Edit: And by the way, this package has ended up with some pretty generally useful routines for raster analysis, all clearly named in src: edge thinning (illustrated in #5), component identification, and boundary tracing. We should also keep those in mind for a repo somewhere to hold these general raster-spatial kind of things - sorta like the beginnings of an spdep for raster. Just a thought.

alternative

Just bundle clipper to replace sf::st_intersects, but that might not speed things up too much (see here), so may still end up needing the above approach?

[mdsumner]

Definitely will be better with pfft and silicate, I'll try your benchmark and see if we can avoid sfdct

[mdsumner]

This is related https://github.com/slu-openGIS/areal