Source code for epygram.myproj
#!/usr/bin/env python
# -*- coding: utf-8 -*-
import math
import numpy
from epygram import config
from epygram.base import RecursiveObject
[docs]class Proj(RecursiveObject):
"""
A class of homemade projections on the sphere, made to handle cases not handled by *pyproj*.
(e.g. sphere geoid for polar stereographic projection...).
"""
[docs] def __init__(self, proj, geoidshape=config.myproj_default_geoid['geoidshape'],
geoidradius=config.myproj_default_geoid['geoidradius'],
x_0=0., y_0=0., **proj_params):
"""
Initializes parameters for projection formulas.
Args:
- *proj*: name of the projection, among ('lambert', 'mercator', 'polar_stereographic')
- *geoidshape*: actually, only 'sphere' is implemented.
- *earth_radius*: can be specified, to use a specific earth radius, in m.
- *x_0*, *y_0* : offset of Origin in (x, y) coordinates of projected map.
- *proj_params*: a set of arguments that depends on the projection.
(all lon/lat are in degrees):
+ lambert: lon_0 = longitude of reference point \n
lat_1 = first secant latitude \n
lat_2 = second secant latitude \n
if tangent, lat_1 = lat_2 = lat_0 = latitude of reference point = tangency latitude
+ mercator: lon_0 = longitude of reference point
lat_ts = tangency latitude (0°) or secant latitude
+ polar_stereographic: lon_0 = longitude of reference point
lat_0 = +/- 90° depending on the projection pole
lat_ts = secant or tangency latitude
"""
deg2rad = math.pi/180.
self._args = proj_params
self._p = {'name':proj}
if geoidshape != 'sphere':
raise NotImplementedError("projected geometry on ellipsoid.")
else:
self._p['earth_radius'] = geoidradius
self._p['FE'] = x_0
self._p['FN'] = y_0
self._p['lambda0'] = proj_params['lon_0']*deg2rad
if proj == 'lambert':
if abs(proj_params['lat_1'] - proj_params['lat_2']) > config.epsilon:
self.secant = True
m1 = math.cos(proj_params['lat_1']*deg2rad)
m2 = math.cos(proj_params['lat_2']*deg2rad)
t1 = math.tan(math.pi/4. - proj_params['lat_1']*deg2rad/2.)
t2 = math.tan(math.pi/4. - proj_params['lat_2']*deg2rad/2.)
lat_0 = (proj_params['lat_1'] + proj_params['lat_1']) / 2.
t0 = math.tan(math.pi/4. - lat_0*deg2rad/2.)
self._p['n'] = (math.log(m1) - math.log(m2)) / (math.log(t1) - math.log(t2))
self._p['F'] = m1 / (self._p['n'] * (t1**self._p['n']))
self._p['r0'] = self._p['earth_radius'] * self._p['F'] * t0**self._p['n']
else:
self.secant = False
lat_0 = proj_params['lat_1']
m0 = math.cos(lat_0*deg2rad)
t0 = math.tan(math.pi/4. - lat_0*deg2rad/2.)
self._p['n'] = math.sin(lat_0*deg2rad)
self._p['F'] = m0 / (self._p['n'] * t0**self._p['n'])
self._p['r0'] = self._p['earth_radius'] * self._p['F'] * t0**self._p['n']
self._p['sign'] = math.copysign(1., lat_0)
elif proj == 'mercator':
if abs(proj_params['lat_ts']) > config.epsilon:
self.secant = True
self._p['k0'] = math.cos(proj_params['lat_ts']*deg2rad)
else:
self.secant = False
self._p['k0'] = 1.
elif proj == 'polar_stereographic':
self._p['sign'] = math.copysign(1., proj_params['lat_0'])
if abs(proj_params['lat_0'] - proj_params['lat_ts']) > config.epsilon:
self.secant = True
mf = math.cos(proj_params['lat_ts']*deg2rad)
tf = math.tan(math.pi/4. - self._p['sign']*proj_params['lat_ts']/2.)
self._p['k0'] = mf / (2.*tf)
else:
self.secant = False
self._p['k0'] = 1.
[docs] def __call__(self, x, y, inverse=False):
"""
Converts lon/lat coordinates (in °) to x/y coord in the projection.
If *inverse* is *True*, makes the inverse conversion, from x/y to lon/lat (in °).
"""
if not isinstance(x, numpy.ndarray):
x = numpy.array(x)
if not isinstance(y, numpy.ndarray):
y = numpy.array(y)
p = self._p
if not inverse:
lon = numpy.pi/180. * x
lat = numpy.pi/180. * y
if self._p['name'] == 'lambert':
r = p['earth_radius'] * p['F'] * numpy.tan(numpy.pi/4. - lat/2.)**p['n']
theta = p['n'] * (lon - p['lambda0'])
E = p['FE'] + r * numpy.sin(theta)
N = p['FN'] + p['r0'] - r * numpy.cos(theta)
elif self._p['name'] == 'mercator':
E = p['FE'] + p['earth_radius'] * p['k0'] * (lon - p['lambda0'])
N = p['FN'] + p['earth_radius'] * p['k0'] * numpy.log(numpy.tan(numpy.pi/4. + lat/2.))
elif self._p['name'] == 'polar_stereographic':
t = numpy.tan(numpy.pi/4. - p['sign']*lat/2.)
rho = 2. * p['earth_radius'] * p['k0'] * t
E = p['FE'] + rho * numpy.sin(lon - p['lambda0'])
N = p['FN'] - p['sign'] * rho * numpy.cos(lon - p['lambda0'])
return (E, N)
elif inverse:
if self._p['name'] == 'lambert':
rp = p['sign'] * numpy.sqrt((x-p['FE'])**2 + (p['r0']-(y-p['FN']))**2)
tp = (rp/(p['earth_radius']*p['F']))**(1./p['n'])
thetap = numpy.arctan((x-p['FE'])/(p['r0']-(y-p['FN'])))
lon = thetap/p['n'] + p['lambda0']
lat = numpy.pi/2. - 2*numpy.arctan(tp)
elif self._p['name'] == 'mercator':
lon = (x - p['FE'])/(p['earth_radius']*p['k0']) + p['lambda0']
lat = 2*numpy.arctan(numpy.exp((y-p['FN'])/(p['earth_radius']*p['k0']))) - numpy.pi/2.
elif self._p['name'] == 'polar_stereographic':
rhop = numpy.sqrt((x - p['FE'])**2 + (y - p['FN'])**2)
tp = rhop / (2 * p['earth_radius'] * p['k0'])
lon = p['lambda0'] -(0.5+0.5*p['sign']*numpy.copysign(1., y-p['FN']))* numpy.pi \
+ numpy.arctan(-p['sign'] * (x-p['FE'])/(y-p['FN']))
lat = p['sign'] * (numpy.pi/2. - 2*numpy.arctan(tp))
return (180./numpy.pi * lon, 180./numpy.pi * lat)