Triangulation wombling

using SpatialBoundaries
using CairoMakie
using StatsBase

Get some points at random

n = 500
x = 2rand(n)
y = rand(n)
z = [(x[i] <= 0.5) & (y[i] <= 0.5) ? rand() : rand() .+ 1.2 for i in eachindex(x)]

scatter(x, y; color = z)

Womble

W = wombling(x, y, z)

TriangulationWomble{Float64}([2.7082315707960385, 12.455620765882022, 2.5730019904789643, 19.837408492786295, 19.478818794497673, 4.427064814926273, 9.115803230152512, 8.52842158565146, 5.05086406606304, 6.5485964083552695  …  21.089326368631642, 10.979910308183655, 12.198642018740028, 6.709806998314583, 3.918757578019098, 12.234301250482654, 19.39516372427419, 35.41822743754425, 4.061021849090529, 4.490113331292793], [337.3078221743234, 203.0141758022063, 257.3845586857974, 186.76908386362453, 335.4164311034109, 98.30928520344224, 347.6681347642833, 329.5927545668783, 230.27260779134613, 174.63622013445294  …  18.172761784986903, 45.6912396037153, 179.37266789939818, 344.4977396000743, 306.4508558014026, 251.52889558478435, 312.62680574541037, 353.21038333337583, 339.28634006776804, 350.2976072454675], [1.6901794717106577, 0.5624508333819352, 1.461575950587589, 1.4988556698479751, 1.532419725888328, 0.4664144606656362, 0.5761386052582695, 1.427143535762774, 1.7305708539801306, 1.4962169649685417  …  1.8079755821585677, 1.8396784104888546, 1.1678191649055332, 1.0989786122755707, 0.5334298107864021, 0.5275329668228416, 1.7491938297102518, 1.7877147019049282, 0.45480741197916036, 0.3926344551778755], [0.3474454668274065, 0.7281847180188837, 0.06407919803115476, 0.8157915520833926, 0.6908704967094019, 0.6381066477121574, 0.5975700654252815, 0.2543872312500289, 0.2980793038389988, 0.5604015911415621  …  0.792920180125248, 0.7699276248786631, 0.14620212406369792, 0.15721800179030096, 0.7204773668562644, 0.7654268457496632, 0.9259364529950318, 0.965232257192644, 0.4560676593113577, 0.46373514643906927])

Get the rate of change

scatter(x, y; color = :lightgrey)
scatter!(W.x, W.y; color = log1p.(W.m))
current_figure()

Angle histogram

f = Figure()
ax = PolarAxis(f[1, 1])
h = fit(Histogram, deg2rad.(values(W.θ)); nbins = 100)
stairs!(ax, h.edges[1], h.weights[[axes(h.weights, 1)..., 1]]; color = :teal, linewidth = 2)
current_figure()

Show the rotation with a color

scatter(W.x, W.y; color = W.θ, colormap = :vikO, colorrange = (0, 360))