Skip to content

Commit fc1121c

Browse files
asinghvi17asinghvi17
committed
make folder
1 parent 5560003 commit fc1121c

File tree

6 files changed

+441
-214
lines changed

6 files changed

+441
-214
lines changed

‎docs/src/experiments/UnitSpherical.jl/UnitSpherical.jl

+79
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,79 @@
1+
# module UnitSpherical
2+
3+
using CoordinateTransformations
4+
using StaticArrays, LinearAlgebra
5+
import GeoInterface as GI, GeoFormatTypes as GFT
6+
7+
import Random
8+
9+
using TestItems # this is a thin package that allows TestItems.@testitem to be parsed.
10+
11+
12+
13+
# Spherical cap implementation
14+
struct SphericalCap{T}
15+
point::UnitSphericalPoint{T}
16+
radius::T
17+
end
18+
19+
SphericalCap(point::UnitSphericalPoint{T}, radius::Number) where T = SphericalCap{T}(point, convert(T, radius))
20+
SphericalCap(point, radius::Number) = SphericalCap(GI.trait(point), point, radius)
21+
function SphericalCap(::GI.PointTrait, point, radius::Number)
22+
return SphericalCap(UnitSphereFromGeographic()(point), radius)
23+
end
24+
25+
SphericalCap(geom) = SphericalCap(GI.trait(geom), geom)
26+
SphericalCap(t::GI.PointTrait, geom) = SphericalCap(t, geom, 0)
27+
# TODO: add implementations for line string and polygon traits
28+
# TODO: add implementations to merge two spherical caps
29+
# TODO: add implementations for multitraits based on this
30+
31+
# TODO: this returns an approximately antipodal point...
32+
33+
34+
# TODO: exact-predicate intersection
35+
# This is all inexact and thus subject to floating point error
36+
function _intersects(x::SphericalCap, y::SphericalCap)
37+
spherical_distance(x.point, y.point) <= max(x.radius, y.radius)
38+
end
39+
40+
_disjoint(x::SphericalCap, y::SphericalCap) = !_intersects(x, y)
41+
42+
function _contains(big::SphericalCap, small::SphericalCap)
43+
dist = spherical_distance(big.point, small.point)
44+
return dist < big.radius #=small.point in big=# &&
45+
dist + small.radius < big.radius # small circle fits in big circle
46+
end
47+
48+
49+
function circumcenter_on_unit_sphere(a::UnitSphericalPoint, b::UnitSphericalPoint, c::UnitSphericalPoint)
50+
LinearAlgebra.normalize(a × b + b × c + c × a)
51+
end
52+
53+
"Get the circumcenter of the triangle (a, b, c) on the unit sphere. Returns a normalized 3-vector."
54+
function SphericalCap(a::UnitSphericalPoint, b::UnitSphericalPoint, c::UnitSphericalPoint)
55+
circumcenter = circumcenter_on_unit_sphere(a, b, c)
56+
circumradius = spherical_distance(a, circumcenter)
57+
return SphericalCap(circumcenter, circumradius)
58+
end
59+
60+
function _is_ccw_unit_sphere(v_0::S, v_c::S, v_i::S) where S <: UnitSphericalPoint
61+
# checks if the smaller interior angle for the great circles connecting u-v and v-w is CCW
62+
return(LinearAlgebra.dot(LinearAlgebra.cross(v_c - v_0,v_i - v_c), v_i) < 0)
63+
end
64+
65+
function angle_between(a::S, b::S, c::S) where S <: UnitSphericalPoint
66+
ab = b - a
67+
bc = c - b
68+
norm_dot = (ab ⋅ bc) / (LinearAlgebra.norm(ab) * LinearAlgebra.norm(bc))
69+
angle = acos(clamp(norm_dot, -1.0, 1.0))
70+
if _is_ccw_unit_sphere(a, b, c)
71+
return angle
72+
else
73+
return 2π - angle
74+
end
75+
end
76+
77+
78+
79+
# end

‎docs/src/experiments/UnitSpherical.jl/cap.jl

+66
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,66 @@
1+
#=
2+
# Spherical caps
3+
=#
4+
# Spherical cap implementation
5+
struct SphericalCap{T}
6+
point::UnitSphericalPoint{T}
7+
radius::T
8+
end
9+
10+
SphericalCap(point::UnitSphericalPoint{T}, radius::Number) where T = SphericalCap{T}(point, convert(T, radius))
11+
SphericalCap(point, radius::Number) = SphericalCap(GI.trait(point), point, radius)
12+
function SphericalCap(::GI.PointTrait, point, radius::Number)
13+
return SphericalCap(UnitSphereFromGeographic()(point), radius)
14+
end
15+
16+
SphericalCap(geom) = SphericalCap(GI.trait(geom), geom)
17+
SphericalCap(t::GI.PointTrait, geom) = SphericalCap(t, geom, 0)
18+
# TODO: add implementations for line string and polygon traits
19+
# TODO: add implementations to merge two spherical caps
20+
# TODO: add implementations for multitraits based on this
21+
22+
# TODO: this returns an approximately antipodal point...
23+
24+
25+
# TODO: exact-predicate intersection
26+
# This is all inexact and thus subject to floating point error
27+
function _intersects(x::SphericalCap, y::SphericalCap)
28+
spherical_distance(x.point, y.point) <= max(x.radius, y.radius)
29+
end
30+
31+
_disjoint(x::SphericalCap, y::SphericalCap) = !_intersects(x, y)
32+
33+
function _contains(big::SphericalCap, small::SphericalCap)
34+
dist = spherical_distance(big.point, small.point)
35+
return dist < big.radius #=small.point in big=# &&
36+
dist + small.radius < big.radius # small circle fits in big circle
37+
end
38+
39+
40+
function circumcenter_on_unit_sphere(a::UnitSphericalPoint, b::UnitSphericalPoint, c::UnitSphericalPoint)
41+
LinearAlgebra.normalize(a × b + b × c + c × a)
42+
end
43+
44+
"Get the circumcenter of the triangle (a, b, c) on the unit sphere. Returns a normalized 3-vector."
45+
function SphericalCap(a::UnitSphericalPoint, b::UnitSphericalPoint, c::UnitSphericalPoint)
46+
circumcenter = circumcenter_on_unit_sphere(a, b, c)
47+
circumradius = spherical_distance(a, circumcenter)
48+
return SphericalCap(circumcenter, circumradius)
49+
end
50+
51+
function _is_ccw_unit_sphere(v_0::S, v_c::S, v_i::S) where S <: UnitSphericalPoint
52+
# checks if the smaller interior angle for the great circles connecting u-v and v-w is CCW
53+
return(LinearAlgebra.dot(LinearAlgebra.cross(v_c - v_0,v_i - v_c), v_i) < 0)
54+
end
55+
56+
function angle_between(a::S, b::S, c::S) where S <: UnitSphericalPoint
57+
ab = b - a
58+
bc = c - b
59+
norm_dot = (ab ⋅ bc) / (LinearAlgebra.norm(ab) * LinearAlgebra.norm(bc))
60+
angle = acos(clamp(norm_dot, -1.0, 1.0))
61+
if _is_ccw_unit_sphere(a, b, c)
62+
return angle
63+
else
64+
return 2π - angle
65+
end
66+
end

‎docs/src/experiments/UnitSpherical.jl/coordinate_transforms.jl

+81
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,81 @@
1+
#=
2+
# Coordinate transformations
3+
=#
4+
# Coordinate transformations from lat/long to geographic and back
5+
"""
6+
UnitSphereFromGeographic()
7+
8+
A transformation that converts a geographic point (latitude, longitude) to a
9+
[`UnitSphericalPoint`] in ℝ³.
10+
11+
Accepts any [GeoInterface-compatible](https://github.com/JuliaGeo/GeoInterface.jl) point.
12+
13+
## Examples
14+
15+
```jldoctest
16+
julia> UnitSphereFromGeographic()(GI.Point(45, 45))
17+
3-element UnitSphericalPoint{Float64} with indices SOneTo(3):
18+
0.5000000000000001
19+
0.5000000000000001
20+
0.7071067811865476
21+
```
22+
23+
```jldoctest
24+
julia> UnitSphereFromGeographic()((45, 45))
25+
3-element UnitSphericalPoint{Float64} with indices SOneTo(3):
26+
0.5000000000000001
27+
0.5000000000000001
28+
0.7071067811865476
29+
```
30+
"""
31+
struct UnitSphereFromGeographic <: CoordinateTransformations.Transformation
32+
end
33+
34+
function (::UnitSphereFromGeographic)(geographic_point)
35+
# Asssume that geographic_point is GeoInterface compatible
36+
# Longitude is directly translatable to a spherical coordinate
37+
# θ (azimuth)
38+
θ = GI.x(geographic_point)
39+
# The polar angle is 90 degrees minus the latitude
40+
# ϕ (polar angle)
41+
ϕ = 90 - GI.y(geographic_point)
42+
# Since this is the unit sphere, the radius is assumed to be 1,
43+
# and we don't need to multiply by it.
44+
sinϕ, cosϕ = sincosd(ϕ)
45+
sinθ, cosθ = sincosd(θ)
46+
47+
return UnitSphericalPoint(
48+
sinϕ * cosθ,
49+
sinϕ * sinθ,
50+
cosϕ
51+
)
52+
end
53+
54+
"""
55+
GeographicFromUnitSphere()
56+
57+
A transformation that converts a [`UnitSphericalPoint`](@ref) in ℝ³ to a
58+
2-tuple geographic point (longitude, latitude), in degrees.
59+
60+
Accepts any 3-element vector, but the input is assumed to be on the unit sphere.
61+
62+
## Examples
63+
64+
```jldoctest
65+
julia> GeographicFromUnitSphere()(UnitSphericalPoint(0.5, 0.5, 1/√(2)))
66+
(45.0, 44.99999999999999)
67+
```
68+
(the inaccuracy is due to the precision of `atan` function)
69+
70+
"""
71+
struct GeographicFromUnitSphere <: CoordinateTransformations.Transformation
72+
end
73+
74+
function (::GeographicFromUnitSphere)(xyz::AbstractVector)
75+
@assert length(xyz) == 3 "GeographicFromUnitCartesian expects a 3D Cartesian vector"
76+
x, y, z = xyz
77+
return (
78+
atand(y, x),
79+
asind(z),
80+
)
81+
end

‎docs/src/experiments/UnitSpherical.jl/point.jl

+163
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,163 @@
1+
#=
2+
3+
4+
# UnitSphericalPoint
5+
6+
This file defines the [`UnitSphericalPoint`](@ref) type, which is
7+
a three-dimensional Cartesian point on the unit 2-sphere (i.e., of radius 1).
8+
9+
This file contains the full implementation of the type as well as a `spherical_distance` function
10+
that computes the great-circle distance between two points on the unit sphere.
11+
12+
```@docs; canonical=false
13+
UnitSphericalPoint
14+
spherical_distance
15+
```
16+
=#
17+
18+
# ## Type definition and constructors
19+
"""
20+
UnitSphericalPoint(v)
21+
22+
A unit spherical point, i.e., point living on the 2-sphere (𝕊²),
23+
represented as Cartesian coordinates in ℝ³.
24+
25+
This currently has no support for heights, only going from lat long to spherical
26+
and back again.
27+
28+
## Examples
29+
30+
```jldoctest
31+
julia> UnitSphericalPoint(1, 0, 0)
32+
UnitSphericalPoint(1.0, 0.0, 0.0)
33+
```
34+
35+
"""
36+
struct UnitSphericalPoint{T} <: StaticArrays.FieldVector{3, T}
37+
x::T
38+
y::T
39+
z::T
40+
end
41+
42+
UnitSphericalPoint{T}(v::SVector{3, T}) where T = UnitSphericalPoint{T}(v...)
43+
UnitSphericalPoint(v::NTuple{3, T}) where T = UnitSphericalPoint{T}(v...)
44+
UnitSphericalPoint{T}(v::NTuple{3, T}) where T = UnitSphericalPoint{T}(v...)
45+
UnitSphericalPoint{T}(v::SVector{3, T}) where T = UnitSphericalPoint{T}(v...)
46+
UnitSphericalPoint(v::SVector{3, T}) where T = UnitSphericalPoint{T}(v...)
47+
## handle the 2-tuple case specifically
48+
UnitSphericalPoint(v::NTuple{2, T}) where T = UnitSphereFromGeographic()(v)
49+
## handle the GeoInterface case, this is the catch-all method
50+
UnitSphericalPoint(v) = UnitSphericalPoint(GI.trait(v), v)
51+
UnitSphericalPoint(::GI.PointTrait, v) = UnitSphereFromGeographic()(v) # since it works on any GI compatible point
52+
## finally, handle the case where a vector is passed in
53+
## we may want it to go to the geographic pipeline _or_ direct materialization
54+
Base.@propagate_inbounds function UnitSphericalPoint(v::AbstractVector{T}) where T
55+
if length(v) == 3
56+
UnitSphericalPoint{T}(v[1], v[2], v[3])
57+
elseif length(v) == 2
58+
UnitSphereFromGeographic()(v)
59+
else
60+
@boundscheck begin
61+
throw(ArgumentError("""
62+
Passed a vector of length `$(length(v))` to the `UnitSphericalPoint` constructor,
63+
which only accepts vectors of lengths:
64+
- **3** (assumed to be on the unit sphere)
65+
- **2** (assumed to be geographic lat/long)
66+
"""))
67+
end
68+
end
69+
end
70+
71+
Base.show(io::IO, p::UnitSphericalPoint) = print(io, "UnitSphericalPoint($(p.x), $(p.y), $(p.z))")
72+
73+
# ## Interface implementations
74+
75+
# StaticArraysCore.jl interface implementation
76+
Base.setindex(p::UnitSphericalPoint, args...) = throw(ArgumentError("`setindex!` on a UnitSphericalPoint is not permitted as it is static."))
77+
StaticArrays.similar_type(::Type{<: UnitSphericalPoint}, ::Type{Eltype}, ::Size{(3,)}) where Eltype = UnitSphericalPoint{Eltype}
78+
# Base math implementation (really just forcing re-wrapping)
79+
# Base.:(*)(a::UnitSphericalPoint, b::UnitSphericalPoint) = a .* b
80+
function Base.broadcasted(f, a::AbstractArray{T}, b::UnitSphericalPoint) where {T <: UnitSphericalPoint}
81+
return Base.broadcasted(f, a, (b,))
82+
end
83+
Base.isnan(p::UnitSphericalPoint) = any(isnan, p)
84+
Base.isinf(p::UnitSphericalPoint) = any(isinf, p)
85+
Base.isfinite(p::UnitSphericalPoint) = all(isfinite, p)
86+
87+
88+
# GeoInterface implementation
89+
## Traits:
90+
GI.trait(::UnitSphericalPoint) = GI.PointTrait()
91+
GI.geomtrait(::UnitSphericalPoint) = GI.PointTrait()
92+
## Coordinate traits:
93+
GI.is3d(::GI.PointTrait, ::UnitSphericalPoint) = true
94+
GI.ismeasured(::GI.PointTrait, ::UnitSphericalPoint) = false
95+
## Accessors:
96+
GI.ncoord(::GI.PointTrait, ::UnitSphericalPoint) = 3
97+
GI.getcoord(::GI.PointTrait, p::UnitSphericalPoint) = p[i]
98+
## Metadata (CRS, extent, etc)
99+
GI.crs(::UnitSphericalPoint) = GFT.ProjString("+proj=cart +R=1 +type=crs") # TODO: make this a full WKT definition
100+
# TODO: extent is a little tricky - do we do a spherical cap or an Extents.Extent?
101+
102+
# ## Spherical distance
103+
"""
104+
spherical_distance(x::UnitSphericalPoint, y::UnitSphericalPoint)
105+
106+
Compute the spherical distance between two points on the unit sphere.
107+
Returns a `Number`, usually Float64 but that depends on the input type.
108+
109+
# Extended help
110+
111+
## Doctests
112+
113+
```jldoctest
114+
julia> spherical_distance(UnitSphericalPoint(1, 0, 0), UnitSphericalPoint(0, 1, 0))
115+
1.5707963267948966
116+
```
117+
118+
```jldoctest
119+
julia> spherical_distance(UnitSphericalPoint(1, 0, 0), UnitSphericalPoint(1, 0, 0))
120+
0.0
121+
```
122+
"""
123+
spherical_distance(x::UnitSphericalPoint, y::UnitSphericalPoint) = acos(clamp(x ⋅ y, -1.0, 1.0))
124+
125+
# ## Random points
126+
Random.rand(rng::Random.AbstractRNG, ::Random.SamplerType{UnitSphericalPoint}) = rand(rng, UnitSphericalPoint{Float64})
127+
function Random.rand(rng::Random.AbstractRNG, ::Random.SamplerType{UnitSphericalPoint{T}}) where T <: Number
128+
θ = 2π * rand(rng, T)
129+
ϕ = acos(2 * rand(rng, T) - 1)
130+
sinθ, cosθ = sincos(θ)
131+
sinϕ, cosϕ = sincos(ϕ)
132+
return UnitSphericalPoint(
133+
sinϕ * cosθ,
134+
sinϕ * sinθ,
135+
cosϕ
136+
)
137+
end
138+
139+
# ## Tests
140+
141+
@testitem "UnitSphericalPoint constructor" begin
142+
using GeometryOps.UnitSpherical
143+
import GeoInterface as GI
144+
145+
northpole = UnitSphericalPoint{Float64}(1, 0, 0)
146+
# test that the constructor works for a vector of length 3
147+
@test UnitSphericalPoint((1, 0, 0)) == northpole
148+
@test UnitSphericalPoint(SVector(1, 0, 0)) == northpole
149+
@test UnitSphericalPoint([1, 0, 0]) == northpole
150+
# test that the constructor works for a tuple of length 2
151+
# and interprets such a thing as a geographic point
152+
@test UnitSphericalPoint((90, 0)) == northpole
153+
@test UnitSphericalPoint([90, 0]) == northpole
154+
@test UnitSphericalPoint(GI.Point((90, 0))) == northpole
155+
156+
157+
end
158+
159+
#=
160+
```@meta
161+
CollapsedDocStrings = true
162+
```
163+
=#

0 commit comments

Comments
 (0)