Adding lat/lon generation for GDT 32769

Adding functions to support generation of lats and lons for GDT 32769.
These functions have been adapted from the NCAR command langauge (ncl).

This commit also adds a function to generate rotation angles for vector
quantities.  The angles should be applied to U/V-based grids in the
following manner:

         formula_u :	Uearth = sin(rot)*Vgrid + cos(rot)*Ugrid
         formula_v :	Vearth = cos(rot)*Vgrid - sin(rot)*Ugrid

This commit references issue #96

grib2io/_grib2io.py

@@ -942,6 +942,8 @@ def grid(self, unrotate=True):
         if self._sha1_section3 in _latlon_datastore.keys():
             return (_latlon_datastore[self._sha1_section3]['latitude'],
                     _latlon_datastore[self._sha1_section3]['longitude'])
+        else:
+            _latlon_datastore[self._sha1_section3] = {}
         gdtn = self.gridDefinitionTemplateNumber.value
         gdtmpl = self.gridDefinitionTemplate
         reggrid = self.gridDefinitionSection[2] == 0 # This means regular 2-d grid
@@ -1016,13 +1018,68 @@ def grid(self, unrotate=True):
             lats = np.where(abslats < 1.e20, lats, 1.e30)
         elif gdtn == 32769:
             # Special NCEP Grid, Rotated Lat/Lon, Arakawa E-Grid (Non-Staggered)
-            # TODO: Work in progress...
-            lons, lats= None, None
+            from grib2io.utils import arakawa_rotated_grid
+            from grib2io.utils.rotated_grid import DEG2RAD
+            di, dj = 0.0, 0.0
+            do_180 = False
+            idir = 1 if self.scanModeFlags[0] == 0 else -1
+            jdir = -1 if self.scanModeFlags[1] == 0 else 1
+            grid_oriented = 0 if self.resolutionAndComponentFlags[4] == 0 else 1
+            do_rot = 1 if grid_oriented == 1 else 0
+            la1 = self.latitudeFirstGridpoint
+            lo1 = self.longitudeFirstGridpoint
+            clon = self.longitudeCenterGridpoint
+            clat = self.latitudeCenterGridpoint
+            lasp = clat - 90.0
+            losp = clon
+            llat, llon = arakawa_rotated_grid.ll2rot(la1,lo1,lasp,losp)
+            la2, lo2 = arakawa_rotated_grid.rot2ll(-llat,-llon,lasp,losp)
+            rlat = -llat
+            rlon = -llon
+            if self.nx == 1:
+                di = 0.0
+            elif idir == 1:
+                ti = rlon
+                while ti < llon:
+                    ti += 360.0
+                di = (ti - llon)/float(self.nx-1)
+            else:
+                ti = llon
+                while ti < rlon:
+                    ti += 360.0
+                di = (ti - rlon)/float(self.nx-1)
+            if self.ny == 1:
+               dj = 0.0
+            else:
+                dj = (rlat - llat)/float(self.ny-1)
+                if dj < 0.0:
+                    dj = -dj
+            if idir == 1:
+                if llon > rlon:
+                    llon -= 360.0
+                if llon < 0 and rlon > 0:
+                    do_180 = True
+            else:
+                if rlon > llon:
+                    rlon -= 360.0
+                if rlon < 0 and llon > 0:
+                    do_180 = True
+            xlat1d = llat + (np.arange(self.ny)*jdir*dj)
+            xlon1d = llon + (np.arange(self.nx)*idir*di)
+            xlons, xlats = np.meshgrid(xlon1d,xlat1d)
+            rot2ll_vectorized = np.vectorize(arakawa_rotated_grid.rot2ll)
+            lats, lons = rot2ll_vectorized(xlats,xlons,lasp,losp)
+            if do_180:
+                lons = np.where(lons>180.0,lons-360.0,lons)
+            vector_rotation_angles_vectorized = np.vectorize(arakawa_rotated_grid.vector_rotation_angles)
+            rots = vector_rotation_angles_vectorized(lats, lons, clat, losp, xlats)
+            _latlon_datastore[self._sha1_section3]['vector_rotation_angles'] = rots
+            del xlat1d, xlon1d, xlats, xlons
         else:
             raise ValueError('Unsupported grid')
 
-        _latlon_datastore[self._sha1_section3] = dict(latitude=lats,
-                                                      longitude=lons)
+        _latlon_datastore[self._sha1_section3]['latitude'] = lats
+        _latlon_datastore[self._sha1_section3]['longitude'] = lons
 
         return lats, lons
 

grib2io/utils/arakawa_rotated_grid.py

@@ -0,0 +1,167 @@
+"""
+Functions for handling an Arakawa Rotated Lat/Lon Grids.
+
+This grid is not often used, but is currently used for the NCEP/RAP using
+[GRIB2 Grid Definition Template 32769](https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_temp3-32769.shtml)
+
+These functions are adapted from the NCAR Command Language (ncl),
+from [NcGRIB2.c](https://github.com/NCAR/ncl/blob/develop/ni/src/ncl/NclGRIB2.c)
+"""
+import math
+import numpy as np
+
+from . import rotated_grid
+
+DEG2RAD = rotated_grid.DEG2RAD
+RAD2DEG = rotated_grid.RAD2DEG
+
+def ll2rot(latin, lonin, latpole, lonpole):
+    """
+    Rotate a latitude/longitude pair.
+
+    Parameters
+    ----------
+    **`latin`**: `float`
+        Latitudes in units of degrees.
+
+    **`lonin`**: `float`
+        Longitudes in units of degrees.
+
+    **`latpole`**: `float`
+        Latitude of Pole.
+
+    **`lonpole`**: `float`
+        Longitude of Pole.
+
+    Returns
+    -------
+    **`tlat, tlons`**
+        Returns two floats of rotated latitude and longitude in units of degrees.
+    """
+    tlon = lonin - lonpole
+
+    # Convert to xyz coordinates
+    x = math.cos(latin * DEG2RAD) * math.cos(tlon * DEG2RAD)
+    y = math.cos(latin * DEG2RAD) * math.sin(tlon * DEG2RAD)
+    z = math.sin(latin * DEG2RAD)
+
+    # Rotate around y axis
+    rotang = (latpole + 90) * DEG2RAD
+    sinrot = math.sin(rotang)
+    cosrot = math.cos(rotang)
+    ry = y
+    rx = x * cosrot + z * sinrot
+    rz = -x * sinrot + z * cosrot
+
+    # Convert back to lat/lon
+    tlat = math.asin(rz) / DEG2RAD
+    if math.fabs(rx) > 0.0001:
+        tlon = math.atan2(ry,rx) / DEG2RAD
+    elif ry > 0:
+        tlon = 90.0
+    else:
+        tlon = -90.0
+
+    if tlon < -180:
+        tlon += 360.0
+    if tlon >= 180:
+        tlon -= 360.0
+
+    return tlat, tlon
+
+
+def rot2ll(latin, lonin, latpole, lonpole):
+    """
+    Unrotate a latitude/longitude pair.
+
+    Parameters
+    ----------
+    **`latin`**: `float`
+        Latitudes in units of degrees.
+
+    **`lonin`**: `float`
+        Longitudes in units of degrees.
+
+    **`latpole`**: `float`
+        Latitude of Pole.
+
+    **`lonpole`**: `float`
+        Longitude of Pole.
+
+    Returns
+    -------
+    **`tlat, tlons`**
+        Returns two floats of unrotated latitude and longitude in units of degrees.
+    """
+    tlon = lonin
+
+    # Convert to xyz coordinates
+    x = math.cos(latin * DEG2RAD) * math.cos(lonin * DEG2RAD)
+    y = math.cos(latin * DEG2RAD) * math.sin(lonin * DEG2RAD)
+    z = math.sin(latin * DEG2RAD)
+
+    # Rotate around y axis
+    rotang = -(latpole + 90) * DEG2RAD
+    sinrot = math.sin(rotang)
+    cosrot = math.cos(rotang)
+    ry = y
+    rx = x * cosrot + z * sinrot
+    rz = -x * sinrot + z * cosrot
+
+    # Convert back to lat/lon
+    tlat = math.asin(rz) / DEG2RAD
+    if math.fabs(rx) > 0.0001:
+        tlon = math.atan2(ry,rx) / DEG2RAD
+    elif ry > 0:
+        tlon = 90.0
+    else:
+        tlon = -90.0
+
+    # Remove the longitude rotation
+    tlon += lonpole
+    if tlon < 0:
+        tlon += 360.0
+    if tlon > 360:
+        tlon -= 360.0
+
+    return tlat, tlon
+
+
+def vector_rotation_angles(tlat, tlon, clat, losp, xlat):
+    """
+    Generate a rotation angle value that can be applied to a vector quantity to
+    make it Earth-oriented.
+
+    Parameters
+    ----------
+    **`tlat`**: `float`
+        True latitude in units of degrees.
+
+    **tlon`**: `float`
+        True longitude in units of degrees..
+
+    **`clat`**: `float`
+        Latitude of center grid point in units of degrees.
+
+    **`losp`**: `float`
+        Longitude of the southern pole in units of degrees.
+
+    **`xlat`**: `float`
+        Latitude of the rotated grid in units of degrees.
+
+    Returns
+    -------
+    **`rot : float`**
+        Rotation angle in units of radians. 
+    """
+    slon = math.sin((tlon-losp)*DEG2RAD)
+    cgridlat = math.cos(xlat*DEG2RAD)
+    if cgridlat <= 0.0:
+        rot = 0.0
+    else:
+        crot = (math.cos(clat*DEG2RAD)*math.cos(tlat*DEG2RAD)+
+                math.sin(clat*DEG2RAD)*math.sin(tlat*DEG2RAD)*
+                math.cos(tlon*DEG2RAD))/cgridlat
+        srot = (-1.0*math.sin(clat*DEG2RAD)*slon)/cgridlat
+        rot = math.atan2(srot,crot)
+    return rot