Albers.py

'''
>>> 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()