Repository: JuliaEarth/CoordRefSystems.jl (julia)

Issue:
Add Sinusoidal projection
https://en.wikipedia.org/wiki/Sinusoidal_projection

Pull request:
closes #192

Patch:
diff --git a/src/CoordRefSystems.jl b/src/CoordRefSystems.jl
index d1fdf773..1fe60d69 100644
--- a/src/CoordRefSystems.jl
+++ b/src/CoordRefSystems.jl
@@ -97,6 +97,7 @@ export
   OrthoSouth,
   TransverseMercator,
   Albers,
+  Sinusoidal,
   datum,
   indomain,
 
diff --git a/src/crs/projected.jl b/src/crs/projected.jl
index 08609965..bf17cb04 100644
--- a/src/crs/projected.jl
+++ b/src/crs/projected.jl
@@ -122,6 +122,7 @@ include("projected/robinson.jl")
 include("projected/orthographic.jl")
 include("projected/transversemercator.jl")
 include("projected/albers.jl")
+include("projected/sinusoidal.jl")
 
 # ----------
 # FALLBACKS
diff --git a/src/crs/projected/sinusoidal.jl b/src/crs/projected/sinusoidal.jl
new file mode 100644
index 00000000..29bc16a8
--- /dev/null
+++ b/src/crs/projected/sinusoidal.jl
@@ -0,0 +1,71 @@
+# ------------------------------------------------------------------
+# Licensed under the MIT License. See LICENSE in the project root.
+# ------------------------------------------------------------------
+
+"""
+    Sinusoidal(x, y)
+    Sinusoidal{Datum}(x, y)
+
+Sinusoidal (Sanson-Flamsteed) coordinates in length units (default to meter)
+with a given `Datum` (default to `WGS84`).
+
+## Examples
+
+```julia
+Sinusoidal(1, 1) # add default units
+Sinusoidal(1m, 1m) # integers are converted converted to floats
+Sinusoidal(1.0km, 1.0km) # length quantities are converted to meters
+Sinusoidal(1.0m, 1.0m)
+Sinusoidal{WGS84Latest}(1.0m, 1.0m)
+```
+
+See [Sinusoidal projection](https://en.wikipedia.org/wiki/Sinusoidal_projection)
+"""
+struct Sinusoidal{Datum,Shift,M<:Met} <: Projected{Datum,Shift}
+  x::M
+  y::M
+end
+
+Sinusoidal{Datum,Shift}(x::M, y::M) where {Datum,Shift,M<:Met} = Sinusoidal{Datum,Shift,float(M)}(x, y)
+Sinusoidal{Datum,Shift}(x::Met, y::Met) where {Datum,Shift} = Sinusoidal{Datum,Shift}(promote(x, y)...)
+Sinusoidal{Datum,Shift}(x::Len, y::Len) where {Datum,Shift} = Sinusoidal{Datum,Shift}(uconvert(m, x), uconvert(m, y))
+Sinusoidal{Datum,Shift}(x::Number, y::Number) where {Datum,Shift} = Sinusoidal{Datum,Shift}(addunit(x, m), addunit(y, m))
+
+Sinusoidal{Datum}(args...) where {Datum} = Sinusoidal{Datum,Shift()}(args...)
+
+Sinusoidal(args...) = Sinusoidal{WGS84Latest}(args...)
+
+Base.convert(::Type{Sinusoidal{Datum,Shift,M}}, coords::Sinusoidal{Datum,Shift}) where {Datum,Shift,M} =
+  Sinusoidal{Datum,Shift,M}(coords.x, coords.y)
+
+constructor(::Type{<:Sinusoidal{Datum,Shift}}) where {Datum,Shift} = Sinusoidal{Datum,Shift}
+
+lentype(::Type{<:Sinusoidal{Datum,Shift,M}}) where {Datum,Shift,M} = M
+
+==(coords₁::Sinusoidal{Datum,Shift}, coords₂::Sinusoidal{Datum,Shift}) where {Datum,Shift} =
+  coords₁.x == coords₂.x && coords₁.y == coords₂.y
+
+# ------------
+# CONVERSIONS
+# ------------
+
+function formulas(::Type{<:Sinusoidal}, ::Type{T}) where {T}
+  fx(λ, ϕ) = λ * cos(ϕ)
+
+  fy(λ, ϕ) = ϕ
+
+  fx, fy
+end
+
+function backward(::Type{<:Sinusoidal}, x, y)
+  ϕ = y
+  λ = x / cos(ϕ)
+
+  λ, ϕ
+end
+
+# ----------
+# FALLBACKS
+# ----------
+
+Base.convert(::Type{Sinusoidal}, coords::CRS{Datum}) where {Datum} = convert(Sinusoidal{Datum}, coords)


Test patch:
diff --git a/test/crs/constructors.jl b/test/crs/constructors.jl
index a03bbe7a..161373a0 100644
--- a/test/crs/constructors.jl
+++ b/test/crs/constructors.jl
@@ -650,5 +650,31 @@
       @test_throws ArgumentError A(T(1) * m, T(1) * s)
       @test_throws ArgumentError A(T(1) * s, T(1) * s)
     end
+
+    @testset "Sinusoidal" begin
+      @test Sinusoidal(T(1), T(1)) == Sinusoidal(T(1) * m, T(1) * m)
+      @test Sinusoidal(T(1) * m, 1 * m) == Sinusoidal(T(1) * m, T(1) * m)
+      @test Sinusoidal(T(1) * km, T(1) * km) == Sinusoidal(T(1000) * m, T(1000) * m)
+
+      c = Sinusoidal(T(1), T(1))
+      @test sprint(show, c) == "Sinusoidal{WGS84Latest}(x: 1.0 m, y: 1.0 m)"
+      if T === Float32
+        @test sprint(show, MIME("text/plain"), c) == """
+        Sinusoidal{WGS84Latest} coordinates
+        ├─ x: 1.0f0 m
+        └─ y: 1.0f0 m"""
+      else
+        @test sprint(show, MIME("text/plain"), c) == """
+        Sinusoidal{WGS84Latest} coordinates
+        ├─ x: 1.0 m
+        └─ y: 1.0 m"""
+      end
+
+      # error: invalid units for coordinates
+      @test_throws ArgumentError Sinusoidal(T(1), T(1) * m)
+      @test_throws ArgumentError Sinusoidal(T(1) * s, T(1) * m)
+      @test_throws ArgumentError Sinusoidal(T(1) * m, T(1) * s)
+      @test_throws ArgumentError Sinusoidal(T(1) * s, T(1) * s)
+    end
   end
 end
diff --git a/test/crs/conversions.jl b/test/crs/conversions.jl
index a160fe52..2c7bb179 100644
--- a/test/crs/conversions.jl
+++ b/test/crs/conversions.jl
@@ -1507,6 +1507,38 @@
       @inferred convert(LatLon{NAD83}, c2)
     end
 
+    @testset "LatLon <> Sinusoidal" begin
+      c1 = LatLon(T(45), T(90))
+      c2 = convert(Sinusoidal, c1)
+      @test allapprox(c2, Sinusoidal(T(7.084329013634147e6), T(5.009377085697311e6)))
+      c3 = convert(LatLon, c2)
+      @test allapprox(c3, c1)
+
+      c1 = LatLon(-T(45), T(90))
+      c2 = convert(Sinusoidal, c1)
+      @test allapprox(c2, Sinusoidal(T(7.084329013634147e6), -T(5.009377085697311e6)))
+      c3 = convert(LatLon, c2)
+      @test allapprox(c3, c1)
+
+      c1 = LatLon(T(45), -T(90))
+      c2 = convert(Sinusoidal, c1)
+      @test allapprox(c2, Sinusoidal(-T(7.084329013634147e6), T(5.009377085697311e6)))
+      c3 = convert(LatLon, c2)
+      @test allapprox(c3, c1)
+
+      c1 = LatLon(-T(45), -T(90))
+      c2 = convert(Sinusoidal, c1)
+      @test allapprox(c2, Sinusoidal(-T(7.084329013634147e6), -T(5.009377085697311e6)))
+      c3 = convert(LatLon, c2)
+      @test allapprox(c3, c1)
+
+      # type stability
+      c1 = LatLon(T(45), T(90))
+      c2 = Sinusoidal(T(7.084329013634147e6), T(5.009377085697311e6))
+      @inferred convert(Sinusoidal, c1)
+      @inferred convert(LatLon, c2)
+    end
+
     @testset "LatLon <> UTM" begin
       UTMNorth32 = utmnorth(32)
       UTMSouth59 = utmsouth(59)
diff --git a/test/crs/crsapi.jl b/test/crs/crsapi.jl
index c71532a6..f15a776c 100644
--- a/test/crs/crsapi.jl
+++ b/test/crs/crsapi.jl
@@ -170,6 +170,8 @@
     @test CoordRefSystems.raw(c) == (T(1), T(2))
     c = AlbersUS(T(1), T(2))
     @test CoordRefSystems.raw(c) == (T(1), T(2))
+    c = Sinusoidal(T(1), T(2))
+    @test CoordRefSystems.raw(c) == (T(1), T(2))
   end
 
   @testset "units" begin
@@ -215,6 +217,8 @@
     @test CoordRefSystems.units(c) == (m, m)
     c = AlbersUS(T(1), T(2))
     @test CoordRefSystems.units(c) == (m, m)
+    c = Sinusoidal(T(1), T(2))
+    @test CoordRefSystems.units(c) == (m, m)
   end
 
   @testset "constructor" begin
@@ -258,6 +262,8 @@
     @test CoordRefSystems.constructor(c) === ShiftedMercator
     c = AlbersUS(T(1), T(2))
     @test CoordRefSystems.constructor(c) === AlbersUS
+    c = Sinusoidal(T(1), T(2))
+    @test CoordRefSystems.constructor(c) === Sinusoidal{WGS84Latest,CoordRefSystems.Shift()}
   end
 
   @testset "reconstruct" begin
@@ -324,6 +330,9 @@
     c = AlbersUS(T(1), T(2))
     rv = CoordRefSystems.raw(c)
     @test CoordRefSystems.reconstruct(typeof(c), rv) == c
+    c = Sinusoidal(T(1), T(2))
+    rv = CoordRefSystems.raw(c)
+    @test CoordRefSystems.reconstruct(typeof(c), rv) == c
   end
 
   @testset "lentype" begin
@@ -367,6 +376,8 @@
     @test CoordRefSystems.lentype(c) == typeof(T(1) * m)
     c = AlbersUS(T(1), T(2))
     @test CoordRefSystems.lentype(c) == typeof(T(1) * m)
+    c = Sinusoidal(T(1), T(2))
+    @test CoordRefSystems.lentype(c) == typeof(T(1) * m)
   end
 
   @testset "mactype" begin
@@ -410,6 +421,8 @@
     @test CoordRefSystems.mactype(c) == T
     c = AlbersUS(T(1), T(2))
     @test CoordRefSystems.mactype(c) == T
+    c = Sinusoidal(T(1), T(2))
+    @test CoordRefSystems.mactype(c) == T
   end
 
   @testset "equality" begin
@@ -434,6 +447,7 @@
     equaltest(ShiftedTM)
     equaltest(ShiftedMercator)
     equaltest(AlbersUS)
+    equaltest(Sinusoidal)
   end
 
   @testset "isapprox" begin
@@ -462,6 +476,7 @@
     ShiftedMercatorWGS = CoordRefSystems.shift(Mercator{WGS84{1762}}, lonₒ=15.0°, xₒ=200.0m, yₒ=200.0m)
     ShiftedMercatorITRF = CoordRefSystems.shift(Mercator{ITRF{2008}}, lonₒ=15.0°, xₒ=200.0m, yₒ=200.0m)
     isapproxtest3D(ShiftedMercatorWGS, ShiftedMercatorITRF)
+    isapproxtest3D(Sinusoidal)
   end
 
   @testset "tol" begin
diff --git a/test/crs/domains.jl b/test/crs/domains.jl
index f61fe867..8e626e24 100644
--- a/test/crs/domains.jl
+++ b/test/crs/domains.jl
@@ -255,4 +255,29 @@
       end
     end
   end
+
+  @testset "Sinusoidal forward" begin
+    for lat in T.(-90:90), lon in T.(-180:180)
+      c1 = LatLon(lat, lon)
+      if indomain(Sinusoidal, c1)
+        c2 = convert(Sinusoidal, c1)
+        @test isfinite(c2.x)
+        @test isfinite(c2.y)
+      else
+        @test_throws ArgumentError convert(Sinusoidal, c1)
+      end
+    end
+  end
+
+  @testset "Sinusoidal inverse" begin
+    # coordinates at the singularity of the projection (lat ≈ ±90) cannot be inverted
+    for lat in T.(-89:89), lon in T.(-180:180)
+      c1 = LatLon(lat, lon)
+      if indomain(Sinusoidal, c1)
+        c2 = convert(Sinusoidal, c1)
+        c3 = convert(LatLon, c2)
+        @test allapprox(c3, c1)
+      end
+    end
+  end
 end
diff --git a/test/crs/mactype.jl b/test/crs/mactype.jl
index 50ae3228..a4c9e2e8 100644
--- a/test/crs/mactype.jl
+++ b/test/crs/mactype.jl
@@ -167,5 +167,13 @@
     c1 = ShiftedMercator(1.0f0, 1.0f0)
     c2 = convert(C, c1)
     @test c2 isa C
+
+    C = Sinusoidal{WGS84Latest,CoordRefSystems.Shift(),Met{T}}
+    c1 = Sinusoidal(1.0, 1.0)
+    c2 = convert(C, c1)
+    @test c2 isa C
+    c1 = Sinusoidal(1.0f0, 1.0f0)
+    c2 = convert(C, c1)
+    @test c2 isa C
   end
 end
diff --git a/test/crs/rand.jl b/test/crs/rand.jl
index aca5698c..9950b286 100644
--- a/test/crs/rand.jl
+++ b/test/crs/rand.jl
@@ -72,5 +72,8 @@
 
     randtest(Albers{23.0°,29.5°,45.5°,WGS84Latest})
     randtest(Albers{23.0°,29.5°,45.5°})
+
+    randtest(Sinusoidal{WGS84Latest})
+    randtest(Sinusoidal)
   end
 end