From 63446b1511b0a507121be3ae3c94514288b7fc8d Mon Sep 17 00:00:00 2001 From: DiracM <32461143+DiracM@users.noreply.github.com> Date: Thu, 16 Mar 2023 19:21:24 +0000 Subject: [PATCH] Added latitude band for the UTMZ point type --- README.md | 15 ++++++----- src/Geodesy.jl | 2 +- src/conversion.jl | 4 +-- src/points.jl | 11 ++++---- src/transformations.jl | 19 ++++++++------ src/utm.jl | 57 +++++++++++++++++++++++++++++++++++++++++ test/conversion.jl | 8 +++--- test/points.jl | 6 ++--- test/transformations.jl | 12 ++++----- 9 files changed, 98 insertions(+), 36 deletions(-) diff --git a/README.md b/README.md index 47f7d2e..0da58eb 100644 --- a/README.md +++ b/README.md @@ -52,9 +52,9 @@ point_enu = ENU(point_enu, point_origin, wgs84) Similarly, we could convert to UTM/UPS coordinates, and two types are provided for this - `UTM` stores 3D coordinates `x`, `y`, and `z` in an unspecified zone, -while `UTMZ` includes the `zone` number and `hemisphere` bool (where `true` = -northern, `false` = southern). To get the canonical zone for your coordinates, -simply use: +while `UTMZ` includes the `zone` number, `hemisphere` bool (where `true` = +northern, `false` = southern) and `latletter` char with the latitude band letter. +To get the canonical zone for your coordinates, simply use: ```julia x_utmz = UTMZ(x_lla, wgs84) ``` @@ -317,11 +317,12 @@ universal polar-stereographic (UPS) coordinates (where the zone is `0`). ##### `UTMZ{T}` - universal transverse-Mercator + zone In addition to the easting `x`, northing `y` and height `z`, the global `UTMZ` type -also encodes the UTM `zone` and `hemisphere`, where `zone` is a `UInt8` and -`hemisphere` is a `Bool` for compact storage. The northern hemisphere is -denoted as `true`, and the southern as `false`. Zone `0` corresponds to the UPS +also encodes the UTM `zone` , `hemisphere` and `latletter`, where `zone` is a `UInt8`, +`hemisphere` is a `Bool` for compact storage and `latletter` a char. The northern hemisphere +is denoted as `true`, and the southern as `false`. Zone `0` corresponds to the UPS projection about the corresponding pole, otherwise `zone` is an integer between -`1` and `60`. +`1` and `60`. The latitude band letter`latletter` char with the latitude band is a char +that can go from `C` to `X`. ##### `ENU{T}` - east-north-up diff --git a/src/Geodesy.jl b/src/Geodesy.jl index e2b73c9..a4cdf82 100644 --- a/src/Geodesy.jl +++ b/src/Geodesy.jl @@ -49,7 +49,7 @@ export ITRF, GDA94, # UTM helpers - utm_zone + utm_zone, utm_lat_letter include("ellipsoids.jl") include("datums.jl") diff --git a/src/conversion.jl b/src/conversion.jl index bf8bc93..86a4f90 100644 --- a/src/conversion.jl +++ b/src/conversion.jl @@ -86,5 +86,5 @@ UTM(enu::ENU, zone::Integer, hemisphere::Bool, origin, datum) = UTMfromENU(origi ### UTMZ <-> UTM coordinates ### ############################### -UTMZ(utm::UTM, zone::Integer, hemisphere::Bool, datum) = UTMZfromUTM(zone, hemisphere, datum)(utm) -UTM(utmz::UTMZ, zone::Integer, hemisphere::Bool, datum) = UTMfromUTMZ(zone, hemisphere, datum)(utmz) +UTMZ(utm::UTM, zone::Integer, hemisphere::Bool, latletter::Char, datum) = UTMZfromUTM(zone, hemisphere, latletter, datum)(utm) +UTM(utmz::UTMZ, zone::Integer, hemisphere::Bool, latletter::Char, datum) = UTMfromUTMZ(zone, hemisphere, latletter, datum)(utmz) diff --git a/src/points.jl b/src/points.jl index 67bc303..1245b94 100644 --- a/src/points.jl +++ b/src/points.jl @@ -118,10 +118,11 @@ struct UTMZ{T <: Number} z::T zone::UInt8 isnorth::Bool # which hemisphere + latletter::Char # which latitude band end -UTMZ(x::T, y::T, zone::Integer, isnorth::Bool) where {T} = UTMZ(x, y, zero(T), UInt8(zone), isnorth) -UTMZ(x, y, z, zone::Integer, isnorth::Bool) = UTMZ(promote(x, y, z)..., UInt8(zone), isnorth) -UTMZ(utm::UTM, zone::Integer, isnorth::Bool) = UTMZ(utm.x, utm.y, utm.z, UInt8(zone), isnorth) +UTMZ(x::T, y::T, zone::Integer, isnorth::Bool, latletter::Char) where {T} = UTMZ(x, y, zero(T), UInt8(zone), isnorth, latletter) +UTMZ(x, y, z, zone::Integer, isnorth::Bool, latletter::Char) = UTMZ(promote(x, y, z)..., UInt8(zone), isnorth, latletter) +UTMZ(utm::UTM, zone::Integer, isnorth::Bool, latletter::Char) = UTMZ(utm.x, utm.y, utm.z, UInt8(zone), isnorth, latletter) UTM(utmz::UTMZ) = UTM(utmz.x, utmz.y, utmz.z) -Base.isapprox(utm1::UTMZ, utm2::UTMZ; atol = 1e-6, kwargs...) = (utm1.zone == utm2.zone) & (utm1.isnorth == utm2.isnorth) & isapprox(utm1.x, utm2.x; atol = atol, kwargs...) & isapprox(utm1.y, utm2.y; atol = atol, kwargs...) & isapprox(utm1.z, utm2.z; atol = atol, kwargs...) # atol in metres (1μm) -Base.show(io::IO, utm::UTMZ) = print(io, "UTMZ($(utm.x), $(utm.y), $(utm.z), zone=$(utm.zone == 0 ? "polar" : Int(utm.zone)) ($(utm.isnorth ? "north" : "south")))") +Base.isapprox(utm1::UTMZ, utm2::UTMZ; atol = 1e-6, kwargs...) = (utm1.zone == utm2.zone) & (utm1.isnorth == utm2.isnorth) & (utm1.latletter == utm2.latletter) & isapprox(utm1.x, utm2.x; atol = atol, kwargs...) & isapprox(utm1.y, utm2.y; atol = atol, kwargs...) & isapprox(utm1.z, utm2.z; atol = atol, kwargs...) # atol in metres (1μm) +Base.show(io::IO, utm::UTMZ) = print(io, "UTMZ($(utm.x), $(utm.y), $(utm.z), zone=$(utm.zone == 0 ? "polar" : Int(utm.zone)) $(utm.latletter) ($(utm.isnorth ? "north" : "south")))") diff --git a/src/transformations.jl b/src/transformations.jl index 5e7e428..b1a0f73 100644 --- a/src/transformations.jl +++ b/src/transformations.jl @@ -476,6 +476,7 @@ Base.show(io::IO, trans::UTMZfromLLA) = print(io, "UTMZfromLLA($(trans.datum))") function (trans::UTMZfromLLA)(lla::LLA) (zone, isnorth) = utm_zone(lla) + latletter = utm_lat_letter(lla.lat) if zone == 0 # Do polar steriographic projection k0 = 0.994 @@ -483,7 +484,7 @@ function (trans::UTMZfromLLA)(lla::LLA) x = x + 2e6 y = y + 2e6 - return UTMZ(x, y, lla.alt, 0, isnorth) + return UTMZ(x, y, lla.alt, 0, isnorth, latletter) else lon_ref = Float64(utm_meridian(zone)) k0 = 0.9996 # Horizontal scaling factor @@ -492,7 +493,7 @@ function (trans::UTMZfromLLA)(lla::LLA) x = 5e5 + x # Convention has 500km offset for easting y = (isnorth ? 0.0 : 1e7) + y # Northing offset for southern hemisphere # also, k = k * k0 - return UTMZ(x, y, lla.alt, zone, isnorth) # Note: scaling not applied to vertical dimension + return UTMZ(x, y, lla.alt, zone, isnorth, latletter) # Note: scaling not applied to vertical dimension end end @@ -511,9 +512,10 @@ Transformation to append the UTM/UPS zone and hemisphere to a raw `UTM` data poi struct UTMZfromUTM{Datum} <: Transformation zone::Int isnorth::Bool + latletter::Char datum::Datum end -UTMZfromUTM(zone::Integer, h, d::D) where {D} = UTMZfromUTM{D}(UInt8(zone), h, d) +UTMZfromUTM(zone::Integer, h, l, d::D) where {D} = UTMZfromUTM{D}(UInt8(zone), h, l, d) Base.show(io::IO, trans::UTMZfromUTM) = print(io, "UTMZfromUTM(zone=$(trans.zone == 0 ? "polar" : trans.zone) ($(trans.isnorth ? "north" : "south")), $(trans.datum))") #Base.show(io::IO, trans::UTMZfromUTM) = print(io, "UTMZfromUTM(zone=$(trans.zone == 0 ? "polar" : trans.zone) ($(trans.isnorth ? "north" : "south")))") @@ -526,15 +528,16 @@ automatically convert the data to the specified zone if necessary. struct UTMfromUTMZ{Datum} <: Transformation zone::Int isnorth::Bool + latletter::Char datum::Datum end -UTMfromUTMZ(zone::Integer, h, d::D) where {D} = UTMfromUTMZ{D}(UInt8(zone), h, d) +UTMfromUTMZ(zone::Integer, h, l, d::D) where {D} = UTMfromUTMZ{D}(UInt8(zone), h, l, d) Base.show(io::IO, trans::UTMfromUTMZ) = print(io, "UTMfromUTMZ(zone=$(trans.zone == 0 ? "polar" : trans.zone) ($(trans.isnorth ? "north" : "south")), $(trans.datum))") #Base.show(io::IO, trans::UTMfromUTMZ) = print(io, "UTMfromUTMZ(zone=$(trans.zone == 0 ? "polar" : trans.zone) ($(trans.isnorth ? "north" : "south")))") -(trans::UTMZfromUTM)(utm::UTM) = UTMZ(utm.x, utm.y, utm.z, trans.zone, trans.isnorth) +(trans::UTMZfromUTM)(utm::UTM) = UTMZ(utm.x, utm.y, utm.z, trans.zone, trans.isnorth, trans.latletter) function (trans::UTMfromUTMZ)(utm::UTMZ) - if trans.zone == utm.zone && trans.isnorth == utm.isnorth + if trans.zone == utm.zone && trans.isnorth == utm.isnorth && trans.latletter == utm.latletter UTM(utm.x, utm.y, utm.z) else # Should this be an error or an automatic transformation to the correct zone? @@ -543,8 +546,8 @@ function (trans::UTMfromUTMZ)(utm::UTMZ) end end -Base.inv(trans::UTMfromUTMZ) = UTMZfromUTM(trans.zone, trans.isnorth, trans.datum) -Base.inv(trans::UTMZfromUTM) = UTMfromUTMZ(trans.zone, trans.isnorth, trans.datum) +Base.inv(trans::UTMfromUTMZ) = UTMZfromUTM(trans.zone, trans.isnorth, trans.latletter, trans.datum) +Base.inv(trans::UTMZfromUTM) = UTMfromUTMZ(trans.zone, trans.isnorth, trans.latletter, trans.datum) ################### ## ECEF <-> UTMZ ## diff --git a/src/utm.jl b/src/utm.jl index 2df9359..9bf1549 100644 --- a/src/utm.jl +++ b/src/utm.jl @@ -54,3 +54,60 @@ function utm_meridian(zone::Integer) end return 6 * zone - 183 end + +""" + utm_lat_letter(lat) + +Calculates the latitude bands from the latitude poistion. +This is part of the military grid reference system (MGRS). +A letter coming before "N" in the alphabet is in the southern +hemisphere and any letter "N" or after is in the northern hemisphere. +""" +function utm_lat_letter(lat::T) where T + + lat_letter::Char = ' ' + if lat <-72 + lat_letter='C' + elseif lat <-64 + lat_letter='D' + elseif lat <-56 + lat_letter='E' + elseif lat <-48 + lat_letter='F' + elseif lat <-40 + lat_letter='G' + elseif lat <-32 + lat_letter='H' + elseif lat <-24 + lat_letter='J' + elseif lat <-16 + lat_letter='K' + elseif lat <-8 + lat_letter='L' + elseif lat <0 + lat_letter='M' + elseif lat <8 + lat_letter='N' + elseif lat <16 + lat_letter='P' + elseif lat <24 + lat_letter='Q' + elseif lat <32 + lat_letter='R' + elseif lat <40 + lat_letter='S' + elseif lat <48 + lat_letter='T' + elseif lat <56 + lat_letter='U' + elseif lat <64 + lat_letter='V' + elseif lat <72 + lat_letter='W' + else + lat_letter='X' + end + + return lat_letter + +end diff --git a/test/conversion.jl b/test/conversion.jl index f025ceb..15042c6 100644 --- a/test/conversion.jl +++ b/test/conversion.jl @@ -22,7 +22,7 @@ @test ECEF(enu, lla_ref, wgs84) ≈ ecef # LLA <-> UTM - (z, h) = (19, true) + (z, h, l) = (19, true, 'T') utm = UTM(327412.48528248386, 4.692686244318043e6, 0.0) @test UTM(lla, z, h, wgs84) ≈ utm @test LLA(utm, z, h, wgs84) ≈ lla @@ -36,7 +36,7 @@ @test ENU(utm, z, h, lla_ref, wgs84) ≈ enu # LLA <-> UTMZ - utmz = UTMZ(327412.48528248386, 4.692686244318043e6, 0, z, h) + utmz = UTMZ(327412.48528248386, 4.692686244318043e6, 0, z, h, l) @test UTMZ(lla, wgs84) ≈ utmz @test LLA(utmz, wgs84) ≈ lla @@ -49,8 +49,8 @@ @test ENU(utmz, lla_ref, wgs84) ≈ enu # UTM <-> UTMZ - @test UTMZ(utm, z, h, wgs84) ≈ utmz - @test UTM(utmz, z, h, wgs84) ≈ utm + @test UTMZ(utm, z, h, l, wgs84) ≈ utmz + @test UTM(utmz, z, h, l, wgs84) ≈ utm # Identity @test ECEF(ecef, wgs84) == ecef diff --git a/test/points.jl b/test/points.jl index 7defa64..86d40aa 100644 --- a/test/points.jl +++ b/test/points.jl @@ -19,9 +19,9 @@ utm = UTM(10., 10., 0.) @test UTM(10., 10.) == utm - utmz = UTMZ(10., 10., 0., 1, true) - @test UTMZ(10., 10., 1, true) == utmz - @test UTMZ(utm, 1, true) == utmz + utmz = UTMZ(10., 10., 0., 1, true, 'D') + @test UTMZ(10., 10., 1, true, 'D') == utmz + @test UTMZ(utm, 1, true, 'D') == utmz @test UTM(utmz) == utm end # @testset diff --git a/test/transformations.jl b/test/transformations.jl index 2654a80..f0db90b 100644 --- a/test/transformations.jl +++ b/test/transformations.jl @@ -30,19 +30,19 @@ utmz = utmz_lla(lla) utmz2 = utmz_lla(lla2) # Test also the poles go to UPS - @test utmz_lla(LLA(89.0, 89.0, 89.0)) ≈ UTMZ(2.111009610242531e6, 1.9980623200455606e6, 89.0, 0, true) - @test lla_utmz(UTMZ(2.111009610242531e6, 1.9980623200455606e6, 89.0, 0, true)) ≈ LLA(89.0, 89.0, 89.0) + @test utmz_lla(LLA(89.0, 89.0, 89.0)) ≈ UTMZ(2.111009610242531e6, 1.9980623200455606e6, 89.0, 0, true, 'X') + @test lla_utmz(UTMZ(2.111009610242531e6, 1.9980623200455606e6, 89.0, 0, true, 'X')) ≈ LLA(89.0, 89.0, 89.0) - @test utmz ≈ UTMZ(239579.67583179142,7.342551466042985e6,1374.7804632852078, UInt8(47), false) + @test utmz ≈ UTMZ(239579.67583179142,7.342551466042985e6,1374.7804632852078, UInt8(47), false, 'J') @test lla ≈ lla_utmz(utmz) - utmz_utm = UTMZfromUTM(47, false, wgs84) - utm_utmz = UTMfromUTMZ(47, false, wgs84) + utmz_utm = UTMZfromUTM(47, false, 'J', wgs84) + utm_utmz = UTMfromUTMZ(47, false, 'J', wgs84) @test utmz == utmz_utm(utm) @test utm == utm_utmz(utmz) # Sometimes this has to do something non-trivial - make sure it does - utm_utmz2 = UTMfromUTMZ(46, false, wgs84) + utm_utmz2 = UTMfromUTMZ(46, false,'X', wgs84) @test utm_utmz2(utmz) ≈ UTM(850009.7418418773, 7.340641097689279e6, 1374.7804632852078) # ENU coordinates