grib2io.utils.arakawa_rotated_grid

  1"""
  2Functions for handling an Arakawa Rotated Lat/Lon Grids.
  3
  4This grid is not often used, but is currently used for the NCEP/RAP using
  5[GRIB2 Grid Definition Template 32769](https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_temp3-32769.shtml)
  6
  7These functions are adapted from the NCAR Command Language (ncl),
  8from [NcGRIB2.c](https://github.com/NCAR/ncl/blob/develop/ni/src/ncl/NclGRIB2.c)
  9"""
 10
 11import math
 12
 13from . import rotated_grid
 14
 15DEG2RAD = rotated_grid.DEG2RAD
 16RAD2DEG = rotated_grid.RAD2DEG
 17
 18
 19def ll2rot(latin: float, lonin: float, latpole: float, lonpole: float) -> tuple[float, float]:
 20    """
 21    Rotate a latitude/longitude pair.
 22
 23    Parameters
 24    ----------
 25    latin
 26        Unrotated latitude in units of degrees.
 27    lonin
 28        Unrotated longitude in units of degrees.
 29    latpole
 30        Latitude of Pole.
 31    lonpole
 32        Longitude of Pole.
 33
 34    Returns
 35    -------
 36    tlat
 37        Rotated latitude in units of degrees.
 38    tlons
 39        Rotated longitude in units of degrees.
 40    """
 41    tlon = lonin - lonpole
 42
 43    # Convert to xyz coordinates
 44    x = math.cos(latin * DEG2RAD) * math.cos(tlon * DEG2RAD)
 45    y = math.cos(latin * DEG2RAD) * math.sin(tlon * DEG2RAD)
 46    z = math.sin(latin * DEG2RAD)
 47
 48    # Rotate around y axis
 49    rotang = (latpole + 90) * DEG2RAD
 50    sinrot = math.sin(rotang)
 51    cosrot = math.cos(rotang)
 52    ry = y
 53    rx = x * cosrot + z * sinrot
 54    rz = -x * sinrot + z * cosrot
 55
 56    # Convert back to lat/lon
 57    tlat = math.asin(rz) / DEG2RAD
 58    if math.fabs(rx) > 0.0001:
 59        tlon = math.atan2(ry, rx) / DEG2RAD
 60    elif ry > 0:
 61        tlon = 90.0
 62    else:
 63        tlon = -90.0
 64
 65    if tlon < -180:
 66        tlon += 360.0
 67    if tlon >= 180:
 68        tlon -= 360.0
 69
 70    return tlat, tlon
 71
 72
 73def rot2ll(latin: float, lonin: float, latpole: float, lonpole: float) -> tuple[float, float]:
 74    """
 75    Unrotate a latitude/longitude pair.
 76
 77    Parameters
 78    ----------
 79    latin
 80        Rotated latitude in units of degrees.
 81    lonin
 82        Rotated longitude in units of degrees.
 83    latpole
 84        Latitude of Pole.
 85    lonpole
 86        Longitude of Pole.
 87
 88    Returns
 89    -------
 90    tlat
 91        Unrotated latitude in units of degrees.
 92    tlons
 93        Unrotated longitude in units of degrees.
 94    """
 95    tlon = lonin
 96
 97    # Convert to xyz coordinates
 98    x = math.cos(latin * DEG2RAD) * math.cos(lonin * DEG2RAD)
 99    y = math.cos(latin * DEG2RAD) * math.sin(lonin * DEG2RAD)
100    z = math.sin(latin * DEG2RAD)
101
102    # Rotate around y axis
103    rotang = -(latpole + 90) * DEG2RAD
104    sinrot = math.sin(rotang)
105    cosrot = math.cos(rotang)
106    ry = y
107    rx = x * cosrot + z * sinrot
108    rz = -x * sinrot + z * cosrot
109
110    # Convert back to lat/lon
111    tlat = math.asin(rz) / DEG2RAD
112    if math.fabs(rx) > 0.0001:
113        tlon = math.atan2(ry, rx) / DEG2RAD
114    elif ry > 0:
115        tlon = 90.0
116    else:
117        tlon = -90.0
118
119    # Remove the longitude rotation
120    tlon += lonpole
121    if tlon < 0:
122        tlon += 360.0
123    if tlon > 360:
124        tlon -= 360.0
125
126    return tlat, tlon
127
128
129def vector_rotation_angles(
130    tlat: float,
131    tlon: float,
132    clat: float,
133    losp: float,
134    xlat: float,
135) -> float:
136    """
137    Generate a rotation angle value.
138
139    The rotation angle value can be applied to a vector quantity to make it
140    Earth-oriented.
141
142    Parameters
143    ----------
144    tlat
145        True latitude in units of degrees.
146    tlon
147        True longitude in units of degrees..
148    clat
149        Latitude of center grid point in units of degrees.
150    losp
151        Longitude of the southern pole in units of degrees.
152    xlat
153        Latitude of the rotated grid in units of degrees.
154
155    Returns
156    -------
157    rot
158        Rotation angle in units of radians.
159    """
160    slon = math.sin((tlon - losp) * DEG2RAD)
161    cgridlat = math.cos(xlat * DEG2RAD)
162    if cgridlat <= 0.0:
163        rot = 0.0
164    else:
165        crot = (
166            math.cos(clat * DEG2RAD) * math.cos(tlat * DEG2RAD) + math.sin(clat * DEG2RAD) * math.sin(tlat * DEG2RAD) * math.cos(tlon * DEG2RAD)
167        ) / cgridlat
168        srot = (-1.0 * math.sin(clat * DEG2RAD) * slon) / cgridlat
169        rot = math.atan2(srot, crot)
170    return rot