io_import_dxf/transverse_mercator.py

   1 # ##### BEGIN GPL LICENSE BLOCK #####
   2 #
   3 #  This program is free software; you can redistribute it and/or
   4 #  modify it under the terms of the GNU General Public License
   5 #  as published by the Free Software Foundation; either version 2
   6 #  of the License, or (at your option) any later version.
   7 #
   8 #  This program is distributed in the hope that it will be useful,
   9 #  but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 #  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  11 #  GNU General Public License for more details.
  12 #
  13 #  You should have received a copy of the GNU General Public License
  14 #  along with this program; if not, write to the Free Software Foundation,
  15 #  Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  16 #
  17 # ##### END GPL LICENSE BLOCK #####
  18
  19 # <pep8 compliant>
  20
  21 from math import sin, cos, atan, atanh, radians, tan, sinh, asin, cosh, degrees
  22
  23 # see conversion formulas at
  24 # http://en.wikipedia.org/wiki/Transverse_Mercator_projection
  25 # http://mathworld.wolfram.com/MercatorProjection.html
  26
  27
  28 class TransverseMercator:
  29     radius = 6378137
  30
  31     def __init__(self, lat=0, lon=0):
  32         self.lat = lat  # in degrees
  33         self.lon = lon  # in degrees
  34         self.lat_rad = radians(self.lat)
  35         self.lon_rad = radians(self.lon)
  36
  37     def fromGeographic(self, lat, lon):
  38         lat_rad = radians(lat)
  39         lon_rad = radians(lon)
  40         B = cos(lat_rad) * sin(lon_rad - self.lon_rad)
  41         x = self.radius * atanh(B)
  42         y = self.radius * (atan(tan(lat_rad) / cos(lon_rad - self.lon_rad)) - self.lat_rad)
  43         return x, y
  44
  45     def toGeographic(self, x, y):
  46         x /= self.radius
  47         y /= self.radius
  48         D = y + self.lat_rad
  49         lon = atan(sinh(x) / cos(D))
  50         lat = asin(sin(D) / cosh(x))
  51
  52         lon = self.lon + degrees(lon)
  53         lat = degrees(lat)
  54         return lat, lon