From 57097531a64783af1fe50323c8b49eb119d54a26 Mon Sep 17 00:00:00 2001 From: mo-AliceLake <69578135+mo-AliceLake@users.noreply.github.com> Date: Thu, 30 Apr 2026 10:33:26 +0100 Subject: [PATCH] Updating wind downscaling code to be more descriptive (#2315) * Updating wind downscaling code to use more descriptive variable names and function names, etc * Adding name to CLA * Adding myself to .mailmap * Adding myself to .mailmap * Adding myself to .mailmap * Adding myself to .mailmap * Trying without space * Fixing unit tests * Updating error log message * Attempting to fix .mailmap by using anon github email * Changing u* -> u_star for clarity * Changing missing spatial coordinates to match old behaviour * Update cli * Pin proj and fix numpy-update-related test failures (#2326) * Pinning proj to 9.7.1 * Adding myself to contributors list * Fix numpy-update test failures * Reverting float() change * Making variable names clearer * Removing duplicate entry * Updating variable names to be more intuitive * Improving variable names and docstrings * Updating variable names in unit tests to match code --- .mailmap | 1 + CONTRIBUTING.md | 1 + envs/conda-forge.yml | 1 + envs/environment_a.yml | 1 + envs/environment_b.yml | 1 + envs/latest.yml | 1 + improver/cli/wind_downscaling.py | 23 +- .../generate_derived_solar_fields.py | 2 +- improver/utilities/spatial.py | 2 +- .../wind_calculations/wind_downscaling.py | 1468 ++++++++++------- .../emos_calibration/test_ApplyEMOS.py | 8 +- .../wind_downscaling/test_FrictionVelocity.py | 2 +- .../test_RoughnessCorrection.py | 316 ++-- 13 files changed, 1042 insertions(+), 785 deletions(-) diff --git a/.mailmap b/.mailmap index 238c5a250d..263025b7e3 100644 --- a/.mailmap +++ b/.mailmap @@ -1,4 +1,5 @@ Aaron Hopkinson +Alice Lake <69578135+mo-AliceLake@users.noreply.github.com> Andrew Creswick AndrewCreswick Anna Booton Anja Schubert <57921950+anja-bom@users.noreply.github.com> diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md index 2ee5f92789..f8f84c05b5 100644 --- a/CONTRIBUTING.md +++ b/CONTRIBUTING.md @@ -62,6 +62,7 @@ below: - Peter Jordan (Met Office, UK) - Anzer Khan (Met Office, UK) - Bruno P. Kinoshita (NIWA, NZ) + - Alice Lake (Met Office, UK) - Phoebe Lambert (Met Office, UK) - Lucy Liu (Bureau of Meteorology, Australia) - Daniel Mentiplay (Bureau of Meteorology, Australia) diff --git a/envs/conda-forge.yml b/envs/conda-forge.yml index c428222482..bfb9beb46b 100644 --- a/envs/conda-forge.yml +++ b/envs/conda-forge.yml @@ -18,6 +18,7 @@ dependencies: - scipy=1.15 - sphinx<9.0.0 # https://github.com/tox-dev/sphinx-autodoc-typehints/issues/586 - pandas<3 # https://github.com/SciTools/iris/issues/6761 + - proj=9.7.1 # https://github.com/metoppv/improver/issues/2325 # Additional libraries to run tests, not included in improver-feedstock - filelock - pytest diff --git a/envs/environment_a.yml b/envs/environment_a.yml index f9ba42c8bf..736daa5b87 100644 --- a/envs/environment_a.yml +++ b/envs/environment_a.yml @@ -10,6 +10,7 @@ dependencies: - iris=3.12 - numpy=2.2 - pandas<3 # https://github.com/SciTools/iris/issues/6761 + - proj=9.7.1 # https://github.com/metoppv/improver/issues/2325 - pytz - scipy=1.15 - sphinx<9.0.0 # https://github.com/tox-dev/sphinx-autodoc-typehints/issues/586 diff --git a/envs/environment_b.yml b/envs/environment_b.yml index c0ca54ff56..3cd5c69315 100644 --- a/envs/environment_b.yml +++ b/envs/environment_b.yml @@ -15,6 +15,7 @@ dependencies: - iris=3.12 - numpy=2.2 - pandas<3 # https://github.com/SciTools/iris/issues/6761 + - proj=9.7.1 # https://github.com/metoppv/improver/issues/2325 - pytz - scipy=1.15 - sphinx<9.0.0 # https://github.com/tox-dev/sphinx-autodoc-typehints/issues/586 diff --git a/envs/latest.yml b/envs/latest.yml index 0ff42cf670..dc6113369e 100644 --- a/envs/latest.yml +++ b/envs/latest.yml @@ -35,6 +35,7 @@ dependencies: - mock - mypy - pandas<3 # https://github.com/SciTools/iris/issues/6761 + - proj=9.7.1 # https://github.com/metoppv/improver/issues/2325 - pytest - pytest-cov - pytest-xdist diff --git a/improver/cli/wind_downscaling.py b/improver/cli/wind_downscaling.py index 6e2e17121d..ff5860251a 100755 --- a/improver/cli/wind_downscaling.py +++ b/improver/cli/wind_downscaling.py @@ -85,6 +85,9 @@ def process( "associated height level has been provided. These units " "will have no effect." ) + + # Attempt to iterate over the wind data by 'realization', and apply corrections to each member. + # If the cube has no realization coordinate (e.g. deterministic model data), treat the entire cube as a single member. try: wind_speed_iterator = wind_speed.slices_over("realization") except CoordinateNotFoundError: @@ -92,20 +95,25 @@ def process( wind_speed_list = iris.cube.CubeList() for wind_speed_slice in wind_speed_iterator: result = wind_downscaling.RoughnessCorrection( - silhouette_roughness, - sigma, - target_orography, - standard_orography, - model_resolution, - z0_cube=vegetative_roughness, + model_silhouette_roughness_cube=silhouette_roughness, + model_orog_stddev_cube=sigma, + target_orog_cube=target_orography, + model_orog_cube=standard_orography, + model_res=model_resolution, + model_z0_cube=vegetative_roughness, height_levels_cube=None, )(wind_speed_slice) wind_speed_list.append(result) - wind_speed = wind_speed_list.merge_cube() + + # Check whether 'realization' exists as a non-dimension coordinate. + # If so, reinsert it as a proper dimension axis so the cube has the expected shape. non_dim_coords = [x.name() for x in wind_speed.coords(dim_coords=False)] if "realization" in non_dim_coords: wind_speed = iris.util.new_axis(wind_speed, "realization") + + # If a specific output height is requested, use apply_extraction to select + # the corresponding height level from the processed wind cube. if output_height_level is not None: constraints = {"height": output_height_level} units = {"height": output_height_level_units} @@ -121,4 +129,5 @@ def process( ) ) wind_speed = single_level + return wind_speed diff --git a/improver/generate_ancillaries/generate_derived_solar_fields.py b/improver/generate_ancillaries/generate_derived_solar_fields.py index f05e486db7..94378da624 100644 --- a/improver/generate_ancillaries/generate_derived_solar_fields.py +++ b/improver/generate_ancillaries/generate_derived_solar_fields.py @@ -404,7 +404,7 @@ def _calc_clearsky_solar_radiation_data( # integrate the irradiance data along the time dimension to get the # accumulated solar irradiance. - solar_radiation_data = np.trapz( + solar_radiation_data = np.trapezoid( irradiance_data, dx=SECONDS_IN_MINUTE * temporal_spacing, axis=0 ) diff --git a/improver/utilities/spatial.py b/improver/utilities/spatial.py index 6a0ed017b9..35586ba72c 100644 --- a/improver/utilities/spatial.py +++ b/improver/utilities/spatial.py @@ -138,7 +138,7 @@ def distance_to_number_of_grid_cells( grid_cells = distance / abs(grid_spacing_metres) if return_int: - grid_cells = int(grid_cells) + grid_cells = int(np.asarray(grid_cells).item()) if grid_cells == 0: zero_distance_error = f"{d_error} gives zero cell extent" raise ValueError(zero_distance_error) diff --git a/improver/wind_calculations/wind_downscaling.py b/improver/wind_calculations/wind_downscaling.py index cd40afab9d..8c59bc9e4a 100644 --- a/improver/wind_calculations/wind_downscaling.py +++ b/improver/wind_calculations/wind_downscaling.py @@ -4,7 +4,6 @@ # See LICENSE in the root of the repository for full licensing details. """Module containing wind downscaling plugins.""" -import copy import itertools from typing import Optional, Tuple, Union @@ -18,597 +17,712 @@ from improver import BasePlugin, PostProcessingPlugin from improver.constants import RMDI -# Scale parameter to determine reference height +# Fractional tolerance used when deciding whether an absolute correction to +# the computed reference height is significant. Corrections smaller than this +# value are treated as noise and ignored. ABSOLUTE_CORRECTION_TOL = 0.04 -# Scaling parameter to determine reference height +# Multiplier used to convert terrain‑related variability (e.g. orographic +# standard deviation or silhouette roughness) into an effective reference +# height for the logarithmic wind‑profile calculation. HREF_SCALE = 2.0 -# Von Karman's constant +# Von Karman's constant, used in the logarithmic wind‑profile equation. VONKARMAN = 0.4 -# Default roughness length for sea points +# Default roughness length (in metres) assigned to sea grid cells when no +# surface roughness information is available. Z0M_SEA = 0.0001 class FrictionVelocity(BasePlugin): - """Class to calculate the friction velocity. - - This holds the function to calculate the friction velocity u_star, - given a reference height h_ref, the velocity at the reference - height u_href and the surface roughness z_0. - + """Compute friction velocity u_star using the logarithmic wind profile: + u_star = K * u(h_ref) / ln(h_ref / z0) + where: + - u(h_ref) is wind speed evaluated at the reference height h_ref, + - z0 is the aerodynamic roughness length, + - K is the von Karman constant. + + Notes: + - h_ref and model_z0 must share the same units. + - The returned u_star has the same velocity units as wspeed_at_h_ref. """ def __init__( - self, u_href: ndarray, h_ref: ndarray, z_0: ndarray, mask: ndarray + self, + wspeed_at_h_ref: ndarray, + h_ref: ndarray, + model_z0: ndarray, + ustar_mask: ndarray, ) -> None: - """Initialise the class. + """Initialize the friction-velocity calculator. Args: - u_href: - A 2D array of float32 for the wind speed at h_ref - h_ref: - A 2D array of float32 for the reference heights - z_0: - A 2D array of float32 for the vegetative roughness lengths - mask: - A 2D array of booleans where True indicates calculate u* + wspeed_at_h_ref: Wind speed evaluated at the reference height h_ref. + h_ref: Effective reference height h_ref. + model_z0: Vegetative roughness length z0. + ustar_mask: 2D boolean array, True where u_star should be computed. - Notes: - * z_0 and h_ref need to have identical units. - * the calculated friction velocity will have the units of the - supplied velocity u_href. + Raises: + ValueError: If any input array has a different size to the others. """ - self.u_href = u_href + self.wspeed_at_h_ref = wspeed_at_h_ref self.h_ref = h_ref - self.z_0 = z_0 - self.mask = mask - - # Check that input cubes are the same size - array_sizes = [np.size(u_href), np.size(h_ref), np.size(z_0), np.size(mask)] - - if not all(x == array_sizes[0] for x in array_sizes): - raise ValueError("Different size input arrays u_href, h_ref, z_0, mask") + self.model_z0 = model_z0 + self.ustar_mask = ustar_mask + + # Check array sizes all the same + sizes = [ + np.size(wspeed_at_h_ref), + np.size(h_ref), + np.size(model_z0), + np.size(ustar_mask), + ] + if not all(s == sizes[0] for s in sizes): + raise ValueError( + "Input arrays must have identical sizes, but sizes are: " + f"wspeed_at_h_ref={sizes[0]}, " + f"h_ref={sizes[1]}, " + f"model_z0={sizes[2]}, " + f"ustar_mask={sizes[3]}" + ) def process(self) -> ndarray: - """Function to calculate the friction velocity. + """Compute friction velocity (u_star) using the logarithmic wind profile. - ustar = K * (u_href / ln(h_ref / z_0)) + The friction velocity is computed according to the neutral-stability + logarithmic wind law: + u_star = K * u(h_ref) / ln(h_ref / z0) + where: + - u(h_ref) is wind speed evaluated at the reference height h_ref, + - z0 is the aerodynamic roughness length, + - K is the von Karman constant. - where ustar is the friction velocity, K is Von Karman's - constant, u_ref is the wind speed at the reference height, - h_ref is the reference height and z_0 is the vegetative - roughness length. + This method applies the calculation only at locations where + ustar_mask is True. All other points are filled with + missing-data indicators. Returns: - A 2D array of float32 friction velocities + 2D float32 array of friction velocity u_star. """ - ustar = np.full(self.u_href.shape, RMDI, dtype=np.float32) - numerator = self.u_href[self.mask] + ustar = np.full(self.wspeed_at_h_ref.shape, RMDI, dtype=np.float32) + + # Values at locations where u_star is to be computed + wind_vals = self.wspeed_at_h_ref[self.ustar_mask] + + # Compute log(h_ref / z0) with np.errstate(invalid="ignore"): - denominator = np.log(self.h_ref[self.mask] / self.z_0[self.mask]) - ustar[self.mask] = VONKARMAN * (numerator / denominator) + log_term = np.log( + self.h_ref[self.ustar_mask] / self.model_z0[self.ustar_mask] + ) + + # Compute u_star following the log-profile relation + ustar[self.ustar_mask] = VONKARMAN * (wind_vals / log_term) + return ustar class RoughnessCorrectionUtilities: - """Class to calculate the height and roughness wind corrections. - - This holds functions to calculate the roughness and height - corrections given the ancil files: - - * standard deviation of height in grid cell as sigma (model grid on - pp grid) - * Silhouette roughness as a_over_s (model grid on pp grid) - * vegetative roughness length z_0 (model grid on pp grid) - * post-processing grid orography pporo - * model grid orography interpolated on post-processing grid modoro - * height level 3D/ 1D grid - * windspeed 3D field on height level 3D grid (from above). + """Utilities for computing wind-speed corrections due to surface roughness + and orographic height differences. + + Provides methods to apply roughness and height adjustments to forecast data + using ancillary inputs: + + - model_silhouette_roughness (ndarray): + Dimensionless measure of sub-grid terrain steepness and associated drag. + - model_orog_stddev (ndarray): + Standard deviation of sub-grid orography height (m), describing terrain variability. + - model_z0 (ndarray): + Vegetative roughness length (m), representing surface drag from land cover. + - target_orog (ndarray): + High-resolution (target) orography (m) to which winds are downscaled. + - model_orog (ndarray): + Model orography (m), representing smoothed terrain. + - height_levels (ndarray): + Heights (m) corresponding to the wind field vertical coordinate. + - wind_field (ndarray): + 3D wind speed field defined on the height levels. + + Note that all input fields must be defined on the same grid as the wind field + (i.e. the target/post-processed grid). In particular, ancillary inputs derived on + the model grid (e.g. model_orog) should be regridded to this target grid before use. """ def __init__( self, - a_over_s: ndarray, - sigma: ndarray, - z_0: ndarray, - pporo: ndarray, - modoro: ndarray, - ppres: float, - modres: float, + model_silhouette_roughness: ndarray, + model_orog_stddev: ndarray, + model_z0: ndarray, + target_orog: ndarray, + model_orog: ndarray, + output_res: float, + model_res: float, ) -> None: - """Set up roughness and height correction. - - This sets up the parameters used for roughness and height - correction given the ancillary file inputs: + """Initialise roughness and height-correction parameters. Args: - a_over_s: - 2D array float32 - Silhouette roughness field, dimensionless - ancillary data, calculated according to Robinson, D. (2008) - - Ancillary file creation for the UM, Unified Model - Documentation Paper 73. - sigma: - 2D array float32 - Standard deviation field of height in the - grid cell, units of length - z_0: - 2D array float32 - Vegetative roughness height field, - units of length - pporo: - 2D array float32 - Post processing grid orography field - modoro: - 2D array float32 - Model orography field interpolated to post - processing grid - ppres: - Float - Grid cell length of post processing grid - modres: - Float - Grid cell length of model grid + model_silhouette_roughness: 2D array of dimensionless silhouette roughness, describing + sub-grid terrain steepness and associated drag. + + model_orog_stddev: 2D array of orographic standard deviation (m), representing + sub-grid terrain height variability. + + model_z0: 2D array of vegetative roughness length (m), controlling + near-surface wind drag from land cover. + + target_orog: 2D array of high-resolution (target) orography (m) to which + winds are downscaled. + + model_orog: 2D array of model orography (m), representing smoothed terrain. + + output_res: Horizontal resolution of the target grid (m). + + model_res: Horizontal resolution of the raw NWP model grid (m). """ - self.a_over_s = a_over_s - self.z_0 = z_0 - self.pporo = pporo - self.modoro = modoro - self.h_over_2 = self.sigma2hover2(sigma) # half peak to trough height - self.hcmask, self.rcmask = self._setmask() # HC mask, RC mask - if self.z_0 is not None: - self.z_0[z_0 <= 0] = Z0M_SEA - self.dx_min = ppres / 2.0 # scales smaller than this not resolved in pp - # the original code had hardcoded 500 - self.dx_max = 3.0 * modres # scales larger than this resolved in model - # the original code had hardcoded 4000 - self.wavenum = self._calc_wav() # k = 2*pi / L + self.model_silhouette_roughness = model_silhouette_roughness + self.model_z0 = model_z0 + self.target_orog = target_orog + self.model_orog = model_orog + + # Half peak‑to‑trough orographic height + self.h_half = self.model_orog_stddev_to_h_half(model_orog_stddev) + + # Height‑correction and roughness‑correction masks + self.hc_mask, self.rc_mask = self._setmask() + + # Replace non‑positive roughness values with the default sea roughness + if self.model_z0 is not None: + self.model_z0[self.model_z0 <= 0] = Z0M_SEA + + # Minimum resolvable scale on the post-processing grid + self.dx_min = output_res / 2.0 + + # Maximum unresolved scale on the model grid + self.dx_max = 3.0 * model_res + + # Wavenumber of terrain variability: k = 2π / L + self.wavenumber = self._calc_wavenumber() + + # Reference height used for roughness correction. self.h_ref = self._calc_h_ref() - self._refinemask() # HC mask needs to be updated for missing orography - self.h_at0 = self._delta_height() # pp orography - model orography + + # Update height correction mask for missing orography + self._refinemask() + + # Height difference between post-processing and model orography + self.h_at0 = self._delta_height() def _refinemask(self) -> None: - """Remask over RMDI and NaN orography. + """Refine the height-correction mask based on invalid orography values. - The mask for HC needs to be False where either of the - orographies (model or pp) has an invalid number. This cannot be - done before because the mask is used to calculate the - wavenumber which can and should be calculated for all points - where h_over_2 and a_over_s is a valid number. + The height-correction mask (hc_mask) must be set to False wherever either + the post-processing or model orography contains invalid values (e.g. + missing data indicators or NaNs). + + This cannot be done earlier because hc_mask is used when calculating + the wavenumber, and the wavenumber should be computed for all points + where both h_half and model_silhouette_roughness are valid (even if the + corresponding orography values are not). """ - self.hcmask[self.pporo == RMDI] = False - self.hcmask[self.modoro == RMDI] = False - self.hcmask[np.isnan(self.pporo)] = False - self.hcmask[np.isnan(self.modoro)] = False + self.hc_mask[ + np.equal(self.target_orog, RMDI) + | np.equal(self.model_orog, RMDI) + | np.isnan(self.target_orog) + | np.isnan(self.model_orog) + ] = False def _setmask(self) -> Tuple[ndarray, ndarray]: - """Create a ~land-sea mask. + """Create the height-correction (hc_mask) and roughness-correction + (rc_mask) masks. - Create a mask that is basically a land-sea mask: - Both, the standard deviation and the silouette roughness, are 0 - over the sea. A standard deviation of 0 results in a RMDI for - h_over_2. + The height-correction mask acts like a land-sea mask: both h_half and + model_silhouette_roughness are zero over the sea, and a standard deviation + of 0 results in a missing data indicator for h_half. - Returns: - - 2D array of booleans- True for land-points, - false for Sea (HC) - - 2D array of booleans- additionally False for - invalid z_0 (RC) + The roughness-correction mask begins as hc_mask but also excludes + points where the vegetative roughness length (model_z0) + is missing or non-positive. """ - hcmask = np.full(self.h_over_2.shape, True, dtype=bool) - hcmask[self.h_over_2 <= 0] = False - hcmask[self.a_over_s <= 0] = False - hcmask[np.isnan(self.h_over_2)] = False - hcmask[np.isnan(self.a_over_s)] = False - rcmask = np.copy(hcmask) - if self.z_0 is not None: - rcmask[self.z_0 <= 0] = False - rcmask[np.isnan(self.z_0)] = False - return hcmask, rcmask + # Height‑correction mask + hc_mask = np.full(self.h_half.shape, True, dtype=bool) + hc_mask[self.h_half <= 0] = False + hc_mask[self.model_silhouette_roughness <= 0] = False + hc_mask[np.isnan(self.h_half)] = False + hc_mask[np.isnan(self.model_silhouette_roughness)] = False + + # Roughness‑correction mask + rc_mask = np.copy(hc_mask) + if self.model_z0 is not None: + rc_mask[self.model_z0 <= 0] = False + rc_mask[np.isnan(self.model_z0)] = False + + return hc_mask, rc_mask @staticmethod - def sigma2hover2(sigma: ndarray) -> ndarray: - """Calculate the half peak-to-trough height. + def model_orog_stddev_to_h_half(model_orog_stddev: ndarray) -> ndarray: + """Convert orography standard deviation into half the peak-to-trough + height. The ancillary data used to estimate the peak to trough height contains the standard deviation of height in a cell. For sine-waves, this relates to the amplitude of the wave as: + amplitude = model_orog_stddev * sqrt(2) - Amplitude = sigma * sqrt(2) - - The amplitude would correspond to half the peak-to-trough - height (h_o_2). + This amplitude corresponds to half the peak-to-trough height (h_half). Args: - sigma: - 2D array of float32 - standard deviation of height in - grid cell. + model_orog_stddev: 2D float32 array containing the standard deviation of terrain + height within each grid cell (metres). Returns: - 2D array of float32 - half peak-to-trough height. - - Comments: - Points that had sigma = 0 (i.e. sea points) are set to - RMDI. + 2D float32 array of half the peak-to-trough height (h_half). + Points with zero or missing input are assigned the missing-data + indicator. """ - h_o_2 = np.full(sigma.shape, RMDI, dtype=np.float32) - h_o_2[sigma > 0] = sigma[sigma > 0] * np.sqrt(2.0) - return h_o_2 - - def _calc_wav(self) -> ndarray: - """Calculate wavenumber k of typical orographic lengthscale. - - Function to calculate wavenumber k of typical orographic - lengthscale L: + h_half = np.full(model_orog_stddev.shape, RMDI, dtype=np.float32) + valid = model_orog_stddev > 0 + h_half[valid] = model_orog_stddev[valid] * np.sqrt(2.0) + return h_half - .. math:: - :label: + def _calc_wavenumber(self) -> ndarray: + """Calculate the wavenumber k associated with the orographic length + scale. - k = 2 * \\pi / L + The orographic length scale L is estimated from the half + peak-to-trough height (h_half) and the silhouette-roughness field + (average of up-slopes per unit length over several cross-sections + through a grid cell) using the relationship + L = 2 * h_half / model_silhouette_roughness - L is approximated from half the peak-to-trough height h_over_2 - and the silhoutte roughness a_over_s (average of up-slopes per - unit length over several cross-sections through a grid cell) - as: + The corresponding wavenumber k is + k = 2π / L = (model_silhouette_roughness * π) / h_half - .. math:: - :label: + h_half is derived from the standard deviation of sub-grid terrain + height (model_orog_stddev) as in model_orog_stddev_to_h_half: + h_half = model_orog_stddev * sqrt(2). - L = 2 * \\rm{h\\_over\\_2} / \\rm{a\\_over\\_s} + Wavenumbers are then limited to the smallest and largest scales that + can be represented by the post-processing grid and the model grid. - a_over_s is dimensionless since it is the sum of up-slopes - measured in the same unit lengths as it is calculated over. + Grid points where h_half is zero or missing are given the missing-data + indicator for the wavenumber. - h_over_2 is calculated from the standard deviation of height in - a grid cell, sigma, as: - - .. math:: - :label: - - \\rm{h\\_over\\_2} = \\sqrt{2} * \\rm{sigma} - - which is based on the assumptions of sine waves, see - sigma2hover2. - - From eq. (1) and (2) it follows that: + Returns: + 2D float32 array of wavenumbers, in units of + inverse units of supplied h_half. + """ + wavenumber = np.full( + self.model_silhouette_roughness.shape, RMDI, dtype=np.float32 + ) - .. math:: - :label: + # Compute wavenumber k for valid height‑correction points + valid = self.hc_mask + wavenumber[valid] = ( + self.model_silhouette_roughness[valid] * np.pi + ) / self.h_half[valid] - k = 2*\\pi / (2 * \\rm{h\\_over\\_2} / \\rm{a\\_over\\_s)} - = \\rm{a\\_over\\_s} * \\pi / \\rm{h\\_over\\_2} + # Apply upper/lower bounds determined by smallest+largest resolvable scales + wavenumber[wavenumber > (np.pi / self.dx_min)] = np.pi / self.dx_min + wavenumber[self.h_half == 0] = RMDI + wavenumber[np.abs(wavenumber) < (np.pi / self.dx_max)] = np.pi / self.dx_max - Returns: - 2D array float32 - wavenumber in units of inverse units of - supplied h_over_2. - """ - wavn = np.full(self.a_over_s.shape, RMDI, dtype=np.float32) - wavn[self.hcmask] = (self.a_over_s[self.hcmask] * np.pi) / self.h_over_2[ - self.hcmask - ] - wavn[wavn > np.pi / self.dx_min] = np.pi / self.dx_min - wavn[self.h_over_2 == 0] = RMDI - wavn[abs(wavn) < np.pi / self.dx_max] = np.pi / self.dx_max - return wavn + return wavenumber def _calc_h_ref(self) -> ndarray: - """Calculate the reference height for roughness correction. + """Calculate the reference height h_ref for roughness correction. - The reference height below which the flow is in equilibrium - with the vegetative roughness is proportional to 1/wavenum - (Howard & Clark, 2007). + The reference height marks the height below which the flow is + considered to be in equilibrium with the vegetative roughness. + This height is proportional to 1 / wavenumber (Howard & Clark, 2007). Vosper (2009) and Clark (2009) argue that at the reference height, the perturbation should have decayed to a fraction - epsilon (ABSOLUTE_CORRECTION_TOL). The factor alpha - implements eq. 1.3 in Clark (2009): UK Climatology - Wind - Screening Tool. See also Vosper (2009) for a motivation. - For a freely available external reference, see the Virtual Met - Mast Version 1 Methodology and Verification paper under - www.thecrownestate.co.uk. + epsilon (ABSOLUTE_CORRECTION_TOL). + + The factor alpha implements eq. 1.3 in Clark (2009): UK Climatology + - Wind Screening Tool. See also Vosper (2009) for a motivation. It is + defined as alpha = -log(ABSOLUTE_CORRECTION_TOL) alpha is the log of scale parameter to determine reference height which is currently set to 0.04 (this corresponds to epsilon in both Vosper and Clark) Returns: - 2D array float32 - reference height for roughness correction + ndarray: 2D array of reference height h_ref for roughness + correction. + + References: + Howard T., Clark P. 2007. Correction and downscaling of NWP wind + speed forecasts. Meteorological Applications 14(2), 105-116. """ alpha = -np.log(ABSOLUTE_CORRECTION_TOL) - tunable_param = np.full(self.wavenum.shape, RMDI, dtype=np.float32) - h_ref = np.full(self.wavenum.shape, RMDI, dtype=np.float32) - tunable_param[self.hcmask] = alpha + np.log( - self.wavenum[self.hcmask] * self.h_over_2[self.hcmask] - ) - tunable_param[tunable_param > 1.0] = 1.0 - tunable_param[tunable_param < 0.0] = 0.0 - h_ref[self.hcmask] = tunable_param[self.hcmask] / self.wavenum[self.hcmask] - h_ref[h_ref < 1.0] = 1.0 - h_ref = np.minimum(h_ref, HREF_SCALE * self.h_over_2) - h_ref[h_ref < 1.0] = 1.0 - h_ref[~self.hcmask] = 0.0 + + tunable = np.full(self.wavenumber.shape, RMDI, dtype=np.float32) + h_ref = np.full(self.wavenumber.shape, RMDI, dtype=np.float32) + + # Compute tunable parameter for valid points + valid = self.hc_mask + tunable[valid] = alpha + np.log(self.wavenumber[valid] * self.h_half[valid]) + tunable = np.clip(tunable, 0.0, 1.0) + + # Compute reference height + h_ref[valid] = tunable[valid] / self.wavenumber[valid] + + # Enforce lower and upper bounds on h_ref + h_ref = np.maximum(h_ref, 1.0) + h_ref = np.minimum(h_ref, HREF_SCALE * self.h_half) + h_ref = np.maximum(h_ref, 1.0) + + # For points outside hc_mask, no roughness correction is applied + h_ref[~self.hc_mask] = 0.0 + return h_ref def calc_roughness_correction( - self, hgrid: ndarray, uold: ndarray, mask: ndarray + self, + height_above_orog: ndarray, + wspeed_original: ndarray, + rc_mask: ndarray, ) -> ndarray: - """Function to perform the roughness correction. + """Apply the roughness correction. Args: - hgrid: - 3D or 1D array float32 - height above orography - uold: - 3D array float32 - original velocities at hgrid. - mask: - 2D array of bools that is True for land-points, False for Sea - and False for invalid z_0. + height_above_orog: 3D or 1D float32 array giving height above orography. + wspeed_original: 3D float32 array containing the original wind speeds. + rc_mask: 2D boolean array. True where roughness correction is valid + (e.g. land points with valid vegetative roughness length), + False elsewhere. Returns: - 3D np.array float32 - Corrected wind speed on hgrid. Above - href, this is equal to uold. - - Comments: - Replace the windspeed profile below the reference height with one - that increases logarithmically with height, bound by the original - velocity uhref at the reference height h_ref and by a 0 velocity at - the vegetative roughness height z_0 + 3D float32 array of roughness-corrected windspeeds. + Above the reference height h_ref, values remain unchanged from + wspeed_original. """ - uhref = self._calc_u_at_h(uold, hgrid, self.h_ref, mask) - if hgrid.ndim == 1: - hgrid = hgrid[np.newaxis, np.newaxis, :] - ustar = FrictionVelocity(uhref, self.h_ref, self.z_0, mask)() - unew = np.copy(uold) - mhref = self.h_ref - mhref[~mask] = RMDI - cond = hgrid < self.h_ref[:, :, np.newaxis] + wspeed_new = np.copy(wspeed_original) + + # Windspeed at the reference height + u_at_href = self._interpolate_wspeed_to_height( + wspeed_original, height_above_orog, self.h_ref, rc_mask + ) + + # Friction velocity u_star + ustar = FrictionVelocity(u_at_href, self.h_ref, self.model_z0, rc_mask)() - # Create array of ones. - arr_ones = np.ones(unew.shape, dtype=np.float32) + # h_ref = 0 where roughness correction does not apply + h_ref = np.copy(self.h_ref) + h_ref[~rc_mask] = 0.0 - first_arg = (ustar[:, :, np.newaxis] * arr_ones)[cond] - sec_arg = np.log( - hgrid / (np.reshape(self.z_0, self.z_0.shape + (1,)) * arr_ones) - )[cond] + # Ensure broadcast correctly (expand 1D to 3D) + if height_above_orog.ndim == 1: + height_above_orog = height_above_orog[np.newaxis, np.newaxis, :] + ustar_3d = ustar[:, :, np.newaxis] * np.ones_like(height_above_orog) + z0_3d = self.model_z0[:, :, np.newaxis] * np.ones_like(height_above_orog) - unew[cond] = (first_arg * sec_arg) / VONKARMAN + # Apply the roughness correction below the reference height + below_href = height_above_orog < h_ref[:, :, np.newaxis] + log_term = np.log(height_above_orog / z0_3d)[below_href] + ustar_term = ustar_3d[below_href] + wspeed_new[below_href] = (ustar_term * log_term) / VONKARMAN - return unew + return wspeed_new - def _calc_u_at_h( + def _interpolate_wspeed_to_height( self, - u_in: ndarray, - h_in: ndarray, - hhere: ndarray, - mask: ndarray, - dolog: bool = False, + wspeed_in: ndarray, + height_levels_in: ndarray, + height_target: ndarray, + valid_mask: ndarray, + use_log_interpolation: bool = False, ) -> ndarray: - """Function to interpolate u_in on h_in at hhere. + """Interpolate wind speed from input height levels to a target height. Args: - u_in: - 3D array float32 - velocity on h_in layer, last dim is height - h_in: - 3D or 1D array float32 - height layer array - hhere: - 2D array float32 - height grid to interpolate at - mask: - 2D array of bools - mask the final result for uath - dolog: - if True, log interpolation, default False + wspeed_in: 3D float32 array of wind speed defined on height_levels_in. + Last dimension is height. + height_levels_in: 3D or 1D float32 array of heights corresponding to wspeed_in. + height_target: 2D float32 array giving the target height at which to interpolate. + valid_mask: 2D boolean array. True where interpolation is permitted. + use_log_interpolation: If True, perform logarithmic interpolation. Otherwise linear. Returns: - 2D array float32 - velocity interpolated at h + 2D float32 array of interpolated wind speed. """ - u_in = np.ma.masked_less(u_in, 0.0) - h_in = np.ma.masked_less(h_in, 0.0) - # h_in.mask = u_in.mask - # If I allow 1D height grids, I think I cannot do the hop over. - - # Ignores the height at the position where u_in is RMDI,"hops over" - hhere = np.ma.masked_less(hhere, 0.0) - upidx = np.argmax(h_in > hhere[:, :, np.newaxis], axis=2) - # loidx = np.maximum(upidx-1, 0) #if RMDI, need below - loidx = np.argmin( - np.ma.masked_less(hhere[:, :, np.newaxis] - h_in, 0.0), axis=2 + + # Mask invalid (negative) heights and speeds + wspeed_in = np.ma.masked_less(wspeed_in, 0.0) + height_levels_in = np.ma.masked_less(height_levels_in, 0.0) + height_target = np.ma.masked_less(height_target, 0.0) + + # Find indices of the first height above and below the target height + above_idx = np.argmax( + height_levels_in > height_target[:, :, np.newaxis], axis=2 ) - if h_in.ndim == 3: - hup = h_in.take( - upidx.flatten() - + np.arange(0, upidx.size * h_in.shape[2], h_in.shape[2]) + # Index of the height just below target + below_idx = np.argmin( + np.ma.masked_less(height_target[:, :, np.newaxis] - height_levels_in, 0.0), + axis=2, + ) + + # Extract bounding heights (upper and lower) + if height_levels_in.ndim == 3: + flat_stride = height_levels_in.shape[2] + h_upper = height_levels_in.take( + above_idx.flatten() + + np.arange(0, above_idx.size * flat_stride, flat_stride) ) - hlow = h_in.take( - loidx.flatten() - + np.arange(0, loidx.size * h_in.shape[2], h_in.shape[2]) + h_lower = height_levels_in.take( + below_idx.flatten() + + np.arange(0, below_idx.size * flat_stride, flat_stride) ) - elif h_in.ndim == 1: - hup = h_in[upidx].flatten() - hlow = h_in[loidx].flatten() - uup = u_in.take( - upidx.flatten() + np.arange(0, upidx.size * u_in.shape[2], u_in.shape[2]) + else: + h_upper = height_levels_in[above_idx].flatten() + h_lower = height_levels_in[below_idx].flatten() + + # Extract bounding wind‑speed values (upper and lower) + flat_stride_u = wspeed_in.shape[2] + u_upper = wspeed_in.take( + above_idx.flatten() + + np.arange(0, above_idx.size * flat_stride_u, flat_stride_u) ) - ulow = u_in.take( - loidx.flatten() + np.arange(0, loidx.size * u_in.shape[2], u_in.shape[2]) + u_lower = wspeed_in.take( + below_idx.flatten() + + np.arange(0, below_idx.size * flat_stride_u, flat_stride_u) ) - mask = mask.flatten() - uath = np.full(mask.shape, RMDI, dtype=np.float32) - if dolog: - uath[mask] = self._interpolate_log( - hup[mask], hlow[mask], hhere.flatten()[mask], uup[mask], ulow[mask] + + # Choose interpolation method + valid_mask_flat = valid_mask.flatten() + u_at_height = np.full(valid_mask_flat.shape, RMDI, dtype=np.float32) + if use_log_interpolation: + u_at_height[valid_mask_flat] = self._interpolate_log( + h_upper[valid_mask_flat], + h_lower[valid_mask_flat], + height_target.flatten()[valid_mask_flat], + u_upper[valid_mask_flat], + u_lower[valid_mask_flat], ) else: - uath[mask] = self._interpolate_1d( - hup[mask], hlow[mask], hhere.flatten()[mask], uup[mask], ulow[mask] + u_at_height[valid_mask_flat] = self._interpolate_1d( + h_upper[valid_mask_flat], + h_lower[valid_mask_flat], + height_target.flatten()[valid_mask_flat], + u_upper[valid_mask_flat], + u_lower[valid_mask_flat], ) - uath = np.reshape(uath, hhere.shape) - return uath + + # Reshape to 2D field + return np.reshape(u_at_height, height_target.shape) @staticmethod def _interpolate_1d( - xup: ndarray, xlow: ndarray, at_x: ndarray, yup: ndarray, ylow: ndarray + x_upper: ndarray, + x_lower: ndarray, + x_target: ndarray, + y_upper: ndarray, + y_lower: ndarray, ) -> ndarray: """Simple 1D linear interpolation for 2D grid inputs level. Args: - xup: - 2D array float32 - upper x-bins - xlow: - 2D array float32 - lower x-bins - at_x: - 2D array float32 - x values to interpolate y at - yup: - 2D array float32 - y(xup) - ylow: - 2D array float32 - y(xlow) + x_upper:Upper x-coordinates (e.g., upper heights). + x_lower:Lower x-coordinates (e.g., lower heights). + x_target:Target x-values to interpolate at. + y_upper:Values at x_upper. + y_lower:Values at x_lower. Returns: - 2D array float32 - y(at_x) assuming a lin function - between xlow and xup + Interpolated y-values at x_target. Missing-data indicator is + returned where interpolation cannot be performed. """ - interp = np.full(xup.shape, RMDI, dtype=np.float32) - diffs = xup - xlow - interp[diffs != 0] = ylow[diffs != 0] + ( - (at_x[diffs != 0] - xlow[diffs != 0]) - / diffs[diffs != 0] - * (yup[diffs != 0] - ylow[diffs != 0]) - ) - interp[diffs == 0] = at_x[diffs == 0] / xup[diffs == 0] * (yup[diffs == 0]) + interp = np.full(x_upper.shape, RMDI, dtype=np.float32) + diff = x_upper - x_lower + + # Standard linear interpolation when x_upper != x_lower + valid = diff != 0 + interp[valid] = y_lower[valid] + (x_target[valid] - x_lower[valid]) / diff[ + valid + ] * (y_upper[valid] - y_lower[valid]) + + # Fallback for x_upper == x_lower + collapse = ~valid + interp[collapse] = x_target[collapse] / x_upper[collapse] * y_upper[collapse] + return interp @staticmethod def _interpolate_log( - xup: ndarray, xlow: ndarray, at_x: ndarray, yup: ndarray, ylow: ndarray + x_upper: ndarray, + x_lower: ndarray, + x_target: ndarray, + y_upper: ndarray, + y_lower: ndarray, ) -> ndarray: """Simple 1D log interpolation y(x), except if lowest layer is ground level. Args: - xup: - 2D array float32 - upper x-bins - xlow: - 2D array float32 - lower x-bins - at_x: - 2D array float32 - x values to interpolate y at - yup: - 2D array float32 - y(xup) - ylow: - 2D array float32 -y(xlow) + x_upper: Upper x-coordinates (e.g., upper heights). + x_lower: Lower x-coordinates (e.g., lower heights). + x_target: Target x-values to interpolate at. + y_upper: Values at x_upper. + y_lower: Values at x_lower. Returns: - 2D array float32 - y(at_x) assuming a log function - between xlow and xup + Interpolated y-values at x_target. """ - ain = np.full(xup.shape, RMDI, dtype=np.float32) - loginterp = np.full(xup.shape, RMDI, dtype=np.float32) - mfrac = xup / xlow - mtest = (xup / xlow != 1) & (at_x != xup) - ain[mtest] = (yup[mtest] - ylow[mtest]) / np.log(mfrac[mtest]) - loginterp[mtest] = ain[mtest] * np.log(at_x[mtest] / xup[mtest]) + yup[mtest] - mtest = xup / xlow == 1 # below lowest layer, make lin interp - loginterp[mtest] = at_x[mtest] / xup[mtest] * (yup[mtest]) - mtest = at_x == xup # just use yup - loginterp[mtest] = yup[mtest] - return loginterp + out = np.full(x_upper.shape, RMDI, dtype=np.float32) + ratio = x_upper / x_lower + + # Case 1: x_upper != x_lower and x_target != x_upper + # Log interpolation + normal = (ratio != 1.0) & (x_target != x_upper) + a = np.full(x_upper.shape, RMDI, dtype=np.float32) + a[normal] = (y_upper[normal] - y_lower[normal]) / np.log(ratio[normal]) + out[normal] = ( + a[normal] * np.log(x_target[normal] / x_upper[normal]) + y_upper[normal] + ) + + # Case 2: x_upper == x_lower + # Collapse to linear scaling + collapse = ratio == 1.0 + out[collapse] = (x_target[collapse] / x_upper[collapse]) * y_upper[collapse] + + # Case 3: x_target == x_upper + same_level = x_target == x_upper + out[same_level] = y_upper[same_level] + + return out def _calc_height_corr( self, - u_a: ndarray, - heightg: ndarray, - mask: ndarray, + wspeed_outer: ndarray, + height_above_orog: ndarray, + valid_mask: ndarray, onemfrac: Union[float, ndarray], ) -> ndarray: - """Function to calculate the additive height correction. + """Calculate the additive height correction. Args: - u_a: - 2D array float32 - outer velocity, e.g. velocity at h_ref_orig - heightg: - 1D or 3D array float32 - heights above orography - mask: - 3D array of bools - Masks the hc_add result - onemfrac: - Currently, scalar = 1. But can be a function of position and - height, e.g. a 3D array (float32) + wspeed_outer: 2D float32 array of wind speed at the reference height. + height_above_orog: 1D or 3D float32 array of heights above orography. + valid_mask: 3D boolean array where the correction should be applied. + onemfrac (float or ndarray): Currently, scalar = 1. But can be a function + of position and height, e.g. a 3D array (float32) Returns: - 3D array float32 - additive height correction to wind speed + 3D float32 array of additive height correction. Comments: The height correction is a disturbance of the flow that decays exponentially with height. The larger the vertical - offset h_at0 (the higher the unresolved hill), the larger - is the disturbance. + offset (the higher the unresolved hill), the larger the + disturbance. The more smooth the disturbance (the larger the horizontal scale of the disturbance), the smaller the height correction (hence, a larger wavenumber results in a larger disturbance). - hc_add = exp(-height*wavenumber)*u(href)*h_at_0*wavenumber A final factor of 1 is assumed and omitted for the Bessel function term. """ - (xdim, ydim) = u_a.shape - if heightg.ndim == 1: - zdim = heightg.shape[0] - heightg = heightg[np.newaxis, np.newaxis, :] - elif heightg.ndim == 3: - zdim = heightg.shape[2] - ml2 = self.h_at0 * self.wavenum - expon = np.ones([xdim, ydim, zdim], dtype=np.float32) - mult = self.wavenum[:, :, np.newaxis] * heightg - expon[mult > 0.0001] = np.exp(-mult[mult > 0.0001]) - hc_add = expon * u_a[:, :, np.newaxis] * ml2[:, :, np.newaxis] * onemfrac - hc_add[~mask, :] = 0 + nx, ny = wspeed_outer.shape + + # Ensure heights are 3D for broadcasting + if height_above_orog.ndim == 1: + nz = height_above_orog.shape[0] + height_above_orog = height_above_orog[np.newaxis, np.newaxis, :] + else: + nz = height_above_orog.shape[2] + + # Amplitude term + amp = self.h_at0 * self.wavenumber + + # Exponential decay factor exp(-k * z) + decay = np.ones((nx, ny, nz), dtype=np.float32) + kz = self.wavenumber[:, :, np.newaxis] * height_above_orog + decay[kz > 1e-4] = np.exp(-kz[kz > 1e-4]) + + # Full additive height correction + hc_add = ( + decay * wspeed_outer[:, :, np.newaxis] * amp[:, :, np.newaxis] * onemfrac + ) + + # Zero correction where mask is False + hc_add[~valid_mask] = 0.0 + return hc_add def _delta_height(self) -> ndarray: - """Function to calculate pp-grid diff from model grid. - - Calculate the difference between pp-grid height and model - grid height. + """Calculate the difference between post-processing grid height and + model grid height. Returns: - 2D array float32 - height difference, ppgrid-model + 2D float32 array of height difference, defined as target_orog - model_orog. """ - delt_z = np.full(self.pporo.shape, RMDI, dtype=np.float32) - delt_z[self.hcmask] = self.pporo[self.hcmask] - self.modoro[self.hcmask] + delt_z = np.full(self.target_orog.shape, RMDI, dtype=np.float32) + valid = self.hc_mask + delt_z[valid] = self.target_orog[valid] - self.model_orog[valid] return delt_z - def do_rc_hc_all(self, hgrid: ndarray, uorig: ndarray) -> ndarray: - """Function to call HC and RC (height and roughness corrections). + def do_rc_hc_all( + self, height_above_orog: ndarray, wspeed_original: ndarray + ) -> ndarray: + """Apply roughness (RC) and height (HC) corrections to the wind field. Args: - hgrid: - 1D or 3D array float32 - height grid of wind input - uorig: - 3D array float32 - wind speed on these levels + height_above_orog: 1D or 3D float32 array of heights above local orography. + wspeed_original: 3D float32 array of wind speed defined on height_above_orog. Returns: - - RC corrected windspeed - - HC additional part - - Friedrich, M. M., 2016 - Wind Downscaling Program (Internal Met Office Report) + 3D float32 array of wind speed after applying RC and HC. """ - if hgrid.ndim == 3: - condition1 = (hgrid == RMDI).any(axis=2) - self.hcmask[condition1] = False - self.rcmask[condition1] = False - mask_rc = np.copy(self.rcmask) - mask_rc[(uorig == RMDI).any(axis=2)] = False - mask_hc = np.copy(self.hcmask) - mask_hc[(uorig == RMDI).any(axis=2)] = False - if self.z_0 is not None: - unew = self.calc_roughness_correction(hgrid, uorig, mask_rc) + # Remove RC/HC where height inputs contain missing values + if height_above_orog.ndim == 3: + missing_h = (height_above_orog == RMDI).any(axis=2) + self.hc_mask[missing_h] = False + self.rc_mask[missing_h] = False + + # Disable RC/HC wherever the vertical wind profile is missing + mask_rc = np.copy(self.rc_mask) + mask_hc = np.copy(self.hc_mask) + missing_w = (wspeed_original == RMDI).any(axis=2) + mask_rc[missing_w] = False + mask_hc[missing_w] = False + + # 1. Roughness correction + if self.model_z0 is not None: + wspeed_rc = self.calc_roughness_correction( + height_above_orog, wspeed_original, mask_rc + ) else: - unew = uorig - uhref_orig = self._calc_u_at_h(uorig, hgrid, 1.0 / self.wavenum, mask_hc) - mask_hc[uhref_orig <= 0] = False + wspeed_rc = wspeed_original + + # 2. Height correction + # Requires wind speed at the reference height, so interpolate first + uhref_orig = self._interpolate_wspeed_to_height( + wspeed_original, height_above_orog, 1.0 / self.wavenumber, mask_hc + ) + # HC only where u(h_ref) is positive + mask_hc[uhref_orig <= 0.0] = False # Setting this value to 1, is equivalent to setting the # Bessel function to 1. (Friedrich, 2016) # Example usage if the Bessel function was not set to 1 is: # onemfrac = 1.0 - BfuncFrac(nx,ny,nz,heightvec,z_0,waveno, Ustar, UI) onemfrac = 1.0 - hc_add = self._calc_height_corr(uhref_orig, hgrid, mask_hc, onemfrac) - result = unew + hc_add - result[result < 0.0] = 0 # HC can be negative if pporo None: - """Initialise the RoughnessCorrection instance. + """Initialise the RoughnessCorrection plugin. Args: - a_over_s_cube: - 2D - model silhouette roughness on pp grid. dimensionless - sigma_cube: - 2D - standard deviation of model orography height on pp grid. - In m. - pporo_cube: - 2D - pp orography. In m - modoro_cube: - 2D - model orography interpolated on pp grid. In m - modres: - original average model resolution in m + model_silhouette_roughness_cube: + 2D model silhouette roughness (dimensionless). Describes how steep + and rugged unresolved terrain is within a model grid box, and hence + the amount of drag and turbulence it introduces. + This is a static model ancillary field. + + model_orog_stddev_cube: + 2D standard deviation of model orography height (m). Represents the + vertical variability of unresolved terrain within a grid box (i.e. + how large the sub-grid hills and valleys are). + This is a static model ancillary field. + + target_orog_cube: + 2D high-resolution (true) orography (m) that winds are downscaled to. + + model_orog_cube: + 2D model orography (m), representing the smoothed terrain used by + the model. + This is a static model ancillary field. + + model_res: + Native horizontal resolution of the model orography (m), prior to + interpolation onto the standard grid. + + model_z0_cube: + 2D vegetative roughness length (m), representing drag from vegetation + and land cover. Controls the near-surface wind profile. + Historically static, but may now be time-varying (e.g. from StaGE). + height_levels_cube: - 3D or 1D - height of input velocity field. - Can be position dependent - z0_cube: - 2D - vegetative roughness length in m. If not given, do not do - any RC + 1D or 3D height levels of the input wind field (m). + + Notes: + All ancillary inputs must be defined on the same grid as the wind field + (the target / post-processed grid). Fields originating on the model grid + must be regridded prior to use. + + References: + Howard T., Clark P. 2007. Correction and downscaling of NWP wind + speed forecasts. Meteorological Applications 14(2), 105-116. """ - # Standard Python 'float' type is either single or double depending on - # system and there is no reliable method of finding which from the - # variable. So force to numpy.float32 by default. - modres = np.float32(modres) - - x_name, y_name, _, _ = self.find_coord_names(pporo_cube) - # Some checking that all the grids match - if not ( - self.check_ancils( - a_over_s_cube, sigma_cube, z0_cube, pporo_cube, modoro_cube - ) + + model_res = np.float32(model_res) + x_name, y_name, _, _ = self.find_coord_names(target_orog_cube) + + # Check grid consistency + if not self.check_ancils( + model_silhouette_roughness_cube, + model_orog_stddev_cube, + model_z0_cube, + target_orog_cube, + model_orog_cube, ): - raise ValueError("ancillary grids are not consistent") - # I want this ordering. Careful: iris.cube.Cube.slices is unreliable. - self.a_over_s = next(a_over_s_cube.slices([y_name, x_name])) - self.sigma = next(sigma_cube.slices([y_name, x_name])) + raise ValueError("Ancillary grids are not consistent.") + + # Extract 2D [y, x] slices + self.model_silhouette_roughness = next( + model_silhouette_roughness_cube.slices([y_name, x_name]) + ) + self.model_orog_stddev = next(model_orog_stddev_cube.slices([y_name, x_name])) + try: - self.z_0 = next(z0_cube.slices([y_name, x_name])) + self.model_z0 = next(model_z0_cube.slices([y_name, x_name])) except AttributeError: - self.z_0 = z0_cube - self.pp_oro = next(pporo_cube.slices([y_name, x_name])) - self.model_oro = next(modoro_cube.slices([y_name, x_name])) - self.ppres = self.calc_av_ppgrid_res(pporo_cube) - self.modres = modres + self.model_z0 = model_z0_cube + + self.target_orog = next(target_orog_cube.slices([y_name, x_name])) + self.model_orog = next(model_orog_cube.slices([y_name, x_name])) + + # Grid resolutions + self.output_res = self.calc_av_output_res(target_orog_cube) + self.model_res = model_res + + # Optional height levels self.height_levels = height_levels_cube - self.x_name = None - self.y_name = None + + # Store coordinate names + self.x_name = x_name + self.y_name = y_name self.z_name = None self.t_name = None @@ -682,94 +830,126 @@ def find_coord_names(self, cube: Cube) -> Tuple[str, str, str, str]: """Extract x, y, z, and time coordinate names. Args: - cube: - some iris cube to find coordinate names from + cube: Cube from which coordinate names will be extracted. Returns: - - name of the axis name in x-direction - - name of the axis name in y-direction - - name of the axis name in z-direction - - name of the axis name in t-direction + (x_name, y_name, z_name, t_name) """ - clist = {cube.coords()[i].name() for i in range(len(cube.coords()))} + coord_names = {coord.name() for coord in cube.coords()} + + # x coordinate try: - xname = cube.coord(axis="x").name() + x_name = cube.coord(axis="x").name() except CoordinateNotFoundError as exc: - print("'{0}' while xname setting. Args: {1}.".format(exc, exc.args)) + print(f"'{exc}' while determining x_name. Args: {exc.args}") + x_name = None + + # y coordinate try: - yname = cube.coord(axis="y").name() + y_name = cube.coord(axis="y").name() except CoordinateNotFoundError as exc: - print("'{0}' while yname setting. Args: {1}.".format(exc, exc.args)) - if clist.intersection(self.zcoordnames): - zname = list(clist.intersection(self.zcoordnames))[0] - else: - zname = None + print(f"'{exc}' while determining y_name. Args: {exc.args}") + y_name = None - if clist.intersection(self.tcoordnames): - tname = list(clist.intersection(self.tcoordnames))[0] - else: - tname = None - return xname, yname, zname, tname + # Check spatial coordinates exist + missing = [ + name for name, value in (("x", x_name), ("y", y_name)) if value is None + ] + if missing: + raise ValueError( + f"Cube is missing required spatial coordinate(s): {', '.join(missing)}" + ) - def calc_av_ppgrid_res(self, a_cube: Cube) -> float: - """Calculate average grid resolution from a cube. + # z coordinate + z_matches = coord_names.intersection(self.zcoordnames) + z_name = next(iter(z_matches), None) + + # time coordinate + t_matches = coord_names.intersection(self.tcoordnames) + t_name = next(iter(t_matches), None) + + return x_name, y_name, z_name, t_name + + def calc_av_output_res(self, input_cube: Cube) -> float: + """Calculate the average horizontal resolution of the given cube. Args: - a_cube: - Cube to calculate average resolution of + input_cube: Cube from which to determine the grid spacing. Returns: - Average grid resolution. + Average horizontal grid resolution (metres). """ - x_name, y_name, _, _ = self.find_coord_names(a_cube) - [exp_xname, exp_yname] = ["projection_x_coordinate", "projection_y_coordinate"] - exp_unit = Unit("m") - if (x_name != exp_xname) or (y_name != exp_yname): - raise ValueError("cannot currently calculate resolution") - - if a_cube.coord(x_name).bounds is None and a_cube.coord(y_name).bounds is None: - xres = (np.diff(a_cube.coord(x_name).points)).mean() - yres = (np.diff(a_cube.coord(y_name).points)).mean() + # Identify horizontal coordinate names + x_name, y_name, _, _ = self.find_coord_names(input_cube) + + # Expected coordinates and units + expected_x = "projection_x_coordinate" + expected_y = "projection_y_coordinate" + expected_units = Unit("m") + + if x_name != expected_x or y_name != expected_y: + raise ValueError("Cannot calculate resolution: unexpected horizontal axes.") + + x_coord = input_cube.coord(x_name) + y_coord = input_cube.coord(y_name) + + # Use bounds if available, else use point spacing + if x_coord.bounds is None and y_coord.bounds is None: + xres = np.diff(x_coord.points).mean() + yres = np.diff(y_coord.points).mean() else: - xres = (np.diff(a_cube.coord(x_name).bounds)).mean() - yres = (np.diff(a_cube.coord(y_name).bounds)).mean() - if (a_cube.coord(x_name).units != exp_unit) or ( - a_cube.coord(y_name).units != exp_unit - ): - raise ValueError("cube axis have units different from m.") + xres = np.diff(x_coord.bounds).mean() + yres = np.diff(y_coord.bounds).mean() + + # Ensure units are metres + if x_coord.units != expected_units or y_coord.units != expected_units: + raise ValueError("Post-processing grid axes must have units of metres.") + + # Mean absolute resolution return (abs(xres) + abs(yres)) / 2.0 @staticmethod def check_ancils( - a_over_s_cube: Cube, - sigma_cube: Cube, - z0_cube: Optional[Cube], - pp_oro_cube: Cube, - model_oro_cube: Cube, - ) -> ndarray: - """Check ancils grid and units. + model_silhouette_roughness_cube: Cube, + model_orog_stddev_cube: Cube, + model_z0_cube: Optional[Cube], + target_orog_cube: Cube, + model_orog_cube: Cube, + ) -> bool: + """ + Check ancillary inputs for grid consistency and expected units. - Check if ancil cubes are on the same grid and if they have the - expected units. The testing for "same grid" might be replaced - if there is a general utils function made for it or so. + Ensures all ancillary cubes are on the same spatial grid and have + appropriate units for wind downscaling. Args: - a_over_s_cube: - holding the silhouette roughness field - sigma_cube: - holding the standard deviation of height in a grid cell - z0_cube: - holding the vegetative roughness field - pp_oro_cube: - holding the post processing grid orography - model_oro_cube: - holding the model orography on post processing grid + model_silhouette_roughness_cube: + Dimensionless field describing sub-grid terrain ruggedness. + model_orog_stddev_cube: + Standard deviation of sub-grid orography height (m). + model_z0_cube: + Vegetative roughness length (m), representing surface drag + from land cover. + target_orog_cube: + High-resolution (target) orography (m) used for downscaling. + model_orog_cube: + Model orography (m) representing the smoothed terrain. Returns: - Containing bools describing whether or not the tests passed + bool: + True if all ancillary fields share the same x/y grid; + False otherwise. """ - ancil_list = [a_over_s_cube, sigma_cube, pp_oro_cube, model_oro_cube] - unwanted_coord_list = [ + required = [ + model_silhouette_roughness_cube, + model_orog_stddev_cube, + target_orog_cube, + model_orog_cube, + ] + required_units = [1, Unit("m"), Unit("m"), Unit("m")] + + # Coords to strip before comparing horizontal grids + drop_coords = [ "time", "height", "model_level_number", @@ -777,174 +957,192 @@ def check_ancils( "forecast_reference_time", "forecast_period", ] - for field, exp_unit in zip(ancil_list, [1, Unit("m"), Unit("m"), Unit("m")]): - for unwanted_coord in unwanted_coord_list: + + # Clean and check required cubes + for field, expected in zip(required, required_units): + for name in drop_coords: try: - field.remove_coord(unwanted_coord) + field.remove_coord(name) except CoordinateNotFoundError: pass - if field.units != exp_unit: - msg = ( - "{} ancil field has unexpected unit:" - " {} (expected) vs. {} (actual)" + if field.units != expected: + raise ValueError( + f"{field.name()} ancillary has unexpected unit: " + f"expected {expected}, got {field.units}" ) - raise ValueError(msg.format(field.name(), exp_unit, field.units)) - if z0_cube is not None: - ancil_list.append(z0_cube) - for unwanted_coord in unwanted_coord_list: + + # Clean and check z0 if supplied + if model_z0_cube is not None: + for name in drop_coords: try: - z0_cube.remove_coord(unwanted_coord) + model_z0_cube.remove_coord(name) except CoordinateNotFoundError: pass - if z0_cube.units != Unit("m"): - msg = "z0 ancil has unexpected unit: should be {} is {}" - raise ValueError(msg.format(Unit("m"), z0_cube.units)) - permutated_ancil_list = list(itertools.permutations(ancil_list, 2)) - oklist = [] - for entry in permutated_ancil_list: - x_axis_flag = entry[0].coord(axis="y") == entry[1].coord(axis="y") - y_axis_flag = entry[0].coord(axis="x") == entry[1].coord(axis="x") - oklist.append(x_axis_flag & y_axis_flag) - # HybridHeightToPhenomOnPressure._cube_compatibility_check(entr[0], - # entr[1]) - return np.array(oklist).all() # replace by a return value of True + if model_z0_cube.units != Unit("m"): + raise ValueError( + f"z0 ancillary has unexpected unit: " + f"expected m, got {model_z0_cube.units}" + ) + required.append(model_z0_cube) + + # Pairwise x/y grid compatibility check across all ancils + ok_pairs: list[bool] = [] + for a, b in itertools.permutations(required, 2): + try: + same_y = a.coord(axis="y") == b.coord(axis="y") + same_x = a.coord(axis="x") == b.coord(axis="x") + ok_pairs.append(bool(same_x & same_y)) + except CoordinateNotFoundError: + ok_pairs.append(False) + + return all(ok_pairs) def find_coord_order(self, mcube: Cube) -> Tuple[int, int, int, int]: - """Extract coordinate ordering within a cube. + """Return the dimension indices of the x, y, z, and time coordinates. Use coord_dims to assess the dimension associated with a particular dimension coordinate. If a coordinate is not a dimension coordinate, then a NaN value will be returned for that coordinate. - Args: - mcube: - cube to check the order of coordinate axis - Returns: - - position of x axis. - - position of y axis. - - position of z axis. - - position of t axis. + (x_dim, y_dim, z_dim, t_dim) with NaN for coordinates not found. """ coord_names = [self.x_name, self.y_name, self.z_name, self.t_name] - positions = len(coord_names) * [np.nan] - for coord_index, coord_name in enumerate(coord_names): - if mcube.coords(coord_name, dim_coords=True): - (positions[coord_index],) = mcube.coord_dims(coord_name) - return positions + positions = [np.nan, np.nan, np.nan, np.nan] + + for idx, coord_name in enumerate(coord_names): + # Skip missing coord names + if coord_name is None: + continue + # Record coordinate axis number + try: + if mcube.coords(coord_name, dim_coords=True): + (positions[idx],) = mcube.coord_dims(coord_name) + except CoordinateNotFoundError: + # Coordinate does not exist, so leave as NaN + pass + + return tuple(positions) def find_heightgrid(self, wind: Cube) -> ndarray: - """Setup the height grid. + """Find the height grid to use for interpolation. - Setup the height grid either from the 1D or 3D height grid - that was supplied to the plugin or from the z-axis information - from the wind grid. + If no height-levels cube is supplied, use the vertical coordinate + from the wind cube. Otherwise use the provided height-levels cube. Args: - wind: - 3D or 4D - representing the wind data. + wind: 3D or 4D wind-speed cube. Returns: - 1D or 3D array - representing the height grid. + 1D or 3D array of heights (metres). """ + # Case 1: No external height-levels cube provided + # -> use wind cube's z‑axis if self.height_levels is None: - hld = wind.coord(self.z_name).points + return wind.coord(self.z_name).points + + # Case 2: Use the height-levels cube + hld = iris.util.squeeze(self.height_levels) + if np.isnan(hld.data).any() or (hld.data == RMDI).any(): + raise ValueError("Height grid contains invalid points.") + if hld.ndim == 3: + try: + x_dim, y_dim, z_dim, _ = self.find_coord_order(hld) + hld = hld.transpose([y_dim, x_dim, z_dim]) + except Exception: + raise ValueError("Height grid does not align with wind grid.") + elif hld.ndim == 1: + try: + hld = next(hld.slices([self.z_name])) + except CoordinateNotFoundError: + raise ValueError("Height grid z‑coordinate differs from wind grid.") else: - hld = iris.util.squeeze(self.height_levels) - if np.isnan(hld.data).any() or (hld.data == RMDI).any(): - raise ValueError("height grid contains invalid points") - if hld.ndim == 3: - try: - xap, yap, zap, _ = self.find_coord_order(hld) - hld.transpose([yap, xap, zap]) - except KeyError: - raise ValueError("height grid different from wind grid") - elif hld.ndim == 1: - try: - hld = next(hld.slices([self.z_name])) - except CoordinateNotFoundError: - raise ValueError("height z coordinate differs from wind z") - else: - raise ValueError( - "hld must have a dimension length of " - "either 3 or 1" - "hld.ndim is {}".format(hld.ndim) - ) - hld = hld.data - return hld + raise ValueError(f"Height grid must be 1D or 3D, got ndim = {hld.ndim}.") - def check_wind_ancil(self, xwp: int, ywp: int) -> None: - """Check wind vs ancillary file grids. + return hld.data - Check if wind and ancillary files are on the same grid and if - they have the same ordering. + def check_wind_ancil(self, xwp: int, ywp: int) -> None: + """Verify that the wind field and ancillary grids share the same + horizontal orientation. Args: - xwp: - representing the position of the x-axis in the wind cube - ywp: - representing the position of the y-axis of the wind cube + xwp: Dimension index of the x-axis in the wind cube. + ywp: Dimension index of the y-axis in the wind cube. + + Raises: + ValueError: If ancillary grids do not share the same x/y dimension + ordering as the wind cube. """ - xap, yap, _, _ = self.find_coord_order(self.pp_oro) + # Dim-order of ancillary post-processing-grid orography + xap, yap, _, _ = self.find_coord_order(self.target_orog) + + # Compare relative ordering of (x,y) dimensions if xwp - ywp != xap - yap: if np.isnan(xap) or np.isnan(yap): - raise ValueError("ancillary grid different from wind grid") - raise ValueError("xy-orientation: ancillary differ from wind") + raise ValueError("Ancillary grid differs from wind grid.") + raise ValueError( + "XY dimension ordering differs between wind and ancillary grids." + ) def process(self, input_cube: Cube) -> Cube: - """Adjust the 4d wind field - cube - (x, y, z including times). + """Apply roughness (RC) and height (HC) corrections to a 4D wind cube. Args: - input_cube: - The wind cube to be operated upon. Should be wind speed on - height_levels for all desired forecast times. + input_cube: Wind-speed cube (x, y, z, time), defined on height levels for all + desired forecast times. Returns: - The 4d wind field with roughness and height correction - applied in the same order as the input cube. + The wind cube with RC and HC applied. Raises: - TypeError: If input_cube is not a cube. + TypeError: If input_cube is not an iris Cube. + ValueError: If any time slice contains invalid wind data. """ if not isinstance(input_cube, iris.cube.Cube): - msg = "wind input is not a cube, but {}" - raise TypeError(msg.format(type(input_cube))) - (self.x_name, self.y_name, self.z_name, self.t_name) = self.find_coord_names( + raise TypeError(f"Wind input is not a cube, but {type(input_cube)}") + + # Determine coordinate names and dimension ordering for the wind cube + self.x_name, self.y_name, self.z_name, self.t_name = self.find_coord_names( input_cube ) xwp, ywp, zwp, twp = self.find_coord_order(input_cube) + + # Reorder wind cube so dimensions are consistently (y, x, z [, t]) if np.isnan(twp): input_cube.transpose([ywp, xwp, zwp]) else: - input_cube.transpose([ywp, xwp, zwp, twp]) # problems with slices - rchc_list = iris.cube.CubeList() - if self.z_0 is None: - z0_data = None - else: - z0_data = self.z_0.data - roughness_correction = RoughnessCorrectionUtilities( - self.a_over_s.data, - self.sigma.data, + input_cube.transpose([ywp, xwp, zwp, twp]) + + z0_data = None if self.model_z0 is None else self.model_z0.data + rc_utils = RoughnessCorrectionUtilities( + self.model_silhouette_roughness.data, + self.model_orog_stddev.data, z0_data, - self.pp_oro.data, - self.model_oro.data, - self.ppres, - self.modres, + self.target_orog.data, + self.model_orog.data, + self.output_res, + self.model_res, ) self.check_wind_ancil(xwp, ywp) - hld = self.find_heightgrid(input_cube) + height_grid = self.find_heightgrid(input_cube) + + corrected_list = iris.cube.CubeList() for time_slice in input_cube.slices_over("time"): + # Validate wind field (e.g. not contain NaNs or negative values) if np.isnan(time_slice.data).any() or (time_slice.data < 0.0).any(): - msg = "{} has invalid wind data" - raise ValueError(msg.format(time_slice.coord(self.t_name))) - rc_hc = copy.deepcopy(time_slice) - rc_hc.data = roughness_correction.do_rc_hc_all(hld, time_slice.data) - rchc_list.append(rc_hc) - output_cube = rchc_list.merge_cube() - # reorder input_cube and output_cube as original + tcoord = time_slice.coord(self.t_name) + raise ValueError(f"{tcoord} has invalid wind data") + # Compute RC + HC result + rc_hc_cube = time_slice.copy() + rc_hc_cube.data = rc_utils.do_rc_hc_all(height_grid, time_slice.data) + corrected_list.append(rc_hc_cube) + output_cube = corrected_list.merge_cube() + + # Restore the original dimension ordering of both input and output if np.isnan(twp): - input_cube.transpose(np.argsort([ywp, xwp, zwp])) - output_cube.transpose(np.argsort([ywp, xwp, zwp])) + order = np.argsort([ywp, xwp, zwp]) + input_cube.transpose(order) + output_cube.transpose(order) else: input_cube.transpose(np.argsort([ywp, xwp, zwp, twp])) output_cube.transpose(np.argsort([twp, ywp, xwp, zwp])) diff --git a/improver_tests/calibration/emos_calibration/test_ApplyEMOS.py b/improver_tests/calibration/emos_calibration/test_ApplyEMOS.py index 2696bd6765..689f0cd161 100644 --- a/improver_tests/calibration/emos_calibration/test_ApplyEMOS.py +++ b/improver_tests/calibration/emos_calibration/test_ApplyEMOS.py @@ -476,8 +476,8 @@ def test_period_percentiles(self): """Test that cell methods are preserved on the calibrated forecast, if present on the input forecast.""" self.percentiles.coord("time").bounds = [ - int(self.percentiles.coord("time").points - 3600), - int(self.percentiles.coord("time").points), + int((self.percentiles.coord("time").points - 3600).item()), + int(self.percentiles.coord("time").points.item()), ] cell_methods = CellMethod("maximum", coords="time") @@ -499,8 +499,8 @@ def test_period_probabilities(self): without duplication of the cell methods, if cell methods are present on the input probability forecast.""" self.probabilities.coord("time").bounds = [ - int(self.probabilities.coord("time").points - 3600), - int(self.probabilities.coord("time").points), + int((self.probabilities.coord("time").points - 3600).item()), + int(self.probabilities.coord("time").points.item()), ] cell_methods = CellMethod("maximum", coords="time") diff --git a/improver_tests/wind_calculations/wind_downscaling/test_FrictionVelocity.py b/improver_tests/wind_calculations/wind_downscaling/test_FrictionVelocity.py index cceda6a58b..2c21b99ce3 100644 --- a/improver_tests/wind_calculations/wind_downscaling/test_FrictionVelocity.py +++ b/improver_tests/wind_calculations/wind_downscaling/test_FrictionVelocity.py @@ -100,7 +100,7 @@ def test_handles_zero_values(self): def test_handles_different_sized_arrays(self): """Test when if different size arrays have been input""" u_href = np.full([3, 3], 10, dtype=float) - msg = "Different size input arrays u_href, h_ref, z_0, mask" + msg = f"Input arrays must have identical sizes, but sizes are: wspeed_at_h_ref={np.size(u_href)}, h_ref={np.size(self.h_ref)}, model_z0={np.size(self.z_0)}, ustar_mask={np.size(self.mask)}" with self.assertRaisesRegex(ValueError, msg): FrictionVelocity(u_href, self.h_ref, self.z_0, self.mask).process() diff --git a/improver_tests/wind_calculations/wind_downscaling/test_RoughnessCorrection.py b/improver_tests/wind_calculations/wind_downscaling/test_RoughnessCorrection.py index a6a473c025..3e87fead18 100644 --- a/improver_tests/wind_calculations/wind_downscaling/test_RoughnessCorrection.py +++ b/improver_tests/wind_calculations/wind_downscaling/test_RoughnessCorrection.py @@ -86,61 +86,79 @@ class TestMultiPoint: """Test (typically) 3x1 or 3x3 point tests. It constructs cubes for the ancillary fields: - Silhouette roughness (AoS), standard deviation of model height grid - cell (Sigma), vegetative roughness (z_0), post-processing grid - orography (pporog) and model orography (modelorog). If no values + Silhouette roughness (silhouette_roughness), standard deviation of model height grid + cell (model_orog_stddev), vegetative roughness (model_z0), post-processing grid + orography (target_orog) and model orography (model_orog). If no values are supplied, the grids that are set up have equal values at all - x-y points: AoS = 0.2, Sigma = 20, z_0 = 0.2, pporog = 250, - modelorog = 230. + x-y points: silhouette_roughness = 0.2, model_orog_stddev = 20, model_z0 = 0.2, target_orog = 250, + model_orog = 230. """ def __init__( - self, shape=(3, 3), AoS=None, Sigma=None, z_0=0.2, pporog=None, modelorog=None + self, + shape=(3, 3), + silhouette_roughness=None, + model_orog_stddev=None, + model_z0=0.2, + target_orog=None, + model_orog=None, ): """Set up multi-point tests. Args: shape (tuple): Required data shape. - AoS (float or 1D or 2D numpy.ndarray): + silhouette_roughness (float or 1D or 2D numpy.ndarray): Silhouette roughness field - Sigma (float or 1D or 2D numpy.ndarray): + model_orog_stddev (float or 1D or 2D numpy.ndarray): Standard deviation field of height in grid cell - z_0 (float or 1D or 2D numpy.ndarray): + model_z0 (float or 1D or 2D numpy.ndarray): Vegetative roughness field - pporog (float or 1D or 2D numpy.ndarray): + target_orog (float or 1D or 2D numpy.ndarray): Unsmoothed orography field on post-processing grid - modelorog (float or 1D or 2D numpy.ndarray): + model_orog (float or 1D or 2D numpy.ndarray): Model orography field on post-processing grid """ self.shape = shape - if AoS is None: - AoS = np.full(shape, 0.2, dtype=np.float32) - if Sigma is None: - Sigma = np.full(shape, 20.0, dtype=np.float32) - if pporog is None: - pporog = np.full(shape, 250.0, dtype=np.float32) - if modelorog is None: - modelorog = np.full(shape, 230.0, dtype=np.float32) + if silhouette_roughness is None: + silhouette_roughness = np.full(shape, 0.2, dtype=np.float32) + if model_orog_stddev is None: + model_orog_stddev = np.full(shape, 20.0, dtype=np.float32) + if target_orog is None: + target_orog = np.full(shape, 250.0, dtype=np.float32) + if model_orog is None: + model_orog = np.full(shape, 230.0, dtype=np.float32) self.w_cube = None - self.aos_cube = make_ancil_cube(AoS, "silhouette_roughness", 1, shape=shape) + self.model_silhouette_roughness_cube = make_ancil_cube( + silhouette_roughness, "silhouette_roughness", 1, shape=shape + ) self.s_cube = make_ancil_cube( - Sigma, "standard_deviation_of_height_in_grid_cell", "m", shape=shape - ) - if z_0 is None: - self.z0_cube = None - elif isinstance(z_0, float): - z_0 = np.full(shape, z_0, dtype=np.float32) - self.z0_cube = make_ancil_cube(z_0, "vegetative_roughness_length", "m") - elif isinstance(z_0, list): - self.z0_cube = make_ancil_cube( - np.array(z_0), "vegetative_roughness_length", "m", shape=shape + model_orog_stddev, + "standard_deviation_of_height_in_grid_cell", + "m", + shape=shape, + ) + if model_z0 is None: + self.model_z0_cube = None + elif isinstance(model_z0, float): + model_z0 = np.full(shape, model_z0, dtype=np.float32) + self.model_z0_cube = make_ancil_cube( + model_z0, "vegetative_roughness_length", "m" ) - self.poro_cube = make_ancil_cube(pporog, "surface_altitude", "m", shape=shape) - self.moro_cube = make_ancil_cube( - modelorog, "surface_altitude", "m", shape=shape + elif isinstance(model_z0, list): + self.model_z0_cube = make_ancil_cube( + np.array(model_z0), + "vegetative_roughness_length", + "m", + shape=shape, + ) + self.target_orog_cube = make_ancil_cube( + target_orog, "surface_altitude", "m", shape=shape + ) + self.model_orog_cube = make_ancil_cube( + model_orog, "surface_altitude", "m", shape=shape ) def run_hc_rc(self, wind, dtime=1, height=None, aslist=False): @@ -193,12 +211,12 @@ def run_hc_rc(self, wind, dtime=1, height=None, aslist=False): ) plugin = RoughnessCorrection( - self.aos_cube, + self.model_silhouette_roughness_cube, self.s_cube, - self.poro_cube, - self.moro_cube, + self.target_orog_cube, + self.model_orog_cube, 1500.0, - self.z0_cube, + self.model_z0_cube, ) return plugin(self.w_cube) @@ -208,11 +226,11 @@ class TestSinglePoint: A cube is a single 1x x 1y grid, however, the z dimension is not 1. It constructs 1x1 cubes for the ancillary fields Silhouette - roughness (AoS) and standard deviation of model height grid cell - (Sigma), vegetative roughness (z_0), post-processing grid orography - (pporog) and model orography(modelorog). If no values are supplied, - the values are: AoS = 0.2, Sigma = 20, z_0 = 0.2, pporog = 250, - modelorog = 230. + roughness (silhouette_roughness) and standard deviation of model height grid cell + (model_orog_stddev), vegetative roughness (model_z0), post-processing grid orography + (target_orog) and model orography(model_orog). If no values are supplied, + the values are: silhouette_roughness = 0.2, model_orog_stddev = 20, model_z0 = 0.2, target_orog = 250, + model_orog = 230. The height level grid (heightlevels) can be supplied as an 1D array. If nothing is supplied, the height level grid is [0.2, 3, @@ -222,44 +240,51 @@ class TestSinglePoint: def __init__( self, - AoS=0.2, - Sigma=20.0, - z_0=0.2, - pporog=250.0, - modelorog=230.0, + silhouette_roughness=0.2, + model_orog_stddev=20.0, + model_z0=0.2, + target_orog=250.0, + model_orog=230.0, heightlevels=np.array([0.2, 3.0, 13.0, 33.0, 133.0, 333.0, 1133.0]), ): """Set up the single point test for RoughnessCorrection. Args: - AoS (float): + silhouette_roughness (float): Silhouette roughness field - Sigma (float): + model_orog_stddev (float): Standard deviation field of height in grid cell - z_0 (float): + model_z0 (float): Vegetative roughness field - pporog (float): + target_orog (float): Unsmoothed orography on post-processing grid - modelorog (float): + model_orog (float): Model orography on post-processing grid heightlevels (1D numpy.ndarray): Height level array """ self.w_cube = None - self.aos_cube = make_ancil_cube(AoS, "silhouette_roughness", 1, shape=(1, 1)) + self.model_silhouette_roughness_cube = make_ancil_cube( + silhouette_roughness, "silhouette_roughness", 1, shape=(1, 1) + ) self.s_cube = make_ancil_cube( - Sigma, "standard_deviation_of_height_in_grid_cell", "m", shape=(1, 1) + model_orog_stddev, + "standard_deviation_of_height_in_grid_cell", + "m", + shape=(1, 1), ) - if z_0 is None: - self.z0_cube = None + if model_z0 is None: + self.model_z0_cube = None else: - self.z0_cube = make_ancil_cube( - z_0, "vegetative_roughness_length", "m", shape=(1, 1) + self.model_z0_cube = make_ancil_cube( + model_z0, "vegetative_roughness_length", "m", shape=(1, 1) ) - self.poro_cube = make_ancil_cube(pporog, "surface_altitude", "m", shape=(1, 1)) - self.moro_cube = make_ancil_cube( - modelorog, "surface_altitude", "m", shape=(1, 1) + self.target_orog_cube = make_ancil_cube( + target_orog, "surface_altitude", "m", shape=(1, 1) + ) + self.model_orog_cube = make_ancil_cube( + model_orog, "surface_altitude", "m", shape=(1, 1) ) if heightlevels is not None: self.hl_cube = make_point_height_ancil_cube(heightlevels) @@ -281,13 +306,13 @@ def run_hc_rc(self, wind, height=None): """ self.w_cube = make_point_data_cube(wind, "wind_speed", "m s-1") plugin = RoughnessCorrection( - self.aos_cube, - self.s_cube, - self.poro_cube, - self.moro_cube, - 1500.0, - self.z0_cube, - self.hl_cube, + model_silhouette_roughness_cube=self.model_silhouette_roughness_cube, + model_orog_stddev_cube=self.s_cube, + target_orog_cube=self.target_orog_cube, + model_orog_cube=self.model_orog_cube, + model_res=1500.0, + model_z0_cube=self.model_z0_cube, + height_levels_cube=self.hl_cube, ) return plugin(self.w_cube) @@ -299,7 +324,7 @@ class Test1D(unittest.TestCase): passed, as well as testing the general behaviour of points that should not have a height corretion (equal height in model and pp orography) and the correct behaviour of doing roughness correction, - depending on whether or not a vegetative roughness (z_0) cube is + depending on whether or not a vegetative roughness (model_z0) cube is provided. Section 0 are tests where RMDI or np.nan values are passed. @@ -311,73 +336,81 @@ class Test1D(unittest.TestCase): hls = np.array([0.2, 3, 13, 33, 133, 333, 1133], dtype=np.float32) def test_section0a(self): - """Test AoS is RMDI, point should not do anything, uin = uout.""" - landpointtests_hc_rc = TestSinglePoint(AoS=RMDI, heightlevels=self.hls) + """Test silhouette_roughness is RMDI, point should not do anything, uin = uout.""" + landpointtests_hc_rc = TestSinglePoint( + silhouette_roughness=RMDI, heightlevels=self.hls + ) land_hc_rc = landpointtests_hc_rc.run_hc_rc(self.uin) np.testing.assert_array_equal(landpointtests_hc_rc.w_cube, land_hc_rc) def test_section0b(self): - """Test AoS is np.nan, point should not do anything, uin = uout.""" - landpointtests_hc_rc = TestSinglePoint(AoS=np.nan, heightlevels=self.hls) + """Test silhouette_roughness is np.nan, point should not do anything, uin = uout.""" + landpointtests_hc_rc = TestSinglePoint( + silhouette_roughness=np.nan, heightlevels=self.hls + ) land_hc_rc = landpointtests_hc_rc.run_hc_rc(self.uin) np.testing.assert_array_equal(landpointtests_hc_rc.w_cube, land_hc_rc) def test_section0c(self): - """Test Sigma is RMDI, point should not do anything, uin = uout.""" - landpointtests_hc_rc = TestSinglePoint(Sigma=RMDI, heightlevels=self.hls) + """Test model_orog_stddev is RMDI, point should not do anything, uin = uout.""" + landpointtests_hc_rc = TestSinglePoint( + model_orog_stddev=RMDI, heightlevels=self.hls + ) land_hc_rc = landpointtests_hc_rc.run_hc_rc(self.uin) np.testing.assert_array_equal(landpointtests_hc_rc.w_cube, land_hc_rc) def test_section0d(self): - """Test Sigma is np.nan, point should not do anything, uin = uout.""" - landpointtests_hc_rc = TestSinglePoint(Sigma=np.nan, heightlevels=self.hls) + """Test model_orog_stddev is np.nan, point should not do anything, uin = uout.""" + landpointtests_hc_rc = TestSinglePoint( + model_orog_stddev=np.nan, heightlevels=self.hls + ) land_hc_rc = landpointtests_hc_rc.run_hc_rc(self.uin) np.testing.assert_array_equal(landpointtests_hc_rc.w_cube, land_hc_rc) def test_section0e(self): - """Test z_0 is RMDI, point should not do RC. + """Test model_z0 is RMDI, point should not do RC. modeloro = pporo, so point should not do HC, uin = uout. """ landpointtests_hc_rc = TestSinglePoint( - z_0=RMDI, pporog=230.0, heightlevels=self.hls + model_z0=RMDI, target_orog=230.0, heightlevels=self.hls ) land_hc_rc = landpointtests_hc_rc.run_hc_rc(self.uin) np.testing.assert_array_equal(landpointtests_hc_rc.w_cube, land_hc_rc) def test_section0f(self): - """Test z_0 is np.nan, point should not do RC. + """Test model_z0 is np.nan, point should not do RC. modeloro = pporo, so point should not do HC, uin = uout. """ landpointtests_hc_rc = TestSinglePoint( - z_0=np.nan, pporog=230.0, heightlevels=self.hls + model_z0=np.nan, target_orog=230.0, heightlevels=self.hls ) land_hc_rc = landpointtests_hc_rc.run_hc_rc(self.uin) np.testing.assert_array_equal(landpointtests_hc_rc.w_cube, land_hc_rc) def test_section0g(self): - """Test z_0 is RMDI, point should not do RC. + """Test model_z0 is RMDI, point should not do RC. modeloro < pporo, so point should do positive HC, uin < uout. """ - landpointtests_hc_rc = TestSinglePoint(z_0=RMDI, heightlevels=self.hls) + landpointtests_hc_rc = TestSinglePoint(model_z0=RMDI, heightlevels=self.hls) land_hc_rc = landpointtests_hc_rc.run_hc_rc(self.uin) self.assertTrue((land_hc_rc.data > landpointtests_hc_rc.w_cube.data).all()) def test_section0h(self): - """Test pporog is RMDI (QUESTION: or should this fail???) + """Test target_orog is RMDI (QUESTION: or should this fail???) RC could be done for this point, HC cannot. uin >= uout - and since z_0=height[0] + and since model_z0=height[0] uout[0] = 0 """ - landpointtests_hc_rc = TestSinglePoint(pporog=RMDI, heightlevels=self.hls) + landpointtests_hc_rc = TestSinglePoint(target_orog=RMDI, heightlevels=self.hls) land_hc_rc = landpointtests_hc_rc.run_hc_rc(self.uin) self.assertTrue( (land_hc_rc.data <= landpointtests_hc_rc.w_cube.data).all() @@ -385,15 +418,17 @@ def test_section0h(self): ) def test_section0i(self): - """Test pporog is np.nan (QUESTION: or should this fail???) + """Test target_orog is np.nan (QUESTION: or should this fail???) RC could be done for this point, HC cannot. uin >= uout - and since z_0=height[0] + and since model_z0=height[0] uout[0] = 0 """ - landpointtests_hc_rc = TestSinglePoint(pporog=np.nan, heightlevels=self.hls) + landpointtests_hc_rc = TestSinglePoint( + target_orog=np.nan, heightlevels=self.hls + ) land_hc_rc = landpointtests_hc_rc.run_hc_rc(self.uin) self.assertTrue( (land_hc_rc.data <= landpointtests_hc_rc.w_cube.data).all() @@ -401,15 +436,15 @@ def test_section0i(self): ) def test_section0j(self): - """Test modelorog is RMDI (QUESTION: or should this fail???). + """Test model_orog is RMDI (QUESTION: or should this fail???). RC could be done for this point, HC cannot. uin >= uout - and since z_0=height[0] + and since model_z0=height[0] uout[0] = 0 """ - landpointtests_hc_rc = TestSinglePoint(modelorog=RMDI, heightlevels=self.hls) + landpointtests_hc_rc = TestSinglePoint(model_orog=RMDI, heightlevels=self.hls) land_hc_rc = landpointtests_hc_rc.run_hc_rc(self.uin) self.assertTrue( (land_hc_rc.data <= landpointtests_hc_rc.w_cube.data).all() @@ -467,24 +502,24 @@ def test_section0n(self): def test_section1a(self): """Test HC only, HC = 0. - z_0 passed as None, hence RC not performed. - modelorg = pporog, hence HC = 0. + model_z0 passed as None, hence RC not performed. + modelorg = target_orog, hence HC = 0. uin = uout """ - landpointtests_hc = TestSinglePoint(z_0=None, modelorog=250.0) + landpointtests_hc = TestSinglePoint(model_z0=None, model_orog=250.0) land_hc_rc = landpointtests_hc.run_hc_rc(self.uin) np.testing.assert_array_equal(landpointtests_hc.w_cube, land_hc_rc) def test_section1b(self): """Test HC only. - z_0 passed as None, hence RC not performed. - modelorg < pporog, hence positive HC. + model_z0 passed as None, hence RC not performed. + modelorg < target_orog, hence positive HC. uin <= uout, at least one height has uin < uout. """ - landpointtests_hc = TestSinglePoint(z_0=None) + landpointtests_hc = TestSinglePoint(model_z0=None) land_hc_rc = landpointtests_hc.run_hc_rc(self.uin) self.assertTrue( (land_hc_rc.data >= landpointtests_hc.w_cube.data).all() @@ -494,12 +529,12 @@ def test_section1b(self): def test_section1c(self): """Test RC and HC, HC=0. - z_0 passed, hence RC performed. - modelorg == pporog, hence no HC. + model_z0 passed, hence RC performed. + modelorg == target_orog, hence no HC. uin >= uout, at least one height has uin > uout, uout[0] = 0. """ - landpointtests_rc = TestSinglePoint(modelorog=250.0) + landpointtests_rc = TestSinglePoint(model_orog=250.0) land_hc_rc = landpointtests_rc.run_hc_rc(self.uin) self.assertTrue( (land_hc_rc.data <= landpointtests_rc.w_cube.data).all() @@ -510,15 +545,15 @@ def test_section1c(self): def test_section1d(self): """Test RC and HC. - z_0 passed, hence RC performed. - modelorg >> pporog, hence negative HC. + model_z0 passed, hence RC performed. + modelorg >> target_orog, hence negative HC. uin >= uout, at least one height has uin > uout - z_0 = height[0] hence RC[0] results in 0. + model_z0 = height[0] hence RC[0] results in 0. uout[0] RC is 0. HC is negative, negative speeds not allowed. Must be 0. """ - landpointtests_hc_rc = TestSinglePoint(pporog=230.0, modelorog=250.0) + landpointtests_hc_rc = TestSinglePoint(target_orog=230.0, model_orog=250.0) land_hc_rc = landpointtests_hc_rc.run_hc_rc(self.uin) self.assertTrue( (land_hc_rc.data <= landpointtests_hc_rc.w_cube.data).all() @@ -528,26 +563,26 @@ def test_section1d(self): ) def test_section1e(self): - """Test RC and HC, but sea point masked out (AoS). + """Test RC and HC, but sea point masked out (silhouette_roughness). - Sea point according to (AoS=0) => masked out. - z_0 passed, hence RC performed in theory. + Sea point according to (silhouette_roughness=0) => masked out. + model_z0 passed, hence RC performed in theory. uin = uout. """ - landpointtests_hc_rc = TestSinglePoint(AoS=0.0) + landpointtests_hc_rc = TestSinglePoint(silhouette_roughness=0.0) land_hc_rc = landpointtests_hc_rc.run_hc_rc(self.uin) np.testing.assert_array_equal(landpointtests_hc_rc.w_cube, land_hc_rc) def test_section1f(self): - """Test RC and HC, but sea point masked out (Sigma). + """Test RC and HC, but sea point masked out (model_orog_stddev). - Sea point according to (Sigma=0) => masked out - z_0 passed, hence RC performed in theory. + Sea point according to (model_orog_stddev=0) => masked out + model_z0 passed, hence RC performed in theory. uin = uout. """ - landpointtests_hc_rc = TestSinglePoint(Sigma=0.0) + landpointtests_hc_rc = TestSinglePoint(model_orog_stddev=0.0) land_hc_rc = landpointtests_hc_rc.run_hc_rc(self.uin) np.testing.assert_array_equal(landpointtests_hc_rc.w_cube, land_hc_rc) @@ -599,9 +634,9 @@ def test_section2b(self): heights = ((np.arange(10) + 1) ** 2.0) * 12 multip_hc_rc = TestMultiPoint( shape=(3, 1), - AoS=[0, 0.2, 0.2], - pporog=[0, 250, 250], - modelorog=[0, 250, 230], + silhouette_roughness=[0, 0.2, 0.2], + target_orog=[0, 250, 250], + model_orog=[0, 250, 230], ) land_hc_rc = multip_hc_rc.run_hc_rc(uin, dtime=2, height=heights) tidx = land_hc_rc.shape.index(2) @@ -638,11 +673,11 @@ def test_section2c(self): heights = ((np.arange(10) + 1) ** 2.0) * 12 multip_hc_rc = TestMultiPoint( shape=(3, 1), - AoS=[0, 0.2, 0.2], - pporog=[0, 250, 250], - modelorog=[0, 250, 230], + silhouette_roughness=[0, 0.2, 0.2], + target_orog=[0, 250, 250], + model_orog=[0, 250, 230], ) - msg = "wind input is not a cube, but " + msg = "Wind input is not a cube, but " with self.assertRaisesRegex(TypeError, msg): _ = multip_hc_rc.run_hc_rc([uin, uin], dtime=2, height=heights, aslist=True) @@ -656,53 +691,62 @@ def test_section2d(self): self.assertEqual(land_hc_rc.dtype, np.float32) def test_section3a(self): - """As test 1c, however with manipulated z_0 cube. + """As test 1c, however with manipulated model_z0 cube. - All ancillary fields have 1x1 dim, z_0 is on a different grid. + All ancillary fields have 1x1 dim, model_z0 is on a different grid. This should fail with ValueError("ancillary grids are not consistent"). """ - landpointtests_rc = TestSinglePoint(z_0=0.2, pporog=250.0, modelorog=250.0) + landpointtests_rc = TestSinglePoint( + model_z0=0.2, target_orog=250.0, model_orog=250.0 + ) z0_data = np.array( - [landpointtests_rc.z0_cube.data, landpointtests_rc.z0_cube.data] + [landpointtests_rc.model_z0_cube.data, landpointtests_rc.model_z0_cube.data] ) - landpointtests_rc.z0_cube = make_ancil_cube( + landpointtests_rc.model_z0_cube = make_ancil_cube( z0_data, "vegetative_roughness_length", "m", shape=(1, 2) ) - msg = "ancillary grids are not consistent" + msg = "Ancillary grids are not consistent." with self.assertRaisesRegex(ValueError, msg): _ = landpointtests_rc.run_hc_rc(self.uin) def test_section3b(self): - """As test 3a, however with manipulated modelorog cube instead. + """As test 3a, however with manipulated model_orog cube instead. This should fail with ValueError("ancillary grids are not consistent"). """ - landpointtests_rc = TestSinglePoint(z_0=0.2, pporog=250.0, modelorog=250.0) + landpointtests_rc = TestSinglePoint( + model_z0=0.2, target_orog=250.0, model_orog=250.0 + ) moro_data = np.array( - [landpointtests_rc.moro_cube.data, landpointtests_rc.moro_cube.data] + [ + landpointtests_rc.model_orog_cube.data, + landpointtests_rc.model_orog_cube.data, + ] ) - landpointtests_rc.moro_cube = make_ancil_cube( + landpointtests_rc.model_orog_cube = make_ancil_cube( moro_data, "surface_altitude", "m", shape=(1, 2) ) - msg = "ancillary grids are not consistent" + msg = "Ancillary grids are not consistent." with self.assertRaisesRegex(ValueError, msg): _ = landpointtests_rc.run_hc_rc(self.uin) def test_section3c(self): - """As test 3a, however with manipulated z_0 units. + """As test 3a, however with manipulated model_z0 units. This should fail with a wrong units error. """ - landpointtests_rc = TestSinglePoint(z_0=0.2, pporog=250.0, modelorog=250.0) - landpointtests_rc.z0_cube.units = Unit("s") - msg = "z0 ancil has unexpected unit: should be {} is {}" + landpointtests_rc = TestSinglePoint( + model_z0=0.2, target_orog=250.0, model_orog=250.0 + ) + landpointtests_rc.model_z0_cube.units = Unit("s") + msg = "z0 ancillary has unexpected unit: expected {}, got {}" with self.assertRaisesRegex( - ValueError, msg.format(Unit("m"), landpointtests_rc.z0_cube.units) + ValueError, msg.format(Unit("m"), landpointtests_rc.model_z0_cube.units) ): _ = landpointtests_rc.run_hc_rc(self.uin)