I have a set of points in 2D that I’m using to define Voronoi cells using NearestNeighbors.jl. I’d like to do some visualization of this. The two things I would ideally like to be able to do:
That should be fairly simple but I don’t think voronoidelaunay has recipes. Could you post an mwe, I guess it wouldn’t be hard to write a recipe (if noone has solved it before I’ll check back here on Tuesday).
But - basically your Vector shoult extract each cell as a Vector of vertices, use the @series macro to create a new series for each cell (with label = ""), and in that block return a vector of xs and a vector of ys, setting seriestype := :shape. You can then use the fill_z keyword to color the cells by a number or function.
The issue is, partly, that I start with the Voronoi cell centers. I don’t know how to turn them into the cell vertices/edges.
We have added that recently to GMT (the C lib), but it’s only available in its dev version (next version, GMT6.1.0, will come out in one or two weeks). With it, one can do
using GMT
G = nearneighbor("lixo.csv", region=(0,10,0,10), N="n", I=0.05);
imshow(G, fmt=:png)
where lixo.csv is a demo file with
2.53857 5.16657 0
2.48365 6.26811 1
8.68883 4.55983 2
4.11104 7.78623 3
1.79704 6.53027 4
7.17493 3.81713 5
3.41052 8.18161 6
8.35062 1.43348 7
8.17060 4.46765 8
5.27815 1.15172 9
and get
I tried extracting the Voronoi edges from VoronoiDelaunay. What. A. Rabbit. Hole. I might get back to you on that, but otherwise the “raster” approach shown by joa-quim is of course very easy. With Plots and NearestNeighbors it would look something like this:
using Plots, NearestNeighbors
# some function of the center point to color by - here the product of the coordinates
function col(i, j, kdtree, data)
id, dist = knn(kdtree, [i,j], 1)
nearest = data[:,id]
nearest[1]*nearest[2]
end
# build the neighbor tree
data = rand(2, 100)
kdtree = KDTree(data)
# color every point by its nearest center
heatmap([col(i, j, kdtree, data) for i in 0:0.001:1, j in 0:0.001:1], c = :viridis)
Voronoi tesselations, and their duals are almost always a huge hardship to contend with. I have yet to see a “clean” implementation, but, that being said, I would LOVE to have first class spatial partitioning and visualizations of such in Julia.
I’m personally a big fan of the line drawings, with the barycenters of the voronoi cells represented with dots.
The really tricky bit is doing this in 3-Space… 
I tried once… Maybe someone from #geometry wou;d be able too, I’m too mortal.
Line drawings are easy though. It’s the polygons that are hard to extract.
x, y = rand(100).+1, rand(100).+1 #yeah the package only accepts input data in the (1,2) range
tess = DelaunayTessellation()
push!(tess, Point2D.(x,y))
scatter(x,y, lims = (1,2), legend = false)
plot!(getplotxy(voronoiedges(tess))...)
This clumsy long code gives me most of the polygons - but the input order is lost, so it’s difficult to color the polygons by the center points. I think the raster approach is likely the best.
using VoronoiDelaunay, Plots
# create a polygon by putting all the edges (point pairs) in order (the order is seemingly random in the package)
function getpoly(eds)
a = geta.(eds)
b = getb.(eds)
ret = Point2D[a[1], b[1]]
ind = 1
for i in 1:length(a)
inds = setdiff(findall(x->x ≈ ret[end], a), [ind])
if !isempty(inds)
ind = inds[1]
push!(ret, b[ind])
else
inds = setdiff(findall(x->x ≈ ret[end], b), [ind])
ind = inds[1]
push!(ret, a[ind])
end
end
push!(ret, a[1])
ret
end
Base.isapprox(pt1::Point2D, pt2::Point2D) = pt1._x ≈ pt2._x && pt1._y ≈ pt2._y
# create polygons by filtering all of the edges for each point
function getpolys(tess::DelaunayTessellation2D{Point2D})
edges = collect(voronoiedges(tess))
ret = Vector{Point2D}[]
for pt in unique([getgenb.(edges); getgena.(edges)])
ed = filter(x-> getgena(x) ≈ pt || getgenb(x) ≈ pt, edges)
try
push!(ret, getpoly(ed))
catch
end
end
ret
end
@recipe f(pts::AbstractVector{<:Point2D}) = getx.(pts), gety.(pts)
@recipe function f(D::DelaunayTessellation2D{Point2D})
label := ""
for poly in getpolys(D)
@series begin
seriestype := :shape
poly
end
end
()
end
# some random data
x, y = rand(100).+1, rand(100).+1
tess = DelaunayTessellation()
push!(tess, Point2D.(x,y))
plot(tess, lims = (1,2), lc = :white)
@mkborregaard you can get the Voronoi polygons with the GMT’s triangulate program
D = triangulate("lixo.csv", region=(0,10,0,10), voronoi="n", network=true);
triangulate [WARNING]: -Qn is experimental and unstable.
plot(D, aspect="equal")
plot!("lixo.csv", marker=:circle, ms=0.2, fill=:black, fmt=:png, show=true)
can you get the barycenters of the cells with GMT?
Sweet! But are those polygons or outlines? Can you colour in the polygon area?
The GMT man page for the gmtspatial module says
``… The centroid is computed using the mean of the 3-D Cartesian vectors making up the polygon vertices, while the area is obtained via an equal-area projection. ```
Is this what you want?
julia> D = triangulate("lixo.csv", region=(0,10,0,10), voronoi="n", network=true);
triangulate [WARNING]: -Qn is experimental and unstable.
julia> Dc = gmtspatial(D, Q=true) # Compute centroid and area of each polygon
Array{GMT.GMTdataset,1} with 1 segments
First segment DATA
10×3 Array{Float64,2}:
3.01992 6.50059 2.43609
2.67102 4.2825 16.4437
5.8361 7.53071 15.0461
1.10002 6.92018 7.91042
2.13613 8.88082 7.76614
6.23657 4.08067 8.62185
3.65584 1.24974 15.8507
8.46981 1.81582 8.66603
7.80101 5.09087 5.20239
9.11705 6.79032 12.0565
julia> plot(D, aspect="equal")![GMTjl_tmp|504x500]
julia> plot!(Dc, marker=:circle, ms=0.2, fill=:black, fmt=:png, show=true)
Yes, they are polygons and yes they can be painted. I’ll prepare an example.
OK, painting the polygons uncovered an issue in the GMT.jl wrapper and a feature in the GMT program triangulate itself. So I have to resort to the hard core GMT syntax to show it. The issue is how to color the polygons with a colormap. Next version of GMT has a more handy way of assigning colors to polygons. But for now, this works and will be made nicer to read when I fix the wrapper issue.
D,G = triangulate("lixo.csv", region=(0,10,0,10), voronoi="n", grid=true, inc=0.05, network=true, edges=true);
triangulate [WARNING]: -Qn is experimental and unstable.
C = makecpt(cmap=:categorical, range=(0,10,1), cptname="categorical.cpt");
plot("-JX12c/0 -Baf -BWSen -Ccategorical.cpt -R-0.4/10.4/-0.4/10.4 -W1p+cf -P > voro_color.ps", D)
OK, for the record. The GMT.jl issue was just an harmless wrong warning and the GMT feature was actually my lack of attention reading the manual. The image above (and with the file posted at the beginning of this thread) can be obtained with
using GMT
D = triangulate("lixo.csv", region=(0,10,0,10), voronoi="n", network=true, xyz=true);
triangulate [WARNING]: -Qn is experimental and unstable.
C = makecpt(cmap=:categorical, range=(0,10,1));
plot(D, color=C, aspect=:equal, show=true, pen=(lw=1, csymbol=true), savefig="voronoi_colored.png")
yes that is what I wanted :). I’ll have to look at this package thank you!
This, and the one below, Voronoi Visualization - #14 by joa-quim, are exactly what I would like to be able to do.
It does look a lot easier with GMT in this case. The issue is not in Plots, but lies in the VoronoiDelaunay package being particularly optimised for edges but ignorant of polygons. So all my code does is that it tries to group and sort all the disparate unsorted edges to form valid polygons. Once you have them plotting is easy. I’ll be honest and say I find the VoronoiDelaunay package interface very little focused on user experience and very much on internal performance. I’m considering doing a wrapper package with the opposite trade-off…
If you are keen to use Plots, a userplot recipe extracting the centers from the tess object again could do the trick. Happy to give it a go, but let me know if you’ll just use the GMT instead, which seems to be simple and work.