I have a simple GIS problem I was not able to find the solution on the internet.

My aim is quite simple: I have a polygons layer that contains several entities, and I would like to compute a new layer that would contain a line that separates the space surrounding each entity based on an equal distance from the closest entity.

To illustrate the question, I have drawn a schematic that sums up the problem.

enter image description here

Figure A represents a given shapefile that contains four polygons. Figure B includes the second layers (dashed blue line) that separate each entity based on equal distance from each other.

Ideally, I am looking for an answer on R or QGIS, as they are the two software I usually use.

  • 1
    Conceptually, convert the original polygons to points along the boundaries, form voronoi polygons from the points, dissolve the voronoi polygons by aggregating by the original polygon the point came from. Will require some data conversion, intersection, and st_voronoi work. Too late to work out the details but that's one procedure.
    – Spacedman
    Jan 19 '21 at 22:20

This function computes "polygon voronoi" areas - it divides the space up into areas that are nearest spatially dispersed polygons.

polyvor <- function(polys){

    pts = st_coordinates(polys)
    pts = pts[duplicated(pts[,"L3"]),]
    pts = st_as_sf(data.frame(pts), coords=1:2)[,"L3"]
    vpts = do.call(c, st_geometry(pts))
    pols = st_collection_extract(st_voronoi(vpts))

    ## which voronoi polygon is each pt in?
    pi = unlist(st_intersects(pts, pols))

    ## get this in the right order
    polsL = st_as_sf(data.frame(geometry=pols[pi,]))
    ## this is the polygon that each point came from
    polsL$polygon = pts$L3
    ## merge
    v = aggregate(polsL, by=list(polsL$polygon), FUN=mean)


Here's a test, we'll use four polygons from nc:

nc = st_read(system.file("shape/nc.shp", package="sf"))
test1 = nc[c(100,96, 86,88),]

test1_vor = polyvor(test1)

and to see what it looks like I'll plot the voronoi regions and add the original polygons:

plot(test1, add=TRUE, col="red")

enter image description here

Note this might fail if your polygons don't have very dense vertices along their boundaries - in this case you'd need to densify the polygons first (the opposite operation to "simplify"). I don't see an st_densify function in sf but I'm sure the capability is there somewhere. Maybe the solution is to turn the polygons into lines and use st_sample to generate lots of points along the borders...

Note this also fails if you have touching polygons because there are then coincident points in the voronoi calculation and bad things happen. For example with test2 = nc[c(100,98,88,82),], note the thick black boundaries show two adjacent polygons but the thin polygons of the output regions don't really know what they are doing inside there...

enter image description here

  • Dear Spacedman, thanks a lot for your help. Your code looks just what I'm looking for, and my layer contains only separate polygons so I should not be confronted to the issue you exposed in your second example. I'll keep you updated! Best regards, Guillaume
    – G. Papuga
    Jan 22 '21 at 8:14

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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