experimental UnitSpherical module that can express geometry on unit sphere

Project.toml

@@ -15,6 +15,7 @@ GeoFormatTypes = "68eda718-8dee-11e9-39e7-89f7f65f511f"
 GeoInterface = "cf35fbd7-0cd7-5166-be24-54bfbe79505f"
 GeometryOpsCore = "05efe853-fabf-41c8-927e-7063c8b9f013"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SortTileRecursiveTree = "746ee33f-1797-42c2-866d-db2fce69d14d"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
@@ -51,7 +52,8 @@ GeometryOpsCore = "=0.1.6"
 LibGEOS = "0.9.2"
 LinearAlgebra = "1"
 Proj = "1"
-SortTileRecursiveTree = "0.1"
+Random = "1"
+SortTileRecursiveTree = "0.1.2"
 StaticArrays = "1"
 Statistics = "1"
 TGGeometry = "0.1"

docs/make.jl

@@ -7,7 +7,7 @@ CairoMakie.activate!(px_per_unit = 2, type = "svg", inline = true) # TODO: make
 # import packages that activate extensions
 import FlexiJoins, LibGEOS, Proj, TGGeometry
 
-DocMeta.setdocmeta!(GeometryOps, :DocTestSetup, :(using GeometryOps; using GeometryOps.GeometryBasics); recursive=true)
+DocMeta.setdocmeta!(GeometryOps, :DocTestSetup, :(using GeometryOps; using GeometryOps.GeometryBasics; using GeometryOps.GeometryOpsCore); recursive=true)
 
 using GeoInterfaceMakie
 

src/GeometryOps.jl

@@ -40,8 +40,10 @@ include("not_implemented_yet.jl")
 include("utils/utils.jl")
 include("utils/LoopStateMachine/LoopStateMachine.jl")
 include("utils/SpatialTreeInterface/SpatialTreeInterface.jl")
+include("utils/UnitSpherical/UnitSpherical.jl")
 
-using .LoopStateMachine, .SpatialTreeInterface
+# Load utility modules in
+using .LoopStateMachine, .SpatialTreeInterface, .UnitSpherical
 
 
 include("methods/angles.jl")

src/methods/clipping/predicates.jl

@@ -12,7 +12,7 @@ module Predicates
     Return 0 if c is on (a, b) or if a == b. =#
     orient(a, b, c; exact) = _orient(booltype(exact), _tuple_point(a, Float64), _tuple_point(b, Float64), _tuple_point(c, Float64))
 
-    # If `exact` is `true`, use `ExactPredicates` to calculate the orientation.
+    # If `exact` is `true`, use `AdaptivePredicates` to calculate the orientation.
     _orient(::True, a, b, c) = AdaptivePredicates.orient2p(_tuple_point(a, Float64), _tuple_point(b, Float64), _tuple_point(c, Float64))
     # _orient(::True, a, b, c) = ExactPredicates.orient(_tuple_point(a, Float64), _tuple_point(b, Float64), _tuple_point(c, Float64))
     # If `exact` is `false`, calculate the orientation without using `ExactPredicates`.

src/utils/UnitSpherical/UnitSpherical.jl

@@ -0,0 +1,24 @@
+module UnitSpherical
+
+using CoordinateTransformations
+using StaticArrays, LinearAlgebra
+import GeoInterface as GI, GeoFormatTypes as GFT
+
+import Random
+
+# using TestItems # this is a thin package that allows TestItems.@testitem to be parsed.
+
+include("point.jl")
+include("coordinate_transforms.jl")
+include("slerp.jl")
+include("cap.jl")
+include("robustcrossproduct/RobustCrossProduct.jl")
+
+export UnitSphericalPoint, UnitSphereFromGeographic, GeographicFromUnitSphere, 
+       slerp, SphericalCap, spherical_distance
+
+# Re-export from RobustCrossProduct
+using .RobustCrossProduct: robust_cross_product
+export robust_cross_product
+
+end

src/utils/UnitSpherical/cap.jl

@@ -0,0 +1,136 @@
+# # Spherical Caps
+
+#=
+```@meta
+CollapsedDocStrings = true
+```
+
+```@docs; canonical=false
+SphericalCap
+circumcenter_on_unit_sphere
+```
+
+## What is SphericalCap?
+
+A spherical cap represents a section of a unit sphere about some point, bounded by a radius. 
+It is defined by a center point on the unit sphere and a radius (in radians).
+
+Spherical caps are used in:
+- Representing circular regions on a spherical surface
+- Approximating and bounding spherical geometries
+- Spatial indexing and filtering on the unit sphere
+- Implementing containment, intersection, and disjoint predicates
+
+The `SphericalCap` type offers multiple constructors to create caps from:
+- UnitSphericalPoint and radius
+- Geographic coordinates and radius
+- Three points on the unit sphere (circumcircle)
+
+## Examples
+
+```@example sphericalcap
+using GeometryOps
+using GeoInterface
+
+# Create a spherical cap from a point and radius
+point = UnitSphericalPoint(1.0, 0.0, 0.0)  # Point on the unit sphere
+cap = SphericalCap(point, 0.5)  # Cap with radius 0.5 radians
+```
+
+```@example sphericalcap
+# Create a spherical cap from geographic coordinates
+lat, lon = 40.0, -74.0  # New York City (approximate)
+point = GeoInterface.Point(lon, lat)
+cap = SphericalCap(point, 0.1)  # Cap with radius ~0.1 radians
+```
+
+```@example sphericalcap
+# Create a spherical cap from three points (circumcircle)
+p1 = UnitSphericalPoint(1.0, 0.0, 0.0)
+p2 = UnitSphericalPoint(0.0, 1.0, 0.0)
+p3 = UnitSphericalPoint(0.0, 0.0, 1.0)
+cap = SphericalCap(p1, p2, p3)
+```
+
+=#
+
+# Spherical cap implementation
+struct SphericalCap{T}
+    point::UnitSphericalPoint{T}
+    radius::T
+end
+
+SphericalCap(point::UnitSphericalPoint{T}, radius::Number) where T = SphericalCap{T}(point, convert(T, radius))
+SphericalCap(point, radius::Number) = SphericalCap(GI.trait(point), point, radius)
+function SphericalCap(::GI.PointTrait, point, radius::Number)
+    return SphericalCap(UnitSphereFromGeographic()(point), radius)
+end
+
+SphericalCap(geom) = SphericalCap(GI.trait(geom), geom)
+SphericalCap(t::GI.PointTrait, geom) = SphericalCap(t, geom, 0)
+# TODO: add implementations for line string and polygon traits
+# TODO: add implementations to merge two spherical caps
+function _merge(x::SphericalCap, y::SphericalCap)
+    d = spherical_distance(x.point, y.point)
+    newradius = (x.radius + y.radius + d) / 2
+    if newradius < x.radius
+        #x contains y
+        x
+    elseif newradius < y.radius
+        #y contains x
+        y
+    else
+        excenter = 0.5 * (1 + (y.radius - x.radius) / d)
+        newcenter = x.point + slerp(x.point, y.point, excenter)
+        SphericalCap(newcenter, d)
+    end
+end
+# TODO: add implementations for multitraits based on this
+
+# TODO: this returns an approximately antipodal point...
+
+# TODO: exact-predicate intersection
+# This is all inexact and thus subject to floating point error
+function _intersects(x::SphericalCap, y::SphericalCap)
+    spherical_distance(x.point, y.point) <= x.radius + y.radius
+end
+
+_disjoint(x::SphericalCap, y::SphericalCap) = !_intersects(x, y)
+
+function _contains(big::SphericalCap, small::SphericalCap)
+    dist = spherical_distance(big.point, small.point)
+    # small circle fits in big circle
+    return dist + small.radius < big.radius 
+end
+function _contains(cap::SphericalCap, point::UnitSphericalPoint)
+    spherical_distance(cap.point, point) <= cap.radius
+end
+
+
+function circumcenter_on_unit_sphere(a::UnitSphericalPoint, b::UnitSphericalPoint, c::UnitSphericalPoint)
+    LinearAlgebra.normalize(a × b + b × c + c × a)
+end
+
+"Get the circumcenter of the triangle (a, b, c) on the unit sphere.  Returns a normalized 3-vector."
+function SphericalCap(a::UnitSphericalPoint, b::UnitSphericalPoint, c::UnitSphericalPoint)
+    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::S, v_c::S, v_i::S) where S <: UnitSphericalPoint
+    # 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::S, b::S, c::S) where S <: UnitSphericalPoint
+    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

src/utils/UnitSpherical/coordinate_transforms.jl

@@ -0,0 +1,87 @@
+#=
+# Coordinate transformations
+=#
+# Coordinate transformations from lat/long to geographic and back
+"""
+    UnitSphereFromGeographic()
+
+A transformation that converts a geographic point (latitude, longitude) to a 
+[`UnitSphericalPoint`] in ℝ³.
+
+Accepts any [GeoInterface-compatible](https://github.com/JuliaGeo/GeoInterface.jl) point.
+
+## Examples
+
+```jldoctest
+julia> import GeoInterface as GI; using GeometryOps.UnitSpherical
+
+julia> UnitSphereFromGeographic()(GI.Point(45, 45))
+3-element UnitSphericalPoint{Float64} with indices SOneTo(3):
+ 0.5000000000000001
+ 0.5000000000000001
+ 0.7071067811865476
+```
+
+```jldoctest
+julia> using GeometryOps.UnitSpherical
+
+julia> UnitSphereFromGeographic()((45, 45))
+3-element UnitSphericalPoint{Float64} with indices SOneTo(3):
+ 0.5000000000000001
+ 0.5000000000000001
+ 0.7071067811865476
+```
+"""
+struct UnitSphereFromGeographic <: CoordinateTransformations.Transformation 
+end
+
+function (::UnitSphereFromGeographic)(geographic_point)
+    # Asssume that geographic_point is GeoInterface compatible
+    # Longitude is directly translatable to a spherical coordinate
+    # θ (azimuth)
+    θ = GI.x(geographic_point)
+    # The polar angle is 90 degrees minus the latitude
+    # ϕ (polar angle)
+    ϕ = 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.
+    sinϕ, cosϕ = sincosd(ϕ)
+    sinθ, cosθ = sincosd(θ)
+
+    return UnitSphericalPoint(
+        sinϕ * cosθ,
+        sinϕ * sinθ,
+        cosϕ
+    )
+end
+
+"""
+    GeographicFromUnitSphere()
+
+A transformation that converts a [`UnitSphericalPoint`](@ref) in ℝ³ to a 
+2-tuple geographic point (longitude, latitude), in degrees.
+
+Accepts any 3-element vector, but the input is assumed to be on the unit sphere.
+
+## Examples
+
+```jldoctest
+julia> using GeometryOps.UnitSpherical
+
+julia> GeographicFromUnitSphere()(UnitSphericalPoint(0.5, 0.5, 1/√(2)))
+(45.0, 44.99999999999999)
+```
+(the inaccuracy is due to the precision of the `atan` function)
+
+"""
+struct GeographicFromUnitSphere <: CoordinateTransformations.Transformation 
+end
+
+function (::GeographicFromUnitSphere)(xyz::AbstractVector)
+    @assert length(xyz) == 3 "GeographicFromUnitCartesian expects a 3D Cartesian vector"
+    x, y, z = xyz
+    return (
+        atand(y, x),
+        asind(z),
+    )
+end