From b7b4a20e7f3cca27e3b9fae201c8119e35a3e312 Mon Sep 17 00:00:00 2001 From: cpelley Date: Tue, 12 Nov 2024 12:00:20 +0000 Subject: [PATCH 1/3] some changes --- improver/api/__init__.py | 2 +- improver/cli/merge.py | 5 +- improver/cli/orographic_enhancement.py | 62 ++--------- improver/orographic_enhancement.py | 130 ++++++++++++++++++++---- improver/utilities/cube_manipulation.py | 18 ++-- 5 files changed, 133 insertions(+), 84 deletions(-) diff --git a/improver/api/__init__.py b/improver/api/__init__.py index 366f4dbdff..556a99ddea 100644 --- a/improver/api/__init__.py +++ b/improver/api/__init__.py @@ -86,7 +86,7 @@ "OccurrenceBetweenThresholds": "improver.between_thresholds", "OccurrenceWithinVicinity": "improver.utilities.spatial", "OpticalFlow": "improver.nowcasting.optical_flow", - "OrographicEnhancement": "improver.orographic_enhancement", + "MetaOrographicEnhancement": "improver.orographic_enhancement", "OrographicSmoothingCoefficients": "improver.generate_ancillaries.generate_orographic_smoothing_coefficients", "PercentileConverter": "improver.percentile", "PhaseChangeLevel": "improver.psychrometric_calculations.psychrometric_calculations", diff --git a/improver/cli/merge.py b/improver/cli/merge.py index eea62178e9..8a1f84f90d 100755 --- a/improver/cli/merge.py +++ b/improver/cli/merge.py @@ -21,8 +21,5 @@ def process(*cubes: cli.inputcube): iris.cube.Cube: A merged cube. """ - from iris.cube import CubeList - from improver.utilities.cube_manipulation import MergeCubes - - return MergeCubes()(CubeList(cubes)) + return MergeCubes()(*cubes) diff --git a/improver/cli/orographic_enhancement.py b/improver/cli/orographic_enhancement.py index ba436fcf48..45c142e6da 100755 --- a/improver/cli/orographic_enhancement.py +++ b/improver/cli/orographic_enhancement.py @@ -4,44 +4,9 @@ # This file is part of IMPROVER and is released under a BSD 3-Clause license. # See LICENSE in the root of the repository for full licensing details. """Script to calculate orographic enhancement.""" - - from improver import cli -def extract_and_check(cube, height_value, units): - """ - Function to attempt to extract a height level. - If no matching level is available an error is raised. - - Args: - cube (cube): - Cube to be extracted from and checked it worked. - height_value (float): - The boundary height to be extracted with the input units. - units (str): - The units of the height level to be extracted. - Returns: - iris.cube.Cube: - A cube containing the extracted height level. - Raises: - ValueError: If height level is not found in the input cube. - """ - from improver.utilities.cube_extraction import extract_subcube - - # Write constraint in this format so a constraint is constructed that - # is suitable for floating point comparison - height_constraint = [ - "height=[{}:{}]".format(height_value - 0.1, height_value + 0.1) - ] - cube = extract_subcube(cube, height_constraint, units=[units]) - - if cube is not None: - return cube - - raise ValueError("No data available at height {}{}".format(height_value, units)) - - @cli.clizefy @cli.with_output def process( @@ -85,20 +50,13 @@ def process( Precipitation enhancement due to orography on the high resolution input orography grid. """ - from improver.orographic_enhancement import OrographicEnhancement - from improver.wind_calculations.wind_components import ResolveWindComponents - - constraint_info = (boundary_height, boundary_height_units) - - temperature = extract_and_check(temperature, *constraint_info) - humidity = extract_and_check(humidity, *constraint_info) - pressure = extract_and_check(pressure, *constraint_info) - wind_speed = extract_and_check(wind_speed, *constraint_info) - wind_direction = extract_and_check(wind_direction, *constraint_info) - - # resolve u and v wind components - u_wind, v_wind = ResolveWindComponents()(wind_speed, wind_direction) - # calculate orographic enhancement - return OrographicEnhancement()( - temperature, humidity, pressure, u_wind, v_wind, orography - ) + from improver.orographic_enhancement import MetaOrographicEnhancement + + return MetaOrographicEnhancement(boundary_height, boundary_height_units)( + temperature, + humidity, + pressure, + wind_speed, + wind_direction, + orography, + ) \ No newline at end of file diff --git a/improver/orographic_enhancement.py b/improver/orographic_enhancement.py index 04c3cc05cb..db23543a64 100644 --- a/improver/orographic_enhancement.py +++ b/improver/orographic_enhancement.py @@ -6,33 +6,129 @@ This module contains a plugin to calculate the enhancement of precipitation over orography. """ -from typing import Tuple - -import iris import numpy as np -from iris.analysis.cartography import rotate_winds -from iris.cube import Cube from numpy import ndarray from scipy.ndimage import uniform_filter1d +import iris from improver import BasePlugin from improver.constants import R_WATER_VAPOUR from improver.metadata.constants.mo_attributes import MOSG_GRID_ATTRIBUTES from improver.metadata.utilities import generate_mandatory_attributes from improver.nbhood.nbhood import NeighbourhoodProcessing -from improver.psychrometric_calculations.psychrometric_calculations import ( - calculate_svp_in_air, -) +from improver.psychrometric_calculations.psychrometric_calculations import \ + calculate_svp_in_air +from improver.utilities.common_input_handle import as_cube from improver.utilities.cube_checker import check_for_x_and_y_axes -from improver.utilities.cube_manipulation import ( - compare_coords, - enforce_coordinate_ordering, - sort_coord_in_cube, -) -from improver.utilities.spatial import ( - GradientBetweenAdjacentGridSquares, - number_of_grid_cells_to_distance, -) +from improver.utilities.cube_manipulation import (compare_coords, + enforce_coordinate_ordering, + sort_coord_in_cube) +from improver.utilities.spatial import (GradientBetweenAdjacentGridSquares, + number_of_grid_cells_to_distance) +from improver.wind_calculations.wind_components import ResolveWindComponents +from iris.analysis.cartography import rotate_winds +from iris.cube import Cube +from typing import Tuple + + +class MetaOrographicEnhancement(BasePlugin): + """Calculate orographic enhancement""" + def __init__(self, boundary_height: float = 1000.0, boundary_height_units="m"): + """ + Initialise the orographic enhancement plugin. + + Args: + boundary_height (float): + Model height level to extract variables for calculating orographic + enhancement, as proxy for the boundary layer. + boundary_height_units (str): + Units of the boundary height specified for extracting model levels. + """ + self._constraint_info = (boundary_height, boundary_height_units) + + @staticmethod + def extract_and_check(cube, height_value, units): + """ + Function to attempt to extract a height level. + If no matching level is available an error is raised. + + Args: + cube (cube): + Cube to be extracted from and checked it worked. + height_value (float): + The boundary height to be extracted with the input units. + units (str): + The units of the height level to be extracted. + Returns: + iris.cube.Cube: + A cube containing the extracted height level. + Raises: + ValueError: If height level is not found in the input cube. + """ + from improver.utilities.cube_extraction import extract_subcube + + # Write constraint in this format so a constraint is constructed that + # is suitable for floating point comparison + height_constraint = [ + "height=[{}:{}]".format(height_value - 0.1, height_value + 0.1) + ] + cube = extract_subcube(cube, height_constraint, units=[units]) + + if cube is not None: + return cube + + raise ValueError("No data available at height {}{}".format(height_value, units)) + + def process(self, + temperature: Cube, + humidity: Cube, + pressure: Cube, + wind_speed: Cube, + wind_direction: Cube, + orography: Cube, + ) -> Cube: + """ + Uses the ResolveWindComponents() and OrographicEnhancement() plugins. + Outputs data on the high resolution orography grid. + + Args: + temperature: + Cube containing temperature at top of boundary layer. + humidity: + Cube containing relative humidity at top of boundary layer. + pressure: + Cube containing pressure at top of boundary layer. + wind_speed: + Cube containing wind speed values. + wind_direction: + Cube containing wind direction values relative to true north. + orography: + Cube containing height of orography above sea level on high + resolution (1 km) UKPP domain grid. + + Returns: + iris.cube.Cube: + Precipitation enhancement due to orography on the high resolution + input orography grid. + """ + temperature = as_cube(temperature) + humidity = as_cube(humidity) + pressure = as_cube(pressure) + wind_speed = as_cube(wind_speed) + wind_direction = as_cube(wind_direction) + orography = as_cube(orography) + temperature = self.extract_and_check(temperature, *self._constraint_info) + humidity = self.extract_and_check(humidity, *self._constraint_info) + pressure = self.extract_and_check(pressure, *self._constraint_info) + wind_speed = self.extract_and_check(wind_speed, *self._constraint_info) + wind_direction = self.extract_and_check(wind_direction, *self._constraint_info) + + # resolve u and v wind components + u_wind, v_wind = ResolveWindComponents()(wind_speed, wind_direction) + # calculate orographic enhancement + return OrographicEnhancement()( + temperature, humidity, pressure, u_wind, v_wind, orography + ) class OrographicEnhancement(BasePlugin): diff --git a/improver/utilities/cube_manipulation.py b/improver/utilities/cube_manipulation.py index 0e61bf828a..ca3da0a7a2 100644 --- a/improver/utilities/cube_manipulation.py +++ b/improver/utilities/cube_manipulation.py @@ -16,7 +16,7 @@ from improver import BasePlugin from improver.metadata.constants import FLOAT_DTYPE, FLOAT_TYPES from improver.metadata.probabilistic import find_threshold_coordinate -from improver.utilities.common_input_handle import as_cube +from improver.utilities.common_input_handle import as_cube, as_cubelist from improver.utilities.cube_checker import check_cube_coordinates @@ -208,7 +208,7 @@ def _check_time_bounds_ranges(cube: Cube) -> None: def process( self, - cubes_in: Union[List[Cube], CubeList], + *cubes: Union[Cube, CubeList], check_time_bounds_ranges: bool = False, slice_over_realization: bool = False, copy: bool = True, @@ -225,7 +225,7 @@ def process( result of premature iris merging on load). Args: - cubes_in: + *cubes: Cubes to be merged. check_time_bounds_ranges: Flag to check whether scalar time bounds ranges match. @@ -244,16 +244,14 @@ def process( Returns: Merged cube. """ - # if input is already a single cube, return unchanged - if isinstance(cubes_in, iris.cube.Cube): - return cubes_in + cubes = as_cubelist(*cubes) - if len(cubes_in) == 1: + if len(cubes) == 1: # iris merges cubelist into shortest list possible on load # - may already have collapsed across invalid time bounds if check_time_bounds_ranges: - self._check_time_bounds_ranges(cubes_in[0]) - return cubes_in[0] + self._check_time_bounds_ranges(cubes[0]) + return cubes[0] if copy: # create copies of input cubes so as not to modify in place @@ -262,7 +260,7 @@ def process( cube_return = lambda cube: cube cubelist = iris.cube.CubeList([]) - for cube in cubes_in: + for cube in cubes: if slice_over_realization: for real_slice in cube.slices_over("realization"): cubelist.append(cube_return(real_slice)) From b02c772e3e5758a01378dca3fb7678a6ee35391f Mon Sep 17 00:00:00 2001 From: cpelley Date: Tue, 26 Nov 2024 14:50:52 +0000 Subject: [PATCH 2/3] further changes --- .../cli/nowcast_optical_flow_from_winds.py | 20 ++++++++----------- improver/nowcasting/optical_flow.py | 6 +++++- 2 files changed, 13 insertions(+), 13 deletions(-) diff --git a/improver/cli/nowcast_optical_flow_from_winds.py b/improver/cli/nowcast_optical_flow_from_winds.py index 90df1adeac..613c1b36c1 100644 --- a/improver/cli/nowcast_optical_flow_from_winds.py +++ b/improver/cli/nowcast_optical_flow_from_winds.py @@ -19,7 +19,8 @@ def process( steering_flow: inputflow, orographic_enhancement: cli.inputcube, - *cubes: cli.inputcube, + radar_precip_1: cli.inputcube, + radar_precip_2: cli.inputcube, ): """Calculate optical flow components as perturbations from the model steering flow. Advects the older of the two input radar observations to @@ -34,24 +35,19 @@ def process( have names: "grid_eastward_wind" and "grid_northward_wind". orographic_enhancement (iris.cube.Cube): Cube containing the orographic enhancement fields. - cubes (tuple of iris.cube.Cube): - Two radar precipitation observation cubes. + radar_precip_1 (iris.cube.Cube): + First radar precipitation observation cube. + radar_precip_2 (iris.cube.Cube): + Second precipitation observation cube. Returns: iris.cube.CubeList: List of u- and v- advection velocities """ - from iris.cube import CubeList - from improver.nowcasting.optical_flow import ( generate_advection_velocities_from_winds, ) - if len(cubes) != 2: - raise ValueError("Expected 2 radar cubes - got {}".format(len(cubes))) - - advection_velocities = generate_advection_velocities_from_winds( - CubeList(cubes), steering_flow, orographic_enhancement + return generate_advection_velocities_from_winds( + radar_precip_1, radar_precip_2, steering_flow, orographic_enhancement ) - - return advection_velocities diff --git a/improver/nowcasting/optical_flow.py b/improver/nowcasting/optical_flow.py index 22b79231eb..e921b136ab 100644 --- a/improver/nowcasting/optical_flow.py +++ b/improver/nowcasting/optical_flow.py @@ -25,6 +25,7 @@ from improver.nowcasting.utilities import ApplyOrographicEnhancement from improver.utilities.cube_checker import check_for_x_and_y_axes from improver.utilities.cube_manipulation import collapsed +from improver.utilities.common_input_handle import as_cubelist from improver.utilities.spatial import ( calculate_grid_spacing, check_if_grid_is_equal_area, @@ -82,7 +83,7 @@ def _calculate_time_average(wind_cubes, time_coord): def generate_advection_velocities_from_winds( - cubes: CubeList, background_flow: CubeList, orographic_enhancement: Cube + radar_precip_1: Cube, radar_precip_2: Cube, background_flow: CubeList, orographic_enhancement: Cube ) -> CubeList: """Generate advection velocities as perturbations from a non-zero background flow @@ -99,6 +100,9 @@ def generate_advection_velocities_from_winds( Returns: u- and v- advection velocities """ + background_flow = background_flow.extract_cubes("grid_eastward_wind", "grid_northward_wind") + + cubes = as_cubelist(radar_precip_1, radar_precip_2) cubes.sort(key=lambda x: x.coord("time").points[0]) lead_time_seconds = ( From e22bf3c5176b9b1efa84a94dc1865626d3042b8a Mon Sep 17 00:00:00 2001 From: cpelley Date: Mon, 2 Dec 2024 16:59:01 +0000 Subject: [PATCH 3/3] weighted_blending update --- .../blending/calculate_weights_and_blend.py | 134 ++++++++++-------- improver/cli/blend_cycles_and_realizations.py | 5 +- improver/cli/weighted_blending.py | 6 +- 3 files changed, 81 insertions(+), 64 deletions(-) diff --git a/improver/blending/calculate_weights_and_blend.py b/improver/blending/calculate_weights_and_blend.py index 527232373f..c55c4de1b4 100644 --- a/improver/blending/calculate_weights_and_blend.py +++ b/improver/blending/calculate_weights_and_blend.py @@ -51,6 +51,12 @@ def __init__( ynval: Optional[float] = None, cval: Optional[float] = None, inverse_ordering: bool = False, + cycletime: Optional[str] = None, + model_id_attr: Optional[str] = None, + record_run_attr: Optional[str] = None, + spatial_weights: bool = False, + fuzzy_length: float = 20000, + attributes_dict: Optional[Dict[str, str]] = None, ) -> None: """ Initialise central parameters @@ -90,6 +96,38 @@ def __init__( Option to invert weighting order for non-linear weights plugin so that higher blend coordinate values get higher weights (eg if cycle blending over forecast reference time). + cycletime: + The forecast reference time to be used after blending has been + applied, in the format YYYYMMDDTHHMMZ. If not provided, the + blended file takes the latest available forecast reference time + from the input datasets supplied. + model_id_attr: + The name of the dataset attribute to be used to identify the source + model when blending data from different models. + record_run_attr: + The name of the dataset attribute to be used to store model and + cycle sources in metadata, e.g. when blending data from different + models. Requires model_id_attr. + spatial_weights: + If True, this option will result in the generation of spatially + varying weights based on the masks of the data we are blending. + The one dimensional weights are first calculated using the chosen + weights calculation method, but the weights will then be adjusted + spatially based on where there is masked data in the data we are + blending. The spatial weights are calculated using the + SpatiallyVaryingWeightsFromMask plugin. + fuzzy_length: + When calculating spatially varying weights we can smooth the + weights so that areas close to areas that are masked have lower + weights than those further away. This fuzzy length controls the + scale over which the weights are smoothed. The fuzzy length is in + terms of m, the default is 20km. This distance is then converted + into a number of grid squares, which does not have to be an + integer. Assumes the grid spacing is the same in the x and y + directions and raises an error if this is not true. See + SpatiallyVaryingWeightsFromMask for more details. + attributes_dict: + Dictionary describing required changes to attributes after blending """ self.blend_coord = blend_coord self.wts_calc_method = wts_calc_method @@ -110,6 +148,30 @@ def __init__( self.wts_calc_method ) ) + + self._cycletime = cycletime + self._model_id_attr = model_id_attr + self._record_run_attr = record_run_attr + self._spatial_weights = spatial_weights + self._fuzzy_length = fuzzy_length + self._attributes_dict = attributes_dict + + if record_run_attr is not None and model_id_attr is None: + raise ValueError( + "record_run_attr can only be used with model_id_attr, which " + "has not been provided." + ) + + if (wts_calc_method == "linear") and cval: + raise RuntimeError("Method: linear does not accept arguments: cval") + if (wts_calc_method == "nonlinear") and any([y0val, ynval]): + raise RuntimeError("Method: non-linear does not accept arguments: y0val, ynval") + if (wts_calc_method == "dict") and wts_dict is None: + raise RuntimeError('Dictionary is required if wts_calc_method="dict"') + if "model" in blend_coord and model_id_attr is None: + raise RuntimeError("model_id_attr must be specified for model blending") + if record_run_attr is not None and model_id_attr is None: + raise RuntimeError("model_id_attr must be specified for blend model recording") def _calculate_blending_weights(self, cube: Cube) -> Cube: """ @@ -206,13 +268,7 @@ def _remove_zero_weighted_slices( def process( self, - cubelist: Union[List[Cube], CubeList], - cycletime: Optional[str] = None, - model_id_attr: Optional[str] = None, - record_run_attr: Optional[str] = None, - spatial_weights: bool = False, - fuzzy_length: float = 20000, - attributes_dict: Optional[Dict[str, str]] = None, + *cubes: Union[Cube, CubeList], ) -> Cube: """ Merge a cubelist, calculate appropriate blend weights and compute the @@ -220,40 +276,8 @@ def process( given by self.blend_coord. Args: - cubelist: - List of cubes to be merged and blended - cycletime: - The forecast reference time to be used after blending has been - applied, in the format YYYYMMDDTHHMMZ. If not provided, the - blended file takes the latest available forecast reference time - from the input datasets supplied. - model_id_attr: - The name of the dataset attribute to be used to identify the source - model when blending data from different models. - record_run_attr: - The name of the dataset attribute to be used to store model and - cycle sources in metadata, e.g. when blending data from different - models. Requires model_id_attr. - spatial_weights: - If True, this option will result in the generation of spatially - varying weights based on the masks of the data we are blending. - The one dimensional weights are first calculated using the chosen - weights calculation method, but the weights will then be adjusted - spatially based on where there is masked data in the data we are - blending. The spatial weights are calculated using the - SpatiallyVaryingWeightsFromMask plugin. - fuzzy_length: - When calculating spatially varying weights we can smooth the - weights so that areas close to areas that are masked have lower - weights than those further away. This fuzzy length controls the - scale over which the weights are smoothed. The fuzzy length is in - terms of m, the default is 20km. This distance is then converted - into a number of grid squares, which does not have to be an - integer. Assumes the grid spacing is the same in the x and y - directions and raises an error if this is not true. See - SpatiallyVaryingWeightsFromMask for more details. - attributes_dict: - Dictionary describing required changes to attributes after blending + *cubes: + One of more cubes to be merged and blended Returns: Cube of blended data. @@ -266,11 +290,7 @@ def process( UserWarning: If blending masked data without spatial weights. This has not been fully tested. """ - if record_run_attr is not None and model_id_attr is None: - raise ValueError( - "record_run_attr can only be used with model_id_attr, which " - "has not been provided." - ) + cubes = as_cubelist(cubes) # Prepare cubes for weighted blending, including creating custom metadata # for multi-model blending. The merged cube has a monotonically ascending @@ -279,10 +299,10 @@ def process( merger = MergeCubesForWeightedBlending( self.blend_coord, weighting_coord=self.weighting_coord, - model_id_attr=model_id_attr, - record_run_attr=record_run_attr, + model_id_attr=self._model_id_attr, + record_run_attr=self._record_run_attr, ) - cube = merger(cubelist, cycletime=cycletime) + cube = merger(cubes, cycletime=self._cycletime) if "model" in self.blend_coord: self.blend_coord = copy(MODEL_BLEND_COORD) @@ -296,15 +316,15 @@ def process( weights = self._calculate_blending_weights(cube) cube, weights = self._remove_zero_weighted_slices(cube, weights) - if record_run_attr is not None and weights is not None: + if self._record_run_attr is not None and weights is not None: cube = update_record_run_weights(cube, weights, self.blend_coord) # Deal with case of only one input cube or non-zero-weighted slice if len(cube.coord(self.blend_coord).points) == 1: result = cube else: - if spatial_weights: - weights = self._update_spatial_weights(cube, weights, fuzzy_length) + if self._spatial_weights: + weights = self._update_spatial_weights(cube, weights, self._fuzzy_length) elif np.ma.is_masked(cube.data): # Raise warning if blending masked arrays using non-spatial weights. warnings.warn( @@ -316,8 +336,8 @@ def process( BlendingPlugin = WeightedBlendAcrossWholeDimension(self.blend_coord) result = BlendingPlugin(cube, weights=weights) - if record_run_attr is not None: - record_run_coord_to_attr(result, cube, record_run_attr) + if self._record_run_attr is not None: + record_run_coord_to_attr(result, cube, self._record_run_attr) # Remove custom metadata and and update time-type coordinates. Remove # non-time-type coordinate that were previously associated with the blend @@ -327,9 +347,9 @@ def process( result, self.blend_coord, coords_to_remove=coords_to_remove, - cycletime=cycletime, - attributes_dict=attributes_dict, - model_id_attr=model_id_attr, + cycletime=self._cycletime, + attributes_dict=self._attributes_dict, + model_id_attr=self._model_id_attr, ) return result diff --git a/improver/cli/blend_cycles_and_realizations.py b/improver/cli/blend_cycles_and_realizations.py index 2dda709c24..ec3b6cb136 100755 --- a/improver/cli/blend_cycles_and_realizations.py +++ b/improver/cli/blend_cycles_and_realizations.py @@ -42,6 +42,5 @@ def process( for cube in cubes: cubelist.append(collapse_realizations(cube)) - plugin = WeightAndBlend("forecast_reference_time", "linear", y0val=0.5, ynval=0.5,) - cube = plugin(cubelist, cycletime=cycletime,) - return cube + plugin = WeightAndBlend("forecast_reference_time", "linear", y0val=0.5, ynval=0.5, cycletime=cycletime) + return plugin(cubelist) diff --git a/improver/cli/weighted_blending.py b/improver/cli/weighted_blending.py index 32abb7ae31..86a3009590 100755 --- a/improver/cli/weighted_blending.py +++ b/improver/cli/weighted_blending.py @@ -135,10 +135,6 @@ def process( y0val=y0val, ynval=ynval, cval=cval, - ) - - return plugin( - cubes, cycletime=cycletime, model_id_attr=model_id_attr, record_run_attr=record_run_attr, @@ -146,3 +142,5 @@ def process( fuzzy_length=fuzzy_length, attributes_dict=attributes_config, ) + + return plugin(cubes)