the personal playground of evan louie; developer, designer, photographer, and breaker of the web.
   (  )   /\   _                 (     
    \ |  (  \ ( \.(               )                      _____
  \  \ \  `  `   ) \             (  ___                 / _   \
 (_`    \+   . x  ( .\            \/   \____-----------/ (o)   \_
- .-               \+  ;          (  O                           \____
                          )        \_____________  `              \  /
(__                +- .( -'.- <. - _  VVVVVVV VV V\                 \/
(_____            ._._: <_ - <- _  (--  _AAAAAAA__A_/                  |
  .    /./.+-  . .- /  +--  - .     \______________//_              \_______
  (__ ' /x  / x _/ (                                  \___'          \     /
 , x / ( '  . / .  /                                      |           \   /
    /  /  _/ /    +                                      /              \/
   '  (__/                                             /                  \

A Clojure(Script) UTM Easting/Northing to Latitude/Longitude Converter

Couldn't find a good one online.

A partial Clojure(Script) port of the UTM package.

Usage

(require '(com.evanlouie.geocoding.utm :refer [utm->latitude-longitude]))

(utm->latitude-longitude 491015.94 5459166.0 10 :zone-letter \\U)
; => {:latitude 49.2851798490208, :longitude -123.12353635192477}

(utm->latitude-longitude 491015.94 5459166.0 10 :zone-letter \"U\")
; => {:latitude 49.2851798490208, :longitude -123.12353635192477}

(utm->latitude-longitude 491015.94 5459166.0 10 :northern-hemisphere? true)
; => {:latitude 49.2851798490208, :longitude -123.12353635192477}

(try
  (utm->longitude-latitude 491015.94 5459166.0 10 :zone-letter \"A\" :validate-inputs? true)
  (catch Exception e
    (prn e)))
; #error {
;   :cause \"Error converint UTM coordinate to longitude/latitude\"
; ....
; }
; => nil

The Code

(ns com.evanlouie.geocoding.utm
  "This is a Clojure(Script) port of the UTM package by Tobias Bieniek
  <Tobias.Bieniek@gmx.de>.

  This package uses host `Math` functions that should operate the same across
  JVM and JS -- however this package has only been thoroughly tested in the
  JVM."
  (:require
   [clojure.spec.alpha :as s]
   [clojure.string]))


;;------------------------------------------------------------------------------
;; Constants


(def ^:private K0 0.9996)
(def ^:private E 0.00669438)
(def ^:private E2 (Math/pow E 2))
(def ^:private E3 (Math/pow E 3))
(def ^:private E_P2 (/ E (- 1 E)))
(def ^:private SQRT_E (Math/sqrt (- 1 E)))
(def ^:private _E (-> 1 (- SQRT_E) (/ (+ 1 SQRT_E))))
(def ^:private _E2 (Math/pow _E 2))
(def ^:private _E3 (Math/pow _E 3))
(def ^:private _E4 (Math/pow _E 4))
(def ^:private _E5 (Math/pow _E 5))
(def ^:private M1 (-> 1
                      (- (-> E (/ 4)))
                      (- (-> 3 (* E2) (/ 64)))
                      (- (-> 5 (* E3) (/ 256)))))
(def ^:private M2 (-> (-> 3 (* E) (/ 8))
                      (+ (-> 3 (* E2) (/ 32)))
                      (+ (-> 45 (* E3) (/ 1024)))))
(def ^:private M3 (-> (-> 15 (* E2) (/ 256))
                      (+ (-> 45 (* E3) (/ 1024)))))
(def ^:private M4 (-> 35 (* E3) (/ 3072)))
(def ^:private P2 (-> (-> 3 (/ 2) (* _E))
                      (- (-> 27 (/ 32) (* _E3)))
                      (+ (-> 269 (/ 512) (* _E5)))))
(def ^:private P3 (-> (-> 21 (/ 16) (* _E2))
                      (- (-> 55 (/ 32) (* _E4)))))
(def ^:private P4 (-> (-> 151 (/ 96) (* _E3))
                      (- (-> 417 (/ 128) (* _E5)))))
(def ^:private P5 (-> 1097 (/ 512) (* _E4)))
(def ^:private R 6378137)
(def ^:private zone-letters (map char "CDEFGHJKLMNPQRSTUVWX"))


;;------------------------------------------------------------------------------
;; Specs
;; NOTE fdefs are located with the functions themselves.


(s/def ::zone-number (s/and pos-int? #(>= % 1) #(<= % 60)))
(s/def ::zone-letter (s/or :string (s/and string? #(= 1 (count %)) #(some #{(first (map char %))} zone-letters))
                           :char (s/and char? #(some #{%} zone-letters))))
(s/def ::easting (s/and number? #(>= % 100000) #(< % 1000000))) ; max == 999999
(s/def ::northing (s/and number? #(>= % 0) #(<= % 10000000)))
(s/def ::latitude (s/and number? #(>= % -90.0) #(<= % 90.0)))
(s/def ::longitude (s/and number? #(>= % -180.0) #(<= % 180.0)))
(s/def ::northern-hemisphere? boolean?)
(s/def ::utm-coordinate (s/keys :req-un [::easting
                                         ::northing
                                         ::zone-number
                                         ::zone-letter]))
(s/def ::validate-inputs? (s/nilable #(some? %)))


;;------------------------------------------------------------------------------
;; Helpers


(defn- longitude->zone-number
  "Reference: https://www.youtube.com/watch?v=HwC5CcWvSFQ

  Example:
  (longitude->utm-zone-number 123.1207) => 10"
  [longitude]
  (when (not (s/valid? ::longitude longitude))
    (throw (ex-info "Invalid longitude"
                    {:failed-spec (s/explain-data ::longitude longitude)})))
  (int (+ 31 (/ longitude 6))))

(defn- latitude->zone-letter
  [latitude]
  (let [zone-letters (clojure.string/split "CDEFGHJKLMNPQRSTUVWXX" #"")
        idx          (int (Math/floor (/ (+ latitude 80) 8)))]
    (get zone-letters idx)))

(defn- zone-number->central-longitude
  [zone-number]
  (-> (dec zone-number)
      (* 6)
      (- 180)
      (+ 3)))

(defn- radians->degrees
  [radians]
  (-> radians
      (/ (Math/PI))
      (* 180)))

(defn- degrees->radians
  [degrees]
  (-> degrees
      (* (Math/PI))
      (/ 180)))

(defn- zone-letter->northern-hemisphere?
  "Returns a `boolean` denoting if the provided `zone-letter` is in the northern
  hemisphere."
  [zone-letter]
  (let [conformed (s/conform ::zone-letter zone-letter)]
    (when (s/invalid? conformed)
      (throw (ex-info "Failed to computer if coordinate in northern-hemisphere based on zone-letter"
                      {:zone-letter zone-letter
                       :cause       (s/explain-data ::zone-letter zone-letter)})))
    (case (first conformed)
      :string (>= (byte (first (map char zone-letter)))
                  (byte \N))
      :char (>= (byte zone-letter) (byte \N))
      (throw (ex-info "Invalid zone-letter"
                      {:zone-letter zone-letter})))))

(defn- latitudinal-identifier->northern-hemisphere?
  [identifier]
  (let [conformed (s/conform ::latitudinal-identifier identifier)]
    (when (s/invalid? conformed)
      (throw (ex-info "Invalid :latitudinal-identifier"
                      (s/explain-data ::latitudinal-identifier identifier))))
    (case (first conformed)
      ;; conformed will [:zone-letter [:char \X]] or [:zone-letter [:string "X"]]
      :zone-letter          (zone-letter->northern-hemisphere? (second (second conformed)))
      ;; conformed will be [:northern-hemisphere true]]
      :northern-hemisphere? (second conformed)
      (throw (ex-info "Invalid :latitudinal-identifier"
                      {:latitudinal-identifier identifier})))))


;;------------------------------------------------------------------------------
;; Public


(s/fdef utm->longitude-latitude
  :args (s/cat :easting ::easting
               :northing ::northing
               :zone-number ::zone-number
               :options (s/keys* :req-un [(or ::northern-hemisphere?
                                              ::zone-letter)]
                                 :opt-un [::validate-inputs?]))
  :ret (s/keys :req-un [::latitude ::longitude]))

(defn utm->longitude-latitude
  "Convert the provided `easting`, `northing`, `zone-number` to a map containing
  `:latitude`and `:longitude`.

  As a zone-letter is not techincally required to do UTM to latitude/longitude
  conversion, an option of EITHER `:zone-letter` or `:northern-hemisphere?` is
  REQUIRED:

  - `:zone-letter` is either a `char` or a `string` of length 1 of the
    UTM zone-letter where the `easting`/`northing` is located.
  - `:northern-hemipshere?` is a `boolean` saying if the `easting`/`northing` is
    in the northern hemisphere.

  When `:validate-inputs?` is non-nil, the passed `easting`, `northing`, and
  `:zone-letter` will be validated prior conversion.
  Unexpected behaviour can occur if any properties of the utm-coordinate are
  not valid (e.g. an `easting` which is greater than 999999).
  This will have a minor performance hit when run against large datasets.
  Defaults to `false`.

  Examples:

  user> (utm->latitude-longitude 491015.94 5459166.0 10 :zone-letter \\U)
  => {:latitude 49.2851798490208, :longitude -123.12353635192477}

  user> (utm->latitude-longitude 491015.94 5459166.0 10 :zone-letter \"U\")
  => {:latitude 49.2851798490208, :longitude -123.12353635192477}

  user> (utm->latitude-longitude 491015.94 5459166.0 10 :northern-hemisphere? true)
  => {:latitude 49.2851798490208, :longitude -123.12353635192477}

  user> (try
          (utm->longitude-latitude 491015.94 5459166.0 10 :zone-letter \"A\" :validate-inputs? true)
          (catch Exception e
            (prn e)))
  #error {
    :cause \"Error converint UTM coordinate to longitude/latitude\"
  ....
  }
  => nil"
  [easting northing zone-number & {:keys [northern-hemisphere? zone-letter validate-inputs?]}]
  (when validate-inputs?
    (doseq [[spec value optional?] [[::easting easting]
                                    [::northing northing]
                                    [::zone-number zone-number]
                                    [::zone-letter zone-letter true]]]
      (when (not (s/valid? (if optional? (s/nilable spec) spec)
                           value))
        (throw (ex-info "Error converting UTM coordinate to longitude/latitude"
                        {:reason (str "Invalid " (name spec))
                         spec value
                         :cause (s/explain-data spec value)})))))

  (let [north?      (or northern-hemisphere?
                        (zone-letter->northern-hemisphere? zone-letter))
        x           (- easting 500000)
        y           (if north? northing (- northing 1e7))
        m           (-> y (/ K0))
        mu          (-> m (/ (-> R (* M1))))
        p-rad       (+ mu
                       (* P2 (Math/sin (* 2 mu)))
                       (* P3 (Math/sin (* 4 mu)))
                       (* P4 (Math/sin (* 6 mu)))
                       (* P5 (Math/sin (* 8 mu))))
        p-sin       (Math/sin p-rad)
        p-sin2      (Math/pow p-sin 2)
        p-cos       (Math/cos p-rad)
        p-tan       (Math/tan p-rad)
        p-tan2      (Math/pow p-tan 2)
        p-tan4      (Math/pow p-tan 4)
        ep-sin      (- 1 (* E p-sin2))
        ep-sin-sqrt (Math/sqrt ep-sin)
        n           (/ R ep-sin-sqrt)
        r           (/ (- 1 E) ep-sin)
        c           (* _E p-cos p-cos)
        c2          (* c c)
        d           (/ x (* n K0))
        d2          (Math/pow d 2)
        d3          (Math/pow d 3)
        d4          (Math/pow d 4)
        d5          (Math/pow d 5)
        d6          (Math/pow d 6)

        lat-rad (-> p-rad
                    (- (-> p-tan
                           (/ r)
                           (* (-> d2
                                  (/ 2)
                                  (- (-> d4
                                         (/ 24)
                                         (* (-> 5
                                                (+ (* 3 p-tan2))
                                                (+ (* 10 c))
                                                (- (* 4 c2))
                                                (- (* 9 E_P2))))))))))
                    (+ (-> d6
                           (/ 720)
                           (* (-> 61
                                  (+ (* 90 p-tan2))
                                  (+ (* 298 c))
                                  (+ (* 45 p-tan4))
                                  (- (* 252 E_P2))
                                  (- (* 3 c2)))))))

        long-rad (-> d
                     (- (-> d3
                            (/ 6)
                            (* (-> 1
                                   (+ (* 2 p-tan2))
                                   (+ c)))))
                     (+ (-> d5
                            (/ 120)
                            (* (-> 5
                                   (- (* 2 c))
                                   (+ (* 28 p-tan2))
                                   (- (* 3 c2))
                                   (+ (* 8 E_P2))
                                   (+ (* 24 p-tan4))))))
                     (/ p-cos))

        lat-deg  (radians->degrees lat-rad)
        long-deg (-> (radians->degrees long-rad)
                     (+ (zone-number->central-longitude zone-number)))]

    {:latitude  lat-deg
     :longitude long-deg}))

(s/fdef longitude-latitude->utm
  :args (s/cat :latitude ::latitude
               :longitude ::longitude
               :options (s/keys* :opt-un [::zone-number
                                          ::validate-inputs?]))
  :ret ::utm-coordinate)

(defn longitude-latitude->utm
  "Converts the provided `longitude`, `latitude` and optional `zone-number` to
  a `::utm-coordinate`.

  If `:zone-number` is not specified, it will be  computed from the provided
  `longitude`.

  If `:validate-inputs?` is non-nil, `longitude`, `latitude`, and `zone-number`
  (if provided) will be validated prior to converstion.
  Unexpected behaviour can occur if any input values are invalid (e.g. a
  `latitude` greater than 90).
  This may have a minor performance hit when run against large datasets.
  Defaults to `false`"
  [latitude longitude & {:keys [zone-number validate-inputs?]}]
  (when validate-inputs?
    (doseq [[spec value optional?] [[::longitude longitude
                                     ::latitude latitude
                                     ::zone-number zone-number true]]]
      (when (not (s/valid? (if optional? (s/nilable spec) spec)
                           value))
        (throw (ex-info "Error converting longitude/latitude to UTM coordinate"
                        {:reason (str "Invalid " (name spec))
                         spec    value
                         :cause  (s/explain-data spec value)})))))

  (let [lat-rad (degrees->radians latitude)
        lat-sin (Math/sin lat-rad)
        lat-cos (Math/cos lat-rad)

        lat-tan  (Math/tan lat-rad)
        lat-tan2 (Math/pow lat-tan 2)
        lat-tan4 (Math/pow lat-tan 4)

        zone-num    (if (nil? zone-number)
                      (longitude->zone-number longitude)
                      zone-number)
        zone-letter (latitude->zone-letter latitude)

        lon-rad         (degrees->radians longitude)
        central-lon     (zone-number->central-longitude zone-num)
        central-lon-rad (degrees->radians central-lon)

        n (-> R (/ (Math/sqrt (-> 1 (- (* E lat-sin lat-sin))))))
        c (* E_P2 lat-cos lat-cos)

        a  (-> lon-rad (- central-lon-rad) (* lat-cos))
        a2 (Math/pow a 2)
        a3 (Math/pow a 3)
        a4 (Math/pow a 4)
        a5 (Math/pow a 5)
        a6 (Math/pow a 6)

        m (-> R
              (* (-> M1
                     (* lat-rad)
                     (- (-> M2 (* (Math/sin (* 2 lat-rad)))))
                     (+ (-> M3 (* (Math/sin (* 4 lat-rad)))))
                     (- (-> M4 (* (Math/sin (* 6 lat-rad))))))))

        easting (-> K0
                    (* n)
                    (* (-> a
                           (+ (-> a3
                                  (/ 6)
                                  (* (-> 1
                                         (- lat-tan2)
                                         (+ c)))))
                           (+ (-> a5
                                  (/ 120)
                                  (* (-> 5
                                         (- (* 18 lat-tan2))
                                         (+ lat-tan4)
                                         (+ (* 72 c))
                                         (- (* 58 E_P2))))))))
                    (+ 500000))

        northing (-> K0
                     (* (-> m
                            (+ (-> n
                                   (* lat-tan)
                                   (* (-> a2
                                          (/ 2)
                                          (+ (-> a4
                                                 (/ 24)
                                                 (* (-> 5
                                                        (- lat-tan2)
                                                        (+ (* 9 c))
                                                        (+ (* 4 c c))))))
                                          (+ (-> a6
                                                 (/ 720)
                                                 (* (-> 61
                                                        (- (* 58 lat-tan2))
                                                        (+ lat-tan4)
                                                        (+ (* 600 c))
                                                        (- (* 330 E_P2)))))))))))))

        northing (if (neg? latitude)
                   (+ northing 1e7)
                   northing)]

    {:easting     easting
     :northing    northing
     :zone-number zone-num
     :zone-letter zone-letter}))