Merge pull request #33 from sebhahn/fix-proj-init

change pyproj.Geod to use as instance variable

pygeogrids/geodetic_datum.py

@@ -1,4 +1,4 @@
-# Copyright (c) 2015
+# Copyright (c) 2016
 # Vienna University of Technology, Department of Geodesy and Geoinformation
 # All rights reserved.
 
@@ -26,18 +26,11 @@
 # POSSIBILITY OF SUCH DAMAGE.
 
 
-"""
-Created on Nov 12, 2015
-
-@author: Christoph Reimer christoph.reimer@geo.tuwien.ac.at
-"""
-
-
 import numpy as np
 import pyproj
 
 
-class GeodeticDatum(pyproj.Geod):
+class GeodeticDatum():
     """
     Class representing a geodetic datum providing transformation and
     calculation functionality
@@ -48,13 +41,12 @@ class GeodeticDatum(pyproj.Geod):
         String of geodetic datum (ellipsoid) as provided in pyproj
 
     """
-    def __new__(self, ellString, **kwargs):
-        kwargs['ellps'] = ellString
-        return super(GeodeticDatum, self).__new__(self, None, **kwargs)
 
-    def __init__(self, ellString, **kwargs):
-        self.name = ellString
-        self.e = np.sqrt(self.es)
+    def __init__(self, ellps, **kwargs):
+        kwargs['ellps'] = ellps
+        self.geod = pyproj.Geod(**kwargs)
+        self.geod.e = np.sqrt(self.geod.es)
+        self.name = ellps
 
     def getParameter(self):
         """
@@ -73,7 +65,7 @@ def getParameter(self):
         x, y, z : np.array
             3D cartesian coordinates
         """
-        return self.a, self.b, self.f, self.e
+        return self.geod.a, self.geod.b, self.geod.f, self.geod.e
 
     def toECEF(self, lon, lat):
         """
@@ -102,7 +94,7 @@ def toECEF(self, lon, lat):
 
         x = N * np.cos(lat) * np.cos(lon)
         y = N * np.cos(lat) * np.sin(lon)
-        z = N * (1 - self.es) * np.sin(lat)
+        z = N * (1 - self.geod.es) * np.sin(lat)
 
         return x, y, z
 
@@ -140,7 +132,7 @@ def GeocentricLat(self, lat):
         if _element_iterable(lat):
             lat = np.array(lat, dtype=np.float64)
 
-        return np.rad2deg(np.arctan((1 - self.e ** 2) *
+        return np.rad2deg(np.arctan((1 - self.geod.e ** 2) *
                                     np.tan(np.deg2rad(lat))))
 
     def GeodeticLat(self, lat):
@@ -160,7 +152,7 @@ def GeodeticLat(self, lat):
         if _element_iterable(lat):
             lat = np.array(lat, dtype=np.float64)
         return np.rad2deg(np.arctan(np.tan(np.deg2rad(lat)) /
-                                    (1 - self.e ** 2)))
+                                    (1 - self.geod.e ** 2)))
 
     def ReducedLat(self, lat):
         """
@@ -178,7 +170,7 @@ def ReducedLat(self, lat):
         """
         if _element_iterable(lat):
             lat = np.array(lat, dtype=np.float64)
-        return np.rad2deg(np.arctan(np.sqrt(1 - self.e ** 2) *
+        return np.rad2deg(np.arctan(np.sqrt(1 - self.geod.e ** 2) *
                                     np.tan(np.deg2rad(lat))))
 
     def GeocentricDistance(self, lon, lat):
@@ -217,8 +209,8 @@ def EllN(self, lat):
 
         if _element_iterable(lat):
             lat = np.array(lat, dtype=np.float64)
-        return self.a / np.sqrt(1 - (self.es) *
-                                (np.sin(np.deg2rad(lat))) ** 2)
+        return self.geod.a / np.sqrt(1 - (self.geod.es) *
+                                     (np.sin(np.deg2rad(lat))) ** 2)
 
     def EllM(self, lat):
         """
@@ -236,8 +228,8 @@ def EllM(self, lat):
         """
         if _element_iterable(lat):
             lat = np.array(lat, dtype=np.float64)
-        return (self.a * (1 - self.es)) / \
-               ((1 - self.es) *
+        return (self.geod.a * (1 - self.geod.es)) / \
+               ((1 - self.geod.es) *
                 (np.sin(np.deg2rad(lat)) ** 2) ** (3. / 2.))
 
     def GaussianRadi(self, lat):
@@ -303,7 +295,7 @@ def MeridianArcDist(self, lat1, lat2):
         if _element_iterable(lat1) and lat1.shape == lat2.shape:
             lat1 = np.array(lat1, dtype=np.float64)
             lat2 = np.array(lat2, dtype=np.float64)
-        fazi, bazi, dist = self.inv(0., lat1, 0., lat2)
+        fazi, bazi, dist = self.geod.inv(0., lat1, 0., lat2)
         return dist
 
 
@@ -328,4 +320,4 @@ def _element_iterable(el):
     except (TypeError, IndexError):
         iterable = False
 
-    return iterable
+    return iterable

pygeogrids/netcdf.py

@@ -103,8 +103,8 @@ def save_lonlat(filename, arrlon, arrlat, geodatum, arrcell=None,
         crs = ncfile.createVariable('crs', np.dtype('int32').char)
         setattr(crs, 'grid_mapping_name', 'latitude_longitude')
         setattr(crs, 'longitude_of_prime_meridian', 0.)
-        setattr(crs, 'semi_major_axis', geodatum.a)
-        setattr(crs, 'inverse_flattening', 1. / geodatum.f)
+        setattr(crs, 'semi_major_axis', geodatum.geod.a)
+        setattr(crs, 'inverse_flattening', 1. / geodatum.geod.f)
         setattr(crs, 'ellipsoid_name', geodatum.name)
 
         gpi = ncfile.createVariable('gpi', np.dtype('int32').char, dim)

tests/test_geodatic_datum.py

@@ -1,4 +1,4 @@
-# Copyright (c) 2015
+# Copyright (c) 2016
 # Vienna University of Technology, Department of Geodesy and Geoinformation
 # All rights reserved.
 
@@ -25,13 +25,6 @@
 # ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 # POSSIBILITY OF SUCH DAMAGE.
 
-
-"""
-Created on Nov 12, 2015
-
-@author: Christoph Reimer christoph.reimer@geo.tuwien.ac.at
-"""
-
 import unittest
 import numpy.testing as nptest
 import numpy as np
@@ -45,35 +38,32 @@ def setUp(self):
 
     def test_toECEF(self):
         x, y, z = self.datum.toECEF(0, 90)
-        nptest.assert_almost_equal(np.array([0, 0, self.datum.b]),
+        nptest.assert_almost_equal(np.array([0, 0, self.datum.geod.b]),
                                    np.array([x, y, z]),
                                    decimal=5)
 
         x, y, z = self.datum.toECEF(0, 0)
-        nptest.assert_almost_equal(np.array([self.datum.a, 0, 0]),
+        nptest.assert_almost_equal(np.array([self.datum.geod.a, 0, 0]),
                                    np.array([x, y, z]),
                                    decimal=5)
 
     def test_ParallelRadi(self):
         r = self.datum.ParallelRadi(0.)
-        nptest.assert_almost_equal(r, self.datum.a, decimal=5)
+        nptest.assert_almost_equal(r, self.datum.geod.a, decimal=5)
 
         r = self.datum.ParallelRadi(90.)
         nptest.assert_almost_equal(r, 0., decimal=5)
 
     def test_GeocentricDistance(self):
         r = self.datum.GeocentricDistance(0., 0.)
-        nptest.assert_almost_equal(r, self.datum.a, decimal=5)
+        nptest.assert_almost_equal(r, self.datum.geod.a, decimal=5)
 
         r = self.datum.GeocentricDistance(0., 90.)
-        nptest.assert_almost_equal(r, self.datum.b, decimal=5)
+        nptest.assert_almost_equal(r, self.datum.geod.b, decimal=5)
 
     def test_ParallelArcDist(self):
         dist = self.datum.ParallelArcDist(0., 0., 360.)
-        nptest.assert_almost_equal(dist, self.datum.a * 2 * np.pi)
-
-
-
+        nptest.assert_almost_equal(dist, self.datum.geod.a * 2 * np.pi)
 
 
 if __name__ == "__main__":