diff --git a/docs/Project.toml b/docs/Project.toml
index fccffcc39..1cc496415 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -37,6 +37,7 @@ MultiFloats = "bdf0d083-296b-4888-a5b6-7498122e68a5"
 NaturalEarth = "436b0209-26ab-4e65-94a9-6526d86fea76"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 Proj = "c94c279d-25a6-4763-9509-64d165bea63e"
+Quickhull = "baca2653-9c88-49d2-a488-12dbf25fc4fb"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Shapefile = "8e980c4a-a4fe-5da2-b3a7-4b4b0353a2f4"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
diff --git a/docs/src/experiments/spherical_delaunay.jl b/docs/src/experiments/spherical_delaunay.jl
new file mode 100644
index 000000000..d8251c4c3
--- /dev/null
+++ b/docs/src/experiments/spherical_delaunay.jl
@@ -0,0 +1,217 @@
+#=
+# Spherical Delaunay triangulation using stereographic projections
+
+This file encodes two approaches to spherical Delaunay triangulation.
+
+1. `StereographicDelaunayTriangulation` projects the coordinates to a rotated stereographic projection, then computes the Delaunay triangulation on the plane.
+    This approach was taken from d3-geo-voronoi.
+    ```
+    Copyright 2018-2021 Philippe Rivière
+
+    Permission to use, copy, modify, and/or distribute this software for any purpose
+    with or without fee is hereby granted, provided that the above copyright notice
+    and this permission notice appear in all copies.
+
+    THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES WITH
+    REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND
+    FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY SPECIAL, DIRECT,
+    INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS
+    OF USE, DATA OR PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER
+    TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
+    THIS SOFTWARE.
+    ```
+2. `SphericalConvexHull` computes the convex hull of the points in 3D Cartesian space, which is by definition the Delaunay triangulation.  This triangulation is currently done using the unreleased Quickhull.jl.
+=#
+
+import GeometryOps as GO, GeoInterface as GI, GeoFormatTypes as GFT
+import Proj # for easy stereographic projection - TODO implement in Julia
+import DelaunayTriangulation as DelTri # Delaunay triangulation on the 2d plane
+import Quickhull # convex hulls in d+1 are Delaunay triangulations in dimension d
+import CoordinateTransformations, Rotations
+
+using Downloads # does what it says on the tin
+using JSON3 # to load data
+using CairoMakie, GeoMakie # for plotting
+import Makie: Point3d
+
+abstract type SphericalTriangulationAlgorithm end
+struct StereographicDelaunayTriangulation <: SphericalTriangulationAlgorithm end
+struct SphericalConvexHull <: SphericalTriangulationAlgorithm end 
+
+spherical_triangulation(input_points; kwargs...) = spherical_triangulation(StereographicDelaunayTriangulation(), input_points; kwargs...)
+
+function spherical_triangulation(::StereographicDelaunayTriangulation, input_points; facetype = CairoMakie.GeometryBasics.TriangleFace)
+    # @assert GI.crstrait(first(input_points)) isa GI.AbstractGeographicCRSTrait
+    points = GO.tuples(input_points)
+    # @assert points isa Vector{GI.Point}
+
+    # In 
+    pivot_ind = findfirst(x -> all(isfinite, x), points)
+    
+    pivot_point = points[pivot_ind]
+    necessary_rotation = #=Rotations.RotY(-π) *=# Rotations.RotY(-deg2rad(90-pivot_point[2])) * Rotations.RotZ(-deg2rad(pivot_point[1]))
+
+    net_transformation_to_corrected_cartesian = CoordinateTransformations.LinearMap(necessary_rotation) ∘ UnitCartesianFromGeographic()
+    stereographic_points = (StereographicFromCartesian() ∘ net_transformation_to_corrected_cartesian).(points)
+    triangulation_points = copy(stereographic_points)
+    
+    # Iterate through the points and find infinite/invalid points
+    max_radius = 1
+    really_far_idxs = Int[]
+    for (i, point) in enumerate(stereographic_points)
+        m = hypot(point[1], point[2])
+        if !isfinite(m) || m > 1e32
+            push!(really_far_idxs, i)
+        elseif m > max_radius
+            max_radius = m
+        end
+    end
+    @debug max_radius length(really_far_idxs)
+    far_value = 1e6 * sqrt(max_radius)
+
+    if !isempty(really_far_idxs)
+        triangulation_points[really_far_idxs] .= (Point2(far_value, 0.0),)
+    end
+
+    boundary_points = reverse([
+        (-far_value, -far_value),
+        (-far_value, far_value),
+        (far_value, far_value),
+        (far_value, -far_value),
+        (-far_value, -far_value),
+    ])
+
+    # boundary_nodes, pts = DelTri.convert_boundary_points_to_indices(boundary_points; existing_points=triangulation_points)
+
+    # triangulation = DelTri.triangulate(triangulation_points; boundary_nodes)
+    triangulation = DelTri.triangulate(triangulation_points)
+    #  for diagnostics, try `fig, ax, sc = triplot(triangulation, show_constrained_edges=true, show_convex_hull=true)`
+
+    # First, get all the "solid" faces, ie, faces not attached to boundary nodes
+    original_triangles = collect(DelTri.each_solid_triangle(triangulation))
+    boundary_face_inds = findall(Base.Fix1(DelTri.is_boundary_triangle, triangulation), original_triangles)
+    # Now we construct the faces.  However, we need to reverse the order of the points 
+    # because the Delaunay triangulation is constructed in the stereographic plane,
+    # and orientation there is the opposite of orientation on the sphere.
+    faces = map(facetype ∘ reverse, view(original_triangles, setdiff(axes(original_triangles, 1), boundary_face_inds)))
+    for boundary_face in view(original_triangles, boundary_face_inds)
+        push!(faces, reverse(facetype(map(boundary_face) do i; first(DelTri.is_boundary_node(triangulation, i)) ? pivot_ind : i end)))
+    end
+    
+    for ghost_face in DelTri.each_ghost_triangle(triangulation)
+        push!(faces, reverse(facetype(map(ghost_face) do i; DelTri.is_ghost_vertex(i) ? pivot_ind : i end)))
+    end
+    # Remove degenerate triangles
+    filter!(faces) do face
+        !(DelTri.geti(face) == DelTri.getj(face) || DelTri.getj(face) == DelTri.getk(face) || DelTri.geti(face) == DelTri.getk(face))
+    end
+
+    return faces
+end
+
+function spherical_triangulation(::SphericalConvexHull, input_points; facetype = CairoMakie.GeometryBasics.TriangleFace{Int32})
+    points = GO.tuples(input_points) # we have to decompose the points into tuples, so they work with Quickhull.jl
+    # @assert points isa Vector{GI.Point}
+    cartesian_points = map(UnitCartesianFromGeographic(), points)
+    # The Delaunay triangulation of points on a sphere is simply the convex hull of those points in 3D Cartesian space.
+    # We can use e.g Quickhull.jl to get us such a convex hull.
+    hull = Quickhull.quickhull(cartesian_points)
+    # We return only the faces from these triangulation methods, so we simply map
+    # the facetype to the returned values from `Quickhull.facets`.
+    return map(facetype, Quickhull.facets(hull))
+end
+
+
+# necessary coordinate transformations
+
+struct StereographicFromCartesian <: CoordinateTransformations.Transformation
+end
+
+function (::StereographicFromCartesian)(xyz::AbstractVector)
+    @assert length(xyz) == 3 "StereographicFromCartesian expects a 3D Cartesian vector"
+    x, y, z = xyz
+    # The Wikipedia definition has the north pole at infinity,
+    # this implementation has the south pole at infinity.
+    return Point2(x/(1-z), y/(1-z))
+end
+
+struct CartesianFromStereographic <: CoordinateTransformations.Transformation
+end
+
+function (::CartesianFromStereographic)(stereographic_point)
+    X, Y = stereographic_point
+    x2y2_1 = X^2 + Y^2 + 1
+    return Point3(2X/x2y2_1, 2Y/x2y2_1, (x2y2_1 - 2)/x2y2_1)
+end
+
+struct UnitCartesianFromGeographic <: CoordinateTransformations.Transformation 
+end
+
+function (::UnitCartesianFromGeographic)(geographic_point)
+    # Longitude is directly translatable to a spherical coordinate
+    # θ (azimuth)
+    θ = deg2rad(GI.x(geographic_point))
+    # The polar angle is 90 degrees minus the latitude
+    # ϕ (polar angle)
+    ϕ = deg2rad(90 - GI.y(geographic_point))
+    # Since this is the unit sphere, the radius is assumed to be 1,
+    # and we don't need to multiply by it.
+    return Point3(
+        sin(ϕ) * cos(θ),
+        sin(ϕ) * sin(θ),
+        cos(ϕ)
+    )
+end
+
+struct GeographicFromUnitCartesian <: CoordinateTransformations.Transformation 
+end
+
+function (::GeographicFromUnitCartesian)(xyz::AbstractVector)
+    @assert length(xyz) == 3 "GeographicFromUnitCartesian expects a 3D Cartesian vector"
+    x, y, z = xyz
+    return Point2(
+        atan(y, x),
+        atan(hypot(x, y), z),
+    )
+end
+
+function randsphericalangles(n)
+    θ = 2π .* rand(n)
+    ϕ = acos.(2 .* rand(n) .- 1)
+    return Point2.(θ, ϕ)
+end
+
+function randsphere(n)
+    θϕs = randsphericalangles(n)
+    return Point3.(
+        sin.(last.(θϕs)) .* cos.(first.(θϕs)),
+        sin.(last.(θϕs)) .* sin.(first.(θϕs)),
+        cos.(last.(θϕs))
+    )
+end
+
+# θϕs = randsphericalangles(50)
+# θs, ϕs = first.(θϕs), last.(θϕs)
+# pts = Point3.(
+#     sin.(ϕs) .* cos.(θs),
+#     sin.(ϕs) .* sin.(θs),
+#     cos.(ϕs)
+# )
+
+# f, a, p = scatter(pts; color = [fill(:blue, 49); :red])
+
+# function Makie.rotate!(t::Makie.Transformable, rotation::Rotations.Rotation)
+#     quat = Rotations.QuatRotation(rotation)
+#     Makie.rotate!(Makie.Absolute, t, Makie.Quaternion(quat.q.v1, quat.q.v2, quat.q.v3, quat.q.s))
+# end
+
+
+# rotate!(p, Rotations.RotX(π/2))
+# rotate!(p, Rotations.RotX(0))
+
+# pivot_point = θϕs[end]
+# necessary_rotation = Rotations.RotY(pi) * Rotations.RotY(-pivot_point[2]) * Rotations.RotZ(-pivot_point[1])
+
+# rotate!(p, necessary_rotation)
+
+# f
\ No newline at end of file
diff --git a/docs/src/experiments/spherical_natural_neighbour_interpolation.jl b/docs/src/experiments/spherical_natural_neighbour_interpolation.jl
new file mode 100644
index 000000000..1a7e340a9
--- /dev/null
+++ b/docs/src/experiments/spherical_natural_neighbour_interpolation.jl
@@ -0,0 +1,307 @@
+#=
+# Spherical natural neighbour interpolatoin
+=#
+
+import GeometryOps as GO, GeoInterface as GI, GeoFormatTypes as GFT
+import Proj # for easy stereographic projection - TODO implement in Julia
+import DelaunayTriangulation as DelTri # Delaunay triangulation on the 2d plane
+import CoordinateTransformations, Rotations
+
+using Downloads # does what it says on the tin
+using JSON3 # to load data
+using CairoMakie, GeoMakie # for plotting
+import Makie: Point3d
+
+include(joinpath(@__DIR__, "spherical_delaunay.jl"))
+
+using LinearAlgebra
+using GeometryBasics
+struct SphericalCap{T}
+    point::Point3{T}
+    radius::T
+end
+
+function SphericalCap(point::Point3{T1}, radius::T2) where {T1, T2}
+    return SphericalCap{promote_type(T1, T2)}(point, radius)
+end
+
+
+function circumcenter_on_unit_sphere(a, b, c)
+    LinearAlgebra.normalize(a × b + b × c + c × a)
+end
+
+spherical_distance(x::Point3, y::Point3) = acos(clamp(x ⋅ y, -1.0, 1.0))
+
+"Get the circumcenter of the triangle (a, b, c) on the unit sphere.  Returns a normalized 3-vector."
+function SphericalCap(a::Point3, b::Point3, c::Point3)
+    circumcenter = circumcenter_on_unit_sphere(a, b, c)
+    circumradius = spherical_distance(a, circumcenter)
+    return SphericalCap(circumcenter, circumradius)
+end
+
+function _is_ccw_unit_sphere(v_0,v_c,v_i)
+    # checks if the smaller interior angle for the great circles connecting u-v and v-w is CCW
+    return(LinearAlgebra.dot(LinearAlgebra.cross(v_c - v_0,v_i - v_c), v_i) < 0)
+end
+
+function angle_between(a, b, c)
+    ab = b - a
+    bc = c - b
+    norm_dot = (ab ⋅ bc) / (LinearAlgebra.norm(ab) * LinearAlgebra.norm(bc))
+    angle =  acos(clamp(norm_dot, -1.0, 1.0))
+    if _is_ccw_unit_sphere(a, b, c)
+        return angle
+    else
+        return 2π - angle
+    end
+end
+
+function bowyer_watson_envelope!(applicable_points, query_point, points, faces, caps = map(splat(SphericalCap), (view(points, face) for face in faces)); applicable_cap_indices = Int64[])
+    # brute force for now, but try the jump and search algorithm later
+    # can use e.g GeometryBasics.decompose(Point3{Float64}, GeometryBasics.Sphere(Point3(0.0), 1.0), 5) 
+    # to get starting points, or similar
+    empty!(applicable_cap_indices)
+    for (i, cap) in enumerate(caps)
+        if cap.radius ≥ spherical_distance(query_point, cap.point)
+            push!(applicable_cap_indices, i)
+        end
+    end
+    # Now that we have the face indices, we need to get the applicable points
+    empty!(applicable_points)
+    # for face_idx in applicable_cap_indices
+    #     append!(applicable_points, faces[face_idx])
+    # end
+    for i in applicable_cap_indices
+        current_face = faces[i]
+        edge_reoccurs = false
+        for current_edge in ((current_face[1], current_face[2]), (current_face[2], current_face[3]), (current_face[3], current_face[1]))
+            for j in applicable_cap_indices
+                if j == i
+                    continue # can't compare a triangle to itself
+                end
+                face_to_compare = faces[j]
+                for edge_to_compare in ((face_to_compare[1], face_to_compare[2]), (face_to_compare[2], face_to_compare[3]), (face_to_compare[3], face_to_compare[1]))
+                    if edge_to_compare == current_edge || reverse(edge_to_compare) == current_edge
+                        edge_reoccurs = true
+                        break
+                    end
+                end
+            end
+            if !edge_reoccurs # edge is unique
+                push!(applicable_points, current_edge[1])
+                push!(applicable_points, current_edge[2])
+            end
+            edge_reoccurs = false
+        end
+    end
+    unique!(applicable_points)
+    # Pick a random point (in this case, the first point) and 
+    # compute the angle subtended by the first point, the query
+    # point, and each individual point.  Sorting these angles
+    # allows us to "wind" the polygon in a clockwise way.
+    three_d_points = view(points, applicable_points)
+    angles = [angle_between(three_d_points[1], query_point, point) for point in three_d_points]
+    pt_inds = sortperm(angles)
+    permute!(applicable_points, pt_inds)
+    # push!(applicable_points, applicable_points[begin])
+    return applicable_points
+end
+
+import NaturalNeighbours: previndex_circular, nextindex_circular
+function laplace_ratio(points, envelope, i #= current vertex index =#, interpolation_point)
+    u = envelope[i]
+    prev_u = envelope[previndex_circular(envelope, i)]
+    next_u = envelope[nextindex_circular(envelope, i)]
+    p, q, r = points[u], points[prev_u], points[next_u]
+    g1 = circumcenter_on_unit_sphere(q, p, interpolation_point)
+    g2 = circumcenter_on_unit_sphere(p, r, interpolation_point)
+    ℓ = spherical_distance(g1, g2)
+    d = spherical_distance(p, interpolation_point)
+    w = ℓ / d
+    return w, u, prev_u, next_u
+end
+
+struct NaturalCoordinates{F, I}
+    coordinates::Vector{F}
+    indices::Vector{I}
+    interpolation_point::Point3{Float64}
+end
+
+function laplace_nearest_neighbour_coords(points, faces, interpolation_point, caps = map(splat(SphericalCap), (view(points, face) for face in faces)); envelope = Int64[], cap_idxs = Int64[])
+    envelope = bowyer_watson_envelope!(envelope, interpolation_point, points, faces, caps; applicable_cap_indices = cap_idxs)
+    coords = NaturalCoordinates(Float64[], Int64[], interpolation_point)
+    for i in eachindex(envelope)
+        w, u, prev_u, next_u = laplace_ratio(points, envelope, i, interpolation_point)
+        push!(coords.coordinates, w)
+        push!(coords.indices, u)
+    end
+    coords.coordinates ./= sum(coords.coordinates)
+    return coords
+end
+
+function eval_laplace_coordinates(points, faces, zs, interpolation_point, caps = map(splat(SphericalCap), (view(points, face) for face in faces)); envelope = Int64[], cap_idxs = Int64[])
+    coords = laplace_nearest_neighbour_coords(points, faces, interpolation_point, caps; envelope, cap_idxs)
+    if isempty(coords.coordinates)
+        return Inf
+    end
+    return sum((coord * zs[idx] for (coord, idx) in zip(coords.coordinates, coords.indices)))
+end
+
+function get_voronoi_polygon(point_index, points, faces, caps)
+    # For a Voronoi polygon, we can simply retrieve all triangles that have
+    # point_index as a vertex, and then filter them together.  
+end
+
+
+
+# These points are known to be good points, i.e., lon, lat, alt
+geographic_points = Point3{Float64}.(JSON3.read(read(Downloads.download("https://gist.githubusercontent.com/Fil/6bc12c535edc3602813a6ef2d1c73891/raw/3ae88bf307e740ddc020303ea95d7d2ecdec0d19/points.json"), String)))
+cartesian_points = UnitCartesianFromGeographic().(geographic_points)
+z_values = last.(geographic_points)
+
+faces = spherical_triangulation(geographic_points)
+
+faces2 = spherical_triangulation(SphericalConvexHull(), geographic_points)
+
+unique!(sort!(reduce(vcat, faces))) # so how am I getting this index?
+
+
+caps = map(splat(SphericalCap), (view(cartesian_points, face) for face in faces))
+
+lons = -180.0:1.0:180.0
+lats = -90.0:1.0:90.0
+
+eval_laplace_coordinates(cartesian_points, faces, z_values, LinearAlgebra.normalize(Point3(1.0, 1.0, 0.0)))
+eval_laplace_coordinates(cartesian_points, faces2, z_values, LinearAlgebra.normalize(Point3(1.0, 1.0, 0.0)))
+
+values = map(UnitCartesianFromGeographic().(Point2.(lons, lats'))) do point
+    eval_laplace_coordinates(cartesian_points, faces, z_values, point)
+end
+
+heatmap(lons, lats, values; axis = (; aspect = DataAspect()))
+
+f = Figure();
+a = LScene(f[1, 1])
+p = meshimage!(a, lons, lats, rotl90(values); npoints = (720, 360))
+p.transformation.transform_func[] = Makie.PointTrans{3}(UnitCartesianFromGeographic())
+# scatter!(a, cartesian_points)
+f # not entirely sure what's going on here
+# diagnostics
+# f, a, p = scatter(reduce(vcat, (view(cartesian_points, face) for face in view(faces, neighbour_inds))))
+# scatter!(query_point; color = :red, markersize = 40)
+
+query_point = LinearAlgebra.normalize(Point3(1.0, 1.0, 0.5))
+pt_inds = bowyer_watson_envelope!(Int64[], query_point, cartesian_points, faces, caps)
+three_d_points = cartesian_points[pt_inds]
+angles = [angle_between(three_d_points[1], query_point, point) for point in three_d_points]
+
+# pushfirst!(angles, 0)
+pt_inds[sortperm(angles)]
+f, a, p = scatter([query_point]; markersize = 30, color = :green, axis = (; type = LScene));
+scatter!(a, cartesian_points)
+scatter!(a, view(cartesian_points, pt_inds[sortperm(angles)]); color = eachindex(pt_inds), colormap = :RdBu, markersize = 20)
+wireframe!(a, GeometryBasics.Mesh(cartesian_points, faces); alpha = 0.3)
+
+f
+
+
+function Makie.convert_arguments(::Type{Makie.Mesh}, cap::SphericalCap)
+    offset_point = LinearAlgebra.normalize(cap.point + LinearAlgebra.normalize(Point3(cos(cap.radius), sin(cap.radius), 0.0)))
+    points = [cap.point + Rotations.AngleAxis(θ, cap.point...) * offset_point for θ in LinRange(0, 2π, 20)]
+    push!(points, cap.point)
+    faces = [GeometryBasics.TriangleFace(i, i+1, 21) for i in 1:19]
+    return (GeometryBasics.normal_mesh(points, faces),)
+end
+
+
+# Animate the natural neighbours of a loxodromic curve
+# To create a loxodromic curve we use a Mercator projection
+import GeometryOps as GO, GeoInterface as GI, GeoFormatTypes as GFT
+using CairoMakie, GeoMakie
+start_point = (-95.7129, 90.0)
+end_point = (-130.0, 0.0)
+
+geographic_path = GI.LineString([start_point, end_point]; crs = GFT.EPSG(4326))
+mercator_path = GO.reproject(geographic_path; target_crs = GFT.EPSG(3857))
+sampled_loxodromic_path = GO.segmentize(GO.LinearSegments(; max_distance = 1e5), mercator_path) # max distance in mercator coordinates!
+sampled_loxodromic_cartesian_path = GO.transform(UnitCartesianFromGeographic(), GO.reproject(sampled_loxodromic_path; target_crs = GFT.EPSG(4326))) |> GI.getpoint
+
+lines(sampled_loxodromic_cartesian_path)
+
+f, a, p = lines(sampled_loxodromic_path; axis = (; type = GeoAxis))
+lines!(a, GeoMakie.coastlines(); color = (:black, 0.5))
+f
+
+
+
+# A minimal test case to debug
+randpoints = randsphere(50)
+push!(randpoints, LinearAlgebra.normalize(Point3d(0.0, 1.0, 1.0)))
+
+zs = [zeros(length(randpoints)-1); 1.0]
+
+faces = spherical_triangulation(SphericalConvexHull(), randpoints)
+caps = map(splat(SphericalCap), (view(randpoints, face) for face in faces))
+values = map(UnitCartesianFromGeographic().(Point2.(lons, lats'))) do query_point
+    eval_laplace_coordinates(randpoints, faces, zs, query_point, caps)
+end
+
+scatter(randpoints; color = zs)
+
+heatmap(lons, lats, values)
+
+
+f = Figure();
+a = LScene(f[1, 1])
+p = meshimage!(a, lons, lats, rotl90(values); npoints = (720, 360))
+p.transformation.transform_func[] = Makie.PointTrans{3}(UnitCartesianFromGeographic())
+# scatter!(a, cartesian_points)
+f # not entirely sure what's going on here
+
+f, a, p = wireframe(GeometryBasics.Mesh(randpoints, faces))
+
+struct UnitSphereSegments <: GO.SegmentizeMethod
+    "The maximum distance between segments, in radians (solid angle units).  A nice value here for the globe is ~0.01." 
+    max_distance::Float64
+end
+
+function GO._fill_linear_kernel!(p, q, npoints = 10)
+    # perform slerp interpolation between p and q
+    # return npoints points on the great circle between p and q
+    # using spherical linear interpolation
+    # return the points as a vector of Point3{Float64}
+    ω = acos(clamp(dot(p, q), -1.0, 1.0))
+    if isapprox(ω, 0, atol=1e-8)
+        return fill(p, (npoints,))
+    end
+    return [
+        (sin((1 - t) * ω) * p + sin(t * ω) * q) / sin(ω)
+        for t in range(0, 1, length=npoints)
+    ]
+end
+
+function _segmentize(method::UnitSphereSegments, geom, T::Union{GI.LineStringTrait, GI.LinearRingTrait})
+    first_coord = GI.getpoint(geom, 1)
+    x1, y1, z1 = GI.x(first_coord), GI.y(first_coord), GI.z(first_coord)
+    new_coords = NTuple{3, Float64}[]
+    sizehint!(new_coords, GI.npoint(geom))
+    push!(new_coords, (x1, y1, z1   ))
+    for coord in Iterators.drop(GI.getpoint(geom), 1)
+        x2, y2, z2 = GI.x(coord), GI.y(coord), GI.z(coord)
+        _fill_linear_kernel!(method, new_coords, x1, y1, z1, x2, y2, z2)
+        x1, y1, z1 = x2, y2, z2
+    end 
+    return rebuild(geom, new_coords)
+end
+
+function get_face_lines(points, faces)
+    return GeometryBasics.MultiLineString(map(faces) do face
+        return GeometryBasics.LineString(vcat(
+            interpolate_line_unit_sphere(points[face[1]], points[face[2]], 10),
+            interpolate_line_unit_sphere(points[face[2]], points[face[3]], 10),
+            interpolate_line_unit_sphere(points[face[3]], points[face[1]], 10),
+        ))
+    end)
+end
+
+lines(get_face_lines(cartesian_points, faces))
\ No newline at end of file
diff --git a/docs/src/experiments/spherical_triangulation_experiments.jl b/docs/src/experiments/spherical_triangulation_experiments.jl
new file mode 100644
index 000000000..adc835bca
--- /dev/null
+++ b/docs/src/experiments/spherical_triangulation_experiments.jl
@@ -0,0 +1,180 @@
+
+include(joinpath(@__DIR__, "spherical_delaunay.jl"))
+
+# These points are known to be good points, i.e., lon, lat, alt
+points = Point3{Float64}.(JSON3.read(read(Downloads.download("https://gist.githubusercontent.com/Fil/6bc12c535edc3602813a6ef2d1c73891/raw/3ae88bf307e740ddc020303ea95d7d2ecdec0d19/points.json"), String)))
+faces = spherical_triangulation(points)
+# or
+# faces = Quickhull.quickhull(map(UnitCartesianFromGeographic(), points)) |> Quickhull.faces |> collect
+# @assert isempty(setdiff(unique!(sort!(reduce(vcat, faces))), 1:length(points)))
+
+# This is the super-cool scrollable 3D globe (though it's a bit deformed...)
+f, a, p = Makie.mesh(map(UnitCartesianFromGeographic(), points), faces; color = last.(points), colormap = Reverse(:RdBu), colorrange = (-20, 40), shading = NoShading)
+
+# We can also replicate the observable notebook almost exactly (just missing ExactPredicates):
+f, a, p = Makie.mesh(points, faces; axis = (; type = GeoAxis, dest = "+proj=bertin1953 +lon_0=-16.5 +lat_0=-42 +x_0=7.93 +y_0=0.09"), color = last.(points), colormap = Reverse(:RdBu), colorrange = (-20, 40), shading = NoShading)
+# Whoops, this doesn't look so good!  Let's try to do this more "manually" instead.
+# We'll use NaturalNeighbours.jl for this, but we first have to reconstruct the triangulation
+# with the same faces, but in a Bertin projection...
+using NaturalNeighbours
+lonlat2bertin = Proj.Transformation(GFT.EPSG(4326), GFT.ProjString("+proj=bertin1953 +type=crs +lon_0=-16.5 +lat_0=-42 +x_0=7.93 +y_0=0.09"); always_xy = true)
+lons = LinRange(-180, 180, 300)
+lats = LinRange(-90, 90, 150)
+bertin_points = lonlat2bertin.(lons, lats')
+
+projected_points = GO.reproject(GO.tuples(points), source_crs = GFT.EPSG(4326), target_crs = "+proj=bertin1953 +lon_0=-16.5 +lat_0=-42 +x_0=7.93 +y_0=0.09")
+
+ch = DelTri.convex_hull(projected_points) # assumes each point is in the triangulation
+boundary_nodes = DelTri.get_vertices(ch) 
+bertin_boundary_poly = GI.Polygon([GI.LineString(DelTri.get_points(ch)[DelTri.get_vertices(ch)])])
+
+tri = DelTri.Triangulation(projected_points, faces, boundary_nodes)
+itp = NaturalNeighbours.interpolate(tri, last.(points); derivatives = true)
+
+mat = [
+    if GO.contains(bertin_boundary_poly, (x, y))
+        itp(x, y; method = Nearest(), project = false)
+    else
+        NaN
+    end
+    for (x, y) in bertin_points
+]
+# TODO: this currently doesn't work, because some points are not inside a triangle and so cannot be located.
+# Options are:
+# 1. Reject all points outside the convex hull of the projected points.  (Tried but failed)
+# 2. ...
+
+
+
+
+
+# Now, we proceed with the script implementation which has quite a bit of debug information + plots
+pivot_ind = findfirst(isfinite, points)
+
+# debug
+point_colors = fill(:blue, length(points))
+point_colors[pivot_ind] = :red
+# end debug
+
+pivot_point = points[pivot_ind]
+necessary_rotation = #=Rotations.RotY(-π) *=# Rotations.RotY(-deg2rad(90-pivot_point[2])) * Rotations.RotZ(-deg2rad(pivot_point[1]))
+#
+net_transformation_to_corrected_cartesian = CoordinateTransformations.LinearMap(necessary_rotation) ∘ UnitCartesianFromGeographic()
+
+scatter(map(net_transformation_to_corrected_cartesian, points); color = point_colors)
+#
+stereographic_points = (StereographicFromCartesian() ∘ net_transformation_to_corrected_cartesian).(points)
+scatter(stereographic_points; color = point_colors)
+#
+triangulation_points = copy(stereographic_points)
+
+
+# Iterate through the points and find infinite/invalid points
+max_radius = 1
+really_far_idxs = Int[]
+for (i, point) in enumerate(stereographic_points)
+    m = hypot(point[1], point[2])
+    if !isfinite(m) || m > 1e32
+        push!(really_far_idxs, i)
+    elseif m > max_radius
+        max_radius = m
+    end
+end
+@show max_radius length(really_far_idxs)
+far_value = 1e6 * sqrt(max_radius)
+
+if !isempty(really_far_idxs)
+    triangulation_points[really_far_idxs] .= (Point2(far_value, 0.0),)
+end
+
+boundary_points = reverse([
+    (-far_value, -far_value),
+    (-far_value, far_value),
+    (far_value, far_value),
+    (far_value, -far_value),
+    (-far_value, -far_value),
+])
+
+boundary_nodes, pts = DelTri.convert_boundary_points_to_indices(boundary_points; existing_points=triangulation_points)
+
+triangulation = DelTri.triangulate(pts; boundary_nodes)
+# triangulation = DelTri.triangulate(triangulation_points)
+triplot(triangulation)
+DelTri.validate_triangulation(triangulation)
+#  for diagnostics, try `fig, ax, sc = triplot(triangulation, show_constrained_edges=true, show_convex_hull=true)`
+
+# First, get all the "solid" faces, ie, faces not attached to boundary nodes
+original_triangles = collect(DelTri.each_solid_triangle(triangulation))
+boundary_face_inds = findall(Base.Fix1(DelTri.is_boundary_triangle, triangulation), original_triangles)
+faces = map(CairoMakie.GeometryBasics.TriangleFace, view(original_triangles, setdiff(axes(original_triangles, 1), boundary_face_inds)))
+
+for boundary_face in view(original_triangles, boundary_face_inds)
+    push!(faces, CairoMakie.GeometryBasics.TriangleFace(map(boundary_face) do i; first(DelTri.is_boundary_node(triangulation, i)) ? pivot_ind : i end))
+end
+
+for ghost_face in DelTri.each_ghost_triangle(triangulation)
+    push!(faces, CairoMakie.GeometryBasics.TriangleFace(map(ghost_face) do i; DelTri.is_ghost_vertex(i) ? pivot_ind : i end))
+end
+# Remove degenerate triangles
+filter!(faces) do face
+    !(face[1] == face[2] || face[2] == face[3] || face[1] == face[3])
+end
+
+wireframe(CairoMakie.GeometryBasics.normal_mesh(map(UnitCartesianFromGeographic(), points), faces))
+
+grid_lons = LinRange(-180, 180, 10)
+grid_lats = LinRange(-90, 90, 10)
+lon_grid = [CairoMakie.GeometryBasics.LineString([Point{2, Float64}((grid_lons[j], -90)), Point{2, Float64}((grid_lons[j], 90))]) for j in 1:10] |> vec
+lat_grid = [CairoMakie.GeometryBasics.LineString([Point{2, Float64}((-180, grid_lats[i])), Point{2, Float64}((180, grid_lats[i]))]) for i in 1:10] |> vec
+
+sampled_grid = GO.segmentize(lat_grid; max_distance = 0.01)
+
+geographic_grid = GO.transform(UnitCartesianFromGeographic(), sampled_grid) .|> x -> GI.LineString{true, false}(x.geom)
+stereographic_grid = GO.transform(sampled_grid) do point
+    (StereographicFromCartesian() ∘ UnitCartesianFromGeographic())(point)
+end .|> x -> GI.LineString{false, false}(x.geom)
+
+
+
+fig = Figure(); ax = LScene(fig[1, 1])
+for (i, line) in enumerate(geographic_grid)
+    lines!(ax, line; linewidth = 3, color = Makie.resample_cmap(:viridis, length(stereographic_grid))[i])
+end
+fig
+
+#
+fig = Figure(); ax = LScene(fig[1, 1])
+for (i, line) in enumerate(stereographic_grid)
+    lines!(ax, line; linewidth = 3, color = Makie.resample_cmap(:viridis, length(stereographic_grid))[i])
+end
+fig
+
+all(reduce(vcat, faces) |> sort! |> unique! .== 1:length(points))
+
+faces[findall(x -> 1 in x, faces)]
+
+# try NaturalNeighbours, fail miserably
+using NaturalNeighbours
+
+itp = NaturalNeighbours.interpolate(triangulation, last.(points); derivatives = true)
+lons = LinRange(-180, 180, 360)
+lats = LinRange(-90, 90, 180)
+_x = vec([x for x in lons, _ in lats])
+_y = vec([y for _ in lons, y in lats])
+
+sibson_1_vals = itp(_x, _y; method=Sibson(1))
+
+transformed_points = GO.transform(points) do point
+    # first, spherical to Cartesian
+    longitude, latitude = deg2rad(GI.x(point)), deg2rad(GI.y(point))
+    radius = 1 # we will operate on the unit sphere
+    xyz = Point3d(
+        radius * sin(latitude) * cos(longitude),
+        radius * sin(latitude) * sin(longitude),
+        radius * cos(latitude)
+    )
+    # then, rotate Cartesian so the pivot point is at the south pole
+    rotated_xyz = Rotations.AngleAxis
+end
+
+
diff --git a/src/transformations/segmentize.jl b/src/transformations/segmentize.jl
index 623637bc5..8f199e4d5 100644
--- a/src/transformations/segmentize.jl
+++ b/src/transformations/segmentize.jl
@@ -182,7 +182,7 @@ _segmentize(method, geom) = _segmentize(method, geom, GI.trait(geom))
 This is a method which performs the common functionality for both linear and geodesic algorithms, 
 and calls out to the "kernel" function which we've defined per linesegment.
 =#
-function _segmentize(method::Union{LinearSegments, GeodesicSegments}, geom, T::Union{GI.LineStringTrait, GI.LinearRingTrait})
+function _segmentize(method::SegmentizeMethod, geom, T::Union{GI.LineStringTrait, GI.LinearRingTrait})
     first_coord = GI.getpoint(geom, 1)
     x1, y1 = GI.x(first_coord), GI.y(first_coord)
     new_coords = NTuple{2, Float64}[]