0

I need to do a pair-wise removal of a set of polygons within a larger set of polygons.

pt<-data.frame("x"= c(0,0),"y"= c(0,5),pair= c("a","b"))
pt<- st_as_sf(pt, coords = c("x","y"), dim="XY", crs=  "+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0")

big.c<-st_buffer(pt, dist=7)
big.c$size<-"x"
lil.c<- st_buffer(pt, dist=1)
lil.c$size<-"y"
circ<-as.data.frame(rbind(big.c,lil.c)) %>% 
  pivot_wider(names_from = "size", values_from = "geometry")

Now I want to remove the inner circle (y[1]) from the corresponding big circle (x[1]) by row. If I just use st_difference it will remove any little circles that overlap with big circles, which I don't want. I want donuts.

Secondly, I will need to do this on a very large dataset. St_difference isn't particularly fast and I'm worried about running into a RAM issue. It would be great to create a lapply() so I can parallelize it. Unfortunately, my brain struggles with lapply() and I can't figure out how to enter 2 columns/lists in a lapply. Clearly this one doesn't work:

diff<- lapply(circ[,c("x","y")], function(x,y) {st_difference(x,y)})

because plot(diff[[1]]) plots the original circles. I've considered doing it in ArcGIS Pro because I thought 'Pairwise Erase' would work quicker, but that gives the same result as st_difference, with multiple holes in each overlapping big circle, and I don't know how to customize it.

How do I get my donuts?

1 Answer 1

0

How about applying a little filter based on pair to exclude the failed donuts and only retain the ones you want to keep?


diff <- st_difference(big.c, lil.c) 

diff
#> Simple feature collection with 4 features and 4 fields
#> Geometry type: POLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -7 ymin: -7 xmax: 7 ymax: 12
#> CRS:           +proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0
#>     pair size pair.1 size.1                       geometry
#> 1      a    x      a      y POLYGON ((6.990407 -0.36635...
#> 2      b    x      a      y POLYGON ((6.990407 4.633648...
#> 1.1    a    x      b      y POLYGON ((6.990407 -0.36635...
#> 2.1    b    x      b      y POLYGON ((6.990407 4.633648...

diff <- diff |> dplyr::filter(pair == pair.1)

diff
#> Simple feature collection with 2 features and 4 fields
#> Geometry type: POLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -7 ymin: -7 xmax: 7 ymax: 12
#> CRS:           +proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0
#>     pair size pair.1 size.1                       geometry
#> 1      a    x      a      y POLYGON ((6.990407 -0.36635...
#> 2.1    b    x      b      y POLYGON ((6.990407 4.633648...

st_geometry(diff[1, ]) |> plot(border = "blue")
st_geometry(diff[2, ]) |> plot(border = "red", add = TRUE)

Are you actually running into issues or are you just worried about? I'd at least try this on your actual dataset.

1
  • Somehow I didn't realize that st_difference was creating a different polygon for every instance of smaller circle within a bigger circle. Yes, filter works! However, I did still have the issues of allocating enough RAM for st_difference. I worked around it by cutting my data into chunks and feeding it through a for loop. It took some time to run but it got the job done. Thank you @falk_env
    – Wanner
    Aug 19, 2022 at 7:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.