grib2io.utils.rotated_grid

  1"""Tools for working with Rotated Lat/Lon Grids."""
  2
  3import numpy as np
  4from numpy.typing import NDArray
  5
  6RAD2DEG = 57.29577951308232087684
  7DEG2RAD = 0.01745329251994329576
  8
  9
 10def rotate(
 11    latin: NDArray[np.float32],
 12    lonin: NDArray[np.float32],
 13    aor: NDArray[np.float32],
 14    splat: NDArray[np.float32],
 15    splon: NDArray[np.float32],
 16) -> tuple[NDArray[np.float32], NDArray[np.float32]]:
 17    """
 18    Perform grid rotation.
 19
 20    This function is adapted from ECMWF's ecCodes library void function,
 21    rotate().
 22
 23    https://github.com/ecmwf/eccodes/blob/develop/src/grib_geography.cc
 24
 25    Parameters
 26    ----------
 27    latin
 28        Latitudes in units of degrees.
 29    lonin
 30        Longitudes in units of degrees.
 31    aor
 32        Angle of rotation as defined in GRIB2 GDTN 4.1.
 33    splat
 34        Latitude of South Pole as defined in GRIB2 GDTN 4.1.
 35    splon
 36        Longitude of South Pole as defined in GRIB2 GDTN 4.1.
 37
 38    Returns
 39    -------
 40    lats
 41        `numpy.ndarrays` with `dtype=numpy.float32` of grid latitudes in units
 42        of degrees.
 43    lons
 44        `numpy.ndarrays` with `dtype=numpy.float32` of grid longitudes in units
 45        of degrees.
 46    """
 47    zsycen = np.sin(DEG2RAD * (splat + 90.0))
 48    zcycen = np.cos(DEG2RAD * (splat + 90.0))
 49    zxmxc = DEG2RAD * (lonin - splon)
 50    zsxmxc = np.sin(zxmxc)
 51    zcxmxc = np.cos(zxmxc)
 52    zsyreg = np.sin(DEG2RAD * latin)
 53    zcyreg = np.cos(DEG2RAD * latin)
 54    zsyrot = zcycen * zsyreg - zsycen * zcyreg * zcxmxc
 55
 56    zsyrot = np.where(zsyrot > 1.0, 1.0, zsyrot)
 57    zsyrot = np.where(zsyrot < -1.0, -1.0, zsyrot)
 58
 59    pyrot = np.arcsin(zsyrot) * RAD2DEG
 60
 61    zcyrot = np.cos(pyrot * DEG2RAD)
 62    zcxrot = (zcycen * zcyreg * zcxmxc + zsycen * zsyreg) / zcyrot
 63    zcxrot = np.where(zcxrot > 1.0, 1.0, zcxrot)
 64    zcxrot = np.where(zcxrot < -1.0, -1.0, zcxrot)
 65    zsxrot = zcyreg * zsxmxc / zcyrot
 66
 67    pxrot = np.arccos(zcxrot) * RAD2DEG
 68
 69    pxrot = np.where(zsxrot < 0.0, -pxrot, pxrot)
 70
 71    return pyrot, pxrot
 72
 73
 74def unrotate(
 75    latin: NDArray[np.float32],
 76    lonin: NDArray[np.float32],
 77    aor: NDArray[np.float32],
 78    splat: NDArray[np.float32],
 79    splon: NDArray[np.float32],
 80) -> tuple[NDArray[np.float32], NDArray[np.float32]]:
 81    """
 82    Perform grid un-rotation.
 83
 84    This function is adapted from ECMWF's ecCodes library void function,
 85    unrotate().
 86
 87    https://github.com/ecmwf/eccodes/blob/develop/src/grib_geography.cc
 88
 89    Parameters
 90    ----------
 91    latin
 92        Latitudes in units of degrees.
 93    lonin
 94        Longitudes in units of degrees.
 95    aor
 96        Angle of rotation as defined in GRIB2 GDTN 4.1.
 97    splat
 98        Latitude of South Pole as defined in GRIB2 GDTN 4.1.
 99    splon
100        Longitude of South Pole as defined in GRIB2 GDTN 4.1.
101
102    Returns
103    -------
104    lats
105        `numpy.ndarrays` with `dtype=numpy.float32` of grid latitudes in units
106        of degrees.
107    lons
108        `numpy.ndarrays` with `dtype=numpy.float32` of grid longitudes in units
109        of degrees.
110    """
111    lon_x = lonin
112    lat_y = latin
113
114    latr = lat_y * DEG2RAD
115    lonr = lon_x * DEG2RAD
116
117    xd = np.cos(lonr) * np.cos(latr)
118    yd = np.sin(lonr) * np.cos(latr)
119    zd = np.sin(latr)
120
121    t = -(90.0 + splat)
122    o = -splon
123
124    sin_t = np.sin(DEG2RAD * t)
125    cos_t = np.cos(DEG2RAD * t)
126    sin_o = np.sin(DEG2RAD * o)
127    cos_o = np.cos(DEG2RAD * o)
128
129    x = cos_t * cos_o * xd + sin_o * yd + sin_t * cos_o * zd
130    y = -cos_t * sin_o * xd + cos_o * yd - sin_t * sin_o * zd
131    z = -sin_t * xd + cos_t * zd
132
133    ret_lat = 0
134    ret_lon = 0
135
136    # Then convert back to 'normal' (lat,lon)
137    # Uses arcsin, to convert back to degrees, put in range -1 to 1 in case of slight rounding error
138    # avoid error on calculating e.g. asin(1.00000001)
139    z = np.where(z > 1.0, 1.0, z)
140    z = np.where(z < -1.0, -1.0, z)
141
142    ret_lat = np.arcsin(z) * RAD2DEG
143    ret_lon = np.arctan2(y, x) * RAD2DEG
144
145    # Still get a very small rounding error, round to 6 decimal places
146    ret_lat = np.round(ret_lat * 1000000.0) / 1000000.0
147    ret_lon = np.round(ret_lon * 1000000.0) / 1000000.0
148
149    ret_lon -= aor
150
151    return ret_lat, ret_lon