0

I compute the Voronoi diagram using certain graphical points with the deldir package and then apply a geographical border to the result much like this question . However, the intersection with the borders I want to apply cut the polygons in several parts and I would like not to have this behavior. The problem is best seen executing the next example.

library(sp)
library(rgeos)
library(deldir)

points <- structure(c(-2.916667, -2.85, -2.8, -2.966667, -2.9997, -2.766667, 
                      -2.75, 47.616667, 47.766667, 47.7, 47.683333, 47.5856, 47.533333, 
                      47.616667), .Dim = c(7L, 2L),
                    .Dimnames = list(c("2", "38", "103", "42", "167", "173", "175"), c("x", "y")))
bound <- structure(c(-2.702027, -2.67545683, -2.74795532, -2.81112671, 
                     -2.83035278, -2.85850525, -2.87222692, -2.89558411, -2.86605835, 
                     -2.85598337, -2.78846741, -2.74177551, -2.73610594, -2.73078918, 
                     -2.81730652, -2.83333365, -2.83333365, -2.88871765, -2.91618347, 
                     -2.95394897, -2.93746948, -2.94190707, -2.98347473, -3.03909302, 
                     -3.06130189, -3.050263, -3.063889, -2.933333, -2.933333, -2.802012, 
                     -2.777703, -2.774573, -2.725, -2.725, -2.719737, -2.69853, -2.702027, 
                     47.563739, 47.49288545, 47.49122489, 47.49122489, 47.49029694, 
                     47.51209932, 47.51999275, 47.53342905, 47.54965239, 47.54923113, 
                     47.54640812, 47.54640812, 47.57055474, 47.59319878, 47.616347, 
                     47.61156744, 47.61156745, 47.59505101, 47.58810479, 47.60014432, 
                     47.55660371, 47.55763166, 47.56726061, 47.56679731, 47.56747229, 
                     47.657139, 47.681667, 47.76, 47.794445, 47.838218, 47.799324, 
                     47.771154, 47.711666, 47.688334, 47.682018, 47.575979, 47.563739
                     ), .Dim = c(37L, 2L), .Dimnames = list(NULL, c("x", "y")))
boundPol <- SpatialPolygons(list(Polygons(list(Polygon(bound)), ID = 'boundId')), proj4string = sp::CRS('+proj=merc'))
rw <- as.numeric(t(bbox(boundPol)))
tl <- deldir::tile.list(deldir::deldir(points[, 1], points[, 2], rw = rw))
ids <- rownames(points)
polys <- vector(mode = 'list', length = length(tl))
for (nbr in seq_along(tl)){
    xy <- tl[[nbr]]
    p <- matrix(c(xy$x, xy$x[1], xy$y, xy$y[1]), length(xy$x) + 1, 2)
    polys[[nbr]] <- sp::Polygons(list(sp::Polygon(p)), ID = as.character(ids[nbr]))
}
spPolys <- sp::SpatialPolygons(polys, proj4string = sp::CRS('+proj=merc'))
v <- rgeos::gIntersection(spPolys, boundPol, byid = TRUE, id = ids)

npolys <- sapply(v@polygons, function(x) length(x@Polygons))
plot(v)
plot(v[1, ], col = 'red', add = TRUE)

Does anybody has the solution so the split polygons can be re factored into their neighbours?

enter image description here

Community
  • 1
  • 1
Usobi
  • 1,816
  • 4
  • 18
  • 25
  • What is your expected output? – MichaelChirico Nov 22 '15 at 18:33
  • I would like the small part of the red polygon to be integrated into its neighbours (single neighbour in this case). – Usobi Nov 22 '15 at 18:48
  • 1
    Based on what criteria, exactly, in general? Just that it's isolated? – MichaelChirico Nov 22 '15 at 19:16
  • Currently I take all the polygon parts in which the original point (before the triangulation) did *not* lie and integrate the 'leftover' polygon in the neighboring polygons by applying Delaunay triangulation to them... But I was wondering if there was a better way. – Usobi Nov 23 '15 at 17:00

0 Answers0