Skip to content

Instantly share code, notes, and snippets.

@migurski

migurski/Lambert.py

Created Jul 5, 2012
Embed
What would you like to do?
Lambert Conformal Conic implementation with US-specific transformation
"""
>>> 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()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.