CODEWITHSUNDEEP
rgisr-sfgeosphere
Dean MacGregor
Dean MacGregor

Reputation: 18691

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.

Upvotes: 1

Views: 269

Answers (1)

Spacedman
Spacedman

Reputation: 94277

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)

Upvotes: 2

Related Questions