source: lambert/lambert.py @ d5dacc74a117f96a578399fc0a885e9b09731390

#-*- encoding: utf-8 -*-
from math import *


# Some constants
X = 0         # Easting
Y = 1         # Northing
PHI = 0       # Latitude
LAMBDA = 1    # Longitude
EPS10 = 1.e-10
HALFPI = pi / 2.0
QUARTERPI = pi / 4.0

def gradesToRadians(grades):
  return grades / 200.0 * pi


def sexagesimalToDecimal(degrees=0.0, minutes=0.0, seconds=0.0):
  """ Convert a sexagesimal (degrees, minutes, seconds) to a decimal representation """
  return degrees + (minutes / 60.0) + (seconds / (60.0*60.0))


def clampPositive(number):
  """ Clamp negative numbers to 0.0 """
  if number < 0.0:
    return 0.0
  else:
    return number


class Ellipsoid:
  """ An ellipsoid ; can be specified using semi minor axis length or inverse flattening """
  
  def __init__(self, semiMajor, semiMinor=None, inverseFlattening=None):
    self.semiMajor = semiMajor
    self.semiMinor = semiMinor
    # Some precalculated properties
    if inverseFlattening is None:
      self.inverseFlattening = self.semiMajor / (self.semiMajor - self.semiMinor)
    else:
      self.inverseFlattening = inverseFlattening
      self.semiMinor = self.semiMajor - (self.semiMajor / self.inverseFlattening)
    self.flattening = 1.0 / self.inverseFlattening
    self.squaredExcentricity = 2.0 * self.flattening - self.flattening * self.flattening
    self.excentricity = sqrt(self.squaredExcentricity)


STD_ELLIPSOIDS = {
    "International 1924": Ellipsoid(6378388.0, inverseFlattening=297.0),
    "Clarke 1866": Ellipsoid(6378206.400, inverseFlattening=294.97870),
    "Clarke 1880 (IGN)": Ellipsoid(6378249.2, semiMinor=6356515.0),
    "GRS 1980": Ellipsoid(6378137.0, inverseFlattening=298.257222101),
    }


class LambertProjectionParameters2SP:
  """ Lambert conic conformal projection (2 standard parallels) parameters """
  
  def __init__(self, ellipsoid, firstParallel, secondParallel, origin, shift):
    self.ellipsoid = ellipsoid
    self.firstParallel = firstParallel
    self.secondParallel = secondParallel
    self.origin = origin
    self.shift = shift
  
  def __str__(self):
    return "first parallel: %f\nsecond_parallel: %f" % \
        (degrees(self.firstParallel), degrees(self.secondParallel))


def msfn(sinphi, cosphi, es):
  return (cosphi / sqrt(1.0 - es * sinphi * sinphi))


def tsfn(phi, sinphi, e):
  sinphi *= e;
  return (tan(0.5 * ((pi / 2.0) - phi)) /
      pow((1. - sinphi) / (1. + sinphi), .5 * e))


class LambertConicConformalProjection2SP:
  """ Lambert conic conformal projection (2 standard parallels) """
  
  _thetaAdjustment = 0.0          # Theta adjustment for belgium modified case
  
  def __init__(self, parameters):
    self.setParameters(parameters)

  def setParameters(self, parameters):
    self.parameters = parameters
    projectionLatitude1 = parameters.firstParallel
    projectionLatitude2 = parameters.secondParallel
    projectionLatitude = parameters.origin[1]
    e = self.parameters.ellipsoid.excentricity
    es = self.parameters.ellipsoid.squaredExcentricity
    self.n = sinphi = sin(projectionLatitude1)
    cosphi = cos(projectionLatitude1)
    secant = abs(projectionLatitude1 - projectionLatitude2) >= 1e-10
    spherical = (es == 0.0)
    if not spherical:
      m1 = msfn(sinphi, cosphi, es)
      ml1 = tsfn(projectionLatitude1, sinphi, e)
      if secant:
        sinphi = sin(projectionLatitude2)
        self.n = log(m1 /
            msfn(sinphi, cos(projectionLatitude2), es))
        self.n /= log(ml1 / tsfn(projectionLatitude2, sinphi, e))
      self.c = self.rho0 = m1 * pow(ml1, -self.n) / self.n
      if abs(abs(projectionLatitude) - HALFPI) < 1e-10:
        self.rho0 *= 0.0
      else:
        self.rho0 *= pow(tsfn(projectionLatitude, sin(projectionLatitude), e), self.n)
    else:
      if secant:
        self.n = log(cosphi / cos(projectionLatitude2)) / \
            log(tan(QUARTERPI + .5 * projectionLatitude2) / \
            tan(QUARTERPI + .5 * projectionLatitude1))
      self.c = cosphi * pow(tan(QUARTERPI + .5 * projectionLatitude1), self.n) / self.n
      if abs(abs(projectionLatitude) - HALFPI) < 1e-10:
        self.rho0 = 0.0
      else:
        self.rho0 = self.c * pow(tan(QUARTERPI + .5 * projectionLatitude), -n)

  def project(self, coordinates):
    """ Projects geographical coordinates to (x, y) using Lambert projection. """
    coordinates = (coordinates[0] - self.parameters.origin[0], coordinates[1])
    if fabs(fabs(coordinates[1]) - (pi / 2.0)) < EPS10:
      if coordinates[1] * self.n <= 0.0:
        raise Exception
      rho = 0.0
    else:
      if self.parameters.ellipsoid.squaredExcentricity == 0.0:
        rho = self.c * pow(tan((pi / 4.0) + 0.5 * coordinates[1]), -n)
      else:
        rho = self.c * pow(tsfn(coordinates[1], sin(coordinates[1]), self.parameters.ellipsoid.excentricity), self.n)
    a = coordinates[0] * self.n
    x = self.parameters.ellipsoid.semiMajor * (rho * sin(a)) + self.parameters.shift[0]
    y = self.parameters.ellipsoid.semiMajor * (self.rho0 - rho * cos(a)) + self.parameters.shift[1]
    return (x, y)


  def inverse_project(self, coordinates):
    """ Projects geographical coordinates to (x, y) using Lambert projection """


class LambertConicConformalProjection2SPBelgium(LambertConicConformalProjection2SP):
  """ Lambert conic conformal projection (2 standard parallels Belgium) """
  
  _thetaAdjustment = radians(sexagesimalToDecimal(seconds=29.2985))


STD_PROJECTIONS_PARAMETERS = {
  "Lambert II Carto": LambertProjectionParameters2SP(
      ellipsoid=STD_ELLIPSOIDS["Clarke 1880 (IGN)"],
      firstParallel=radians(sexagesimalToDecimal(degrees=45.0, minutes=53.0, seconds=56.108)),   # 45°53'56.108"
      secondParallel=radians(sexagesimalToDecimal(degrees=47.0, minutes=41.0, seconds=45.652)),  # 47°41'45.652"
      origin=(radians(sexagesimalToDecimal(2.0, 20.0, 14.025)), radians(sexagesimalToDecimal(46.0, 48.0))),
      shift=(600000.0, 2200000.0)),
  }

STD_PROJECTIONS = {
  "Lambert II Carto": LambertConicConformalProjection2SP(STD_PROJECTIONS_PARAMETERS["Lambert II Carto"]),
  }
      

def toDegrees(vector):
  return [degrees(x) for x in vector]

if __name__ == "__main__":
  projection = STD_PROJECTIONS["Lambert II Carto"]
  print projection.parameters
  coords = (radians(5.82), radians(43.95))
  result = projection.project(coords)
  #result = ((result[0] - 840005) / 16.666 , (result[1] - 1889995) / 16.666)
  print "coords: %s,\nprojected: %s" % (coords, result)
  #print "...", [degrees(x) for x in projection.unProject((840005, 1889995))]