updates to layout core module packages in a more logical way, including renaming methods from objects to nodes and nodes.objid to nodes.id

daemon/core/location/corelocation.py

@@ -0,0 +1,273 @@
+"""
+location.py: definition of CoreLocation class that is a member of the
+Session object. Provides conversions between Cartesian and geographic coordinate
+systems. Depends on utm contributed module, from
+https://pypi.python.org/pypi/utm (version 0.3.0).
+"""
+
+import logging
+
+from core.emulator.enumerations import RegisterTlvs
+from core.location import utm
+
+
+class CoreLocation(object):
+    """
+    Member of session class for handling global location data. This keeps
+    track of a latitude/longitude/altitude reference point and scale in
+    order to convert between X,Y and geo coordinates.
+    """
+    name = "location"
+    config_type = RegisterTlvs.UTILITY.value
+
+    def __init__(self):
+        """
+        Creates a MobilityManager instance.
+
+        :return: nothing
+        """
+        # ConfigurableManager.__init__(self)
+        self.reset()
+        self.zonemap = {}
+        self.refxyz = (0.0, 0.0, 0.0)
+        self.refscale = 1.0
+        self.zoneshifts = {}
+        self.refgeo = (0.0, 0.0, 0.0)
+        for n, l in utm.ZONE_LETTERS:
+            self.zonemap[l] = n
+
+    def reset(self):
+        """
+        Reset to initial state.
+        """
+        # (x, y, z) coordinates of the point given by self.refgeo
+        self.refxyz = (0.0, 0.0, 0.0)
+        # decimal latitude, longitude, and altitude at the point (x, y, z)
+        self.setrefgeo(0.0, 0.0, 0.0)
+        # 100 pixels equals this many meters
+        self.refscale = 1.0
+        # cached distance to refpt in other zones
+        self.zoneshifts = {}
+
+    def px2m(self, val):
+        """
+        Convert the specified value in pixels to meters using the
+        configured scale. The scale is given as s, where
+        100 pixels = s meters.
+
+        :param val: value to use in converting to meters
+        :return: value converted to meters
+        """
+        return (val / 100.0) * self.refscale
+
+    def m2px(self, val):
+        """
+        Convert the specified value in meters to pixels using the
+        configured scale. The scale is given as s, where
+        100 pixels = s meters.
+
+        :param val: value to convert to pixels
+        :return: value converted to pixels
+        """
+        if self.refscale == 0.0:
+            return 0.0
+        return 100.0 * (val / self.refscale)
+
+    def setrefgeo(self, lat, lon, alt):
+        """
+        Record the geographical reference point decimal (lat, lon, alt)
+        and convert and store its UTM equivalent for later use.
+
+        :param lat: latitude
+        :param lon: longitude
+        :param alt: altitude
+        :return: nothing
+        """
+        self.refgeo = (lat, lon, alt)
+        # easting, northing, zone
+        e, n, zonen, zonel = utm.from_latlon(lat, lon)
+        self.refutm = ((zonen, zonel), e, n, alt)
+
+    def getgeo(self, x, y, z):
+        """
+        Given (x, y, z) Cartesian coordinates, convert them to latitude,
+        longitude, and altitude based on the configured reference point
+        and scale.
+
+        :param x: x value
+        :param y: y value
+        :param z: z value
+        :return: lat, lon, alt values for provided coordinates
+        :rtype: tuple
+        """
+        # shift (x,y,z) over to reference point (x,y,z)
+        x -= self.refxyz[0]
+        y = -(y - self.refxyz[1])
+        if z is None:
+            z = self.refxyz[2]
+        else:
+            z -= self.refxyz[2]
+        # use UTM coordinates since unit is meters
+        zone = self.refutm[0]
+        if zone == "":
+            raise ValueError("reference point not configured")
+        e = self.refutm[1] + self.px2m(x)
+        n = self.refutm[2] + self.px2m(y)
+        alt = self.refutm[3] + self.px2m(z)
+        (e, n, zone) = self.getutmzoneshift(e, n)
+        try:
+            lat, lon = utm.to_latlon(e, n, zone[0], zone[1])
+        except utm.OutOfRangeError:
+            logging.exception("UTM out of range error for n=%s zone=%s xyz=(%s,%s,%s)", n, zone, x, y, z)
+            lat, lon = self.refgeo[:2]
+        # self.info("getgeo(%s,%s,%s) e=%s n=%s zone=%s  lat,lon,alt=" \
+        #          "%.3f,%.3f,%.3f" % (x, y, z, e, n, zone, lat, lon, alt))
+        return lat, lon, alt
+
+    def getxyz(self, lat, lon, alt):
+        """
+        Given latitude, longitude, and altitude location data, convert them
+        to (x, y, z) Cartesian coordinates based on the configured
+        reference point and scale. Lat/lon is converted to UTM meter
+        coordinates, UTM zones are accounted for, and the scale turns
+        meters to pixels.
+
+        :param lat: latitude
+        :param lon: longitude
+        :param alt: altitude
+        :return: converted x, y, z coordinates
+        :rtype: tuple
+        """
+        # convert lat/lon to UTM coordinates in meters
+        e, n, zonen, zonel = utm.from_latlon(lat, lon)
+        _rlat, _rlon, ralt = self.refgeo
+        xshift = self.geteastingshift(zonen, zonel)
+        if xshift is None:
+            xm = e - self.refutm[1]
+        else:
+            xm = e + xshift
+        yshift = self.getnorthingshift(zonen, zonel)
+        if yshift is None:
+            ym = n - self.refutm[2]
+        else:
+            ym = n + yshift
+        zm = alt - ralt
+
+        # shift (x,y,z) over to reference point (x,y,z)
+        x = self.m2px(xm) + self.refxyz[0]
+        y = -(self.m2px(ym) + self.refxyz[1])
+        z = self.m2px(zm) + self.refxyz[2]
+        return x, y, z
+
+    def geteastingshift(self, zonen, zonel):
+        """
+        If the lat, lon coordinates being converted are located in a
+        different UTM zone than the canvas reference point, the UTM meters
+        may need to be shifted.
+        This picks a reference point in the same longitudinal band
+        (UTM zone number) as the provided zone, to calculate the shift in
+        meters for the x coordinate.
+
+        :param zonen: zonen
+        :param zonel: zone1
+        :return: the x shift value
+        """
+        rzonen = int(self.refutm[0][0])
+        # same zone number, no x shift required
+        if zonen == rzonen:
+            return None
+        z = (zonen, zonel)
+        # x shift already calculated, cached
+        if z in self.zoneshifts and self.zoneshifts[z][0] is not None:
+            return self.zoneshifts[z][0]
+
+        rlat, rlon, _ralt = self.refgeo
+        # ea. zone is 6deg band
+        lon2 = rlon + 6 * (zonen - rzonen)
+        # ignore northing
+        e2, _n2, _zonen2, _zonel2 = utm.from_latlon(rlat, lon2)
+        # NOTE: great circle distance used here, not reference ellipsoid!
+        xshift = utm.haversine(rlon, rlat, lon2, rlat) - e2
+        # cache the return value
+        yshift = None
+        if z in self.zoneshifts:
+            yshift = self.zoneshifts[z][1]
+        self.zoneshifts[z] = (xshift, yshift)
+        return xshift
+
+    def getnorthingshift(self, zonen, zonel):
+        """
+        If the lat, lon coordinates being converted are located in a
+        different UTM zone than the canvas reference point, the UTM meters
+        may need to be shifted.
+        This picks a reference point in the same latitude band (UTM zone letter)
+        as the provided zone, to calculate the shift in meters for the
+        y coordinate.
+
+        :param zonen: zonen
+        :param zonel:  zone1
+        :return: calculated y shift
+        """
+        rzonel = self.refutm[0][1]
+        # same zone letter, no y shift required
+        if zonel == rzonel:
+            return None
+        z = (zonen, zonel)
+        # y shift already calculated, cached
+        if z in self.zoneshifts and self.zoneshifts[z][1] is not None:
+            return self.zoneshifts[z][1]
+
+        rlat, rlon, _ralt = self.refgeo
+        # zonemap is used to calculate degrees difference between zone letters
+        latshift = self.zonemap[zonel] - self.zonemap[rzonel]
+        # ea. latitude band is 8deg high
+        lat2 = rlat + latshift
+        _e2, n2, _zonen2, _zonel2 = utm.from_latlon(lat2, rlon)
+        # NOTE: great circle distance used here, not reference ellipsoid
+        yshift = -(utm.haversine(rlon, rlat, rlon, lat2) + n2)
+        # cache the return value
+        xshift = None
+        if z in self.zoneshifts:
+            xshift = self.zoneshifts[z][0]
+        self.zoneshifts[z] = (xshift, yshift)
+        return yshift
+
+    def getutmzoneshift(self, e, n):
+        """
+        Given UTM easting and northing values, check if they fall outside
+        the reference point's zone boundary. Return the UTM coordinates in a
+        different zone and the new zone if they do. Zone lettering is only
+        changed when the reference point is in the opposite hemisphere.
+
+        :param e: easting value
+        :param n: northing value
+        :return: modified easting, northing, and zone values
+        :rtype: tuple
+        """
+        zone = self.refutm[0]
+        rlat, rlon, _ralt = self.refgeo
+        if e > 834000 or e < 166000:
+            num_zones = (int(e) - 166000) / (utm.R / 10)
+            # estimate number of zones to shift, E (positive) or W (negative)
+            rlon2 = self.refgeo[1] + (num_zones * 6)
+            _e2, _n2, zonen2, zonel2 = utm.from_latlon(rlat, rlon2)
+            xshift = utm.haversine(rlon, rlat, rlon2, rlat)
+            # after >3 zones away from refpt, the above estimate won't work
+            # (the above estimate could be improved)
+            if not 100000 <= (e - xshift) < 1000000:
+                # move one more zone away
+                num_zones = (abs(num_zones) + 1) * (abs(num_zones) / num_zones)
+                rlon2 = self.refgeo[1] + (num_zones * 6)
+                _e2, _n2, zonen2, zonel2 = utm.from_latlon(rlat, rlon2)
+                xshift = utm.haversine(rlon, rlat, rlon2, rlat)
+            e = e - xshift
+            zone = (zonen2, zonel2)
+        if n < 0:
+            # refpt in northern hemisphere and we crossed south of equator
+            n += 10000000
+            zone = (zone[0], 'M')
+        elif n > 10000000:
+            # refpt in southern hemisphere and we crossed north of equator
+            n -= 10000000
+            zone = (zone[0], 'N')
+        return e, n, zone

daemon/core/location/event.py

@@ -0,0 +1,248 @@
+"""
+event.py: event loop implementation using a heap queue and threads.
+"""
+
+import heapq
+import threading
+import time
+
+
+class Timer(threading.Thread):
+    """
+    Based on threading.Timer but cancel() returns if the timer was
+    already running.
+    """
+
+    def __init__(self, interval, function, args=None, kwargs=None):
+        """
+        Create a Timer instance.
+
+        :param interval: time interval
+        :param function: function to call when timer finishes
+        :param args: function arguments
+        :param kwargs: function keyword arguments
+        """
+        super(Timer, self).__init__()
+        self.interval = interval
+        self.function = function
+
+        self.finished = threading.Event()
+        self._running = threading.Lock()
+
+        # validate arguments were provided
+        if args:
+            self.args = args
+        else:
+            self.args = []
+
+        # validate keyword arguments were provided
+        if kwargs:
+            self.kwargs = kwargs
+        else:
+            self.kwargs = {}
+
+    def cancel(self):
+        """
+        Stop the timer if it hasn't finished yet.  Return False if
+        the timer was already running.
+
+        :return: True if canceled, False otherwise
+        :rtype: bool
+        """
+        locked = self._running.acquire(False)
+        if locked:
+            self.finished.set()
+            self._running.release()
+        return locked
+
+    def run(self):
+        """
+        Run the timer.
+
+        :return: nothing
+        """
+        self.finished.wait(self.interval)
+        with self._running:
+            if not self.finished.is_set():
+                self.function(*self.args, **self.kwargs)
+            self.finished.set()
+
+
+class Event(object):
+    """
+    Provides event objects that can be used within the EventLoop class.
+    """
+
+    def __init__(self, eventnum, event_time, func, *args, **kwds):
+        """
+        Create an Event instance.
+
+        :param eventnum: event number
+        :param event_time: event time
+        :param func: event function
+        :param args: function arguments
+        :param kwds: function keyword arguments
+        """
+        self.eventnum = eventnum
+        self.time = event_time
+        self.func = func
+        self.args = args
+        self.kwds = kwds
+        self.canceled = False
+
+    def __cmp__(self, other):
+        """
+        Comparison function.
+
+        :param Event other: event to compare with
+        :return: comparison result
+        :rtype: int
+        """
+        tmp = cmp(self.time, other.time)
+        if tmp == 0:
+            tmp = cmp(self.eventnum, other.eventnum)
+        return tmp
+
+    def run(self):
+        """
+        Run an event.
+
+        :return: nothing
+        """
+        if self.canceled:
+            return
+        self.func(*self.args, **self.kwds)
+
+    def cancel(self):
+        """
+        Cancel event.
+
+        :return: nothing
+        """
+        # XXX not thread-safe
+        self.canceled = True
+
+
+class EventLoop(object):
+    """
+    Provides an event loop for running events.
+    """
+
+    def __init__(self):
+        """
+        Creates a EventLoop instance.
+        """
+        self.lock = threading.RLock()
+        self.queue = []
+        self.eventnum = 0
+        self.timer = None
+        self.running = False
+        self.start = None
+
+    def __run_events(self):
+        """
+        Run events.
+
+        :return: nothing
+        """
+        schedule = False
+        while True:
+            with self.lock:
+                if not self.running or not self.queue:
+                    break
+                now = time.time()
+                if self.queue[0].time > now:
+                    schedule = True
+                    break
+                event = heapq.heappop(self.queue)
+            if event.time > now:
+                raise ValueError("invalid event time: %s > %s", event.time, now)
+            event.run()
+
+        with self.lock:
+            self.timer = None
+            if schedule:
+                self.__schedule_event()
+
+    def __schedule_event(self):
+        """
+        Schedule event.
+
+        :return: nothing
+        """
+        with self.lock:
+            if not self.running:
+                raise ValueError("scheduling event while not running")
+            if not self.queue:
+                return
+            delay = self.queue[0].time - time.time()
+            if self.timer:
+                raise ValueError("timer was already set")
+            self.timer = Timer(delay, self.__run_events)
+            self.timer.daemon = True
+            self.timer.start()
+
+    def run(self):
+        """
+        Start event loop.
+
+        :return: nothing
+        """
+        with self.lock:
+            if self.running:
+                return
+            self.running = True
+            self.start = time.time()
+            for event in self.queue:
+                event.time += self.start
+            self.__schedule_event()
+
+    def stop(self):
+        """
+        Stop event loop.
+
+        :return: nothing
+        """
+        with self.lock:
+            if not self.running:
+                return
+            self.queue = []
+            self.eventnum = 0
+            if self.timer is not None:
+                self.timer.cancel()
+                self.timer = None
+            self.running = False
+            self.start = None
+
+    def add_event(self, delaysec, func, *args, **kwds):
+        """
+        Add an event to the event loop.
+
+        :param float delaysec: delay in seconds for event
+        :param func: event function
+        :param args: event arguments
+        :param kwds: event keyword arguments
+        :return: created event
+        :rtype: Event
+        """
+        with self.lock:
+            eventnum = self.eventnum
+            self.eventnum += 1
+            evtime = float(delaysec)
+            if self.running:
+                evtime += time.time()
+            event = Event(eventnum, evtime, func, *args, **kwds)
+
+            if self.queue:
+                prevhead = self.queue[0]
+            else:
+                prevhead = None
+
+            heapq.heappush(self.queue, event)
+            head = self.queue[0]
+            if prevhead is not None and prevhead != head:
+                if self.timer is not None and self.timer.cancel():
+                    self.timer = None
+            if self.running and self.timer is None:
+                self.__schedule_event()
+        return event

daemon/core/location/utm.py

@@ -0,0 +1,259 @@
+"""
+utm
+===
+
+.. image:: https://travis-ci.org/Turbo87/utm.png
+
+Bidirectional UTM-WGS84 converter for python
+
+Usage
+-----
+
+::
+
+  import utm
+
+Convert a (latitude, longitude) tuple into an UTM coordinate::
+
+  utm.from_latlon(51.2, 7.5)
+  >>> (395201.3103811303, 5673135.241182375, 32, 'U')
+
+Convert an UTM coordinate into a (latitude, longitude) tuple::
+
+  utm.to_latlon(340000, 5710000, 32, 'U')
+  >>> (51.51852098408468, 6.693872395145327)
+
+Speed
+-----
+
+The library has been compared to the more generic pyproj library by running the
+unit test suite through pyproj instead of utm. These are the results:
+
+* with pyproj (without projection cache): 4.0 - 4.5 sec
+* with pyproj (with projection cache): 0.9 - 1.0 sec
+* with utm: 0.4 - 0.5 sec
+
+Authors
+-------
+
+* Tobias Bieniek <Tobias.Bieniek@gmx.de>
+
+License
+-------
+
+Copyright (C) 2012 Tobias Bieniek <Tobias.Bieniek@gmx.de>
+
+Permission is hereby granted, free of charge, to any person obtaining a copy of
+this software and associated documentation files (the "Software"), to deal in
+the Software without restriction, including without limitation the rights to
+use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies
+of the Software, and to permit persons to whom the Software is furnished to do
+so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.
+"""
+
+import math
+
+__all__ = ['to_latlon', 'from_latlon']
+
+
+class OutOfRangeError(ValueError):
+    pass
+
+
+K0 = 0.9996
+
+E = 0.00669438
+E2 = E * E
+E3 = E2 * E
+E_P2 = E / (1.0 - E)
+
+SQRT_E = math.sqrt(1 - E)
+_E = (1 - SQRT_E) / (1 + SQRT_E)
+_E3 = _E * _E * _E
+_E4 = _E3 * _E
+
+M1 = (1 - E / 4 - 3 * E2 / 64 - 5 * E3 / 256)
+M2 = (3 * E / 8 + 3 * E2 / 32 + 45 * E3 / 1024)
+M3 = (15 * E2 / 256 + 45 * E3 / 1024)
+M4 = (35 * E3 / 3072)
+
+P2 = (3 * _E / 2 - 27 * _E3 / 32)
+P3 = (21 * _E3 / 16 - 55 * _E4 / 32)
+P4 = (151 * _E3 / 96)
+
+R = 6378137
+
+ZONE_LETTERS = [
+    (84, None), (72, 'X'), (64, 'W'), (56, 'V'), (48, 'U'), (40, 'T'),
+    (32, 'S'), (24, 'R'), (16, 'Q'), (8, 'P'), (0, 'N'), (-8, 'M'), (-16, 'L'),
+    (-24, 'K'), (-32, 'J'), (-40, 'H'), (-48, 'G'), (-56, 'F'), (-64, 'E'),
+    (-72, 'D'), (-80, 'C')
+]
+
+
+def to_latlon(easting, northing, zone_number, zone_letter):
+    zone_letter = zone_letter.upper()
+
+    if not 100000 <= easting < 1000000:
+        raise OutOfRangeError('easting out of range (must be between 100.000 m and 999.999 m)')
+    if not 0 <= northing <= 10000000:
+        raise OutOfRangeError('northing out of range (must be between 0 m and 10.000.000 m)')
+    if not 1 <= zone_number <= 60:
+        raise OutOfRangeError('zone number out of range (must be between 1 and 60)')
+    if not 'C' <= zone_letter <= 'X' or zone_letter in ['I', 'O']:
+        raise OutOfRangeError('zone letter out of range (must be between C and X)')
+
+    x = easting - 500000
+    y = northing
+
+    if zone_letter < 'N':
+        y -= 10000000
+
+    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))
+
+    p_sin = math.sin(p_rad)
+    p_sin2 = p_sin * p_sin
+
+    p_cos = math.cos(p_rad)
+
+    p_tan = p_sin / p_cos
+    p_tan2 = p_tan * p_tan
+    p_tan4 = p_tan2 * p_tan2
+
+    ep_sin = 1 - E * p_sin2
+    ep_sin_sqrt = math.sqrt(1 - E * p_sin2)
+
+    n = R / ep_sin_sqrt
+    r = (1 - E) / ep_sin
+
+    c = _E * p_cos ** 2
+    c2 = c * c
+
+    d = x / (n * K0)
+    d2 = d * d
+    d3 = d2 * d
+    d4 = d3 * d
+    d5 = d4 * d
+    d6 = d5 * d
+
+    latitude = (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))
+
+    longitude = (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
+
+    return (math.degrees(latitude),
+            math.degrees(longitude) + zone_number_to_central_longitude(zone_number))
+
+
+def from_latlon(latitude, longitude):
+    if not -80.0 <= latitude <= 84.0:
+        raise OutOfRangeError('latitude out of range (must be between 80 deg S and 84 deg N)')
+    if not -180.0 <= longitude <= 180.0:
+        raise OutOfRangeError('northing out of range (must be between 180 deg W and 180 deg E)')
+
+    lat_rad = math.radians(latitude)
+    lat_sin = math.sin(lat_rad)
+    lat_cos = math.cos(lat_rad)
+
+    lat_tan = lat_sin / lat_cos
+    lat_tan2 = lat_tan * lat_tan
+    lat_tan4 = lat_tan2 * lat_tan2
+
+    lon_rad = math.radians(longitude)
+
+    zone_number = latlon_to_zone_number(latitude, longitude)
+    central_lon = zone_number_to_central_longitude(zone_number)
+    central_lon_rad = math.radians(central_lon)
+
+    zone_letter = latitude_to_zone_letter(latitude)
+
+    n = R / math.sqrt(1 - E * lat_sin ** 2)
+    c = E_P2 * lat_cos ** 2
+
+    a = lat_cos * (lon_rad - central_lon_rad)
+    a2 = a * a
+    a3 = a2 * a
+    a4 = a3 * a
+    a5 = a4 * a
+    a6 = a5 * a
+
+    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 ** 2) +
+                                        a6 / 720 * (61 - 58 * lat_tan2 + lat_tan4 + 600 * c - 330 * E_P2)))
+
+    if latitude < 0:
+        northing += 10000000
+
+    return easting, northing, zone_number, zone_letter
+
+
+def latitude_to_zone_letter(latitude):
+    for lat_min, zone_letter in ZONE_LETTERS:
+        if latitude >= lat_min:
+            return zone_letter
+
+    return None
+
+
+def latlon_to_zone_number(latitude, longitude):
+    if 56 <= latitude <= 64 and 3 <= longitude <= 12:
+        return 32
+
+    if 72 <= latitude <= 84 and longitude >= 0:
+        if longitude <= 9:
+            return 31
+        elif longitude <= 21:
+            return 33
+        elif longitude <= 33:
+            return 35
+        elif longitude <= 42:
+            return 37
+
+    return int((longitude + 180) / 6) + 1
+
+
+def zone_number_to_central_longitude(zone_number):
+    return (zone_number - 1) * 6 - 180 + 3
+
+
+def haversine(lon1, lat1, lon2, lat2):
+    """
+    Calculate the great circle distance between two points
+    on the earth (specified in decimal degrees)
+    """
+    # convert decimal degrees to radians
+    lon1, lat1, lon2, lat2 = map(math.radians, [lon1, lat1, lon2, lat2])
+    # haversine formula
+    dlon = lon2 - lon1
+    dlat = lat2 - lat1
+    a = math.sin(dlat / 2) ** 2 + math.cos(lat1) * math.cos(lat2) * math.sin(dlon / 2) ** 2
+    c = 2 * math.asin(math.sqrt(a))
+    m = 6367000 * c
+    return m