diff --git a/src/Cartography.jl b/src/Cartography.jl index 9ef52fb..407e63a 100644 --- a/src/Cartography.jl +++ b/src/Cartography.jl @@ -68,9 +68,16 @@ export Robinson, OrthoNorth, OrthoSouth, + UTM, + UTMNorth, + UTMSouth, datum, indomain, + # hemispheres + North, + South, + # codes EPSG, ESRI diff --git a/src/crs/projected.jl b/src/crs/projected.jl index 61a197d..8a7e04b 100644 --- a/src/crs/projected.jl +++ b/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 diff --git a/src/crs/projected/utm.jl b/src/crs/projected/utm.jl new file mode 100644 index 0000000..04f8f8c --- /dev/null +++ b/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" diff --git a/test/crs.jl b/test/crs.jl index 861aefa..a187908 100644 --- a/test/crs.jl +++ b/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))