[WIP] Spherical triangulation and interpolation

docs/src/experiments/spherical_delaunay_stereographic.jl

@@ -0,0 +1,332 @@
+#=
+# Spherical Delaunay triangulation using stereographic projections
+
+This is the approach which d3-geo-voronoi and friends use. The alternative is STRIPACK which computes in 3D on the sphere.
+The 3D approach is basically that the 3d convex hull of a set of points on the sphere is its Delaunay triangulation.
+=#
+
+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
+
+function delaunay_triangulate_spherical(input_points; facetype = CairoMakie.GeometryBasics.TriangleFace)
+ # @assert GI.crstrait(first(input_points)) isa GI.AbstractGeographicCRSTrait
+ points = GO.tuples(input_points)
+
+ 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
+ boundary_faces = Iterators.filter(Base.Fix1(DelTri.is_boundary_triangle, triangulation), DelTri.each_solid_triangle(triangulation)) |> collect
+ faces = map(facetype, DelTri.each_solid_triangle(triangulation))
+
+ for ghost_face in DelTri.each_ghost_triangle(triangulation)
+ push!(faces, facetype(map(ghost_face) do i; i ≤ 0 ? pivot_ind : i end))
+ end
+
+ return faces
+end
+
+# 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 = delaunay_triangulate_spherical(points)
+
+# This is the super-cool scrollable 3D globe (though it's a bit deformed... :D)
+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())
+ 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(triangulation_points; boundary_nodes)
+triangulation = DelTri.triangulate(triangulation_points)
+triplot(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
+boundary_faces = Iterators.filter(Base.Fix1(DelTri.is_boundary_triangle, triangulation), DelTri.each_solid_triangle(triangulation)) |> collect
+faces = map(CairoMakie.GeometryBasics.TriangleFace, DelTri.each_solid_triangle(triangulation))
+
+for ghost_face in DelTri.each_ghost_triangle(triangulation)
+ push!(faces, CairoMakie.GeometryBasics.TriangleFace(map(ghost_face) do i; i ≤ 0 ? pivot_ind : i end))
+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))
+
+# 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
+
+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
+
+
+
+
+function randsphericalangles(n)
+ θ = 2π .* rand(n)
+ ϕ = acos.(2 .* rand(n) .- 1)
+ return Point2.(θ, ϕ)
+end
+
+function randsphere(n)
+ θϕ = 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