diff --git a/lib/iris/aux_factory.py b/lib/iris/aux_factory.py index 6953fc9a1f..8a17ef1641 100644 --- a/lib/iris/aux_factory.py +++ b/lib/iris/aux_factory.py @@ -1,4 +1,4 @@ -# (C) British Crown Copyright 2010 - 2013, Met Office +# (C) British Crown Copyright 2010 - 2014, Met Office # # This file is part of Iris. # @@ -366,7 +366,8 @@ def _remap_with_bounds(self, dependency_dims, derived_dims): return nd_values_by_key def _shape(self, nd_values_by_key): - nd_values = nd_values_by_key.values() + nd_values = sorted(nd_values_by_key.values(), + key=lambda value: value.ndim) shape = list(nd_values.pop().shape) for array in nd_values: for i, size in enumerate(array.shape): @@ -446,8 +447,8 @@ def dependencies(self): return {'delta': self.delta, 'sigma': self.sigma, 'orography': self.orography} - def _derive(self, delta, sigma, surface_pressure): - temp = delta + sigma * surface_pressure + def _derive(self, delta, sigma, orography): + temp = delta + sigma * orography return temp def make_coord(self, coord_dims_func): @@ -557,12 +558,12 @@ class HybridPressureFactory(AuxCoordFactory): """ - def __init__(self, delta=None, sigma=None, surface_pressure=None): + def __init__(self, delta=None, sigma=None, surface_air_pressure=None): """ Creates a hybrid-height coordinate factory with the formula: p = ap + b * ps - At least one of `delta` or `surface_pressure` must be provided. + At least one of `delta` or `surface_air_pressure` must be provided. Args: @@ -570,41 +571,70 @@ def __init__(self, delta=None, sigma=None, surface_pressure=None): The coordinate providing the `ap` term. * sigma: Coord The coordinate providing the `b` term. - * surface_pressure: Coord + * surface_air_pressure: Coord The coordinate providing the `ps` term. """ super(HybridPressureFactory, self).__init__() + # Check that provided coords meet necessary conditions. + self._check_dependencies(delta, sigma, surface_air_pressure) + + self.delta = delta + self.sigma = sigma + self.surface_air_pressure = surface_air_pressure + + self.standard_name = 'air_pressure' + self.attributes = {} + + @property + def units(self): + if self.delta is not None: + units = self.delta.units + else: + units = self.surface_air_pressure.units + return units + + @staticmethod + def _check_dependencies(delta, sigma, + surface_air_pressure): + # Check for sufficient coordinates. + if (delta is None and (sigma is None or + surface_air_pressure is None)): + msg = 'Unable to contruct hybrid pressure coordinate factory ' \ + 'due to insufficient source coordinates.' + raise ValueError(msg) + + # Check bounds. if delta and delta.nbounds not in (0, 2): raise ValueError('Invalid delta coordinate: must have either 0 or' ' 2 bounds.') if sigma and sigma.nbounds not in (0, 2): raise ValueError('Invalid sigma coordinate: must have either 0 or' ' 2 bounds.') - if surface_pressure and surface_pressure.nbounds: - msg = 'Surface pressure coordinate {!r} has bounds.' \ - ' These will be disregarded.'.format(surface_pressure.name()) + if surface_air_pressure and surface_air_pressure.nbounds: + msg = 'Surface pressure coordinate {!r} has bounds. These will' \ + ' be disregarded.'.format(surface_air_pressure.name()) warnings.warn(msg, UserWarning, stacklevel=2) - self.delta = delta - self.sigma = sigma - self.surface_pressure = surface_pressure + # Check units. + if sigma is not None and not sigma.units.is_dimensionless(): + raise ValueError('Invalid units: sigma must be dimensionless.') + if delta is not None and surface_air_pressure is not None and \ + delta.units != surface_air_pressure.units: + msg = 'Incompatible units: delta and ' \ + 'surface_air_pressure must have the same units.' + raise ValueError(msg) + + if delta is not None: + units = delta.units + else: + units = surface_air_pressure.units - self.standard_name = 'air_pressure' - if delta is None and surface_pressure is None: - raise ValueError('Unable to determine units: no delta or' - ' surface_pressure available.') - incompatible = (delta and surface_pressure and - delta.units != surface_pressure.units) - if incompatible: - raise ValueError('Incompatible units: delta and surface_pressure' - ' must have the same units.') - self.units = (delta and delta.units) or surface_pressure.units - if not self.units.is_convertible('Pa'): - raise ValueError('Invalid units: delta and/or surface_pressure' - ' must be expressed in pressure units.') - self.attributes = {} + if not units.is_convertible('Pa'): + msg = 'Invalid units: delta and ' \ + 'surface_air_pressure must have units of pressure.' + raise ValueError(msg) @property def dependencies(self): @@ -614,10 +644,10 @@ def dependencies(self): """ return {'delta': self.delta, 'sigma': self.sigma, - 'surface_pressure': self.surface_pressure} + 'surface_air_pressure': self.surface_air_pressure} - def _derive(self, delta, sigma, surface_pressure): - temp = delta + sigma * surface_pressure + def _derive(self, delta, sigma, surface_air_pressure): + temp = delta + sigma * surface_air_pressure return temp def make_coord(self, coord_dims_func): @@ -645,7 +675,7 @@ def make_coord(self, coord_dims_func): def calc_points(): return self._derive(nd_points_by_key['delta'], nd_points_by_key['sigma'], - nd_points_by_key['surface_pressure']) + nd_points_by_key['surface_air_pressure']) shape = self._shape(nd_points_by_key) points = LazyArray(shape, calc_points) @@ -660,21 +690,22 @@ def calc_points(): def calc_bounds(): delta = nd_values_by_key['delta'] sigma = nd_values_by_key['sigma'] - surface_pressure = nd_values_by_key['surface_pressure'] + surface_air_pressure = nd_values_by_key['surface_air_pressure'] ok_bound_shapes = [(), (1,), (2,)] if delta.shape[-1:] not in ok_bound_shapes: raise ValueError('Invalid delta coordinate bounds.') if sigma.shape[-1:] not in ok_bound_shapes: raise ValueError('Invalid sigma coordinate bounds.') - if surface_pressure.shape[-1:] not in [(), (1,)]: + if surface_air_pressure.shape[-1:] not in [(), (1,)]: warnings.warn('Surface pressure coordinate has bounds. ' 'These are being disregarded.') - surface_pressure_pts = nd_points_by_key['surface_pressure'] - surface_pressure_pts_shape = list( - surface_pressure_pts.shape) - surface_pressure = surface_pressure_pts.reshape( - surface_pressure_pts_shape.append(1)) - return self._derive(delta, sigma, surface_pressure) + surface_air_pressure_pts = nd_points_by_key[ + 'surface_air_pressure'] + surface_air_pressure_pts_shape = list( + surface_air_pressure_pts.shape) + surface_air_pressure = surface_air_pressure_pts.reshape( + surface_air_pressure_pts_shape.append(1)) + return self._derive(delta, sigma, surface_air_pressure) b_shape = self._shape(nd_values_by_key) bounds = LazyArray(b_shape, calc_bounds) @@ -698,19 +729,14 @@ def update(self, old_coord, new_coord=None): any dependency using old_coord is updated to use new_coord. """ - if self.delta is old_coord: - if new_coord and new_coord.nbounds not in (0, 2): - raise ValueError('Invalid delta coordinate:' - ' must have either 0 or 2 bounds.') - self.delta = new_coord - elif self.sigma is old_coord: - if new_coord and new_coord.nbounds not in (0, 2): - raise ValueError('Invalid sigma coordinate:' - ' must have either 0 or 2 bounds.') - self.sigma = new_coord - elif self.surface_pressure is old_coord: - if new_coord and new_coord.nbounds: - msg = 'Surface pressure coordinate {!r} has bounds. ' \ - 'These will be disregarded.'.format(new_coord.name()) - warnings.warn(msg, UserWarning, stacklevel=2) - self.surface_pressure = new_coord + new_dependencies = self.dependencies + for name, coord in self.dependencies.items(): + if old_coord is coord: + new_dependencies[name] = new_coord + try: + self._check_dependencies(**new_dependencies) + except ValueError as e: + msg = 'Failed to update dependencies. ' + e.message + raise ValueError(msg) + else: + setattr(self, name, new_coord) diff --git a/lib/iris/cube.py b/lib/iris/cube.py index c3e06c96ba..2a19062a8e 100644 --- a/lib/iris/cube.py +++ b/lib/iris/cube.py @@ -750,8 +750,8 @@ def _add_unique_aux_coord(self, coord, data_dims): coord.name(), i, coord.shape[i])) elif coord.shape != (1,): - raise ValueError('You must supply the data-dimensions for ' - 'multi-valued coordinates.') + raise ValueError('Missing data dimensions for multi-valued' + ' coordinate {!r}'.format(coord.name())) self._aux_coords_and_dims.append([coord, data_dims]) diff --git a/lib/iris/etc/pp_save_rules.txt b/lib/iris/etc/pp_save_rules.txt index 83e75deaa4..bb9ff57296 100644 --- a/lib/iris/etc/pp_save_rules.txt +++ b/lib/iris/etc/pp_save_rules.txt @@ -651,23 +651,46 @@ THEN pp.lblev = scalar_coord(cm, 'air_potential_temperature').points[0] pp.blev = scalar_coord(cm, 'air_potential_temperature').points[0] -# single hybrid_height level (reversed from the equivalent load rule) +# single hybrid_height level IF + has_aux_factory(cm, iris.aux_factory.HybridHeightFactory) scalar_coord(cm, 'model_level_number') is not None scalar_coord(cm, 'model_level_number').bounds is None - scalar_coord(cm, 'level_height') is not None - scalar_coord(cm, 'level_height').bounds is not None - scalar_coord(cm, 'sigma') is not None - scalar_coord(cm, 'sigma').bounds is not None + aux_factory(cm, iris.aux_factory.HybridHeightFactory).dependencies['delta'] is not None + aux_factory(cm, iris.aux_factory.HybridHeightFactory).dependencies['delta'].bounds is not None + aux_factory(cm, iris.aux_factory.HybridHeightFactory).dependencies['sigma'] is not None + aux_factory(cm, iris.aux_factory.HybridHeightFactory).dependencies['sigma'].bounds is not None THEN pp.lbvc = 65 pp.lblev = scalar_coord(cm, 'model_level_number').points[0] - pp.blev = scalar_coord(cm, 'level_height').points[0] - pp.brlev = scalar_coord(cm, 'level_height').bounds[0, 0] - pp.brsvd[0] = scalar_coord(cm, 'level_height').bounds[0, 1] - pp.bhlev = scalar_coord(cm, 'sigma').points[0] - pp.bhrlev = scalar_coord(cm, 'sigma').bounds[0, 0] - pp.brsvd[1] = scalar_coord(cm, 'sigma').bounds[0, 1] + pp.blev = aux_factory(cm, iris.aux_factory.HybridHeightFactory).dependencies['delta'].points[0] + pp.brlev = aux_factory(cm, iris.aux_factory.HybridHeightFactory).dependencies['delta'].bounds[0, 0] + pp.brsvd[0] = aux_factory(cm, iris.aux_factory.HybridHeightFactory).dependencies['delta'].bounds[0, 1] + pp.bhlev = aux_factory(cm, iris.aux_factory.HybridHeightFactory).dependencies['sigma'].points[0] + pp.bhrlev = aux_factory(cm, iris.aux_factory.HybridHeightFactory).dependencies['sigma'].bounds[0, 0] + pp.brsvd[1] = aux_factory(cm, iris.aux_factory.HybridHeightFactory).dependencies['sigma'].bounds[0, 1] + +# single hybrid pressure level +IF + has_aux_factory(cm, iris.aux_factory.HybridPressureFactory) + scalar_coord(cm, 'model_level_number') is not None + scalar_coord(cm, 'model_level_number').bounds is None + aux_factory(cm, iris.aux_factory.HybridPressureFactory).dependencies['delta'] is not None + aux_factory(cm, iris.aux_factory.HybridPressureFactory).dependencies['delta'].bounds is not None + aux_factory(cm, iris.aux_factory.HybridPressureFactory).dependencies['sigma'] is not None + aux_factory(cm, iris.aux_factory.HybridPressureFactory).dependencies['sigma'].bounds is not None +THEN + pp.lbvc = 9 + pp.lblev = scalar_coord(cm, 'model_level_number').points[0] + + # Note that sigma and delta are swapped around from the hybrid height rules above. + pp.blev = aux_factory(cm, iris.aux_factory.HybridPressureFactory).dependencies['sigma'].points[0] + pp.brlev = aux_factory(cm, iris.aux_factory.HybridPressureFactory).dependencies['sigma'].bounds[0, 0] + pp.brsvd[0] = aux_factory(cm, iris.aux_factory.HybridPressureFactory).dependencies['sigma'].bounds[0, 1] + + pp.bhlev = aux_factory(cm, iris.aux_factory.HybridPressureFactory).dependencies['delta'].points[0] + pp.bhrlev = aux_factory(cm, iris.aux_factory.HybridPressureFactory).dependencies['delta'].bounds[0, 0] + pp.brsvd[1] = aux_factory(cm, iris.aux_factory.HybridPressureFactory).dependencies['delta'].bounds[0, 1] #MDI diff --git a/lib/iris/fileformats/grib/load_rules.py b/lib/iris/fileformats/grib/load_rules.py index ffd8aeefd2..8fe3f6031c 100644 --- a/lib/iris/fileformats/grib/load_rules.py +++ b/lib/iris/fileformats/grib/load_rules.py @@ -1,4 +1,4 @@ -# (C) British Crown Copyright 2013, Met Office +# (C) British Crown Copyright 2013 - 2014, Met Office # # This file is part of Iris. # @@ -15,8 +15,8 @@ # You should have received a copy of the GNU Lesser General Public License # along with Iris. If not, see . -# DO NOT EDIT DIRECTLY -# Auto-generated from SciTools/iris-code-generators:tools/gen_rules.py +# Historically this was auto-generated from +# SciTools/iris-code-generators:tools/gen_rules.py import warnings @@ -352,7 +352,7 @@ def convert(grib): aux_coords_and_dims.append((AuxCoord(grib.scaledValueOfFirstFixedSurface, standard_name='model_level_number', attributes={'positive': 'up'}), None)) aux_coords_and_dims.append((DimCoord(grib.pv[grib.scaledValueOfFirstFixedSurface], long_name='level_pressure', units='Pa'), None)) aux_coords_and_dims.append((AuxCoord(grib.pv[grib.numberOfCoordinatesValues/2 + grib.scaledValueOfFirstFixedSurface], long_name='sigma'), None)) - factories.append(Factory(HybridPressureFactory, [{'long_name': 'level_pressure'}, {'long_name': 'sigma'}, Reference('surface_pressure')])) + factories.append(Factory(HybridPressureFactory, [{'long_name': 'level_pressure'}, {'long_name': 'sigma'}, Reference('surface_air_pressure')])) if grib._originatingCentre != 'unknown': aux_coords_and_dims.append((AuxCoord(points=grib._originatingCentre, long_name='originating_centre', units='no_unit'), None)) @@ -372,7 +372,7 @@ def convert(grib): (grib.parameterCategory == 3) and \ (grib.parameterNumber == 25) and \ (grib.typeOfFirstFixedSurface == 105): - references.append(ReferenceTarget('surface_pressure', lambda cube: {'standard_name': 'surface_air_pressure', 'units': 'Pa', 'data': np.exp(cube.data)})) + references.append(ReferenceTarget('surface_air_pressure', lambda cube: {'standard_name': 'surface_air_pressure', 'units': 'Pa', 'data': np.exp(cube.data)})) return (factories, references, standard_name, long_name, units, attributes, cell_methods, dim_coords_and_dims, aux_coords_and_dims) diff --git a/lib/iris/fileformats/pp_rules.py b/lib/iris/fileformats/pp_rules.py index 9ec1383a5d..ab29c18956 100644 --- a/lib/iris/fileformats/pp_rules.py +++ b/lib/iris/fileformats/pp_rules.py @@ -15,8 +15,8 @@ # You should have received a copy of the GNU Lesser General Public License # along with Iris. If not, see . -# DO NOT EDIT DIRECTLY -# Auto-generated from SciTools/iris-code-generators:tools/gen_rules.py +# Historically this was auto-generated from +# SciTools/iris-code-generators:tools/gen_rules.py import warnings @@ -307,6 +307,35 @@ def convert(f): (f.lbvc == 19): aux_coords_and_dims.append((DimCoord(f.blev, standard_name='air_potential_temperature', units='K', attributes={'positive': 'up'}), None)) + # Hybrid pressure coordinate + if f.lbvc == 9: + model_level_number = DimCoord(f.lblev, + standard_name='model_level_number', + attributes={'positive': 'up'}) + # The following match the hybrid height scheme, but data has the + # blev and bhlev values the other way around. + #level_pressure = DimCoord(f.blev, + # long_name='level_pressure', + # units='Pa', + # bounds=[f.brlev, f.brsvd[0]]) + #sigma = AuxCoord(f.bhlev, + # long_name='sigma', + # bounds=[f.bhrlev, f.brsvd[1]]) + level_pressure = DimCoord(f.bhlev, + long_name='level_pressure', + units='Pa', + bounds=[f.bhrlev, f.brsvd[1]]) + sigma = AuxCoord(f.blev, + long_name='sigma', + bounds=[f.brlev, f.brsvd[0]]) + aux_coords_and_dims.extend([(model_level_number, None), + (level_pressure, None), + (sigma, None)]) + factories.append(Factory(HybridPressureFactory, + [{'long_name': 'level_pressure'}, + {'long_name': 'sigma'}, + Reference('surface_air_pressure')])) + if f.lbvc == 65: aux_coords_and_dims.append((DimCoord(f.lblev, standard_name='model_level_number', attributes={'positive': 'up'}), None)) aux_coords_and_dims.append((DimCoord(f.blev, long_name='level_height', units='m', bounds=[f.brlev, f.brsvd[0]], attributes={'positive': 'up'}), None)) @@ -381,5 +410,8 @@ def convert(f): if f.lbuser[3] == 33: references.append(ReferenceTarget('orography', None)) + if f.lbuser[3] == 409 or f.lbuser[3] == 1: + references.append(ReferenceTarget('surface_air_pressure', None)) + return (factories, references, standard_name, long_name, units, attributes, cell_methods, dim_coords_and_dims, aux_coords_and_dims) diff --git a/lib/iris/fileformats/rules.py b/lib/iris/fileformats/rules.py index dd648e17ca..979bf89dbb 100644 --- a/lib/iris/fileformats/rules.py +++ b/lib/iris/fileformats/rules.py @@ -594,6 +594,35 @@ def scalar_cell_method(cube, method, coord_name): return found_cell_method +def has_aux_factory(cube, aux_factory_class): + """ + Try to find an class:`~iris.aux_factory.AuxCoordFactory` instance of the + specified type on the cube. + + """ + for factory in cube.aux_factories: + if isinstance(factory, aux_factory_class): + return True + return False + + +def aux_factory(cube, aux_factory_class): + """ + Return the class:`~iris.aux_factory.AuxCoordFactory` instance of the + specified type from a cube. + + """ + aux_factories = [aux_factory for aux_factory in cube.aux_factories if + isinstance(aux_factory, aux_factory_class)] + if not aux_factories: + raise ValueError('Cube does not have an aux factory of ' + 'type {!r}.'.format(aux_factory_class)) + elif len(aux_factories) > 1: + raise ValueError('Cube has more than one aux factory of ' + 'type {!r}.'.format(aux_factory_class)) + return aux_factories[0] + + class _ReferenceError(Exception): """Signals an invalid/missing reference field.""" pass diff --git a/lib/iris/tests/integration/test_pp.py b/lib/iris/tests/integration/test_pp.py index 4affefd923..bfd3f73674 100644 --- a/lib/iris/tests/integration/test_pp.py +++ b/lib/iris/tests/integration/test_pp.py @@ -23,10 +23,22 @@ import mock import numpy as np +from iris.aux_factory import HybridHeightFactory, HybridPressureFactory +from iris.coords import AuxCoord +from iris.cube import Cube import iris.fileformats.pp +import iris.fileformats.pp_rules class TestVertical(tests.IrisTest): + def _test_coord(self, cube, point, bounds=None, **kwargs): + coords = cube.coords(**kwargs) + self.assertEqual(len(coords), 1, 'failed to find exactly one coord' + ' using: {}'.format(kwargs)) + self.assertEqual(coords[0].points, point) + if bounds is not None: + self.assertArrayEqual(coords[0].bounds, [bounds]) + def test_soil_level_round_trip(self): # Use pp.load_cubes() to convert a fake PPField into a Cube. # NB. Use MagicMock so that SplittableInt header items, such as @@ -40,9 +52,7 @@ def test_soil_level_round_trip(self): cube = next(iris.fileformats.pp.load_cubes('DUMMY')) self.assertIn('soil', cube.standard_name) - self.assertEqual(len(cube.coords('soil_model_level_number')), 1) - self.assertEqual(cube.coord('soil_model_level_number').points, - soil_level) + self._test_coord(cube, soil_level, long_name='soil_model_level_number') # Now use the save rules to convert the Cube back into a PPField. field = iris.fileformats.pp.PPField3() @@ -56,7 +66,8 @@ def test_soil_level_round_trip(self): self.assertEqual(field.lblev, soil_level) def test_potential_temperature_level_round_trip(self): - """Check save+load for data on 'potential temperature' levels.""" + # Check save+load for data on 'potential temperature' levels. + # Use pp.load_cubes() to convert a fake PPField into a Cube. # NB. Use MagicMock so that SplittableInt header items, such as # LBCODE, support len(). @@ -67,9 +78,8 @@ def test_potential_temperature_level_round_trip(self): with mock.patch('iris.fileformats.pp.load', new=load) as load: cube = next(iris.fileformats.pp.load_cubes('DUMMY')) - self.assertEqual(len(cube.coords('air_potential_temperature')), 1) - self.assertEqual(cube.coord('air_potential_temperature').points, - potm_value) + self._test_coord(cube, potm_value, + standard_name='air_potential_temperature') # Now use the save rules to convert the Cube back into a PPField. field = iris.fileformats.pp.PPField3() @@ -82,6 +92,161 @@ def test_potential_temperature_level_round_trip(self): self.assertEqual(field.lbvc, 19) self.assertEqual(field.blev, potm_value) + def test_hybrid_pressure_round_trip(self): + # Use pp.load_cubes() to convert fake PPFields into Cubes. + # NB. Use MagicMock so that SplittableInt header items, such as + # LBCODE, support len(). + def field_with_data(scale=1): + x, y = 40, 30 + field = mock.MagicMock(_data=np.arange(1200).reshape(y, x) * scale, + lbcode=[1], lbnpt=x, lbrow=y, + bzx=350, bdx=1.5, bzy=40, bdy=1.5, + lbuser=[0] * 7, lbrsvd=[0] * 4) + field._x_coord_name = lambda: 'longitude' + field._y_coord_name = lambda: 'latitude' + field.coord_system = lambda: None + return field + + # Make a fake reference surface field. + pressure_field = field_with_data(10) + pressure_field.stash = iris.fileformats.pp.STASH(1, 0, 409) + pressure_field.lbuser[3] = 409 + + # Make a fake data field which needs the reference surface. + model_level = 5678 + sigma_lower, sigma, sigma_upper = 0.85, 0.9, 0.95 + delta_lower, delta, delta_upper = 0.05, 0.1, 0.15 + data_field = field_with_data() + data_field.configure_mock(lbvc=9, lblev=model_level, + bhlev=delta, bhrlev=delta_lower, + blev=sigma, brlev=sigma_lower, + brsvd=[sigma_upper, delta_upper]) + + # Convert both fields to cubes. + load = mock.Mock(return_value=iter([pressure_field, data_field])) + with mock.patch('iris.fileformats.pp.load', new=load) as load: + pressure_cube, data_cube = iris.fileformats.pp.load_cubes('DUMMY') + + # Check the reference surface cube looks OK. + self.assertEqual(pressure_cube.standard_name, 'surface_air_pressure') + self.assertEqual(pressure_cube.units, 'Pa') + + # Check the data cube is set up to use hybrid-pressure. + self._test_coord(data_cube, model_level, + standard_name='model_level_number') + self._test_coord(data_cube, delta, [delta_lower, delta_upper], + long_name='level_pressure') + self._test_coord(data_cube, sigma, [sigma_lower, sigma_upper], + long_name='sigma') + aux_factories = data_cube.aux_factories + self.assertEqual(len(aux_factories), 1) + surface_coord = aux_factories[0].dependencies['surface_air_pressure'] + self.assertArrayEqual(surface_coord.points, + np.arange(12000, step=10).reshape(30, 40)) + + # Now use the save rules to convert the Cubes back into PPFields. + pressure_field = iris.fileformats.pp.PPField3() + pressure_field.lbfc = 0 + pressure_field.lbvc = 0 + pressure_field.brsvd = [None, None] + pressure_field.lbuser = [None] * 7 + iris.fileformats.pp._ensure_save_rules_loaded() + iris.fileformats.pp._save_rules.verify(pressure_cube, pressure_field) + + data_field = iris.fileformats.pp.PPField3() + data_field.lbfc = 0 + data_field.lbvc = 0 + data_field.brsvd = [None, None] + data_field.lbuser = [None] * 7 + iris.fileformats.pp._save_rules.verify(data_cube, data_field) + + # The reference surface field should have STASH=409 + self.assertArrayEqual(pressure_field.lbuser, + [None, None, None, 409, None, None, 1]) + + # Check the data field has the vertical coordinate as originally + # specified. + self.assertEqual(data_field.lbvc, 9) + self.assertEqual(data_field.lblev, model_level) + self.assertEqual(data_field.bhlev, delta) + self.assertEqual(data_field.bhrlev, delta_lower) + self.assertEqual(data_field.blev, sigma) + self.assertEqual(data_field.brlev, sigma_lower) + self.assertEqual(data_field.brsvd, [sigma_upper, delta_upper]) + + def test_hybrid_height_with_non_standard_coords(self): + # Check the save rules are using the AuxFactory to find the + # hybrid height coordinates and not relying on their names. + ny, nx = 30, 40 + sigma_lower, sigma, sigma_upper = 0.75, 0.8, 0.75 + delta_lower, delta, delta_upper = 150, 200, 250 + + cube = Cube(np.zeros((ny, nx)), 'air_temperature') + level_coord = AuxCoord(0, 'model_level_number') + cube.add_aux_coord(level_coord) + delta_coord = AuxCoord(delta, bounds=[[delta_lower, delta_upper]], + long_name='moog', units='m') + sigma_coord = AuxCoord(sigma, bounds=[[sigma_lower, sigma_upper]], + long_name='mavis') + surface_altitude_coord = AuxCoord(np.zeros((ny, nx)), + 'surface_altitude', units='m') + cube.add_aux_coord(delta_coord) + cube.add_aux_coord(sigma_coord) + cube.add_aux_coord(surface_altitude_coord, (0, 1)) + cube.add_aux_factory(HybridHeightFactory(delta_coord, sigma_coord, + surface_altitude_coord)) + + field = iris.fileformats.pp.PPField3() + field.lbfc = 0 + field.lbvc = 0 + field.brsvd = [None, None] + field.lbuser = [None] * 7 + iris.fileformats.pp._ensure_save_rules_loaded() + iris.fileformats.pp._save_rules.verify(cube, field) + + self.assertEqual(field.blev, delta) + self.assertEqual(field.brlev, delta_lower) + self.assertEqual(field.bhlev, sigma) + self.assertEqual(field.bhrlev, sigma_lower) + self.assertEqual(field.brsvd, [delta_upper, sigma_upper]) + + def test_hybrid_pressure_with_non_standard_coords(self): + # Check the save rules are using the AuxFactory to find the + # hybrid pressure coordinates and not relying on their names. + ny, nx = 30, 40 + sigma_lower, sigma, sigma_upper = 0.75, 0.8, 0.75 + delta_lower, delta, delta_upper = 0.15, 0.2, 0.25 + + cube = Cube(np.zeros((ny, nx)), 'air_temperature') + level_coord = AuxCoord(0, 'model_level_number') + cube.add_aux_coord(level_coord) + delta_coord = AuxCoord(delta, bounds=[[delta_lower, delta_upper]], + long_name='moog', units='Pa') + sigma_coord = AuxCoord(sigma, bounds=[[sigma_lower, sigma_upper]], + long_name='mavis') + surface_air_pressure_coord = AuxCoord(np.zeros((ny, nx)), + 'surface_air_pressure', + units='Pa') + cube.add_aux_coord(delta_coord) + cube.add_aux_coord(sigma_coord) + cube.add_aux_coord(surface_air_pressure_coord, (0, 1)) + cube.add_aux_factory(HybridPressureFactory( + delta_coord, sigma_coord, surface_air_pressure_coord)) + + field = iris.fileformats.pp.PPField3() + field.lbfc = 0 + field.lbvc = 0 + field.brsvd = [None, None] + field.lbuser = [None] * 7 + iris.fileformats.pp._ensure_save_rules_loaded() + iris.fileformats.pp._save_rules.verify(cube, field) + + self.assertEqual(field.bhlev, delta) + self.assertEqual(field.bhrlev, delta_lower) + self.assertEqual(field.blev, sigma) + self.assertEqual(field.brlev, sigma_lower) + self.assertEqual(field.brsvd, [sigma_upper, delta_upper]) + class TestCoordinateForms(tests.IrisTest): def _common(self, x_coord): diff --git a/lib/iris/tests/test_hybrid.py b/lib/iris/tests/test_hybrid.py index e12b9bff6a..c9ff75d33e 100644 --- a/lib/iris/tests/test_hybrid.py +++ b/lib/iris/tests/test_hybrid.py @@ -1,4 +1,4 @@ -# (C) British Crown Copyright 2010 - 2013, Met Office +# (C) British Crown Copyright 2010 - 2014, Met Office # # This file is part of Iris. # @@ -207,7 +207,8 @@ def test_invalid_dependencies(self): # Cause all warnings to raise Exceptions warnings.simplefilter("error") with self.assertRaises(UserWarning): - factory = HybridPressureFactory(surface_pressure=sigma) + factory = HybridPressureFactory( + sigma=sigma, surface_air_pressure=sigma) def test_bounded_surface_pressure(self): # Start with everything normal diff --git a/lib/iris/tests/unit/aux_factory/__init__.py b/lib/iris/tests/unit/aux_factory/__init__.py new file mode 100644 index 0000000000..8b6504eeea --- /dev/null +++ b/lib/iris/tests/unit/aux_factory/__init__.py @@ -0,0 +1,17 @@ +# (C) British Crown Copyright 2014, Met Office +# +# This file is part of Iris. +# +# Iris is free software: you can redistribute it and/or modify it under +# the terms of the GNU Lesser General Public License as published by the +# Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# Iris is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU Lesser General Public License for more details. +# +# You should have received a copy of the GNU Lesser General Public License +# along with Iris. If not, see . +"""Unit tests for the :mod:`iris.aux_factory` module.""" diff --git a/lib/iris/tests/unit/aux_factory/test_HybridPressureFactory.py b/lib/iris/tests/unit/aux_factory/test_HybridPressureFactory.py new file mode 100644 index 0000000000..0d86cbb8b4 --- /dev/null +++ b/lib/iris/tests/unit/aux_factory/test_HybridPressureFactory.py @@ -0,0 +1,276 @@ +# (C) British Crown Copyright 2014, Met Office +# +# This file is part of Iris. +# +# Iris is free software: you can redistribute it and/or modify it under +# the terms of the GNU Lesser General Public License as published by the +# Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# Iris is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU Lesser General Public License for more details. +# +# You should have received a copy of the GNU Lesser General Public License +# along with Iris. If not, see . +""" +Unit tests for the +`iris.aux_factory.HybridPressureFactory` class. + +""" + +# Import iris.tests first so that some things can be initialised before +# importing anything else. +import iris.tests as tests + +import mock + +import numpy as np + +import iris +from iris.aux_factory import HybridPressureFactory + + +class Test___init__(tests.IrisTest): + def setUp(self): + self.delta = mock.Mock(units=iris.unit.Unit('Pa'), nbounds=0) + self.sigma = mock.Mock(units=iris.unit.Unit('1'), nbounds=0) + self.surface_air_pressure = mock.Mock(units=iris.unit.Unit('Pa'), + nbounds=0) + + def test_insufficient_coords(self): + with self.assertRaises(ValueError): + HybridPressureFactory() + with self.assertRaises(ValueError): + HybridPressureFactory( + delta=None, sigma=self.sigma, + surface_air_pressure=None) + with self.assertRaises(ValueError): + HybridPressureFactory( + delta=None, sigma=None, + surface_air_pressure=self.surface_air_pressure) + + def test_incompatible_delta_units(self): + self.delta.units = iris.unit.Unit('m') + with self.assertRaises(ValueError): + HybridPressureFactory( + delta=self.delta, sigma=self.sigma, + surface_air_pressure=self.surface_air_pressure) + + def test_incompatible_sigma_units(self): + self.sigma.units = iris.unit.Unit('Pa') + with self.assertRaises(ValueError): + HybridPressureFactory( + delta=self.delta, sigma=self.sigma, + surface_air_pressure=self.surface_air_pressure) + + def test_incompatible_surface_air_pressure_units(self): + self.surface_air_pressure.units = iris.unit.Unit('unknown') + with self.assertRaises(ValueError): + HybridPressureFactory( + delta=self.delta, sigma=self.sigma, + surface_air_pressure=self.surface_air_pressure) + + def test_different_pressure_units(self): + self.delta.units = iris.unit.Unit('hPa') + self.surface_air_pressure.units = iris.unit.Unit('Pa') + with self.assertRaises(ValueError): + HybridPressureFactory( + delta=self.delta, sigma=self.sigma, + surface_air_pressure=self.surface_air_pressure) + + def test_too_many_delta_bounds(self): + self.delta.nbounds = 4 + with self.assertRaises(ValueError): + HybridPressureFactory( + delta=self.delta, sigma=self.sigma, + surface_air_pressure=self.surface_air_pressure) + + def test_too_many_sigma_bounds(self): + self.sigma.nbounds = 4 + with self.assertRaises(ValueError): + HybridPressureFactory( + delta=self.delta, sigma=self.sigma, + surface_air_pressure=self.surface_air_pressure) + + def test_factory_metadata(self): + factory = HybridPressureFactory( + delta=self.delta, sigma=self.sigma, + surface_air_pressure=self.surface_air_pressure) + self.assertEqual(factory.standard_name, 'air_pressure') + self.assertIsNone(factory.long_name) + self.assertIsNone(factory.var_name) + self.assertEqual(factory.units, self.delta.units) + self.assertEqual(factory.units, self.surface_air_pressure.units) + self.assertIsNone(factory.coord_system) + self.assertEqual(factory.attributes, {}) + + +class Test_dependencies(tests.IrisTest): + def setUp(self): + self.delta = mock.Mock(units=iris.unit.Unit('Pa'), nbounds=0) + self.sigma = mock.Mock(units=iris.unit.Unit('1'), nbounds=0) + self.surface_air_pressure = mock.Mock(units=iris.unit.Unit('Pa'), + nbounds=0) + + def test_value(self): + kwargs = dict(delta=self.delta, + sigma=self.sigma, + surface_air_pressure=self.surface_air_pressure) + factory = HybridPressureFactory(**kwargs) + self.assertEqual(factory.dependencies, kwargs) + + +class Test_make_coord(tests.IrisTest): + @staticmethod + def coords_dims_func(coord): + mapping = dict(level_pressure=(0,), sigma=(0,), + surface_air_pressure=(1, 2)) + return mapping[coord.name()] + + def setUp(self): + self.delta = iris.coords.DimCoord( + [0.0, 1.0, 2.0], long_name='level_pressure', units='Pa') + self.sigma = iris.coords.DimCoord( + [1.0, 0.9, 0.8], long_name='sigma') + self.surface_air_pressure = iris.coords.AuxCoord( + np.arange(4).reshape(2, 2), 'surface_air_pressure', + units='Pa') + + def test_points_only(self): + # Determine expected coord by manually broadcasting coord points + # knowing the dimension mapping. + delta_pts = self.delta.points[..., np.newaxis, np.newaxis] + sigma_pts = self.sigma.points[..., np.newaxis, np.newaxis] + surf_pts = self.surface_air_pressure.points[np.newaxis, ...] + expected_points = delta_pts + sigma_pts * surf_pts + expected_coord = iris.coords.AuxCoord(expected_points, + standard_name='air_pressure', + units='Pa') + factory = HybridPressureFactory( + delta=self.delta, sigma=self.sigma, + surface_air_pressure=self.surface_air_pressure) + derived_coord = factory.make_coord(self.coords_dims_func) + self.assertEqual(expected_coord, derived_coord) + + def test_none_delta(self): + delta_pts = 0 + sigma_pts = self.sigma.points[..., np.newaxis, np.newaxis] + surf_pts = self.surface_air_pressure.points[np.newaxis, ...] + expected_points = delta_pts + sigma_pts * surf_pts + expected_coord = iris.coords.AuxCoord(expected_points, + standard_name='air_pressure', + units='Pa') + factory = HybridPressureFactory( + sigma=self.sigma, surface_air_pressure=self.surface_air_pressure) + derived_coord = factory.make_coord(self.coords_dims_func) + self.assertEqual(expected_coord, derived_coord) + + def test_none_sigma(self): + delta_pts = self.delta.points[..., np.newaxis, np.newaxis] + sigma_pts = 0 + surf_pts = self.surface_air_pressure.points[np.newaxis, ...] + expected_points = delta_pts + sigma_pts * surf_pts + expected_coord = iris.coords.AuxCoord(expected_points, + standard_name='air_pressure', + units='Pa') + factory = HybridPressureFactory( + delta=self.delta, surface_air_pressure=self.surface_air_pressure) + derived_coord = factory.make_coord(self.coords_dims_func) + self.assertEqual(expected_coord, derived_coord) + + def test_none_surface_air_pressure(self): + # Note absence of broadcasting as multidimensional coord + # is not present. + expected_points = self.delta.points + expected_coord = iris.coords.AuxCoord(expected_points, + standard_name='air_pressure', + units='Pa') + factory = HybridPressureFactory(delta=self.delta, sigma=self.sigma) + derived_coord = factory.make_coord(self.coords_dims_func) + self.assertEqual(expected_coord, derived_coord) + + def test_with_bounds(self): + self.delta.guess_bounds(0) + self.sigma.guess_bounds(0.5) + # Determine expected coord by manually broadcasting coord points + # and bounds based on the dimension mapping. + delta_pts = self.delta.points[..., np.newaxis, np.newaxis] + sigma_pts = self.sigma.points[..., np.newaxis, np.newaxis] + surf_pts = self.surface_air_pressure.points[np.newaxis, ...] + expected_points = delta_pts + sigma_pts * surf_pts + delta_vals = self.delta.bounds.reshape(3, 1, 1, 2) + sigma_vals = self.sigma.bounds.reshape(3, 1, 1, 2) + surf_vals = self.surface_air_pressure.points.reshape(1, 2, 2, 1) + expected_bounds = delta_vals + sigma_vals * surf_vals + expected_coord = iris.coords.AuxCoord(expected_points, + standard_name='air_pressure', + units='Pa', + bounds=expected_bounds) + factory = HybridPressureFactory( + delta=self.delta, sigma=self.sigma, + surface_air_pressure=self.surface_air_pressure) + derived_coord = factory.make_coord(self.coords_dims_func) + self.assertEqual(expected_coord, derived_coord) + + +class Test_update(tests.IrisTest): + def setUp(self): + self.delta = mock.Mock(units=iris.unit.Unit('Pa'), nbounds=0) + self.sigma = mock.Mock(units=iris.unit.Unit('1'), nbounds=0) + self.surface_air_pressure = mock.Mock(units=iris.unit.Unit('Pa'), + nbounds=0) + + self.factory = HybridPressureFactory( + delta=self.delta, sigma=self.sigma, + surface_air_pressure=self.surface_air_pressure) + + def test_good_delta(self): + new_delta_coord = mock.Mock(units=iris.unit.Unit('Pa'), nbounds=0) + self.factory.update(self.delta, new_delta_coord) + self.assertIs(self.factory.delta, new_delta_coord) + + def test_bad_delta(self): + new_delta_coord = mock.Mock(units=iris.unit.Unit('1'), nbounds=0) + with self.assertRaises(ValueError): + self.factory.update(self.delta, new_delta_coord) + + def test_alternative_bad_delta(self): + new_delta_coord = mock.Mock(units=iris.unit.Unit('Pa'), nbounds=4) + with self.assertRaises(ValueError): + self.factory.update(self.delta, new_delta_coord) + + def test_good_surface_air_pressure(self): + new_surface_p_coord = mock.Mock(units=iris.unit.Unit('Pa'), nbounds=0) + self.factory.update(self.surface_air_pressure, new_surface_p_coord) + self.assertIs(self.factory.surface_air_pressure, new_surface_p_coord) + + def test_bad_surface_air_pressure(self): + new_surface_p_coord = mock.Mock(units=iris.unit.Unit('km'), nbounds=0) + with self.assertRaises(ValueError): + self.factory.update(self.surface_air_pressure, new_surface_p_coord) + + def test_non_dependency(self): + old_coord = mock.Mock() + new_coord = mock.Mock() + orig_dependencies = self.factory.dependencies + self.factory.update(old_coord, new_coord) + self.assertEqual(orig_dependencies, self.factory.dependencies) + + def test_none_delta(self): + self.factory.update(self.delta, None) + self.assertIsNone(self.factory.delta) + + def test_none_sigma(self): + self.factory.update(self.sigma, None) + self.assertIsNone(self.factory.sigma) + + def test_insufficient_coords(self): + self.factory.update(self.delta, None) + with self.assertRaises(ValueError): + self.factory.update(self.surface_air_pressure, None) + + +if __name__ == "__main__": + tests.main() diff --git a/lib/iris/tests/unit/fileformats/pp_rules/test_convert.py b/lib/iris/tests/unit/fileformats/pp_rules/test_convert.py index 0ee486b854..9ea3f789b2 100644 --- a/lib/iris/tests/unit/fileformats/pp_rules/test_convert.py +++ b/lib/iris/tests/unit/fileformats/pp_rules/test_convert.py @@ -23,6 +23,8 @@ import mock import types +import numpy as np + from iris.fileformats.pp_rules import convert from iris.util import guess_coord_axis from iris.fileformats.pp import SplittableInt @@ -31,41 +33,102 @@ class TestLBVC(tests.IrisTest): - def test_soil_levels(self): - field = mock.MagicMock(lbvc=6, lblev=1234) + @staticmethod + def _is_potm_level_coord(coord): + return (coord.standard_name == 'air_potential_temperature' and + coord.attributes['positive'] == 'up') + + @staticmethod + def _is_model_level_number_coord(coord): + return (coord.standard_name == 'model_level_number' and + coord.units.is_dimensionless() and + coord.attributes['positive'] == 'up') + + @staticmethod + def _is_level_pressure_coord(coord): + return (coord.name() == 'level_pressure' and + coord.units == 'Pa') + + @staticmethod + def _is_sigma_coord(coord): + return (coord.name() == 'sigma' and + coord.units.is_dimensionless()) + + @staticmethod + def _is_soil_model_level_number_coord(coord): + return (coord.long_name == 'soil_model_level_number' and + coord.units.is_dimensionless() and + coord.attributes['positive'] == 'down') + + def _test_for_coord(self, field, coord_predicate, expected_points, + expected_bounds): (factories, references, standard_name, long_name, units, attributes, cell_methods, dim_coords_and_dims, aux_coords_and_dims) = convert(field) - def is_soil_model_level_coord(coord_and_dims): - coord, dims = coord_and_dims - return coord.long_name == 'soil_model_level_number' + # Check for one and only one matching coordinate. + matching_coords = [coord for coord, _ in aux_coords_and_dims if + coord_predicate(coord)] + self.assertEqual(len(matching_coords), 1) + coord = matching_coords[0] - coords_and_dims = filter(is_soil_model_level_coord, - aux_coords_and_dims) - self.assertEqual(len(coords_and_dims), 1) - coord, dims = coords_and_dims[0] - self.assertEqual(coord.points, 1234) - self.assertIsNone(coord.bounds) - self.assertEqual(guess_coord_axis(coord), 'Z') + # Check points and bounds. + if expected_points is not None: + self.assertArrayEqual(coord.points, expected_points) + + if expected_bounds is None: + self.assertIsNone(coord.bounds) + else: + self.assertArrayEqual(coord.bounds, expected_bounds) + + def test_soil_levels(self): + level = 1234 + field = mock.MagicMock(lbvc=6, lblev=level) + self._test_for_coord(field, TestLBVC._is_soil_model_level_number_coord, + expected_points=np.array([level]), + expected_bounds=None) + + def test_hybrid_pressure_model_level_number(self): + level = 5678 + field = mock.MagicMock(lbvc=9, lblev=level, + blev=20, brlev=23, bhlev=42, + bhrlev=45, brsvd=[17, 40]) + self._test_for_coord(field, TestLBVC._is_model_level_number_coord, + expected_points=np.array([level]), + expected_bounds=None) + + def test_hybrid_pressure_delta(self): + delta_point = 12.0 + delta_lower_bound = 11.0 + delta_upper_bound = 13.0 + field = mock.MagicMock(lbvc=9, lblev=5678, + blev=20, brlev=23, bhlev=delta_point, + bhrlev=delta_lower_bound, + brsvd=[17, delta_upper_bound]) + self._test_for_coord(field, TestLBVC._is_level_pressure_coord, + expected_points=np.array([delta_point]), + expected_bounds=np.array([[delta_lower_bound, + delta_upper_bound]])) + + def test_hybrid_pressure_sigma(self): + sigma_point = 0.5 + sigma_lower_bound = 0.6 + sigma_upper_bound = 0.4 + field = mock.MagicMock(lbvc=9, lblev=5678, + blev=sigma_point, brlev=sigma_lower_bound, + bhlev=12, bhrlev=11, + brsvd=[sigma_upper_bound, 13]) + self._test_for_coord(field, TestLBVC._is_sigma_coord, + expected_points=np.array([sigma_point]), + expected_bounds=np.array([[sigma_lower_bound, + sigma_upper_bound]])) def test_potential_temperature_levels(self): potm_value = 27.32 field = mock.MagicMock(lbvc=19, blev=potm_value) - (factories, references, standard_name, long_name, units, - attributes, cell_methods, dim_coords_and_dims, - aux_coords_and_dims) = convert(field) - - def is_potm_level_coord(coord_and_dims): - coord, dims = coord_and_dims - return coord.standard_name == 'air_potential_temperature' - - coords_and_dims = filter(is_potm_level_coord, aux_coords_and_dims) - self.assertEqual(len(coords_and_dims), 1) - coord, dims = coords_and_dims[0] - self.assertArrayEqual(coord.points, [potm_value]) - self.assertIsNone(coord.bounds) - self.assertEqual(guess_coord_axis(coord), 'Z') + self._test_for_coord(field, TestLBVC._is_potm_level_coord, + expected_points=np.array([potm_value]), + expected_bounds=None) class TestLBTIM(tests.IrisTest):