Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
6 changes: 2 additions & 4 deletions src/CSET/operators/plot.py
Original file line number Diff line number Diff line change
Expand Up @@ -1493,10 +1493,9 @@ def _plot_and_save_postage_stamp_power_spectrum_series(
plt.subplot(grid_size, grid_size, subplot)

frequency = member.coord("frequency").points
power_spectrum = member.data

ax = plt.gca()
ax.plot(frequency, power_spectrum[0])
ax.plot(frequency, member.data)
ax.set_title(f"Member #{member.coord(stamp_coordinate).points[0]}")

# Overall figure title.
Expand All @@ -1521,10 +1520,9 @@ def _plot_and_save_postage_stamps_in_single_plot_power_spectrum_series(
# Loop over all slices along the stamp_coordinate
for member in cube.slices_over(stamp_coordinate):
frequency = member.coord("frequency").points
power_spectrum = member.data
ax.plot(
frequency,
power_spectrum[0],
member.data,
label=f"Member #{member.coord(stamp_coordinate).points[0]}",
)

Expand Down
144 changes: 42 additions & 102 deletions src/CSET/operators/read.py
Original file line number Diff line number Diff line change
Expand Up @@ -34,7 +34,11 @@

from CSET._common import iter_maybe
from CSET.operators._stash_to_lfric import STASH_TO_LFRIC
from CSET.operators._utils import get_cube_yxcoordname
from CSET.operators._utils import (
get_cube_coordindex,
get_cube_yxcoordname,
is_spatialdim,
)


class NoDataError(FileNotFoundError):
Expand Down Expand Up @@ -187,10 +191,14 @@ def read_cubes(

# Select sub region.
cubes = _cutout_cubes(cubes, subarea_type, subarea_extent)

# Merge and concatenate cubes now metadata has been fixed.
cubes = cubes.merge()
cubes = cubes.concatenate()

# Squeeze single valued coordinates into scalar coordinates.
cubes = iris.cube.CubeList(iris.util.squeeze(cube) for cube in cubes)

# Ensure dimension coordinates are bounded.
for cube in cubes:
for dim_coord in cube.coords(dim_coords=True):
Expand All @@ -213,18 +221,10 @@ def _load_model(
input_files = _check_input_files(paths)
# If unset, a constraint of None lets everything be loaded.
logging.debug("Constraint: %s", constraint)
cubes = iris.load(
input_files, constraint, callback=_create_callback(is_ensemble=False)
)
cubes = iris.load(input_files, constraint, callback=_loading_callback)
# Make the UM's winds consistent with LFRic.
_fix_um_winds(cubes)

# Reload with ensemble handling if needed.
if _is_ensemble(cubes):
cubes = iris.load(
input_files, constraint, callback=_create_callback(is_ensemble=True)
)

# Add model_name attribute to each cube to make it available at any further
# step without needing to pass it as function parameter.
if model_name is not None:
Expand Down Expand Up @@ -356,84 +356,28 @@ def _cutout_cubes(
return cutout_cubes


def _is_ensemble(cubelist: iris.cube.CubeList) -> bool:
"""Test if a CubeList is likely to be ensemble data.

If cubes either have a realization dimension, or there are multiple files
for the same time-step, we can assume it is ensemble data.
"""
unique_cubes = set()
for cube in cubelist:
# Ignore realization of 0, as that is given to deterministic data.
if cube.coords("realization") and any(cube.coord("realization").points != 0):
return True
# Compare XML representation of cube structure check for duplicates.
cube_content = cube.xml()
if cube_content in unique_cubes:
logging.info("Ensemble data loaded.")
return True
else:
unique_cubes.add(cube_content)
logging.info("Deterministic data loaded.")
return False


def _create_callback(is_ensemble: bool) -> callable:
def _loading_callback(cube: iris.cube.Cube, field, filename: str) -> iris.cube.Cube:
"""Compose together the needed callbacks into a single function."""

def callback(cube: iris.cube.Cube, field, filename: str):
if is_ensemble:
_ensemble_callback(cube, field, filename)
else:
_deterministic_callback(cube, field, filename)

_um_normalise_callback(cube, field, filename)
_lfric_normalise_callback(cube, field, filename)
_lfric_time_coord_fix_callback(cube, field, filename)
_normalise_var0_varname(cube)
_fix_spatial_coords_callback(cube)
_fix_pressure_coord_callback(cube)
_fix_um_radtime(cube)
_fix_cell_methods(cube)
_convert_cube_units_callback(cube)
_grid_longitude_fix_callback(cube)
_fix_lfric_cloud_base_altitude(cube)
_proleptic_gregorian_fix(cube)
_lfric_time_callback(cube)
_lfric_forecast_period_standard_name_callback(cube)

return callback


def _ensemble_callback(cube, field, filename):
"""Add a realization coordinate to a cube.

Uses the filename to add an ensemble member ('realization') to each cube.
Assumes data is formatted enuk_um_0XX/enukaa_pd0HH.pp where XX is the
ensemble member.

Arguments
---------
cube: Cube
ensemble member cube
field
Raw data variable, unused.
filename: str
filename of ensemble member data
"""
if not cube.coords("realization"):
if "em" in filename:
# Assuming format is *emXX*
loc = filename.find("em") + 2
member = np.int32(filename[loc : loc + 2])
else:
# Assuming raw fields files format is enuk_um_0XX/enukaa_pd0HH
member = np.int32(filename[-15:-13])

cube.add_aux_coord(iris.coords.AuxCoord(member, standard_name="realization"))
# Most callbacks operate in-place, but save the cube when returned!
_realization_callback(cube, field, filename)
_um_normalise_callback(cube, field, filename)
_lfric_normalise_callback(cube, field, filename)
cube = _lfric_time_coord_fix_callback(cube, field, filename)
_normalise_var0_varname(cube)
_fix_spatial_coords_callback(cube)
_fix_pressure_coord_callback(cube)
_fix_um_radtime(cube)
_fix_cell_methods(cube)
cube = _convert_cube_units_callback(cube)
cube = _grid_longitude_fix_callback(cube)
_fix_lfric_cloud_base_altitude(cube)
_proleptic_gregorian_fix(cube)
_lfric_time_callback(cube)
_lfric_forecast_period_standard_name_callback(cube)
return cube


def _deterministic_callback(cube, field, filename):
def _realization_callback(cube, field, filename):
"""Give deterministic cubes a realization of 0.

This means they can be handled in the same way as ensembles through the rest
Expand All @@ -442,7 +386,7 @@ def _deterministic_callback(cube, field, filename):
# Only add if realization coordinate does not exist.
if not cube.coords("realization"):
cube.add_aux_coord(
iris.coords.AuxCoord(np.int32(0), standard_name="realization", units="1")
iris.coords.DimCoord(0, standard_name="realization", units="1")
)


Expand Down Expand Up @@ -494,7 +438,9 @@ def _lfric_normalise_callback(cube: iris.cube.Cube, field, filename):
cube.attributes["um_stash_source"] = str(sorted(ast.literal_eval(stash_list)))


def _lfric_time_coord_fix_callback(cube: iris.cube.Cube, field, filename):
def _lfric_time_coord_fix_callback(
cube: iris.cube.Cube, field, filename
) -> iris.cube.Cube:
"""Ensure the time coordinate is a DimCoord rather than an AuxCoord.

The coordinate is converted and replaced if not. SLAMed LFRic data has this
Expand All @@ -520,12 +466,10 @@ def _lfric_time_coord_fix_callback(cube: iris.cube.Cube, field, filename):
for i in range(len(time_coord.bounds))
]
iris.util.promote_aux_coord_to_dim_coord(cube, time_coord)

# Force single-valued coordinates to be scalar coordinates.
return iris.util.squeeze(cube)
return cube


def _grid_longitude_fix_callback(cube: iris.cube.Cube):
def _grid_longitude_fix_callback(cube: iris.cube.Cube) -> iris.cube.Cube:
"""Check grid_longitude coordinates are in the range -180 deg to 180 deg.

This is necessary if comparing two models with different conventions --
Expand All @@ -535,10 +479,8 @@ def _grid_longitude_fix_callback(cube: iris.cube.Cube):
model data bounds may not extend exactly to 0. or 360.
Input cubes on non-rotated grid coordinates are not impacted.
"""
import CSET.operators._utils as utils

try:
y, x = utils.get_cube_yxcoordname(cube)
y, x = get_cube_yxcoordname(cube)
except ValueError:
# Don't modify non-spatial cubes.
return cube
Expand Down Expand Up @@ -577,23 +519,21 @@ def _fix_spatial_coords_callback(cube: iris.cube.Cube):
particularly where comparing multiple input models with differing spatial
coordinates.
"""
import CSET.operators._utils as utils

# Check if cube is spatial.
if not utils.is_spatialdim(cube):
if not is_spatialdim(cube):
# Don't modify non-spatial cubes.
return cube
return

# Get spatial coords and dimension index.
y_name, x_name = utils.get_cube_yxcoordname(cube)
ny = utils.get_cube_coordindex(cube, y_name)
nx = utils.get_cube_coordindex(cube, x_name)
y_name, x_name = get_cube_yxcoordname(cube)
ny = get_cube_coordindex(cube, y_name)
nx = get_cube_coordindex(cube, x_name)

# Translate [grid_latitude, grid_longitude] to an unrotated 1-d DimCoord
# [latitude, longitude] for instances where rotated_pole=90.0
if "grid_latitude" in [coord.name() for coord in cube.coords(dim_coords=True)]:
coord_system = cube.coord("grid_latitude").coord_system
pole_lat = coord_system.grid_north_pole_latitude
pole_lat = getattr(coord_system, "grid_north_pole_latitude", None)
if pole_lat == 90.0:
lats = cube.coord("grid_latitude").points
lons = cube.coord("grid_longitude").points
Expand Down
11 changes: 1 addition & 10 deletions tests/operators/test_read.py
Original file line number Diff line number Diff line change
Expand Up @@ -91,9 +91,8 @@ def test_read_cubes_ensemble_separate_files():
)
# Check ensemble members have been merged into a single cube.
assert len(cubes) == 1
cube = cubes[0]
# Check realization is an integer.
for point in cube.coord("realization").points:
for point in cubes[0].coord("realization").points:
assert isinstance(int(point), int)


Expand Down Expand Up @@ -121,14 +120,6 @@ def test_read_cubes_incorrect_number_of_model_names():
)


def test_fieldsfile_ensemble_naming():
"""Extracting the realization from the fields file naming convention."""
cube = iris.cube.Cube([0])
filename = "myfieldsfile_enuk_um_001/enukaa_pd000"
read._ensemble_callback(cube, None, filename)
assert cube.coord("realization").points[0] == 1


def test_read_cube():
"""Returns a Cube rather than CubeList."""
from CSET.operators import constraints
Expand Down
Binary file modified tests/test_data/exeter_em01.nc
Binary file not shown.
Binary file modified tests/test_data/exeter_em02.nc
Binary file not shown.
2 changes: 1 addition & 1 deletion tests/workflow_utils/test_fetch_data.py
Original file line number Diff line number Diff line change
Expand Up @@ -242,7 +242,7 @@ def test_FilesystemFileRetriever(tmp_path):
# Check symlink points to correct file.
with open((tmp_path / "exeter_em01.nc"), "rb") as fp:
digest = hashlib.file_digest(fp, "sha256").hexdigest()
assert digest == "67899970eeca75b9378f0275ce86db3d1d613f2bc7a178540912848dc8a69ca7"
assert digest == "23761fd2456b2dbbb18f35fa4909561569af0851ebda6e3705e70e366c86ac09"


def test_FilesystemFileRetriever_no_files(tmp_path, caplog):
Expand Down