CODEWITHSUNDEEP
rvoronoiareas
AlwaysInTheDark
AlwaysInTheDark

Reputation: 67

Overlapping circle areas on a Voronoi diagram

I've generated a Voronoi diagram from a series of points on an excel file. I would like to create a series of circles of the same radii each centred around a second data set of coordinates and determine the corresponding overlapping areas for each point.

I have used a function defined in https://gis.stackexchange.com/questions/190917/r-voronoi-tesselation-from-long-lat-data to create the Voronoi tessellation as a SpatialPolygonsDataFrame. However, I'm unsure how to generate the circles from the second data set. I used st_buffer at a defined radius, and tried to convert the Voronoi results using st_as_sf, but this didn't seem to work. When trying use st_intersection to determine the areas, I received the error st_crs(x) == st_crs(y), which seems to come from the Voronoi set not having defined coordinates.

I'm unsure which function is the best to use and whether sf has any additional tools that could be used to output the area that each segment of the Voronoi tessellation is overlapped by these radii.

I have found similar issues from https://gis.stackexchange.com/questions/229453/create-a-circle-of-defined-radius-around-a-point-and-then-find-the-overlapping-a?noredirect=1&lq=1 and https://gis.stackexchange.com/questions/140504/extracting-intersection-areas-in-r, but I cannot find a way to make these methods compatible with the Voronoi function.

A sample test data set that I've been using is included here:

Hubs <- cbind('X Coord' = c(52.37639999,52.36989975,52.86299896,52.01011658,51.67409897,50.84980011,51.88669968,52.1048943,52.0746994), 
          'Y Coord' = c(4.894589901,4.876679897,6.493730068,4.703330517,4.548630238,5.687580109,4.491449833,5.0528496,4.310130119),
          'Info 1'  = c(13,15,62,24,9,46,73,97,69))

References <- cbind('X Coord' = c(51.88529968,52.3360136,52.37440109,51.92269897,51.9192276,51.43019867,52.0780896,51.90299988,52.04100037,50.98810196),
                'Y Coord' = c(4.493549824,4.8750552,4.895800114,4.47382021,4.481627464,5.451980114,5.1470653,4.458670139,4.318089962,5.78059721))

Upvotes: 0

Views: 610

Answers (2)

Robert Hijmans
Robert Hijmans

Reputation: 47401

Example data:

library(dismo)
p <- cbind(x=c(17, 42, 85, 70, 19, 53, 26, 84, 84, 46, 48, 85, 4, 95, 48, 54, 66, 74, 50, 48), 
      y=c(28, 73, 38, 56, 43, 29, 63, 22, 46, 45, 7, 60, 46, 34, 14, 51, 70, 31, 39, 26))

p1 <- p[11:20,]
p2 <- p[1:10,]

Create voronoi tesselation and circles

vr <- voronoi(p1)
names(vr) <- "vID"
cr <- circles(p2, 5, lonlat=FALSE, dissolve=FALSE)
cr <- SpatialPolygonsDataFrame(polygons(cr), data.frame(cID=1:10))

Intersect and compute area

x <- intersect(vr, cr)
crs(x) <- "+proj=utm +zone=1 +datum=WGS84" 
x$area <- area(x)
data.frame(x)

That was using Spatial objects. Now with sf following AlwaysInTheDark, but without the loop

library(sf)
sf2 <- st_as_sf(data.frame(p2), coords = c("x", "y"), crs = 3035)
circles <- st_buffer(sf2, dist = 5)
circles$cID <- 1:nrow(circles)

sf1 <- st_as_sf(data.frame(p1), coords = c("x", "y"), crs = 3035)
voroi <- st_collection_extract(st_voronoi(st_combine(sf1)))

d <- st_intersection(circles, voroi)
d$area <- st_area(d)

d

Upvotes: 1

AlwaysInTheDark
AlwaysInTheDark

Reputation: 67

Whilst I was able to come up with this solution, this may not be the most efficient solution. 

library(raster)
library(tidyverse)
library(sf)
library(sp)
library(rgdal)
library(rgeos)
library(maptools)
#> Linking to GEOS 3.5.0, GDAL 2.1.1, proj.4 4.9.3

voronoipolygons = function(layer) {
require(deldir)
crds = layer@coords
z = deldir(crds[,1], crds[,2])
w = tile.list(z)
polys = vector(mode='list', length=length(w))
require(sp)
for (i in seq(along=polys)) {
    pcrds = cbind(w[[i]]$x, w[[i]]$y)
    pcrds = rbind(pcrds, pcrds[1,])
    polys[[i]] = Polygons(list(Polygon(pcrds)), ID=as.character(i))
}
SP = SpatialPolygons(polys)
voronoi = SpatialPolygonsDataFrame(SP, data=data.frame(dummy = 
seq(length(SP)), row.names=sapply(slot(SP, 'polygons'), 

function(x) slot(x, 'ID'))))
}

Hubs = read.csv("Test_Hub_Data.csv", header=TRUE)
coordinates(Hubs) <- c("X.Coord", "Y.Coord")

vp <- voronoipolygons(Hubs)

#Voronoi Setup
crs(vp) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
vp_test <- sf::st_as_sf(vp) %>% st_transform(3035)

References = read.csv("Test_Reference_Points.csv", header = TRUE)

# Convert to sf, set the crs to EPSG:4326 (lat/long), 
# and transform to EPSG:3035
Ref_sf <- st_as_sf(References, coords = c("X.Coord", "Y.Coord"), crs = 4326) %>% st_transform(3035)

#Initialise matrix of areas x: Voronoi Segments, y: Nodes
areas <- matrix(0,length(References[,1]),length(Hubs[,1]))

for (i in 1:length(References[,1])){
   #New circles need to be generated for each point
   # Buffer circles by 100m
   Ref_circle <- st_buffer(Ref_sf[i,1], dist = 100)

   for(j in 1:length(Hubs[,1])){
   # Intersect the circles with the polygons
        vp_int_circle <- st_intersection(vp_test[j,1], Ref_circle)

        #Append array
        areas[i,j] = sum(st_area(vp_int_circle))
   }
}

Upvotes: 0

Related Questions