FIX: Rotate HRRR grid-relative winds to true north before use as constraint

pydda/__init__.py

@@ -12,7 +12,8 @@
 from . import constraints
 from . import io
 
-__version__ = "2.4.1"
+
+__version__ = "2.5.0"
 
 print("Welcome to PyDDA %s" % __version__)
 print("If you are using PyDDA in your publications, please cite:")

pydda/constraints/model_data.py

@@ -533,10 +533,73 @@ def make_constraint_from_wrf(Grid, file_path, wrf_time, radar_loc, vel_field=Non
     return Grid
 
 
-def add_hrrr_constraint_to_grid(Grid, file_path, method="nearest"):
+def _uv_grid_relative_to_true_north(
+    u, v, lon, lat, lov=-97.5, truelat1=38.5, truelat2=38.5
+):
+    """
+    HRRR (and other WRF-based models) output u and v components that are
+    relative to the model's Lambert Conformal Conic grid rather than true
+    north. This rotates them to earth-relative (true north) components
+    using the standard WRF/NCEP grid-rotation formula (equivalent to the
+    DCOMPUTEUVMET routine used by wrf-python's ``uvmet`` diagnostic).
+
+    Parameters
+    ----------
+    u, v: array_like
+        Grid-relative u and v wind components. Must be broadcastable
+        against lon/lat.
+    lon, lat: array_like
+        Longitude (in the range [-180, 180]) and latitude of each point,
+        broadcastable against u and v.
+    lov: float
+        Standard longitude (line of view/center longitude) of the Lambert
+        Conformal Conic projection, in degrees. Defaults to -97.5, the
+        value used by the operational CONUS HRRR grid.
+    truelat1, truelat2: float
+        The two standard (true) latitudes of the Lambert Conformal Conic
+        projection, in degrees. Default to 38.5, the values used by the
+        operational CONUS HRRR grid (a tangent cone).
+
+    Returns
+    -------
+    u_true, v_true: array_like
+        The u and v wind components relative to true north.
+    """
+    rad_per_deg = np.pi / 180.0
+    if np.isclose(truelat1, truelat2):
+        cone = np.sin(np.abs(truelat1) * rad_per_deg)
+    else:
+        cone = (
+            np.log(np.cos(truelat1 * rad_per_deg))
+            - np.log(np.cos(truelat2 * rad_per_deg))
+        ) / (
+            np.log(np.tan(np.pi / 4.0 + truelat2 * rad_per_deg / 2.0))
+            - np.log(np.tan(np.pi / 4.0 + truelat1 * rad_per_deg / 2.0))
+        )
+
+    diff = lon - lov
+    diff = np.where(diff > 180.0, diff - 360.0, diff)
+    diff = np.where(diff < -180.0, diff + 360.0, diff)
+    alpha = diff * cone * rad_per_deg * np.sign(lat)
+    cos_alpha = np.cos(alpha)
+    sin_alpha = np.sin(alpha)
+
+    u_true = v * sin_alpha + u * cos_alpha
+    v_true = v * cos_alpha - u * sin_alpha
+    return u_true, v_true
+
+
+def add_hrrr_constraint_to_grid(
+    Grid, file_path, method="nearest", lov=-97.5, truelat1=38.5, truelat2=38.5
+):
     """
     This function will read an HRRR GRIB2 file and create the constraining
-    u, v, and w fields for the model constraint
+    u, v, and w fields for the model constraint.
+
+    The u and v winds in the HRRR GRIB2 file are relative to the HRRR's
+    Lambert Conformal Conic grid rather than true north. This function
+    rotates them to true north (earth-relative) components before
+    interpolating them onto the analysis grid.
 
     Parameters
     ----------
@@ -548,6 +611,16 @@ def add_hrrr_constraint_to_grid(Grid, file_path, method="nearest"):
     method: str
             Interpolation method: 'nearest' for nearest neighbor,
             'linear' for linear.
+    lov: float
+        The standard longitude (center longitude) of the HRRR's Lambert
+        Conformal Conic projection, in degrees. Only change this if the
+        input file uses a different Lambert Conformal Conic grid than the
+        operational CONUS HRRR.
+    truelat1, truelat2: float
+        The two standard (true) latitudes of the HRRR's Lambert Conformal
+        Conic projection, in degrees. Only change these if the input file
+        uses a different Lambert Conformal Conic grid than the
+        operational CONUS HRRR.
 
     Returns
     -------
@@ -577,6 +650,19 @@ def add_hrrr_constraint_to_grid(Grid, file_path, method="nearest"):
     lon = the_grib.variables["longitude"].data[:, :]
     lon[lon > 180] = lon[lon > 180] - 360
 
+    # HRRR's u and v are grid relative (relative to its Lambert Conformal
+    # Conic grid). Rotate them to true north before interpolating onto
+    # the analysis grid.
+    u_true, v_true = _uv_grid_relative_to_true_north(
+        grb_u.data[:, :, :],
+        grb_v.data[:, :, :],
+        lon[np.newaxis, :, :],
+        lat[np.newaxis, :, :],
+        lov=lov,
+        truelat1=truelat1,
+        truelat2=truelat2,
+    )
+
     # Convert geometric height to geopotential height
     EARTH_MEAN_RADIUS = 6.3781e6
     gh = gh.data[:, :, :]
@@ -610,7 +696,7 @@ def add_hrrr_constraint_to_grid(Grid, file_path, method="nearest"):
     lat_flattened = lat_flattened[the_box]
     height_flattened = height_flattened[the_box]
 
-    u_flattened = grb_u.data[:, :, :].flatten()
+    u_flattened = u_true.flatten()
     u_flattened = u_flattened[the_box]
     if method == "nearest":
         u_interp = NearestNDInterpolator(
@@ -627,7 +713,7 @@ def add_hrrr_constraint_to_grid(Grid, file_path, method="nearest"):
 
     u_new = u_interp(radar_grid_alt, radar_grid_lat, radar_grid_lon)
 
-    v_flattened = grb_v.data[:, :, :].flatten()
+    v_flattened = v_true.flatten()
     v_flattened = v_flattened[the_box]
     if method == "nearest":
         v_interp = NearestNDInterpolator(
@@ -668,7 +754,7 @@ def add_hrrr_constraint_to_grid(Grid, file_path, method="nearest"):
     )
 
     # Free up memory
-    del grb_u, grb_v, grb_w, lat, lon
+    del grb_u, grb_v, grb_w, lat, lon, u_true, v_true
     del the_grib
     gc.collect()
     return new_grid

pydda/tests/test_initialization.py

@@ -204,3 +204,53 @@ def test_hrrr_data():
     assert Grid["U_hrrr"].max() > 15
     assert Grid["V_hrrr"].max() > 15
     assert Grid["W_hrrr"].max() > 0
+
+
+def test_hrrr_uv_rotated_to_true_north():
+    # HRRR's u and v are relative to its Lambert Conformal Conic grid, not
+    # true north. add_hrrr_constraint_to_grid must rotate them before they
+    # are interpolated onto the analysis grid.
+    grid_shape = (4, 4, 4)
+    grid_limits = ((0, 5000.0), (-50000.0, 50000.0), (-50000.0, 50000.0))
+    file_path = pydda.tests.get_sample_file("ruc2anl_130_20110520_0800_001.grb2")
+
+    def make_grid(origin_lat, origin_lon):
+        Grid = pyart.testing.make_empty_grid(grid_shape, grid_limits)
+        for field in ("origin_latitude", "radar_latitude"):
+            getattr(Grid, field)["data"] = np.array([origin_lat])
+        for field in ("origin_longitude", "radar_longitude"):
+            getattr(Grid, field)["data"] = np.array([origin_lon])
+        Grid.init_point_longitude_latitude()
+        fdata3 = np.zeros(grid_shape)
+        Grid.add_field("zero_field", {"data": fdata3, "_FillValue": -9999.0})
+        return pydda.io.read_from_pyart_grid(Grid)
+
+    # Place the analysis domain well away from HRRR's -97.5 degree central
+    # meridian, where the grid-rotation correction has a large, easily
+    # measurable effect.
+    Grid = make_grid(38.5, -85.0)
+    Grid = pydda.constraints.add_hrrr_constraint_to_grid(Grid, file_path)
+    u_true_north = Grid["U_hrrr"].values
+    v_true_north = Grid["V_hrrr"].values
+
+    # Setting both true latitudes to the equator collapses the Lambert
+    # Conformal cone factor to zero, which makes the rotation a no-op.
+    # This gives an otherwise identical baseline of the raw, grid-relative
+    # winds to compare against.
+    Grid_grid_relative = make_grid(38.5, -85.0)
+    Grid_grid_relative = pydda.constraints.add_hrrr_constraint_to_grid(
+        Grid_grid_relative, file_path, truelat1=0.0, truelat2=0.0
+    )
+    u_grid_relative = Grid_grid_relative["U_hrrr"].values
+    v_grid_relative = Grid_grid_relative["V_hrrr"].values
+
+    # The rotation should meaningfully change the wind components this far
+    # from the central meridian...
+    assert np.abs(u_true_north - u_grid_relative).max() > 0.5
+    assert np.abs(v_true_north - v_grid_relative).max() > 0.5
+
+    # ...while preserving wind speed, since a coordinate rotation cannot
+    # change the magnitude of the wind vector.
+    speed_true_north = np.sqrt(u_true_north**2 + v_true_north**2)
+    speed_grid_relative = np.sqrt(u_grid_relative**2 + v_grid_relative**2)
+    np.testing.assert_allclose(speed_true_north, speed_grid_relative, rtol=1e-4)

setup.py

@@ -41,8 +41,9 @@
 LICENSE = "BSD"
 PLATFORMS = "Linux, Windows, OSX"
 MAJOR = 2
-MINOR = 4
-MICRO = 1
+MINOR = 5
+MICRO = 0
+
 
 # SCRIPTS = glob.glob('scripts/*')
 # TEST_SUITE = 'nose.collector'