diff --git a/src/CSET/cset_workflow/app/fetch_fcst/bin/fetch_data.py b/src/CSET/cset_workflow/app/fetch_fcst/bin/fetch_data.py index 39b4c92f2..692b83c67 100644 --- a/src/CSET/cset_workflow/app/fetch_fcst/bin/fetch_data.py +++ b/src/CSET/cset_workflow/app/fetch_fcst/bin/fetch_data.py @@ -19,6 +19,7 @@ import itertools import logging import os +import shutil import ssl import urllib.parse import urllib.request @@ -177,6 +178,25 @@ def _get_needed_environment_variables() -> dict: return variables +def _get_needed_environment_variables_nimrod() -> dict: + """Load the needed variables from the environment to retrieve UK Nimrod data.""" + variables = { + "field": [ + os.environ["NIMROD_FIELDS_COMP_HOUR"], + os.environ["NIMROD_FIELDS_1KM"], + os.environ["NIMROD_FIELDS_2KM"], + ], + "raw_path": "/data/users/radar/UKnimrod", + "date_type": "initiation", + "data_time": datetime.fromisoformat(os.environ["CYLC_TASK_CYCLE_POINT"]), + "forecast_length": isodate.parse_duration(os.environ["ANALYSIS_LENGTH"]), + "model_identifier": "Nimrod", + "rose_datac": os.environ["ROSE_DATAC"], + } + logging.debug("Environment variables loaded for Nimrod: %s", variables) + return variables + + def _template_file_path( raw_path: str, date_type: Literal["validity", "initiation"], @@ -270,3 +290,126 @@ def fetch_data(file_retriever: FileRetrieverABC): # exiting the with block. if not files_found: raise FileNotFoundError("No files found for model!") + + +def _template_nimrod_files(field, use_date) -> dict: + """Write the dictionary to retrieve Nimrod obs.""" + files = { + "rainaccum_comp_hour": { + "basedir": "/data/users/radar/UKnimrod", + "dir": "rainaccum_comp_hour", + "fname": "_nimrod_ng_radar_rainaccum_comp_hour", + "weights_dir": "rainaccwt_comp_hour", + "weights_fname": "_nimrod_ng_radar_rainaccwt_comp_hour", + "freq": "1hour", + }, + "rainaccum_comp_hour_1km_cutout": { + "basedir": "/data/users/radar/UKnimrod", + "dir": "rainaccum_comp_hour_1km_cutout", + "fname": "_nimrod_ng_radar_rainaccum_comp_hour_1km_cutout_775X640_eng", + "weights_dir": "rainaccwt_comp_hour_1km_cutout", + "weights_fname": "_nimrod_ng_radar_rainaccwt_comp_hour_1km_cutout_775X640_eng", + "freq": "1hour", + }, + "rainaccum_comp_hour_2km": { + "basedir": "/data/users/radar/UKnimrod", + "dir": "rainaccum_comp_hour_2km", + "fname": "_nimrod_ng_radar_rainaccum_comp_hour_2km", + "weights_dir": "rainaccwt_comp_hour_2km", + "weights_fname": "_nimrod_ng_radar_rainaccwt_comp_hour_2km", + "freq": "1hour", + }, + } + + # initialise the dictionary to return, putting in null values of "None" + filepath = {"obs": "None", "weights": "None"} + + if field in files.keys(): + # form the filename of the Nimrod obs file to retrieve from the RMED archive + nextfile = ( + files[field]["basedir"] + + "/" + + files[field]["dir"] + + "/" + + f"{use_date.year}/{use_date.strftime('%Y%m%d%H')}00" + + files[field]["fname"] + ) + filepath["obs"] = nextfile + + # form the filename of the corresponding Nimrod obs weighting file to retrieve + nextfile = ( + files[field]["basedir"] + + "/" + + files[field]["weights_dir"] + + "/" + + f"{use_date.year}/{use_date.strftime('%Y%m%d%H')}00" + + files[field]["weights_fname"] + ) + filepath["weights"] = nextfile + + return filepath + + +def _get_nimrod_file_paths(start_date, end_date, fields): + filepaths: set[str] = set() + use_date = start_date + while use_date <= end_date: + for next_field in fields: + # form the Nimrod file names + files = _template_nimrod_files(next_field, use_date) + if files["obs"] != "None": + filepaths.add(files["obs"]) + filepaths.add(files["weights"]) + + # advance the date/time counter on a hour + use_date = use_date + timedelta(hours=1) + + return filepaths + + +# def fetch_nimrod(nimrod_retriever: FileRetrieverABC): +def fetch_nimrod(): + # def fetch_nimrod(file_retriever: FileRetrieverABC): + """Fetch the observations corresponding to a model run. + + The following environment variables need to be set: + * NIMROD_FIELDs + #* ANALYSIS_OFFSET + #* ANALYSIS_LENGTH + #* CYLC_TASK_CYCLE_POINT + * DATA_PATH + #* DATA_PERIOD + #* DATE_TYPE + #* MODEL_IDENTIFIER + * ROSE_DATAC + + Parameters + ---------- + nimrod_retriever: ObsRetriever + ObsRetriever implementation to use. Defaults to FilesystemFileRetriever. + + Raises + ------ + FileNotFound: + If no observations are available. + """ + v = _get_needed_environment_variables_nimrod() + + # Prepare output directory. + cycle_nimrod_dir = f"{v['rose_datac']}/data/Nimrod" + os.makedirs(cycle_nimrod_dir, exist_ok=True) + logging.debug("Nimrod directory: %s", cycle_nimrod_dir) + + # form the Nimrod files to use + start_date = v["data_time"] + end_date = v["data_time"] + v["forecast_length"] + filepaths: set[str] = set() + filepaths = _get_nimrod_file_paths(start_date, end_date, v["field"]) + + # copy the Nimrod files from the archive to the local directory + print(" Copying Nimrod files to directory: ", cycle_nimrod_dir) + for i in filepaths: + print(" Grabbing Nimrod file: ", i) + shutil.copy(i, cycle_nimrod_dir) + + return diff --git a/src/CSET/cset_workflow/app/fetch_nimrod/bin/constants.py b/src/CSET/cset_workflow/app/fetch_nimrod/bin/constants.py new file mode 100644 index 000000000..86f525308 --- /dev/null +++ b/src/CSET/cset_workflow/app/fetch_nimrod/bin/constants.py @@ -0,0 +1,9 @@ +"""Global constants.""" + +#: Number of seconds in an hour. +HOUR_IN_SECONDS = 3600 +#: Number of metres in a kilometre. +M_IN_KM = 1000 +#: Tolerance used when selecting by particular real valued +#: quantities (e.g. level number, latitude, longitude) +REAL_MATCH_TOL = 1.0e-4 diff --git a/src/CSET/cset_workflow/app/fetch_nimrod/bin/cube_utils.py b/src/CSET/cset_workflow/app/fetch_nimrod/bin/cube_utils.py new file mode 100644 index 000000000..24dcd126a --- /dev/null +++ b/src/CSET/cset_workflow/app/fetch_nimrod/bin/cube_utils.py @@ -0,0 +1,797 @@ +""" +General utility functions for use with Iris :class:`iris.cube.Cube` or \ +:class:`iris.cube.CubeList` objects. +""" + +import collections +import time +import warnings + +import agg_regrid +import constants +import dask +import iris +import iris.analysis +import matplotlib +import numpy as np +import shapely.geometry as sgeom +from ants.regrid.rectilinear import TwoStage +from iris.experimental.regrid_conservative import regrid_conservative_via_esmpy + +warnings.filterwarnings("ignore") + + +def get_temporal_range(cube): + """ + Returns the time range spanned by the time coordinate of an input \ + :class:`iris.cube.Cube`, in the form [first time, last time]. + Time coordinate bounds are taken into account if present. + """ + # Get a list of coordinate names for the cube + coord_names = [coord.name() for coord in cube.coords()] + # Check there is a time coordinate present + if "time" not in coord_names: + raise ValueError("Time coordinate not found in cube") + time_coord = cube.coord("time") + time_unit = time_coord.units + if time_coord.has_bounds(): + # Use time bounds if present... + times = time_coord.bounds + else: + # ...else just use the time points + times = time_coord.points + # Convert from numeric times to datetimes + times = time_unit.num2date(times) + return [np.min(times), np.max(times)] + + +def get_spatial_coords(cube): + """Returns the x, y coordinates of an input :class:`iris.cube.Cube`.""" + # Usual names for spatial coordinates + X_COORD_NAMES = ["longitude", "grid_longitude", "projection_x_coordinate", "x"] + Y_COORD_NAMES = ["latitude", "grid_latitude", "projection_y_coordinate", "y"] + + # Get a list of coordinate names for the cube + coord_names = [coord.name() for coord in cube.coords()] + + # Check which x-coordinate we have, if any + x_coords = [coord for coord in coord_names if coord in X_COORD_NAMES] + if len(x_coords) != 1: + raise ValueError("Could not identify a unique x-coordinate in cube") + x_coord = cube.coord(x_coords[0]) + + # Check which y-coordinate we have, if any + y_coords = [coord for coord in coord_names if coord in Y_COORD_NAMES] + if len(y_coords) != 1: + raise ValueError("Could not identify a unique y-coordinate in cube") + y_coord = cube.coord(y_coords[0]) + + return [x_coord, y_coord] + + +def get_spatial_coord_dims(cube): + """ + Returns a tuple of the dimensions corresponding to the x, y coordinates \ + of an input :class:`iris.cube.Cube`. + """ + # Get the spatial coordinates of the cube... + x_coord, y_coord = get_spatial_coords(cube) + # ...then the cube dimensions corresponding to these coordinates + x_dim = cube.coord_dims(x_coord)[0] + y_dim = cube.coord_dims(y_coord)[0] + return (x_dim, y_dim) + + +def get_spatial_extent(cube): + """ + Takes an :class:`iris.cube.Cube` and returns its spatial extent in the \ + form of a list [min x, max x, min y, max y]. + """ + # Get the spatial coordinates of the cube + x_coord, y_coord = get_spatial_coords(cube) + return [ + min(x_coord.points), + max(x_coord.points), + min(y_coord.points), + max(y_coord.points), + ] + + +def get_non_spatial_coords(cube): + """ + Returns a list of the non-spatial coordinates of an input \ + :class:`iris.cube.Cube`. + """ + # Get a list of the cube coordinates + coords = cube.coords() + # Get the spatial coordinates of the cube + x_coord, y_coord = get_spatial_coords(cube) + # Remove the spatial coordinates from the list of coordinates + coords.remove(x_coord) + coords.remove(y_coord) + return coords + + +def guess_bounds(cube): + """ + Takes an input :class:`iris.cube.Cube`, guesses bounds on the x, y \ + coordinates, then returns the cube. Such bounds are often required by \ + regridding algorithms. + """ + # Loop over spatial coordinates + for axis in ["x", "y"]: + coord = cube.coord(axis=axis) + # Check that this is not a variable resolution grid + # TODO: Does bounds really not work with variable resolution? AVD + # seem to think it won't... + try: + _ = iris.util.regular_step(coord) + except: + raise ValueError("Cannot guess bounds for a variable resolution grid") + # Guess bounds if there aren't any + if coord.bounds is None: + coord.guess_bounds() + return cube + + +def remove_bounds(cube): + """ + Takes an input :class:`iris.cube.Cube`, removes any bounds on the x, y \ + coordinates, then returns the cube. + """ + # Loop over spatial coordinates + for axis in ["x", "y"]: + cube.coord(axis=axis).bounds = None + return cube + + +def remove_cell_method(cube, cell_method): + """ + Removes the supplied :class:`iris.coords.CellMethod` from an input + :class:`iris.cube.Cube`, then returns the cube. + """ + cell_methods = [cm for cm in cube.cell_methods if cm != cell_method] + cube.cell_methods = () + for cm in cell_methods: + cube.add_cell_method(cm) + return cube + + +def get_grid(cube): + """ + Takes an :class:`iris.cube.Cube` and returns its spatial grid as a 2D \ + :class:`iris.cube.Cube`. + """ + # Take a copy of the input cube + grid_cube = cube.copy() + # Get the spatial coordinates of the cube + x_coord, y_coord = get_spatial_coords(grid_cube) + + # Take an arbitrary spatial slice of the cube + grid_cube = next(grid_cube.slices([y_coord, x_coord])) + # We don't care about the data of this slice so just zero it + grid_cube.data = np.zeros(grid_cube.data.shape) + # Add a mask to show which points in the spatial grid are to be + # considered valid. By default, all of them (False means not masked in + # masked arrays) + mask = np.zeros(grid_cube.data.shape, dtype=bool) + grid_cube.data = np.ma.array(grid_cube.data, mask=mask) + # Remove all non-spatial coordinates from the cube + coords = grid_cube.coords() + + coords.remove(x_coord) + coords.remove(y_coord) + + for coord in coords: + grid_cube.remove_coord(coord) + + # Set some metadata to signify this is just a grid cube + grid_cube.standard_name = None + grid_cube.long_name = "Grid" + grid_cube.units = "no_unit" + grid_cube.cell_methods = None + grid_cube.attributes = {} + + return grid_cube + + +def get_grid_spacing(cube): + """ + Takes a :class:`iris.cube.Cube` and returns a tuple (dx, dy) of the grid \ + spacings in the x, y directions. + """ + # Get the spatial coordinates of the cube + x_coord, y_coord = get_spatial_coords(cube) + # Note that regular_step will fail for variable resolution grids + dx = iris.util.regular_step(x_coord) + dy = iris.util.regular_step(y_coord) + return (dx, dy) + + +def get_min_max_grid_spacing(cube): + """ + Takes a :class:`iris.cube.Cube` and returns a tuple of the minimum and \ + maximum grid spacings in either of the x, y directions. + """ + # Get the grid spacings in the x, y directions + dx, dy = get_grid_spacing(cube) + return (np.min([dx, dy]), np.max([dx, dy])) + + +def trim_cube(cube, xmin=None, xmax=None, ymin=None, ymax=None, ignore_bounds=True): + """ + Trims the spatial extent of an :class:`iris.cube.Cube` to lie in a \ + specified range. + + Arguments: + + * **cube** - an :class:`iris.cube.Cube` object. + + Keyword arguments: + + * **xmin** - minimum value of the cube x-coordinate to select. + * **xmax** - maximum value of the cube x-coordinate to select. + * **ymin** - minimum value of the cube y-coordinate to select. + * **ymax** - maximum value of the cube y-coordinate to select. + * **ignore_bounds** - set to True to ignore any bounds on the coordinates \ + and just consider their points. + + Returns + ------- + * **trimmed_cube** - a copy of the input :class:`iris.cube.Cube` with \ + x, y coordinate values restricted to the specified \ + ranges. + """ + # Take a copy of the input cube + trimmed_cube = cube.copy() + + # Get the spatial coordinates of the cube + x_coord, y_coord = get_spatial_coords(trimmed_cube) + + # First trim the cube in the x direction... + if xmin is None: + xmin = np.min(x_coord.points) + xmin -= constants.REAL_MATCH_TOL + if xmax is None: + xmax = np.max(x_coord.points) + xmax += constants.REAL_MATCH_TOL + x_extent = iris.coords.CoordExtent(x_coord, xmin, xmax) + + trimmed_cube = trimmed_cube.intersection(x_extent, ignore_bounds=ignore_bounds) + + # ...now trim the cube in the y direction + if ymin is None: + ymin = np.min(y_coord.points) + ymin -= constants.REAL_MATCH_TOL + if ymax is None: + ymax = np.max(y_coord.points) + ymax += constants.REAL_MATCH_TOL + y_extent = iris.coords.CoordExtent(y_coord, ymin, ymax) + trimmed_cube = trimmed_cube.intersection(y_extent, ignore_bounds=ignore_bounds) + + return trimmed_cube + + +def trim_variable_grid(cube): + """ + Removes any variable resolution part of the spatial grid of an input \ + :class:`iris.cube.Cube`. + """ + # Take a copy of the input cube + trimmed_cube = cube.copy() + # Loop over spatial coordinates + for axis in ["x", "y"]: + delta = None + # Keep looping until the variable resolution part of the grid has been + # fully removed + while delta is None: + coord = trimmed_cube.coord(axis=axis) + try: + # This will fail for variable resolution grids + delta = iris.util.regular_step(coord) + except: + # Remove one grid cell from either end of the axis + trimmed_coord = coord.points[1:-1] + trimmed_cube = trimmed_cube.intersection( + iris.coords.CoordExtent( + coord, np.min(trimmed_coord), np.max(trimmed_coord) + ) + ) + return trimmed_cube + + +def get_boundary(cube, coord_system_out=None): + """ + Determines the boundary of the spatial grid of an :class:`iris.cube.Cube` \ + in a given coordinate system. + + Arguments: + + * **cube** - an :class:`iris.cube.Cube` object. + + Keyword arguments: + + * **coord_system_out** - an instance of \ + :class:`iris.coord_systems.CoordSystem` that \ + specifies the coordinate system the boundary is \ + to be computed in. If not supplied, the boundary \ + is determined in the same coordinate system as \ + the input cube. + + Returns + ------- + * **boundary** - a :class:`shapely.geometry.Polygon` object describing the \ + boundary of the cube spatial grid in the specified \ + coordinate system. + """ + coord_system_in = cube.coord_system() + if coord_system_out is None: + coord_system_out = coord_system_in + + # Check desired coordinate system is either GeogCS or RotatedGeogCS + if not isinstance( + coord_system_out, (iris.coord_systems.GeogCS, iris.coord_systems.RotatedGeogCS) + ): + msg = ("Unsupported coordinate system supplied {0:s}").format(coord_system_out) + raise NotImplementedError(msg) + + # Get x, y coordinates of cube spatial grid + x, y = iris.analysis.cartography.get_xy_grids(cube) + + if coord_system_out == coord_system_in: + # Nothing to do + pass + else: + # Transform the x, y coordinates of the cube spatial grid to the + # new coordinate system, coord_system_out + xy = coord_system_out.as_cartopy_crs().transform_points( + coord_system_in.as_cartopy_crs(), x, y + ) + x = xy[:, :, 0] + y = xy[:, :, 1] + + # Construct a Shapely Polygon object describing the boundary of the cube + # spatial grid by taking the edges of the x, y coordinate arrays + left = list(map(list, zip(*[x[:, 0], y[:, 0]], strict=False))) + top = list(map(list, zip(*[x[-1, :], y[-1, :]], strict=False))) + right = list( + map(list, zip(*[np.flipud(x[:, -1]), np.flipud(y[:, -1])], strict=False)) + ) + bottom = list( + map(list, zip(*[np.flipud(x[0, :]), np.flipud(y[0, :])], strict=False)) + ) + + boundary = left + top + right + bottom + boundary = sgeom.Polygon(boundary) + + return boundary + + +def mask_outside_boundary(cube, boundary, distance=None): + """ + Masks cube data outside a given spatial boundary. + + Arguments: + + * **cube** - an :class:`iris.cube.Cube` object. + * **boundary** - a :class:`shapely.geometry.Polygon` object specifying the \ + boundary of a spatial region. Data at points outside this \ + boundary will be masked. + + Keyword arguments: + + * **distance** - a factor to grow (distance > 0) or shrink (distance < 0) \ + the specified boundary by. + + Returns + ------- + * **masked_cube** - a new cube that is the same as the input \ + :class:`iris.cube.Cube` but with data masked at points \ + outside the supplied boundary. + """ + # Get spatial grid of input cube + grid_cube = get_grid(cube) + + # Get cube x, y coordinates of grid cube as a list of points + x, y = iris.analysis.cartography.get_xy_grids(grid_cube) + points = list(map(list, zip(x.ravel(), y.ravel(), strict=False))) + + # Grow or shrink boundary polygon by specified factor + if distance is not None: + boundary = boundary.buffer(distance) + + # Turn the boundary into a Path object + x, y = boundary.exterior.xy + x = np.asarray(x) + y = np.asarray(y) + boundary = list(map(list, zip(x.ravel(), y.ravel(), strict=False))) + + boundary = matplotlib.path.Path(boundary) + + # Find the points inside the boundary + points_inside_boundary = boundary.contains_points(points) + points_inside_boundary = points_inside_boundary.reshape(grid_cube.shape) + + # The array points_inside_boundary is boolean, with points inside the + # boundary set to True and those outside set to False. To use this as + # an array mask (see below), flip True <-> False (masked arrays have + # masked elements set to True) + points_inside_boundary = np.invert(points_inside_boundary) + + # Use points_inside_boundary as a mask on the data of the input cube + x_dim, y_dim = get_spatial_coord_dims(cube) + # Make a mask the same size as the input cube + mask = iris.util.broadcast_to_shape( + points_inside_boundary, cube.shape, [y_dim, x_dim] + ) + # Take a copy of the input cube for masking + masked_cube = cube.copy() + if isinstance(masked_cube.data, np.ma.MaskedArray): + # If the input cube already has a mask, combine it with the new mask + mask = np.ma.mask_or(mask, masked_cube.data.mask) + masked_cube.data.mask = mask + else: + # Mask the cube data + masked_cube.data = np.ma.array(masked_cube.data, mask=mask) + return masked_cube + + +def check_regrid_required(cube, grid_cube): + """ + Given two :class:`iris.cube.Cube` objects - a data cube and a cube \ + defining a desired spatial grid - this tests if their coordinate systems \ + and grid spacings are the same. If not, a regrid of the data cube would \ + be required to get it onto the spatial grid defined by the grid cube. + + Arguments: + + * **cube** - an :class:`iris.cube.Cube` object holding the data that \ + may need to be regridded. + * **grid_cube** - an :class:`iris.cube.Cube` object defining the desired \ + spatial grid. + + Returns + ------- + * **regrid_required** - False if the coordinate systems and grid \ + spacings of the data cube and grid cube are the \ + same. True otherwise. + """ + # Get the coordinate system and grid spacing of the input cube + cube_coord_system = cube.coord_system() + cube_grid_spacing = get_grid_spacing(cube) + + # Get the coordinate system and grid spacing of the cube defining the + # desired spatial grid + grid_coord_system = grid_cube.coord_system() + grid_spacing = get_grid_spacing(grid_cube) + + # Start by assuming no regrid would be required + regrid_required = False + + # Check if the coordinate systems of the cube and the desired grid + # differ. If so, a regrid would be required. + if cube_coord_system != grid_coord_system: + regrid_required = True + + # Check if the resolution of the cube and the desired grid differ. + # If so, a regrid would be required. + if not all(np.isclose(np.asarray(cube_grid_spacing), np.asarray(grid_spacing))): + regrid_required = True + + return regrid_required + + +def regrid_cube(cube, grid_cube, method="linear", mdtol=1): + """ + Regrids a cube onto a specified spatial grid. + + Arguments: + + * **cube** - an input :class:`iris.cube.Cube` to be regridded. + * **grid_cube** - an :class:`iris.cube.Cube` object defining the desired \ + spatial grid. + + Keyword arguments: + + * **method** - string specifying the regridding method to use. Can be \ + any of: + + * nearest - nearest-neighbour scheme :class:`iris.analysis.Nearest` \ + (non-conservative). + * linear - linear scheme :class:`iris.analysis.Linear` (non-conservative). + * area_weighted - area-weighted scheme \ + :class:`iris.analysis.AreaWeighted` (conservative). + * esmpy - regrid using ESMPy :func:`iris.experimental.regrid_conservative\ +.regrid_conservative_via_esmpy` (conservative). + * agg - regrid using Anti-Grain Geometry (conservative). + * two_stage - ANTS two-stage regridder (linear followed by \ + area-weighted; approximately conservative). + + * **mdtol** - Tolerance of missing data. The value returned in each \ + element of the returned array will be masked if the \ + fraction of missing data exceeds mdtol. This fraction \ + is calculated based on the area of masked cells within \ + each target cell. mdtol=0 means no masked data is \ + tolerated while mdtol=1 will mean the resulting \ + element will be masked if and only if all the overlapping \ + elements of the source grid are masked. Defaults to 1. \ + This is only used for method="AREA_WEIGHTED". + + Returns + ------- + * **result** - the input :class:`iris.cube.Cube` regridded onto the \ + desired spatial grid. + """ + # Guess spatial bounds on input cubes + cube = guess_bounds(cube) + grid_cube = guess_bounds(grid_cube) + + # Do regridding + t_start = time.time() + + if method == "nearest": + # Nearest-neighbour + print("Nearest-neighbour regridding...") + result = cube.regrid( + grid_cube, iris.analysis.Nearest(extrapolation_mode="mask") + ) + elif method == "linear": + # Linear + print("Linear regridding...") + result = cube.regrid(grid_cube, iris.analysis.Linear(extrapolation_mode="mask")) + elif method == "area_weighted": + # Area-weighted + # Requires cube and grid_cube to have the same coordinate system and + # all spatial coordinates must be rectilinear (1D) + print("Area-weighted regridding...") + result = cube.regrid(grid_cube, iris.analysis.AreaWeighted(mdtol=mdtol)) + elif method == "esmpy": + # ESMPy + print("Regridding via ESMPy...") + result = regrid_conservative_via_esmpy(cube, grid_cube) + elif method == "agg": + # AGG regrid + print("Regridding via agg regrid...") + result = cube.regrid(grid_cube, agg_regrid.AreaWeighted()) + elif method == "two_stage": + # ANTS two_stage (linear followed by area-weighted) regrid + with dask.config.set({"multiprocessing.context": "spawn"}): + print("Regridding via ANTS two-stage...") + result = cube.regrid(grid_cube, TwoStage(mdtol=mdtol)) + else: + raise ValueError("Unrecognised regrid method {0:s}".format(method)) + + # Timing information + t_end = time.time() + time_taken = t_end - t_start + print("Time taken {0:.2f} s".format(time_taken)) + + return result + + +def unrotate_pole(cube): + """ + Transforms the spatial coordinates of an input :class:`iris.cube.Cube` \ + object from a rotated pole coordinate system to standard \ + latitude/longitude coordinates. + """ + # Take a copy of the input cube + unrot_cube = cube.copy() + + # If the cube has a coordinate system with a rotated pole, do the + # unrotation + coord_system = unrot_cube.coord_system() + if isinstance(coord_system, iris.coord_systems.RotatedGeogCS): + # Get the latitude/longitude of the rotated pole + pole_lat = coord_system.grid_north_pole_latitude + pole_lon = coord_system.grid_north_pole_longitude + + # x, y coordinates in the rotated pole coordinate system + x, y = iris.analysis.cartography.get_xy_grids(unrot_cube) + + # Unrotate the x, y coordinates to standard lat/lon + x, y = iris.analysis.cartography.unrotate_pole(x, y, pole_lon, pole_lat) + + # Remove the existing spatial coordinates of the cube (in the + # rotated pole coordinate system)... + x_coord, y_coord = get_spatial_coords(unrot_cube) + x_dim, y_dim = get_spatial_coord_dims(unrot_cube) + unrot_cube.remove_coord(x_coord) + unrot_cube.remove_coord(y_coord) + + # ...and replace them with new standard lat/lon coordinates (2D) + target_coord_system = iris.coord_systems.GeogCS( + iris.fileformats.pp.EARTH_RADIUS + ) + x_coord = iris.coords.AuxCoord( + x, + standard_name="longitude", + units="degrees_east", + coord_system=target_coord_system, + ) + y_coord = iris.coords.AuxCoord( + y, + standard_name="latitude", + units="degrees_north", + coord_system=target_coord_system, + ) + unrot_cube.add_aux_coord(x_coord, (y_dim, x_dim)) + unrot_cube.add_aux_coord(y_coord, (y_dim, x_dim)) + + return unrot_cube + + +def rotate_pole(cube, pole_lat, pole_lon): + """ + Transforms the spatial coordinates of an :class:`iris.cube.Cube` object \ + from standard latitude/longitude coordinates to a rotated pole coordinate \ + system. + + Arguments: + + * **cube** - an input :class:`iris.cube.Cube` object. + * **pole_lat** - latitude of the rotated North pole. + * **pole_lon** - longitude of the rotated North pole. + """ + # Take a copy of the input cube + rot_cube = cube.copy() + + # If the cube has a standard latitude/longitude coordinate system, do the + # rotation + coord_system = rot_cube.coord_system() + if isinstance(coord_system, iris.coord_systems.GeogCS): + # Get the standard lat/lon x, y coordinates + x, y = iris.analysis.cartography.get_xy_grids(rot_cube) + + # Rotate the x, y coordinates to the rotated pole coordinate system + x, y = iris.analysis.cartography.rotate_pole(x, y, pole_lon, pole_lat) + + # Remove the existing spatial coordinates of the cube... + x_coord, y_coord = get_spatial_coords(rot_cube) + x_dim, y_dim = get_spatial_coord_dims(rot_cube) + rot_cube.remove_coord(x_coord) + rot_cube.remove_coord(y_coord) + + # ...and replace them with new (2D) coordinates in the rotated + # pole coordinate system + target_coord_system = iris.coord_systems.RotatedGeogCS( + pole_lat, pole_lon, ellipsoid=coord_system + ) + x_coord = iris.coords.AuxCoord( + x, + standard_name="grid_longitude", + units="degrees", + coord_system=target_coord_system, + ) + y_coord = iris.coords.AuxCoord( + y, + standard_name="grid_latitude", + units="degrees", + coord_system=target_coord_system, + ) + rot_cube.add_aux_coord(x_coord, (y_dim, x_dim)) + rot_cube.add_aux_coord(y_coord, (y_dim, x_dim)) + + return rot_cube + + +def remove_duplicates(cubelist): + """ + Removes any duplicate :class:`iris.cube.Cube` objects from an \ + :class:`iris.cube.CubeList`. + """ + # Nothing to do if the cubelist is empty + if not cubelist: + return cubelist + # Build up a list of indices of the cubes to remove because they are + # duplicated + indices_to_remove = [] + for i in range(len(cubelist) - 1): + cube_i = cubelist[i] + for j in range(i + 1, len(cubelist)): + cube_j = cubelist[j] + if cube_i == cube_j: + if j not in indices_to_remove: + indices_to_remove.append(j) + # Only keep unique cubes + cubelist = iris.cube.CubeList( + [cube for index, cube in enumerate(cubelist) if index not in indices_to_remove] + ) + return cubelist + + +def repeat_scalar_coord_along_dim_coord(cubelist, scalar_coord_name, dim_coord_name): + """ + For each :class:`iris.cube.Cube` in a given :class:`iris.cube.CubeList`, \ + this extends (by repetition) a specified scalar coordinate along the \ + dimension corresponding to a specified dimension coordinate. + """ + for cube in cubelist: + scalar_coord = cube.coord(scalar_coord_name) + # Check the coordinate referenced by scalar_coord_name is indeed + # a scalar coordinate. Otherwise nothing to do. + if scalar_coord.points.size == 1: + # Get the data value held by the scalar coordinate... + scalar_coord_val = scalar_coord.points[0] + # ...then remove it from the cube + cube.remove_coord(scalar_coord) + # Extract the dimension coordinate matching dim_coord_name + dim_coord = cube.coord(dim_coord_name) + # Get the dimension spanned by this dimension coordinate... + dim = cube.coord_dims(dim_coord_name)[0] + # ...and its length + dim_size = dim_coord.points.size + # Construct an auxillary coordinate by replicating the data + # value from the scalar coordinate to match the size of the + # specified dimension coordinate + scalar_coord = iris.coords.AuxCoord( + np.repeat(scalar_coord_val, dim_size), long_name=scalar_coord_name + ) + + # Add the new auxillary coordinate to the cube + cube.add_aux_coord(scalar_coord, dim) + + return cubelist + + +def extract_overlapping(cubelist, coord_name): + """ + Extracts regions from cubes in a :class:`iris.cube.CubeList` such that \ + the specified coordinate is the same across all cubes. + + Arguments: + + * **cubelist** - an input :class:`iris.cube.CubeList`. + * **coord_name** - a string specifying the name of the coordinate \ + over which to perform the extraction. + + Returns a :class:`iris.cube.CubeList` where the coordinate corresponding \ + to coord_name is the same for all cubes. + """ + # Build a list of all Cell instances for this coordinate by + # looping through all cubes in the supplied cubelist + all_cells = [] + for cube in cubelist: + for cell in cube.coord(coord_name).cells(): + all_cells.append(cell) + + # Work out which coordinate Cell instances are common across + # all cubes in the cubelist... + cell_counts = collections.Counter(all_cells) + # unique_cells = cell_counts.keys() + unique_cells = list(cell_counts.keys()) + unique_cell_counts = list(cell_counts.values()) + num_cubes = len(cubelist) + common_cells = [ + unique_cells[i] + for i, count in enumerate(unique_cell_counts) + if count == num_cubes + ] + # ...and use these to subset the cubes in the cubelist + constraint = iris.Constraint( + coord_values={coord_name: lambda cell: cell in common_cells} + ) + + cubelist = iris.cube.CubeList([cube.extract(constraint) for cube in cubelist]) + return cubelist + + +def cube_data_func(cube_data_1, cube_data_2): + """ + Takes two :class:`numpy.ndarray` arrays and combines them in quadrature. + """ + return np.sqrt(cube_data_1**2 + cube_data_2**2) + + +def cube_units_func(cube_1, cube_2): + """ + Compares the units of two :class:`iris.cube.Cube` objects to see if \ + they are the same. + """ + if cube_1.units != getattr(cube_2, "units", cube_1.units): + raise ValueError("Cube units do not match") + return cube_1.units + + +#: Function to compute the vector magnitude of two cubes (e.g. wind speed) +vector_magnitude = iris.analysis.maths.IFunc(cube_data_func, cube_units_func) diff --git a/src/CSET/cset_workflow/app/fetch_nimrod/bin/fetch-nimrod.py b/src/CSET/cset_workflow/app/fetch_nimrod/bin/fetch-nimrod.py new file mode 100755 index 000000000..34765e5dd --- /dev/null +++ b/src/CSET/cset_workflow/app/fetch_nimrod/bin/fetch-nimrod.py @@ -0,0 +1,768 @@ +#! /usr/bin/env python3 + +"""Retrieve UK radar observations. Specific to the Met Office.""" + +import datetime + +# import cube_utils +import time + +import dateutil +import iris +import iris.coords +import iris.cube + +# import agg_regrid +import numpy as np +from iris.experimental.regrid_conservative import regrid_conservative_via_esmpy + +# from ants.regrid.rectilinear import TwoStage +# from CSET._workflow_utils.fetch_data import FileRetrieverABC, fetch_nimrod +from CSET.cset_workflow.app.fetch_fcst.bin.fetch_data import ( + FileRetrieverABC, + fetch_nimrod, +) + + +def insert_datetime(filename, date_time): + """ + Insert a datetime. + + Inserts a datetime into a file name containing date formatting characters. + + Arguments: + + * **filename** - the name of a file. If this contains any of the special \ + date formatting characters + + * %Y - 4-digit year + * %m - 2-digit month + * %d - 2-digit day + * %H - 2-digit hour + * %M - 2-digit minute + + then these are replaced with numeric values derived from the components \ + of the supplied :class:`datetime.datetime` object. + * **date_time** - a :class:`datetime.datetime` object specifying the \ + datetime to insert in the given filename. + + Returns the input filename with date formatting characters replaced by \ + the appropriate components of date_time. + """ + filename = filename.replace("%Y", "{0:04d}".format(date_time.year)) + filename = filename.replace("%m", "{0:02d}".format(date_time.month)) + filename = filename.replace("%d", "{0:02d}".format(date_time.day)) + filename = filename.replace("%H", "{0:02d}".format(date_time.hour)) + filename = filename.replace("%M", "{0:02d}".format(date_time.minute)) + + return filename + + +def get_mean_period(time_coord): + """ + Get mean period. + + Use the bounds information of an input :class:`iris.coords.Coord` object \ + representing a time coordinate to determine a unique time period over \ + which to compute time means. + """ + # Check that the input time coordinate has bounds + if not time_coord.has_bounds(): + msg = ( + "Supplied time coordinate must have bounds to define time " + "periods over which to compute time means" + ) + raise ValueError(msg) + + # Check the input time intervals to make sure they specify a sensible + # meaning period + time_unit = time_coord.units + mean_period = [ + (end_time - start_time) + for start_time, end_time in time_unit.num2date(time_coord.bounds) + ] + mean_period = list(set(mean_period)) + if len(mean_period) != 1: + msg = ( + "Could not determine a unique time meaning period from " + "input time coordinate" + ) + raise ValueError(msg) + #: Number of seconds in an hour. + HOUR_IN_SECONDS = 3600 + mean_period = mean_period[0].total_seconds() / float(HOUR_IN_SECONDS) + + return mean_period + + +def guess_bounds(cube): + """ + Guess bounds. + + Takes an input :class:`iris.cube.Cube`, guesses bounds on the x, y \ + coordinates, then returns the cube. Such bounds are often required by \ + regridding algorithms. + """ + # Loop over spatial coordinates + for axis in ["x", "y"]: + coord = cube.coord(axis=axis) + # Check that this is not a variable resolution grid + # TODO: Does bounds really not work with variable resolution? AVD + # seem to think it won't... + try: + _ = iris.util.regular_step(coord) + except ValueError as err: + raise ValueError( + "Cannot guess bounds for a variable resolution grid" + ) from err + # Guess bounds if there aren't any + if coord.bounds is None: + coord.guess_bounds() + return cube + + +def get_spatial_coords(cube): + """ + Get spatial coords. + + Returns the x, y coordinates of an input :class:`iris.cube.Cube`. + """ + # Usual names for spatial coordinates + X_COORD_NAMES = ["longitude", "grid_longitude", "projection_x_coordinate", "x"] + Y_COORD_NAMES = ["latitude", "grid_latitude", "projection_y_coordinate", "y"] + + # Get a list of coordinate names for the cube + coord_names = [coord.name() for coord in cube.coords()] + + # Check which x-coordinate we have, if any + x_coords = [coord for coord in coord_names if coord in X_COORD_NAMES] + if len(x_coords) != 1: + raise ValueError("Could not identify a unique x-coordinate in cube") + x_coord = cube.coord(x_coords[0]) + + # Check which y-coordinate we have, if any + y_coords = [coord for coord in coord_names if coord in Y_COORD_NAMES] + if len(y_coords) != 1: + raise ValueError("Could not identify a unique y-coordinate in cube") + y_coord = cube.coord(y_coords[0]) + + return [x_coord, y_coord] + + +def get_spatial_coord_dims(cube): + """ + Get spatial coords dims. + + Returns a tuple of the dimensions corresponding to the x, y coordinates \ + of an input :class:`iris.cube.Cube`. + """ + # Get the spatial coordinates of the cube... + x_coord, y_coord = get_spatial_coords(cube) + # ...then the cube dimensions corresponding to these coordinates + x_dim = cube.coord_dims(x_coord)[0] + y_dim = cube.coord_dims(y_coord)[0] + return (x_dim, y_dim) + + +def get_grid_spacing(cube): + """ + Grid spacing. + + Takes a :class:`iris.cube.Cube` and returns a tuple (dx, dy) of the grid \ + spacings in the x, y directions. + """ + # Get the spatial coordinates of the cube + x_coord, y_coord = get_spatial_coords(cube) + # Note that regular_step will fail for variable resolution grids + dx = iris.util.regular_step(x_coord) + dy = iris.util.regular_step(y_coord) + return (dx, dy) + + +def check_regrid_required(cube, grid_cube): + """ + Check grid. + + Given two :class:`iris.cube.Cube` objects - a data cube and a cube \ + defining a desired spatial grid - this tests if their coordinate systems \ + and grid spacings are the same. If not, a regrid of the data cube would \ + be required to get it onto the spatial grid defined by the grid cube. + + Arguments: + + * **cube** - an :class:`iris.cube.Cube` object holding the data that \ + may need to be regridded. + * **grid_cube** - an :class:`iris.cube.Cube` object defining the desired \ + spatial grid. + + Returns + ------- + * **regrid_required** - False if the coordinate systems and grid \ + spacings of the data cube and grid cube are the \ + same. True otherwise. + """ + # Get the coordinate system and grid spacing of the input cube + cube_coord_system = cube.coord_system() + cube_grid_spacing = get_grid_spacing(cube) + + # Get the coordinate system and grid spacing of the cube defining the + # desired spatial grid + grid_coord_system = grid_cube.coord_system() + grid_spacing = get_grid_spacing(grid_cube) + + # Start by assuming no regrid would be required + regrid_required = False + + # Check if the coordinate systems of the cube and the desired grid + # differ. If so, a regrid would be required. + if cube_coord_system != grid_coord_system: + regrid_required = True + + # Check if the resolution of the cube and the desired grid differ. + # If so, a regrid would be required. + if not all(np.isclose(np.asarray(cube_grid_spacing), np.asarray(grid_spacing))): + regrid_required = True + + return regrid_required + + +def regrid_cube(cube, grid_cube, method="linear", mdtol=1): + """ + Regrid a cube onto a specified spatial grid. + + Arguments: + + * **cube** - an input :class:`iris.cube.Cube` to be regridded. + * **grid_cube** - an :class:`iris.cube.Cube` object defining the desired \ + spatial grid. + + Keyword arguments: + + * **method** - string specifying the regridding method to use. Can be \ + any of: + + * nearest - nearest-neighbour scheme :class:`iris.analysis.Nearest` \ + (non-conservative). + * linear - linear scheme :class:`iris.analysis.Linear` (non-conservative). + * area_weighted - area-weighted scheme \ + :class:`iris.analysis.AreaWeighted` (conservative). + * esmpy - regrid using ESMPy :func:`iris.experimental.regrid_conservative\ +.regrid_conservative_via_esmpy` (conservative). + * agg - regrid using Anti-Grain Geometry (conservative). + * two_stage - ANTS two-stage regridder (linear followed by \ + area-weighted; approximately conservative). + + * **mdtol** - Tolerance of missing data. The value returned in each \ + element of the returned array will be masked if the \ + fraction of missing data exceeds mdtol. This fraction \ + is calculated based on the area of masked cells within \ + each target cell. mdtol=0 means no masked data is \ + tolerated while mdtol=1 will mean the resulting \ + element will be masked if and only if all the overlapping \ + elements of the source grid are masked. Defaults to 1. \ + This is only used for method="AREA_WEIGHTED". + + Returns + ------- + * **result** - the input :class:`iris.cube.Cube` regridded onto the \ + desired spatial grid. + """ + # Guess spatial bounds on input cubes + cube = guess_bounds(cube) + grid_cube = guess_bounds(grid_cube) + + # Do regridding + t_start = time.time() + + if method == "nearest": + # Nearest-neighbour + print("Nearest-neighbour regridding...") + result = cube.regrid( + grid_cube, iris.analysis.Nearest(extrapolation_mode="mask") + ) + elif method == "linear": + # Linear + print("Linear regridding...") + result = cube.regrid(grid_cube, iris.analysis.Linear(extrapolation_mode="mask")) + elif method == "area_weighted": + # Area-weighted + # Requires cube and grid_cube to have the same coordinate system and + # all spatial coordinates must be rectilinear (1D) + print("Area-weighted regridding...") + result = cube.regrid(grid_cube, iris.analysis.AreaWeighted(mdtol=mdtol)) + elif method == "esmpy": + # ESMPy + print("Regridding via ESMPy...") + result = regrid_conservative_via_esmpy(cube, grid_cube) + # elif method == "agg": + # # AGG regrid + # print("Regridding via agg regrid...") + # result = cube.regrid(grid_cube, agg_regrid.AreaWeighted()) + # elif method == "two_stage": + # # ANTS two_stage (linear followed by area-weighted) regrid + # with dask.config.set({"multiprocessing.context": "spawn"}): + # print("Regridding via ANTS two-stage...") + # result = cube.regrid(grid_cube, TwoStage(mdtol=mdtol)) + else: + raise ValueError("Unrecognised regrid method {0:s}".format(method)) + + # Timing information + t_end = time.time() + time_taken = t_end - t_start + print("Time taken {0:.2f} s".format(time_taken)) + + return result + + +def get_spatial_extent(cube): + """ + Get spatial extent. + + Takes an :class:`iris.cube.Cube` and returns its spatial extent in the \ + form of a list [min x, max x, min y, max y]. + """ + # Get the spatial coordinates of the cube + x_coord, y_coord = get_spatial_coords(cube) + return [ + min(x_coord.points), + max(x_coord.points), + min(y_coord.points), + max(y_coord.points), + ] + + +def trim_cube(cube, xmin=None, xmax=None, ymin=None, ymax=None, ignore_bounds=True): + """ + Trim a cube. + + Trims the spatial extent of an :class:`iris.cube.Cube` to lie in a \ + specified range. + + Arguments: + + * **cube** - an :class:`iris.cube.Cube` object. + + Keyword arguments: + + * **xmin** - minimum value of the cube x-coordinate to select. + * **xmax** - maximum value of the cube x-coordinate to select. + * **ymin** - minimum value of the cube y-coordinate to select. + * **ymax** - maximum value of the cube y-coordinate to select. + * **ignore_bounds** - set to True to ignore any bounds on the coordinates \ + and just consider their points. + + Returns + ------- + * **trimmed_cube** - a copy of the input :class:`iris.cube.Cube` with \ + x, y coordinate values restricted to the specified \ + ranges. + """ + REAL_MATCH_TOL = 1.0e-4 + + # Take a copy of the input cube + trimmed_cube = cube.copy() + + # Get the spatial coordinates of the cube + x_coord, y_coord = get_spatial_coords(trimmed_cube) + + # First trim the cube in the x direction... + if xmin is None: + xmin = np.min(x_coord.points) + xmin -= REAL_MATCH_TOL + if xmax is None: + xmax = np.max(x_coord.points) + xmax += REAL_MATCH_TOL + x_extent = iris.coords.CoordExtent(x_coord, xmin, xmax) + + trimmed_cube = trimmed_cube.intersection(x_extent, ignore_bounds=ignore_bounds) + + # ...now trim the cube in the y direction + if ymin is None: + ymin = np.min(y_coord.points) + ymin -= REAL_MATCH_TOL + if ymax is None: + ymax = np.max(y_coord.points) + ymax += REAL_MATCH_TOL + y_extent = iris.coords.CoordExtent(y_coord, ymin, ymax) + trimmed_cube = trimmed_cube.intersection(y_extent, ignore_bounds=ignore_bounds) + + return trimmed_cube + + +def get_grid(cube): + """ + Get grid. + + Takes an :class:`iris.cube.Cube` and returns its spatial grid as a 2D \ + :class:`iris.cube.Cube`. + """ + # Take a copy of the input cube + grid_cube = cube.copy() + # Get the spatial coordinates of the cube + x_coord, y_coord = get_spatial_coords(grid_cube) + + # Take an arbitrary spatial slice of the cube + grid_cube = next(grid_cube.slices([y_coord, x_coord])) + # We don't care about the data of this slice so just zero it + grid_cube.data = np.zeros(grid_cube.data.shape) + # Add a mask to show which points in the spatial grid are to be + # considered valid. By default, all of them (False means not masked in + # masked arrays) + mask = np.zeros(grid_cube.data.shape, dtype=bool) + grid_cube.data = np.ma.array(grid_cube.data, mask=mask) + # Remove all non-spatial coordinates from the cube + coords = grid_cube.coords() + + coords.remove(x_coord) + coords.remove(y_coord) + + for coord in coords: + grid_cube.remove_coord(coord) + + # Set some metadata to signify this is just a grid cube + grid_cube.standard_name = None + grid_cube.long_name = "Grid" + grid_cube.units = "no_unit" + grid_cube.cell_methods = None + grid_cube.attributes = {} + + return grid_cube + + +class NimrodRetriever(FileRetrieverABC): + """ + Class for Nimrod. + + Class for UK precipitation radar observations, a subclass of \ + :class:`FileRetrieverABC`. + """ + + def __init__( + self, + name, + data_dir, + filename, + weights_filename, + temporal_range, + regrid_method="linear", + field_name=None, + ): + """ + Initialise class. + + Arguments: + + * **name** - name for this observational dataset. + * **data_dir** - directory where UK radar 1-hourly accumulation data \ + files are located. + * **filename** - names of the UK radar 1-hourly accumulation data \ + files. Can contain any of the following special date \ + formatting characters: + + * %Y - 4-digit year + * %m - 2-digit month + * %d - 2-digit day + * %H - 2-digit hour + * %M - 2-digit minute + + * **weights_filename** - names of the files holding the corresponding \ + weights for the 1-hour radar accumulations. \ + Can contain any of the following special date \ + formatting characters: + + * %Y - 4-digit year + * %m - 2-digit month + * %d - 2-digit day + * %H - 2-digit hour + * %M - 2-digit minute + + * **temporal_range** - a list of the form [first_time, last_time] \ + specifying the time range covered by available \ + UK radar data. first_time and last_time can be \ + either :class:`datetime.datetime` objects or \ + None for an open-ended range. + + Keyword arguments: + + * **field_name** - name of the field to extract from data files for \ + this dataset. + """ + FileRetrieverABC.__init__( + self, + name, + data_dir, + filename, + temporal_range, + regrid_method=regrid_method, + field_name=field_name, + ) + self.weights_filename = weights_filename + # UK radar accumulations have an hourly frequency + self.__frequency_in_hours = 1.0 + # The minimum number of valid rainrate values there must be in an + # hourly accumulation for it to be used (the maximum possible number + # is 13) + self.__min_weight = 11 + + def get_native_grid(self, grid_shift=None): + """ + Return the spatial grid. + + Returns the spatial grid this observational dataset is defined on \ + as an :class:`iris.cube.Cube` object. + + Keyword arguments: + + * **grid_shift** * - Which grid to shift to + + .. warning :: + grid_shift is not implemented for this class and must be set to None + + Raises + ------ + NotImplementedError: Raised if grid_shift is set to anything other than None + """ + if grid_shift is not None: + raise NotImplementedError("grid_shift is not implemented for %s") + + # Choose a random data file to extract the grid from + sample_file = self.get_sample_file() + if sample_file is None: + raise ValueError("Cannot determine grid without a sample data file") + + # Load sample data + if self.field_name is None: + cube = iris.load_cube(sample_file) + else: + cube = iris.load_cube(sample_file, self.field_name) + + # Extract spatial grid + + grid_cube = get_grid(cube) + + return grid_cube + + def get_data(self, time_coord, grid_cube=None, cutout=None): + """ + Compute mean UK radar precipitation rates over a given time period. + + Arguments: + + * **time_coord** - a :class:`iris.coords.Coord` instance defining the \ + time periods to compute mean precipitation rates \ + for. The coordinate must have bounds as these are \ + used to compute the time meaning period to use. + + Keyword arguments: + + * **grid_cube** - an :class:`iris.cube.Cube` object defining a spatial \ + grid to regrid UK radar data onto. + * **cutout** - a list specifying the longitude/latitude extent of a \ + subregion to extract from UK radar data, of the form \ + [min lon, max lon, min lat, max lat]. + + Returns + ------- + * **uk_radar_mean_cube** - an :class:`iris.cube.Cube` object holding \ + time-meaned UK radar precipitation rates, \ + with a meaning period derived from the input \ + time coordinate. + """ + # Work out the period over which to compute time means + mean_period = get_mean_period(time_coord) + + # The meaning period must be an integer multiple of an hour to use + # UK radar hourly accumulations + if mean_period % self.__frequency_in_hours != 0: + msg = ( + "Meaning period must be an integer multiple of an hour " + "to use UK radar data. Currently {0:f} hours".format(mean_period) + ) + raise ValueError(msg) + # Number of hourly periods in this meaning period + num_uk_radar_periods = int(mean_period / self.__frequency_in_hours) + + # Names of UK radar data files... + uk_radar_filename = "/".join([self.data_dir, self.filename]) + # ...and associated weights + uk_radar_weights_filename = "/".join([self.data_dir, self.weights_filename]) + + # Determine the start and end points of the desired time range + time_unit = time_coord.units + first_time = time_unit.num2date( + np.min(time_coord.bounds) + self.__frequency_in_hours + ) + last_time = time_unit.num2date(np.max(time_coord.bounds)) + # Load all required data + uk_radar_cubes = iris.cube.CubeList() + uk_radar_weights_cubes = iris.cube.CubeList() + for time_loop in dateutil.rrule.rrule( + dateutil.rrule.HOURLY, + interval=1, + dtstart=datetime.datetime.fromisoformat(first_time.isoformat()), + until=datetime.datetime.fromisoformat(last_time.isoformat()), + ): + # Load hourly accumulation data + uk_radar_file = insert_datetime(uk_radar_filename, time_loop) + uk_radar_cube = iris.load_cube(uk_radar_file) + # Adjust the time coordinate to be the midpoint of the hourly-mean + # period, which is how times are typically represented in Iris + shifted_time_coord = uk_radar_cube.coord("time").copy( + points=(uk_radar_cube.coord("time").points - 0.5), + bounds=uk_radar_cube.coord("time").bounds, + ) + uk_radar_cube.replace_coord(shifted_time_coord) + + # Load corresponding accumulation weights + uk_radar_weights_file = insert_datetime( + uk_radar_weights_filename, time_loop + ) + uk_radar_weights_cube = iris.load_cube(uk_radar_weights_file) + # Adjust the time coordinate to be the midpoint of the hourly-mean + # period, which is how times are represented in Iris + shifted_time_coord = uk_radar_weights_cube.coord("time").copy( + points=(uk_radar_weights_cube.coord("time").points - 0.5), + bounds=uk_radar_weights_cube.coord("time").bounds, + ) + uk_radar_weights_cube.replace_coord(shifted_time_coord) + + # TODO: Do we want to mask out values below self.__min_weight + # for each hourly accumulation, rather than after they have been + # combined to make an X-hourly mean as below? + # For example, imagine if we are producing 6-hourly accumulations + # and the weights are 12, 12, 12, 12, 12, 7. The average weight + # over the 6 hour period is then 11.17 which exceeds + # self.__min_weight=11. However, the last rainrate value was + # clearly dodgy (weight=7) so we may want to mask out this pixel + # despite having 5 good rainrate values (weight=12). + # Mask out data where the accumulation weights are below the + # acceptable threshold + # uk_radar_cube.data = np.ma.masked_where( + # uk_radar_weights_cube.data < self.__min_weight, + # uk_radar_cube.data) + + # Remove all attributes from the radar cubes to avoid future merge problems. + # See https://code.metoffice.gov.uk/trac/rmedtoolbox/ticket/152#comment:24 + uk_radar_cube.attributes = None + uk_radar_weights_cube.attributes = None + + # Store data for this hour + uk_radar_cubes.append(uk_radar_cube) + uk_radar_weights_cubes.append(uk_radar_weights_cube) + + # Merge all hourly radar accumulations and accumulation weights + # into single cubes + uk_radar_cube = uk_radar_cubes.merge_cube() + uk_radar_weights_cube = uk_radar_weights_cubes.merge_cube() + + # Check that the number of times in the cube is a multiple of + # num_uk_radar_periods, otherwise we cannot compute the mean over + # the desired period + num_times = uk_radar_cube.coord("time").points.size + if num_times % num_uk_radar_periods != 0: + msg = "Cannot compute {0:.1f}-hour mean due to missing times".format( + mean_period + ) + raise ValueError(msg) + + # Convert units from mm*32 to mm + uk_radar_cube.data = uk_radar_cube.data.astype(np.float64) / 32.0 + uk_radar_cube.units = "mm" + + # Now compute the mean over the desired meaning period using a + # rolling window (in mm/h) + uk_radar_mean_data = iris.util.rolling_window( + uk_radar_cube.data, + window=num_uk_radar_periods, + step=num_uk_radar_periods, + axis=0, + ) + uk_radar_mean_data = np.sum( + uk_radar_mean_data, axis=1, dtype=np.float64 + ) / float(mean_period) + + # Do the same for the accumulation weights + uk_radar_mean_weights = iris.util.rolling_window( + uk_radar_weights_cube.data, + window=num_uk_radar_periods, + step=num_uk_radar_periods, + axis=0, + ) + uk_radar_mean_weights = np.sum( + uk_radar_mean_weights, axis=1, dtype=np.float64 + ) / float(mean_period) + + # Mask out data where the accumulation weights, averaged over the + # desired meaning period, are below the threshold + mask = uk_radar_mean_weights < self.__min_weight + uk_radar_mean_data = np.ma.array(uk_radar_mean_data, mask=mask) + + # Create an Iris cube to hold the UK radar data meaned over the + # desired period + x_coord, y_coord = get_spatial_coords(uk_radar_cube) + x_dim, y_dim = get_spatial_coord_dims(uk_radar_cube) + cell_method = iris.coords.CellMethod( + method="mean", coords="time", intervals="{0:.1f} hour".format(mean_period) + ) + uk_radar_mean_cube = iris.cube.Cube( + uk_radar_mean_data, + standard_name=uk_radar_cube.standard_name, + units="mm h-1", + dim_coords_and_dims=[(time_coord, 0), (y_coord, y_dim), (x_coord, x_dim)], + cell_methods=(cell_method,), + ) + + # Regrid UK radar data from its native grid onto the specified grid + if grid_cube is not None: + # Check if regridding is actually required + regrid_required = check_regrid_required(uk_radar_mean_cube, grid_cube) + if regrid_required: + # Do the regrid + uk_radar_mean_cube = regrid_cube( + uk_radar_mean_cube, grid_cube, method=self.regrid_method + ) + else: + # If both the coordinate system and resolution of the UK + # radar data cube and the desired grid are the same, no need + # to do any regridding, just take a cutout of the radar data + # cube instead + xmin, xmax, ymin, ymax = get_spatial_extent(grid_cube) + uk_radar_mean_cube = trim_cube( + uk_radar_mean_cube, xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax + ) + + # Use any mask associated with the grid cube to mask the UK + # radar data + if isinstance(grid_cube.data, np.ma.MaskedArray): + mask = grid_cube.data.mask + # Make a mask the same size as the UK radar data cube + mask = iris.util.broadcast_to_shape( + mask, uk_radar_mean_cube.shape, [y_dim, x_dim] + ) + if isinstance(uk_radar_mean_cube.data, np.ma.MaskedArray): + # If the UK radar data cube already has a mask, combine it + # with the mask from the grid cube + mask = np.ma.mask_or(mask, uk_radar_mean_cube.data.mask) + uk_radar_mean_cube.data.mask = mask + else: + # Mask the UK radar data using the mask from the grid cube + uk_radar_mean_cube.data = np.ma.array( + uk_radar_mean_cube.data, mask=mask + ) + + # If required, take a spatial cutout + # This can be helpful otherwise we can end up with large cubes! + if cutout is not None: + xmin, xmax, ymin, ymax = cutout + uk_radar_mean_cube = trim_cube( + uk_radar_mean_cube, xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax + ) + return uk_radar_mean_cube + + def get_file(self, file_path: str, output_dir: str) -> bool: + """Save a file from the filesystem to the output directory.""" + return True + + +# fetch_nimrod(NimrodRetriever()) +fetch_nimrod() diff --git a/src/CSET/cset_workflow/app/fetch_nimrod/rose-app.conf b/src/CSET/cset_workflow/app/fetch_nimrod/rose-app.conf new file mode 100644 index 000000000..187112a4d --- /dev/null +++ b/src/CSET/cset_workflow/app/fetch_nimrod/rose-app.conf @@ -0,0 +1,3 @@ +[command] +default=echo "Please set ROSE_APP_COMMAND_KEY to your storage system."; false +fetchnimrod=app_env_wrapper fetch-nimrod.py diff --git a/src/CSET/cset_workflow/flow.cylc b/src/CSET/cset_workflow/flow.cylc index 69b6a1f4f..b0bb4edc9 100644 --- a/src/CSET/cset_workflow/flow.cylc +++ b/src/CSET/cset_workflow/flow.cylc @@ -38,14 +38,14 @@ final cycle point = {{CSET_TRIAL_END_DATE}} # Runs for every forecast initiation time to process the data in parallel. {% for date in CSET_CASE_DATES %} R1/{{date}} = """ - setup_complete[^] => FETCH_DATA:succeed-all => fetch_complete + setup_complete[^] => FETCH_DATA:succeed-all => FETCH_NIMROD:succeed-all => fetch_complete fetch_complete & parbake_recipes => bake_recipes => cycle_complete """ {% endfor %} {% elif CSET_CYCLING_MODE == "trial" %} # Analysis from each forecast. {{CSET_TRIAL_CYCLE_PERIOD}} = """ - setup_complete[^] => FETCH_DATA:succeed-all => fetch_complete + setup_complete[^] => FETCH_DATA:succeed-all => FETCH_NIMROD:succeed-all => fetch_complete fetch_complete & parbake_recipes => bake_recipes => cycle_complete """ {% endif %} @@ -97,6 +97,12 @@ final cycle point = {{CSET_TRIAL_END_DATE}} [[[environment]]] ANALYSIS_LENGTH = {{ANALYSIS_LENGTH}} + [[FETCH_NIMROD]] + script = rose task-run -v --app-key=fetch_nimrod + execution time limit = PT1H + [[[environment]]] + ANALYSIS_LENGTH = {{ANALYSIS_LENGTH}} + [[METPLUS]] [[[environment]]] {% if METPLUS_GRID_STAT|default(False) %} @@ -118,6 +124,9 @@ final cycle point = {{CSET_TRIAL_END_DATE}} [[fetch_complete]] inherit = DUMMY_TASK + [[dummy_fetch_nimrod]] + inherit = DUMMY_TASK, FETCH_NIMROD + [[cycle_complete]] inherit = DUMMY_TASK @@ -179,6 +188,17 @@ final cycle point = {{CSET_TRIAL_END_DATE}} [[[environment]]] DO_CASE_AGGREGATION = True + {% if NIMROD_RADAR_OBS %} + [[fetch_nimrod]] + # Fetch Nimrod radar observations. + inherit = FETCH_NIMROD + [[[environment]]] + ROSE_APP_COMMAND_KEY = "fetchnimrod" + NIMROD_FIELDS_COMP_HOUR = {{NIMROD_FIELDS_COMP_HOUR}} + NIMROD_FIELDS_1KM = {{NIMROD_FIELDS_1KM}} + NIMROD_FIELDS_2KM = {{NIMROD_FIELDS_2KM}} + {% endif %} + [[housekeeping]] # Housekeep input data files. [[[environment]]] diff --git a/src/CSET/cset_workflow/meta/rose-meta.conf b/src/CSET/cset_workflow/meta/rose-meta.conf index 99a67b25a..c4c30957a 100644 --- a/src/CSET/cset_workflow/meta/rose-meta.conf +++ b/src/CSET/cset_workflow/meta/rose-meta.conf @@ -1,6 +1,91 @@ # Diagnostics settings are split into a separate file. import=meta/diagnostics +################################################################################ +# Observations +################################################################################ + +[Observations] +ns=Observations +sort-key=sec-d +title=Observational Comparisons + +################################################################################ +# Spatial (2D map) fields +################################################################################ + +[Observations/Types] +ns=Observations/Types +sort-key=sec-d1 +title=UK radar observations + +# Surface Nimrod radar rainfall fields. +[template variables=NIMROD_RADAR_OBS] +ns=Observations/Types +description=Use UK Nimrod radar observations +type=python_boolean +compulsory=true +sort-key=nimrod-radar-1 +trigger=template variables=NIMROD_FIELDS_COMP_HOUR: True; + template variables=NIMROD_FIELDS_1KM: True; + template variables=NIMROD_FIELDS_2KM: True; +# template variables=USE_WMO_STATION_NUMBERS: True; +# template variables=SURFACE_SYNOP_INTERVAL: True; +# template variables=SURFACE_SYNOP_OFFSET: True; + +#[Observations/Types] +#ns=Observations/Types +#sort-key=sec-d2 +#title=UKNimrod observations + +#[template variables=NIMROD_FIELDS] +#ns=Observations/Types +#title=UKNimrod radar rainfall products +#description=UK Nimrod rainfall products +#help=Name of field of the UKNimrod rainfall product. +#compulsory=true +#type=python_list +#values="rainaccum_comp_hour", "rainaccum_comp_hour_1km_cutout", "rainaccum_comp_hour_2km" +#value-titles=rainaccum_comp_hour, rainaccum_comp_hour_1km_cutout, rainaccum_comp_hour_2km +#sort-key=nimrod-radar-2 + +[template variables=NIMROD_FIELDS_COMP_HOUR] +ns=Observations/Types +title=UK Nimrod radar accumulated rainfall +description=Default UK Nimrod rainfall accumulation +help=Switch to select default gridded UK Nimrod rainfall accumulation. +compulsory=true +#type=python_list +values="rainaccum_comp_hour", "" +value-titles=rainaccum_comp_hour, off +sort-key=nimrod-radar-3 + +[template variables=NIMROD_FIELDS_1KM] +ns=Observations/Types +title=UK Nimrod radar accumulated rainfall 1km gridding +description=1km gridded UK Nimrod rainfall accumulation +help=Switch to select 1km gridded UK Nimrod rainfall accumulation. +compulsory=true +#type=python_list +values="rainaccum_comp_hour_1km_cutout", "" +value-titles=rainaccum_comp_hour_1km_cutout, off +sort-key=nimrod-radar-4 + +[template variables=NIMROD_FIELDS_2KM] +ns=Observations/Types +title=UK Nimrod radar accumulated rainfall 2km gridding +description=2km gridded UK Nimrod rainfall accumulation +help=Switch to select 2km gridded UK Nimrod rainfall accumulation. +compulsory=true +#type=python_list +values="rainaccum_comp_hour_2km", "" +value-titles=rainaccum_comp_hour_2km, off +sort-key=nimrod-radar-5 + + +################################################################ + + [ns=Setup] sort-key=sec-a title=General setup options diff --git a/src/CSET/loaders/spatial_field.py b/src/CSET/loaders/spatial_field.py index b248b982b..bec0f1655 100644 --- a/src/CSET/loaders/spatial_field.py +++ b/src/CSET/loaders/spatial_field.py @@ -24,6 +24,27 @@ def load(conf: Config): # Load a list of model detail dictionaries. models = get_models(conf.asdict()) + # Surface (2D) fields for radar. + if conf.SPATIAL_SURFACE_FIELD: + # for model, field, method in itertools.product( + for model, field in itertools.product( + models, conf.SURFACE_FIELDS, conf.SPATIAL_SURFACE_FIELD_METHOD + ): + yield RawRecipe( + recipe="generic_surface_spatial_plot_sequence_radar.yaml", + model_ids=model["id"], + variables={ + "VARNAME": field, + # "MODEL_NAME": model["name"], + # "METHOD": method, + # "SUBAREA_TYPE": conf.SUBAREA_TYPE if conf.SELECT_SUBAREA else None, + # "SUBAREA_EXTENT": conf.SUBAREA_EXTENT + # if conf.SELECT_SUBAREA + # else None, + }, + aggregation=False, + ) + # Surface (2D) fields. if conf.SPATIAL_SURFACE_FIELD: for model, field, method in itertools.product( diff --git a/src/CSET/operators/_colorbar_definition.json b/src/CSET/operators/_colorbar_definition.json index 86331e92e..f79dd9803 100644 --- a/src/CSET/operators/_colorbar_definition.json +++ b/src/CSET/operators/_colorbar_definition.json @@ -1,4 +1,11 @@ { + "Hourly rain accumulation": { + "cmap": "cividis", + "max": 256, + "min": 0, + "ymax": 1.0, + "ymin": 0.0 + }, "air_potential_temperature": { "cmap": "RdYlBu_r", "max": 350, diff --git a/src/CSET/recipes/surface_fields/generic_surface_spatial_plot_sequence_radar.yaml b/src/CSET/recipes/surface_fields/generic_surface_spatial_plot_sequence_radar.yaml new file mode 100644 index 000000000..1dd27733d --- /dev/null +++ b/src/CSET/recipes/surface_fields/generic_surface_spatial_plot_sequence_radar.yaml @@ -0,0 +1,35 @@ +category: Surface Spatial Plot +#title: $MODEL_NAME $VARNAME $METHOD +#description: Extracts and plots the $METHOD of $VARNAME from all times in $MODEL_NAME. +title: $VARNAME +description: Extracts and plots radar data. + +steps: + - operator: read.read_cube + file_paths: $INPUT_PATHS + + constraint: + operator: constraints.combine_constraints + varname_constraint: + operator: constraints.generate_var_constraint + varname: $VARNAME +# cell_methods_constraint: +# operator: constraints.generate_cell_methods_constraint +# cell_methods: [] +# varname: $VARNAME +# level_constraint: +# operator: constraints.generate_level_constraint +# coordinate: "pressure" +# levels: [] +# subarea_type: $SUBAREA_TYPE +# subarea_extent: $SUBAREA_EXTENT + +# - operator: collapse.collapse +# coordinate: [time] +# method: $METHOD + + - operator: plot.spatial_pcolormesh_plot + sequence_coordinate: time + + - operator: write.write_cube_to_nc + overwrite: True