rgisr-sfgeosphere

How to create a geometry between two sets of points


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...

enter image description here

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

enter image description here

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:

enter image description here

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')

enter image description here

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.


Solution

  • 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 
    

    enter image description here

    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:

    enter image description here

    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)