From e02529b0447ef668edcb8841581cb060c91ec21d Mon Sep 17 00:00:00 2001 From: Robert Jackson Date: Fri, 31 Jul 2026 13:29:09 -0500 Subject: [PATCH 1/3] FIX: Rotate HRRR grid-relative winds to true north before use as constraint HRRR u/v are defined relative to its Lambert Conformal Conic grid, not true north. Apply the standard WRF/NCEP grid-rotation correction using the projection's cone factor before interpolating onto the analysis grid. Co-Authored-By: Claude Sonnet 5 --- pydda/constraints/model_data.py | 96 +++++++++++++++++++++++++++++++-- 1 file changed, 91 insertions(+), 5 deletions(-) diff --git a/pydda/constraints/model_data.py b/pydda/constraints/model_data.py index 8fc421a1..9d9c9b37 100644 --- a/pydda/constraints/model_data.py +++ b/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 From 45439c29878f274ada6175a5771a1a630485e5fa Mon Sep 17 00:00:00 2001 From: Robert Jackson Date: Fri, 31 Jul 2026 13:56:25 -0500 Subject: [PATCH 2/3] TST: Add unit test for HRRR grid-relative to true-north wind rotation Uses the bundled HRRR/RUC sample GRIB2 file to verify add_hrrr_constraint_to_grid's true-north correction meaningfully changes u/v away from the central meridian while conserving wind speed, as a rotation must. Co-Authored-By: Claude Sonnet 5 --- pydda/tests/test_initialization.py | 50 ++++++++++++++++++++++++++++++ 1 file changed, 50 insertions(+) diff --git a/pydda/tests/test_initialization.py b/pydda/tests/test_initialization.py index 6baf317d..ee6aa6f0 100644 --- a/pydda/tests/test_initialization.py +++ b/pydda/tests/test_initialization.py @@ -85,3 +85,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) From ef442d88295767b06d6ab78f052c2b415847d7fe Mon Sep 17 00:00:00 2001 From: Robert Jackson Date: Fri, 31 Jul 2026 14:00:13 -0500 Subject: [PATCH 3/3] Bump version to 2.5.0 Co-Authored-By: Claude Sonnet 5 --- pydda/__init__.py | 2 +- setup.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/pydda/__init__.py b/pydda/__init__.py index 1b6ca089..bba227dc 100644 --- a/pydda/__init__.py +++ b/pydda/__init__.py @@ -12,7 +12,7 @@ from . import constraints from . import io -__version__ = "2.3.0" +__version__ = "2.5.0" print("Welcome to PyDDA %s" % __version__) print("If you are using PyDDA in your publications, please cite:") diff --git a/setup.py b/setup.py index afa6d956..30278162 100644 --- a/setup.py +++ b/setup.py @@ -41,7 +41,7 @@ LICENSE = "BSD" PLATFORMS = "Linux, Windows, OSX" MAJOR = 2 -MINOR = 3 +MINOR = 5 MICRO = 0 # SCRIPTS = glob.glob('scripts/*')