# Weighted Voronoi polygons in R

I want to build weighted Voronoi diagram in R, as in How to draw weighted voronoi polygons (thiessen) from points in QGIS? or Create weighted Thiessen polygons?, but I could not find any function for that. I have average coding skills for R, but not enough for this task. Is there any function available?

Eventually I managed to solve my problem, building from the excellent answer https://gis.stackexchange.com/a/17290/103010 from user "whuber" in the already cited Create weighted Thiessen polygons?.

This is meant to be a simple summary of the code, not of the science behind the weighted Voronoi polygons:

1. Calculate the standard Euclidean distance of every pixel in a raster to each point.
2. Apply a weighting function that change the distance value: in the example provided the value increases exponentially with the distance but decreases with higher size of the points (this works for my specific application). Bigger trees will have larger areas of influence but will not override completely smaller ones.
3. Identify which tree claims influence on a pixel: the one with the lower weighted distance.
4. With this info, build a raster, polygon or simply a data.table for your need
``````    # necessary packages
library(data.table)
library(raster)
library(rgeos)

# create base raster
rez= 1
r <- raster(xmn=0, xmx=40, ymn=0, ymx=40,
crs = '+proj=utm +zone=01 +datum=WGS84', resolution=rez)

# create example data
treez <- data.table(treeID=c("a","b","c","d"),
x=c(5,35,25,10),
y=c(35,35,15,10),
size=c(5,20,10,1))

# create an empty datatable with the length equal to the raster size
dt <- data.table(V1=seq(1,length.out = 40^2*(1/rez)^2))

# calculate distance grid for each tree, apply weigth and store the data
for (i in treez[,unique(treeID)]){
# 1) standard Euclidean distance
d <- distanceFromPoints(r, treez[treeID==i,.(x,y)])
# 2) apply your own custom weight
d <- ((d@data@values)^2)/treez[treeID==i,size]
# 3) store the results with the correct name
dt <- cbind(dt,d)
setnames(dt, "d", i)
}

# carry out comparison
# find the column with the lowest value (+1 to account to column V1 at the beginning)
dt[,minval:=which.min(.SD)+1,by=V1,.SDcols=2:length(names(dt))]
# find the name of that column, that is the id of the point
dt[,minlab:=names(dt)[minval]]

# create resulting raster
out <- raster(xmn=0, xmx=40, ymn=0, ymx=40,
crs = '+proj=utm +zone=01 +datum=WGS84', resolution=rez,
vals=dt[,as.factor(minlab)])

# Transform into polygons
outpoly <- rasterToPolygons(out, dissolve=T)
# Plot results with starting points
plot(outpoly);points(treez[,.(x,y)],col="red",pch=16, cex=sqrt(treez[,size]))
`````` • Could you post the plots and describe a little more how this works? – Simbamangu Sep 10 at 12:53