How to create a geometry between two sets of points

Question

I have two sets of points (actually I have many hundreds of pairs of two sets but for example, here's one) that look like this...

What I want to do is create a linestring or multilinestring that, when drawn, looks loosely like this...

The data looks like

onecons<- st_read("https://raw.githubusercontent.com/deanm0000/SOexamples/e5591130b1d0e9c5c5fa8b6b5fc05b3958d13099/example.geojson")
mapdeck(style=mapdeck_style("light")) %>% 
  add_scatterplot(onecons, radius_min_pixels=5,
                   tooltip='nodename', fill_colour = "pos", layer_id = 'nodes')

I think I'm close.

I did

polys<-onecons %>% group_by(pos) %>% summarise() %>% st_convex_hull()
mapdeck(style=mapdeck_style("light")) %>% 
    add_polygon(st_as_sf(st_difference(st_geometry(polys[1,]), st_geometry(polys[2,]))), fill_opacity=.5) %>%
      add_scatterplot(onecons, radius_min_pixels=5,
                       tooltip='nodeest', fill_colour = "pos", layer_id = 'nodes')

which gives me:

A couple issues, if I flip the row reference in st_difference then what I get doesn't make sense ...

mapdeck(style=mapdeck_style("light")) %>% 
  add_polygon(st_as_sf(st_difference(st_geometry(polys[2,]), st_geometry(polys[1,]))), fill_opacity=.5) %>%
  add_scatterplot(onecons, radius_min_pixels=5,
                  tooltip='nodeest', fill_colour = "pos", layer_id = 'nodes')

Aside from that, I want the line I generate to be in between the points, not on them. I can live with it if it does land on the points but I'd rather it be in between.

Lastly, I don't know how to get the relevant border of the polygon to represent my new line. I can use st_boundary to get the full boundary but I don't know how to just get the part that is in between.

Construct voronoi polygons based on the points, dissolve based on the point colour to get two polygons. For two sets of points nicely spaced as above this should get a boundary that bisects the sets. If there's a yellow dot in Houston though, it'll be a separate polygon for yellow in the middle of the blue polygon. — Spacedman, Sep 10 '22 at 08:39

score 2 · Accepted Answer · answered Sep 12 '22 at 08:57

First, define three useful functions. One computes voronoi polygons for a set of points and re-matches them to the points because the order of returned voronoi polygons doesn't match the point order you started with - this mostly taken from ?st_voronoi:

st_pvoronoi <- function(pts){
    ## compute voronoi polys and match back to points
    pols = st_collection_extract(st_voronoi(do.call(c, st_geometry(pts))))
    st_crs(pols) = st_crs(pts)
    st_geometry(pts) = pols[unlist(st_intersects(pts, pols))]
    pts
}

Then a function to dissolve polygons grouped by some variable. This is an alternative to using group_by without the unnecessary complication of having to handle non-standard evaluation of passing a value for splitting:

st_dissolve <- function(polys, v){
    ## dissolve polygons with common value `v`
    splits = split(polys, v)
    dissol = lapply(splits, st_union)
    st_as_sf(
        data.frame(
            ID = names(dissol),
            geometry = do.call(c, dissol)
        )
    )
}

And thirdly a function that returns only lines that exist more than once in a geometry set:

st_common  <- function(x){
    ## convert polys to boundary lines, intersect
    do.call(st_intersection, st_boundary(x)$geometry)
}

Now how to use these. I've roughly digitised your blue and yellow example points:

library(sf)
pts = st_read("pts.gpkg")
summary(pts$type)
plot(st_geometry(pts), col=pts$type+1, pch=19)

giving:

## pts = st_read("pts.gpkg")
Reading layer `pts' from data source 
  `/nobackup/rowlings/Downloads/SO/bord/pts.gpkg' using driver `GPKG'
Simple feature collection with 120 features and 1 field
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -11576280 ymin: 2996930 xmax: -10502040 ymax: 4110600
Projected CRS: WGS 84 / Pseudo-Mercator
## summary(pts$type)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  0.000   0.000   0.000   0.475   1.000   1.000

Then I construct the voronoi polygons, dissolve them on the type, then find the common boundary:

pols = st_pvoronoi(pts)
pp = st_dissolve(pols, pols$type)
pps = st_common(pp)
plot(pps, add=TRUE)

You can see what each of those stages do by inspecting the object at each stage, or plotting. The final plot looks like this:

Note that if you have isolated points within a region of different type then those points will get "pinched off" and a ring line will be returned.

You should put this into a function:

st_borderline <- function(pts, v){
    pols = st_pvoronoi(pts)
    pp = st_dissolve(pols, v)
    pps = st_common(pp)
    return(pps)
}

Then all you need do is:

pts = st_read("pts.gpkg")
my_borderline = st_borderline(pts, pts$type)

How to create a geometry between two sets of points

1 Answers1