Regression in radians calculations

[chrrrisw]

In moving from 1.8.9 to 1.9.5.1 I've noticed that Geod no longer performs calculations in radians correctly.

The following prints different numbers for the two calculations in Python 3.4.3:

import math
from pyproj import Geod

g = Geod(ellps='clrk66')  # Use Clarke 1966 ellipsoid.

# specify the lat/lons of some cities.
boston_lat = 42. + (15. / 60.)
boston_lon = -71. - (7. / 60.)

portland_lat = 45. + (31. / 60.)
portland_lon = -123. - (41. / 60.)

# compute forward and back azimuths, plus distance between Boston and Portland.

# First in degrees
az12_d, az21_d, dist_d = g.inv(
    boston_lon,
    boston_lat,
    portland_lon,
    portland_lat,
    radians=False)

# Now in radians
az12_r, az21_r, dist_r = g.inv(
    math.radians(boston_lon),
    math.radians(boston_lat),
    math.radians(portland_lon),
    math.radians(portland_lat),
    radians=True)

print(az12_d, az21_d, dist_d)
print(math.degrees(az12_r), math.degrees(az21_r), dist_r)

The same occurs for g.fwd() and g.npts().

Regards,
Chris.

[micahcochran]

Results from pyproj 1.9.5.1 Python version 3.4.3 using xonsh.

Results of the "First in degrees" (leaving the variable names out):
(-66.53059478766228, 75.65363415556968, 4164192.7080994677)

Results of the "Now in radians":
(-3811.9222898281223, 4334.633941940914, 4164192.7080994677)

You basically did this--which seems correct to me--which gives a rather large numbers.

>>> math.degrees(-3811.9222898281223)
-218407.059038996
>>> math.degrees(4334.633941940914)
248356.2306073695

If I do this, it matches the degrees. Seems a bit strange.

>>> math.radians(-3811.9222898281223)
-66.53059478766228
>>> math.radians(4334.633941940914)
75.65363415556968

Perhaps there is wrong conversion is being performed in the code. I really don't know. I'm not as familiar with the Geod class. Here is the documentation for the geodesic.h library, which the functions only deals with degrees. So, the radians conversions are pyproj feature.

[chrrrisw]

I've only verified the inv() changes, but the following diff seems to help - I'm working on unit test cases at the moment.

diff --git a/_proj.pyx b/_proj.pyx
index e0c0775..b331471 100644
--- a/_proj.pyx
+++ b/_proj.pyx
@@ -541,9 +541,9 @@ cdef class Geod:
                 az1 = azdata[i]
                 s12 = distdata[i]
             else:
-                lon1 = _dg2rad*lonsdata[i]
-                lat1 = _dg2rad*latsdata[i]
-                az1 = _dg2rad*azdata[i]
+                lon1 = _rad2dg*lonsdata[i]
+                lat1 = _rad2dg*latsdata[i]
+                az1 = _rad2dg*azdata[i]
                 s12 = distdata[i]
             geod_direct(&self._geod_geodesic, lat1, lon1, az1, s12,\
                    &plat2, &plon2, &pazi2)
@@ -561,9 +561,9 @@ cdef class Geod:
                 latsdata[i] = plat2
                 azdata[i] = pazi2
             else:
-                lonsdata[i] = _rad2dg*plon2
-                latsdata[i] = _rad2dg*plat2
-                azdata[i] = _rad2dg*pazi2
+                lonsdata[i] = _dg2rad*plon2
+                latsdata[i] = _dg2rad*plat2
+                azdata[i] = _dg2rad*pazi2

     def _inv(self, object lons1, object lats1, object lons2, object lats2, radians=False):
         """
@@ -621,8 +621,8 @@ cdef class Geod:
             if ps12 != ps12: # check for NaN
                 raise ValueError('undefined inverse geodesic (may be an antipodal point)')
             if radians:
-                lonsdata[i] = _rad2dg*pazi1
-                latsdata[i] = _rad2dg*pazi2
+                lonsdata[i] = _dg2rad*pazi1
+                latsdata[i] = _dg2rad*pazi2
             else:
                 lonsdata[i] = pazi1
                 latsdata[i] = pazi2

[micahcochran]

I did a little bit of comparing the code between 1.8.9 and 1.9.5.1. The differences included the the geodesic library is used in 1.9.5.1, where the Proj.4 had its own native pj_fwd and pj_inv functions (which are now "obsolete"). Proj.4 native functions used radians as their units.

When code was converted, the conversion factors were left alone instead of being switched.

[chrrrisw]

Thanks Micah,
that seems to match the fix, which essentially converts from radians to degrees and then back to radians (instead of the other way around).