Add UTM projection (#28)

* Add 'UTM' projection

* Fix typo

* Apply suggestions

* Format code

* Apply suggestions

* Apply suggestions

* Apply suggestions from code review

Co-authored-by: Júlio Hoffimann <[email protected]>

* Add tests

---------

Co-authored-by: Júlio Hoffimann <[email protected]>

src/Cartography.jl

@@ -68,9 +68,16 @@ export
   Robinson,
   OrthoNorth,
   OrthoSouth,
+  UTM,
+  UTMNorth,
+  UTMSouth,
   datum,
   indomain,
 
+  # hemispheres
+  North,
+  South,
+
   # codes
   EPSG,
   ESRI

src/crs/projected.jl

@@ -43,6 +43,7 @@ include("projected/winkeltripel.jl")
 include("projected/robinson.jl")
 include("projected/orthographic.jl")
 include("projected/transversemercator.jl")
+include("projected/utm.jl")
 
 # ----------
 # FALLBACKS

src/crs/projected/utm.jl

@@ -0,0 +1,133 @@
+# ------------------------------------------------------------------
+# Licensed under the MIT License. See LICENSE in the project root.
+# ------------------------------------------------------------------
+
+"""
+    Hemisphere
+
+Hemisphere of the Earth.
+"""
+abstract type Hemisphere end
+
+"""
+    North
+
+Northern hemisphere.
+"""
+abstract type North <: Hemisphere end
+
+"""
+    South
+
+Southern hemisphere.
+"""
+abstract type South <: Hemisphere end
+
+"""
+    UTM{Hemisphere,Zone,Datum}
+
+UTM (Universal Transverse Mercator) CRS in `Hemisphere` with `Zone` (1 ≤ Zone ≤ 60) and a given `Datum`.
+"""
+struct UTM{Hemisphere,Zone,Datum,M<:Met} <: Projected{Datum}
+  x::M
+  y::M
+  function UTM{Hemisphere,Zone,Datum}(x::M, y::M) where {Hemisphere,Zone,Datum,M<:Met}
+    if !(1 ≤ Zone ≤ 60)
+      throw(ArgumentError("the UTM zone must be an integer between 1 and 60"))
+    end
+    new{Hemisphere,Zone,Datum,float(M)}(x, y)
+  end
+end
+
+UTM{Hemisphere,Zone,Datum}(x::Met, y::Met) where {Hemisphere,Zone,Datum} = UTM{Hemisphere,Zone,Datum}(promote(x, y)...)
+UTM{Hemisphere,Zone,Datum}(x::Len, y::Len) where {Hemisphere,Zone,Datum} =
+  UTM{Hemisphere,Zone,Datum}(uconvert(u"m", x), uconvert(u"m", y))
+UTM{Hemisphere,Zone,Datum}(x::Number, y::Number) where {Hemisphere,Zone,Datum} =
+  UTM{Hemisphere,Zone,Datum}(addunit(x, u"m"), addunit(y, u"m"))
+
+UTM{Hemisphere,Zone}(args...) where {Hemisphere,Zone} = UTM{Hemisphere,Zone,WGS84Latest}(args...)
+
+"""
+    UTMNorth{zone}(x, y)
+    UTMNorth{zone,Datum}(x, y)
+
+UTM (Universal Transverse Mercator) North coordinates in length units (default to meter)
+with `zone` (1 ≤ zone ≤ 60) and a given `Datum` (default to `WGS84`).
+
+## Examples
+
+```julia
+UTMNorth{1}(1, 1) # add default units
+UTMNorth{1}(1u"m", 1u"m") # integers are converted converted to floats
+UTMNorth{1}(1.0u"km", 1.0u"km") # length quantities are converted to meters
+UTMNorth{1}(1.0u"m", 1.0u"m")
+UTMNorth{1,WGS84Latest}(1.0u"m", 1.0u"m")
+```
+"""
+const UTMNorth{Zone,Datum} = UTM{North,Zone,Datum}
+
+"""
+    UTMSouth{zone}(x, y)
+    UTMSouth{zone,Datum}(x, y)
+
+UTM (Universal Transverse Mercator) South coordinates in length units (default to meter)
+with `zone` (1 ≤ zone ≤ 60) and a given `Datum` (default to `WGS84`).
+
+## Examples
+
+```julia
+UTMSouth{1}(1, 1) # add default units
+UTMSouth{1}(1u"m", 1u"m") # integers are converted converted to floats
+UTMSouth{1}(1.0u"km", 1.0u"km") # length quantities are converted to meters
+UTMSouth{1}(1.0u"m", 1.0u"m")
+UTMSouth{1,WGS84Latest}(1.0u"m", 1.0u"m")
+```
+"""
+const UTMSouth{Zone,Datum} = UTM{South,Zone,Datum}
+
+function formulas(C::Type{<:UTM{Hemisphere,Zone,Datum}}, ::Type{T}) where {Hemisphere,Zone,Datum,T}
+  a = numconvert(T, majoraxis(ellipsoid(Datum)))
+  xₒ = falseeasting(T, C) / a
+  yₒ = falsenorthing(T, C) / a
+
+  tmfx, tmfy = formulas(astm(C), T)
+
+  fx(λ, ϕ) = tmfx(λ, ϕ) + xₒ
+  fy(λ, ϕ) = tmfy(λ, ϕ) + yₒ
+
+  fx, fy
+end
+
+function Base.convert(C::Type{UTM{Hemisphere,Zone,Datum}}, coords::LatLon{Datum}) where {Hemisphere,Zone,Datum}
+  T = numtype(coords.lon)
+  xₒ = falseeasting(T, C)
+  yₒ = falsenorthing(T, C)
+
+  tm = convert(astm(C), coords)
+
+  C(tm.x + xₒ, tm.y + yₒ)
+end
+
+function Base.convert(::Type{LatLon{Datum}}, coords::C) where {Hemisphere,Zone,Datum,C<:UTM{Hemisphere,Zone,Datum}}
+  T = numtype(coords.x)
+  xₒ = falseeasting(T, C)
+  yₒ = falsenorthing(T, C)
+
+  tm = astm(C)(coords.x - xₒ, coords.y - yₒ)
+
+  convert(LatLon{Datum}, tm)
+end
+
+# -----------------
+# HELPER FUNCTIONS
+# -----------------
+
+function astm(::Type{<:UTM{Hemisphere,Zone,Datum}}) where {Hemisphere,Zone,Datum}
+  lonₒ = (6 * Zone - 183) * u"°"
+  TransverseMercator{0.9996,0.0u"°",lonₒ,Datum}
+end
+
+falseeasting(::Type{T}, ::Type{<:UTM}) where {T} = T(500000) * u"m"
+
+falsenorthing(::Type{T}, ::Type{<:UTM{North}}) where {T} = T(0) * u"m"
+falsenorthing(::Type{T}, ::Type{<:UTM{South}}) where {T} = T(10000000) * u"m"

test/crs.jl

@@ -1369,6 +1369,30 @@
       @inferred convert(LatLon{WGS84Latest}, c2)
     end
 
+    @testset "LatLon <> UTM" begin
+      c1 = LatLon(T(56), T(12))
+      c2 = convert(UTMNorth{32,WGS84Latest}, c1)
+      @test c2 ≈ UTMNorth{32}(T(687071.439107327), T(6210141.326872105))
+      c3 = convert(LatLon{WGS84Latest}, c2)
+      @test c3 ≈ c1
+
+      c1 = LatLon(-T(44), T(174))
+      c2 = convert(UTMSouth{59,WGS84Latest}, c1)
+      @test c2 ≈ UTMSouth{59}(T(740526.3210524899), T(5123750.873037999))
+      c3 = convert(LatLon{WGS84Latest}, c2)
+      @test c3 ≈ c1
+
+      # type stability
+      c1 = LatLon(T(56), T(12))
+      c2 = LatLon(-T(44), T(174))
+      c3 = UTMNorth{32}(T(687071.439107327), T(6210141.326872105))
+      c4 = UTMSouth{59}(T(740526.3210524899), T(5123750.873037999))
+      @inferred convert(UTMNorth{32,WGS84Latest}, c1)
+      @inferred convert(UTMSouth{59,WGS84Latest}, c2)
+      @inferred convert(LatLon{WGS84Latest}, c3)
+      @inferred convert(LatLon{WGS84Latest}, c4)
+    end
+
     @testset "Projection conversion" begin
       # same datum
       c1 = Lambert(T(10018754.171394622), T(4489858.8869480025))