[Matplotlib-users] pyproj.Geod() azimuth problem

[Knut-Frode Dagestad]

Hi,

I have a question related to use of pyproj.
I am trying to calculate the azimuth orientation of vectors in a general coordinate system, defined by a proj4 string.
My approach is to calculate the azimuth angle of the y-axis of the given projection, and add this to the angle relative to the given coordinate system.

The illustration code below uses pyproj.Geod to calculate azimuth angle along the projection y-axis, but the output (bottom) does not make sense; the azimuth angle is found to be nearly the same for points at very different latitudes. This is the same for other projections such as e.g.:
proj4 = '+proj=lcc +lon_0=15 +lat_0=63 +lat_1=63 +lat_2=63 +R=6.371e+06 +units=m +no_defs'

However, using simply proj4 = '+proj=stere', I get results as expected. But as soon as I add '+lat_0' to the proj4-string, output is clearly wrong.

Does anyone have a clue what is going on here?

Best regards from Knut-Frode

--------------------------------------

import numpy as np
from mpl_toolkits.basemap import pyproj

proj4 = '+proj=stere +lat_0=90 +lon_0=70 +lat_ts=60 +units=m +a=6.371e+06 +e=0 +no_defs'
p = pyproj.Proj(proj4)

# Create three points for testing
lon_1 = np.array([30., 30., 30.])
lat_1 = np.array([60., 80., 89.])

x_1, y_1 = p(lon_1, lat_1) # Calculate projection (x,y) coordinates
lon_2, lat_2 = p(x_1, y_1 + 1000.0, inverse=True, errchk=True) # Find (lon,lat) of a point 1000 m away along y-axis

# Calculate azimuth from p1 to p2 (i.e. along projection y-axis)
g = pyproj.Geod(ellps='WGS84')
az = g.inv(lon_1, lat_1, lon_2, lat_2, radians=False)

# Output: [-40.04826805 -40.00599836 -40.00008458] - does not make sense
print az[0]