Lambert Conformal Conic implementation with US-specific transformation

Lambert.py

"""

>>> usa = USA()

>>> loc = Location(48.38544, -124.72916)
>>> pt = usa.locationProj(loc)
>>> abs(-2109470 - pt.x) < 1
True
>>> abs( 1532790 - pt.y) < 1
True

>>> loc = Location(25.28754, -80.36911)
>>> pt = usa.locationProj(loc)
>>> abs( 1585668 - pt.x) < 1
True
>>> abs(-1224099 - pt.y) < 1
True

>>> pt = Point(-2109470, 1532790)
>>> loc = usa.projLocation(pt)
>>> abs(  48.38544 - loc.lat) < 0.00001
True
>>> abs(-124.72916 - loc.lon) < 0.00001
True

>>> pt = Point(1585668, -1224099)
>>> loc = usa.projLocation(pt)
>>> abs( 25.28754 - loc.lat) < 0.00001
True
>>> abs(-80.36911 - loc.lon) < 0.00001
True

"""
from sys import stderr
from math import pi, log as ln, sin, cos, tan, atan, sqrt, copysign as sgn

def sec(t): return 1 / cos(t)
def cot(t): return 1 / tan(t)

def deg2rad(deg): return pi * deg / 180
def rad2deg(rad): return 180 * rad / pi

from ModestMaps.Geo import Transformation, IProjection, Location
from ModestMaps.Core import Point, Coordinate

try:
    from pyproj import Proj
except ImportError:
    Proj = False

class ConformalConic (IProjection):
    """
    """
    def __init__(self, lat0, lon0, lat1, lat2, xmin, ymin, xmax, ymax):
        """
        """
        self.srs = '+proj=lcc +lat_0=%(lat0)f +lon_0=%(lon0)f +lat_1=%(lat1)f +lat_2=%(lat2)f +a=6378137 +b=6378137 +k=1.0 +units=m +nadgrids=@null +no_defs' % locals()
        self.proj = Proj and Proj(self.srs)

        self.phi0 = deg2rad(lat0)
        self.lam0 = deg2rad(lon0)
        self.phi1 = deg2rad(lat1)
        self.phi2 = deg2rad(lat2)
        
        #
        # http://mathworld.wolfram.com/LambertConformalConicProjection.html
        #
        
        # Equation 4
        self.n = ln(cos(self.phi1) * sec(self.phi2))
        self.n /= ln(tan(pi/4 + self.phi2/2) * cot(pi/4 + self.phi1/2))
        
        # Equation 3
        self.F = cos(self.phi1) * pow(tan(pi/4 + self.phi1/2), self.n)
        self.F /= self.n
        
        # Equation 6
        self.P0 = self.F * pow(cot(pi/4 + self.phi0/2), self.n)
        
        #
        #
        #

        tx = 1 / (xmax - xmin)
        dx = .5 - (xmin/2 + xmax/2) * tx
        
        ty = -tx
        dy = .5 - (ymin/2 + ymax/2) * ty
        
        t = Transformation(tx, 0, dx, 0, ty, dy)
        
        IProjection.__init__(self, 0, t)
    
    def coordinateProj(self, coord):
        """ Convert from Coordinate object to a Point object in +proj=lcc
        """
        coord = coord.zoomTo(self.zoom)
        point = Point(coord.column, coord.row)
        return self.transformation.untransform(point)

    def projCoordinate(self, point):
        """ Convert from Point object in +proj=lcc to a Coordinate object
        """
        point = self.transformation.transform(point)
        coord = Coordinate(point.y, point.x, self.zoom)
        return coord

    def locationProj(self, location):
        """ Convert from Location object to a Point object in +proj=lcc
        """
        if Proj:
            x, y = self.proj(location.lon, location.lat)
            return Point(x, y)
        
        lam = deg2rad(location.lon)
        phi = deg2rad(location.lat)
        
        # Equation 5
        P = self.F * pow(cot(pi/4 + phi/2), self.n)
        
        # Equations 1 & 2
        x = P * sin(self.n * (lam - self.lam0))
        y = self.P0 - P * cos(self.n * (lam - self.lam0))
        
        return Point(x * 6378137., y * 6378137.)

    def projLocation(self, point):
        """ Convert from Point object in +proj=lcc to a Location object
        """
        if Proj:
            lon, lat = self.proj(point.x, point.y, inverse=True)
            return Location(lat, lon)
        
        x = point.x / 6378137.
        y = point.y / 6378137.
        
        # Equation 9
        P = sgn(1, self.n) * sqrt(x**2 + (self.P0 - y)**2)
        
        # Equation 10
        theta = atan(x / (self.P0 - y))
        
        # Equation 7
        phi = 2 * atan(pow(self.F/P, 1/self.n)) - pi/2
        
        # Equation 8
        lam = self.lam0 + theta/self.n
        
        lat = rad2deg(phi)
        lon = rad2deg(lam)
        
        return Location(lat, lon)

class USA (ConformalConic):
    
    def __init__(self):
        """
        """
        # adjust for latitude.
        scale = cos(deg2rad(37.5))
        
        # echo the usual +/- 20m extent of spherical mercator.
        xmin = scale * -6378137 * pi
        ymin = scale * -6378137 * pi

        xmax = scale * 6378137 * pi
        ymax = scale * 6378137 * pi
        
        ConformalConic.__init__(self, 37.5, -96, 29.5, 45.5, xmin, ymin, xmax, ymax)

if __name__ == '__main__':
    import doctest
    doctest.testmod()