Triangulation wombling

using SpatialBoundaries
using StatsPlots

Plot defaults

default(; dpi=500, size=(600, 600), aspectratio=1, c=:davos, frame=:box)

Get some points at random

n = 500
x = rand(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, marker_z = z, lab="")

Womble

W = wombling(x, y, z)

TriangulationWomble{Float64}([0.17692234530713924, 0.8226392375685523, 32.10484408405045, 12.353302414872358, 4.819839047252428, 1.395122847533027, 3.219720019240268, 3.2095180719073295, 11.628943267416647, 9.3439218954165  …  9.651030914849331, 5.350405438116565, 11.779126865173586, 3.104087955291764, 40.984122015465545, 2.8746868739286784, 13.29066106511594, 7.299196253355385, 11.928303657358905, 10.662360749837756], [276.3930642199775, 259.76247745065996, 58.483248577567196, 118.63466978907199, 122.31873642258111, 200.52360007873875, 116.04191129822642, 305.6711539903525, 95.99830991581321, 269.67894768589946  …  197.29761307494897, 347.69332564956323, 198.8030675612142, 258.0554717975009, 265.84403485384723, 129.68039212427647, 85.58846575739412, 179.4196539990786, 345.79470187415166, 336.7579858873496], [0.26424441449185493, 0.132635027845171, 0.46588239398815734, 0.4763763239420615, 0.0799401223391432, 0.08483038695947469, 0.04656837697189679, 0.18154239116004964, 0.07847193748926284, 0.09674968967515  …  0.4821560814971255, 0.45724902215886704, 0.26833594141149136, 0.2146642936003904, 0.07879912404535827, 0.04628887990714733, 0.4845917756014604, 0.4405352826651822, 0.7007360908403121, 0.71378036845711], [0.34257163686644176, 0.14634057917627544, 0.7781799992959533, 0.9550101255378575, 0.699051392106659, 0.29304134820272937, 0.265148593638694, 0.389327934393626, 0.4902284174356338, 0.3728662495671034  …  0.34659527559009295, 0.33087080662758434, 0.17003818970577866, 0.21220292619560396, 0.9353817200452452, 0.883537245776801, 0.6592285799904075, 0.6426853953493149, 0.2648141364435504, 0.287989282293272])

Get the rate of change

scatter(x, y, c=:lightgrey, msw=0.0, lab="", m=:square, ms=3)
scatter!(W.x, W.y, marker_z = log1p.(W.m), lab="")

Angle histogram

stephist(
    deg2rad.(sort(vec(W.θ)));
    proj=:polar,
    lab="",
    c=:teal,
    fill=(0, 0.2, :teal),
    nbins=100,
)

Show the rotation with a color

scatter(W.x, W.y, marker_z = W.θ, c=:vik, clim=(0, 360))