improved st_intersection(polygons) suggestion

[jepol77]

This function seems to work with my data, returning the n-ary intersections of my polygons. Feel free to adapt it.

st_stack <- function(x) {  
  require(dplyr)
  require(sf)
 sf_use_s2(F)
  
  #Function does not work when there are duplicate geometries
  if(length(unique(x$geometry)) != nrow(x))
    stop("Error: Duplicate geometries in input data")
  
  #Faster st_erase function
  st_erase_big = function(x, y){
    x$ref <- 1:nrow(x)
    x_inter = x[apply(st_intersects(x,y), 1, any),] #x in y
    y_inter = y[apply(st_intersects(y,x), 1, any),] #y in x
    difference <- st_difference(x_inter, st_union(y_inter)) #x_inter not in y_inter
    nointer_x = x[!x$ref %in% x_inter$ref,] # x not intersecting
    dplyr::bind_rows(difference,nointer_x)
  } 
  
  #Progress bar
  progress <- function (x, max = 100) {
    percent <- x / max * 100
    cat(sprintf('\r[%-50s] %d%%',
                paste(rep('=', percent / 2), collapse = ''),
                floor(percent)))
    if (x == max)
      cat('\n')
  }
  
a <- st_as_sf(x[1,"geometry"])
a$n_overlaps <- 1

d <- x[,"geometry"]
e <- x[,"geometry"]
e$n <- 1
i <- 2

for(i in 2:nrow(x)){
  progress(match(i, 2:nrow(x)),length(2:nrow(x)))
  b <- st_intersection(a, e[i,]) 
  
  #If no overlaps - just bind them together
  if(nrow(b) == 0){
    a <- plyr::rbind.fill(list(a, d[i,])) %>% 
      mutate(n_overlaps= ifelse(is.na(n_overlaps), 1, n_overlaps)) %>% 
      st_as_sf() %>% 
      st_zm()
   
    next
  }
  
  b$n_overlaps <- b$n_overlaps+b$n
  b$n <- NULL
  
  #Get geometry type
  ty <- st_geometry_type(a, by_geometry = TRUE)
  
  #if some are a collection, only take polygons from those. 
  if(!all(unique(ty) %in% c("POLYGON", "MULTIPOLYGON"))){
    zz <-  a[ty == "GEOMETRYCOLLECTION",] %>% 
      st_collection_extract("POLYGON") 
    a <- rbind(a[ty != "GEOMETRYCOLLECTION",], zz)
  }
  
  #Remove areas from new polygon 
  zz <- st_erase_big(d[i,], a[,"geometry"]) %>% 
    mutate(ref=NULL)

  #Remove areas from old polygons 
  a <- st_erase_big(a, d[i,])%>% 
    mutate(ref=NULL)
 
  #Bind old polygons, new polygon and intersected area
  a <- plyr::rbind.fill(list(a, b, zz)) %>% 
    mutate(n_overlaps= ifelse(is.na(n_overlaps), 1, n_overlaps)) %>% 
    st_as_sf() %>% 
    st_zm()
}
return(a)
}

[edzer]

Thanks; I guess you need this because st_intersection() doesn't do it for your case. Could you share a dataset that fails with st_intersection but succeeds with your code?

[jepol77]

#1668

My solution seems to work without setting precision.

df <- structure(list(lat = c(53.218461547511, 53.2354543814792, 53.2219368310548,
                             53.2142084420887, 53.2051969857393), lon = c(6.57022923569802,
                                                                          6.55005158685787, 6.57015231400986, 6.5633788874466, 6.5692156138519
                             )), row.names = c(NA, -5L), class = c("tbl_df", "tbl", "data.frame"
                             ))

circle_df <- df %>%
  st_as_sf(coords = c("lon", "lat"), crs = 4326) %>%
  st_transform(3035) %>%
  st_buffer(dist = units::set_units(50, "kilometers"))

# Applying st_intersection() gives an error:

circle_df %>%
  st_intersection()

# But the new function does not
c_stack <- st_stack(circle_df)

mapview(c_stack)

[arthurgailes]

Another possible solution. If of interest I can submit a pr:

count_polygon_overlaps <- function(
    shp{
  union_orig <- sf::st_union(shp$geometry)
  
  bounds <- shp$geometry |> sf::st_boundary() |> sf::st_union()
  new_polys <- bounds |> sf::st_polygonize() |> sf:: st_collection_extract('POLYGON')
  new_polys <- new_polys[which(sf::st_geometry_type(new_polys) == 'POLYGON')]
  new_polys <- new_polys[which(
    sf::st_covers(union_orig, new_polys) |> t() |> lengths() > 0)]
  
  new_shp <- sf::st_sf(geometry = new_polys) 
  
  # covers faster than covered_by; better to have index on polygons
  overlap_df <- sf::st_covers(shp$geometry, sf::st_centroid(new_polys)) |> t() |>
    as.data.frame() |>
    dplyr::group_by(row.id) |>
    dplyr::summarize(origins = list(col.id)) |> 
    subset(!is.null(origins))
  
  new_shp[["origins"]] <- rep_len(list(), nrow(new_shp))
  new_shp[["origins"]][overlap_df$row.id] <- overlap_df['origins']
  
  new_shp$overlaps <- lengths(new_shp$origins)
  
  new_shp <- subset(new_shp, overlaps > 0)

  return(new_shp)

}

From the previous example, circle_df2 <- count_polygon_overlaps(circle_df)