Created
June 28, 2013 17:57
-
-
Save migurski/5886670 to your computer and use it in GitHub Desktop.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
''' | |
>>> 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