Skip to content

Instantly share code, notes, and snippets.

@migurski

migurski/Albers.py

Created Jun 28, 2013
Embed
What would you like to do?
'''
>>> arc = Projection()
Geographic to Albers:
>>> p = arc.locationProj(Location(23, -96)); '%.2f %.2f' % (abs(p.x), abs(p.y))
'0.00 0.00'
>>> p = arc.locationProj(Location(33, -106)); '%.2f %.2f' % (p.x, p.y)
'-926031.64 1156991.11'
>>> p = arc.locationProj(Location(33, -86)); '%.2f %.2f' % (p.x, p.y)
'926031.64 1156991.11'
>>> p = arc.locationProj(Location(13, -86)); '%.2f %.2f' % (p.x, p.y)
'1154734.78 -1008658.24'
>>> p = arc.locationProj(Location(13, -96)); '%.2f %.2f' % (p.x, p.y)
'0.00 -1069461.99'
Albers to geographic:
>>> l = arc.projLocation(Point(-2000000, 2000000)); '%.6f %.6f' % (l.lon, l.lat)
'-119.492821 38.728671'
>>> l = arc.projLocation(Point(2000000, 2000000)); '%.6f %.6f' % (l.lon, l.lat)
'-72.507179 38.728671'
>>> l = arc.projLocation(Point(2000000, -2000000)); '%.6f %.6f' % (l.lon, l.lat)
'-80.207835 2.136423'
>>> l = arc.projLocation(Point(-2000000, -2000000)); '%.6f %.6f' % (l.lon, l.lat)
'-111.792165 2.136423'
Albers to tiles:
>>> c = arc.projCoordinate(Point(0, 0)); '%d %.3f %.3f' % (c.zoom, c.column, c.row)
'0 0.500 0.500'
>>> c = arc.projCoordinate(Point(-4000000, 4000000)); '%d %.6f %.6f' % (c.zoom, c.column, c.row)
'0 0.000000 0.000000'
>>> c = arc.projCoordinate(Point(4000000, 4000000)); '%d %.6f %.6f' % (c.zoom, c.column, c.row)
'0 1.000000 0.000000'
>>> c = arc.projCoordinate(Point(4000000, -4000000)); '%d %.6f %.6f' % (c.zoom, c.column, c.row)
'0 1.000000 1.000000'
>>> c = arc.projCoordinate(Point(-4000000, -4000000)); '%d %.6f %.6f' % (c.zoom, c.column, c.row)
'0 0.000000 1.000000'
Tiles to Albers:
>>> p = arc.coordinateProj(Coordinate(.5, .5, 0)); '%.3f %.3f' % (abs(p.x), abs(p.y))
'0.000 0.000'
>>> p = arc.coordinateProj(Coordinate(0, 0, 0)); '%.3f %.3f' % (p.x, p.y)
'-4000000.000 4000000.000'
>>> p = arc.coordinateProj(Coordinate(0, 1, 0)); '%.3f %.3f' % (p.x, p.y)
'4000000.000 4000000.000'
>>> p = arc.coordinateProj(Coordinate(1, 1, 0)); '%.3f %.3f' % (p.x, p.y)
'4000000.000 -4000000.000'
>>> p = arc.coordinateProj(Coordinate(1, 0, 0)); '%.3f %.3f' % (p.x, p.y)
'-4000000.000 -4000000.000'
Tiles to geographic:
>>> l = arc.coordinateLocation(Coordinate(0, 0, 0)); '%.6f %.6f' % (l.lon, l.lat)
'-152.432839 47.900775'
>>> l = arc.coordinateLocation(Coordinate(1, 1, 0)); '%.6f %.6f' % (l.lon, l.lat)
'-69.415736 -25.767692'
Geographic to tiles:
>>> c = arc.locationCoordinate(Location(47.900775, -152.432839)); '%d %.6f %.6f' % (c.zoom, abs(c.column), abs(c.row))
'0 0.000000 0.000000'
>>> c = arc.locationCoordinate(Location(-25.767692, -69.415736)); '%d %.6f %.6f' % (c.zoom, c.column, c.row)
'0 1.000000 1.000000'
'''
from math import pi
from ModestMaps.Geo import IProjection, Location, deriveTransformation
from ModestMaps.Core import Point, Coordinate
from pyproj import Proj
class Projection (IProjection):
'''
'''
srs = '+proj=aea +lat_0=23 +lon_0=-96 +lat_1=29.5 +lat_2=45.5 ' \
+ '+nadgrids=@null +no_defs +a=6378137 +b=6378137'
def __init__(self):
# Transform from raw Albers projection to tile coordinates
t = deriveTransformation(-4000000, 4000000, 0, 0,
4000000, 4000000, 1, 0,
-4000000, -4000000, 0, 1)
IProjection.__init__(self, 0, t)
self.proj = Proj(self.srs)
def rawProject(self, point):
''' ? http://mathworld.wolfram.com/AlbersEqual-AreaConicProjection.html
'''
return Point(*self.proj(point.x, point.y, radians=True))
def rawUnproject(self, point):
''' ? http://mathworld.wolfram.com/AlbersEqual-AreaConicProjection.html
'''
return Point(*self.proj(point.x, point.y, inverse=True, radians=True))
def coordinateProj(self, coord):
""" Convert from Coordinate object to a Point object in Albers
"""
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 Albers to a Coordinate object
"""
point = self.transformation.transform(point)
return Coordinate(point.y, point.x, self.zoom)
def locationProj(self, location):
""" Convert from Location object to a Point object in Albers
"""
point = Point(pi * location.lon / 180, pi * location.lat / 180)
return self.rawProject(point)
def projLocation(self, point):
""" Convert from Point object in Albers to a Location object
"""
point = self.rawUnproject(point)
return Location(180 * point.y / pi, 180 * point.x / pi)
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.