From 2ad90cd41cb3ce4e78049f94f44473f04b9d13bb Mon Sep 17 00:00:00 2001 From: Eric Engle Date: Sat, 11 Nov 2023 17:41:26 -0500 Subject: [PATCH] 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 | 65 +++++++++- grib2io/utils/arakawa_rotated_grid.py | 167 ++++++++++++++++++++++++++ 2 files changed, 228 insertions(+), 4 deletions(-) create mode 100644 grib2io/utils/arakawa_rotated_grid.py diff --git a/grib2io/_grib2io.py b/grib2io/_grib2io.py index 59ded75..e8b160e 100644 --- a/grib2io/_grib2io.py +++ b/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 diff --git a/grib2io/utils/arakawa_rotated_grid.py b/grib2io/utils/arakawa_rotated_grid.py new file mode 100644 index 0000000..b8a54b9 --- /dev/null +++ b/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