Skip to content

Hybrid Pressure load and save support for PP #921

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
bjlittle merged 12 commits into from
Jan 29, 2014
Merged

Hybrid Pressure load and save support for PP #921

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
12 commits
Select commit Hold shift + click to select a range
f0973b4
WIP
pelson Dec 12, 2013
2451fe8
CF aligned alternative implementation of hybrid pressure
esc24 Dec 17, 2013
a922e4a
Add round-trip tests
rhattersley Dec 19, 2013
254f24b
Saving of HybridPressure to pp.
pelson Dec 19, 2013
35b52d2
More units tests.
esc24 Dec 19, 2013
016de97
Find hybrid vertical coords via factory dependencies
rhattersley Dec 20, 2013
662f893
More tests, removed standard name pp saving workaround and added supp…
esc24 Dec 20, 2013
934c531
Removed new factory and swapped sigma and delta coords.
esc24 Jan 21, 2014
c34ae24
Updating test to swap delta and sigma from field vals
esc24 Jan 21, 2014
9431dcb
Review actions
esc24 Jan 27, 2014
e793b99
Fixes following rebase.
esc24 Jan 28, 2014
5a04978
biggus fallout
esc24 Jan 29, 2014
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
136 changes: 81 additions & 55 deletions lib/iris/aux_factory.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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.
#
Expand Down Expand Up @@ -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):
Expand Down Expand Up @@ -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):
Expand Down Expand Up @@ -557,54 +558,83 @@ 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:

* delta: Coord
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)):
Copy link
Member

@bjlittle bjlittle Jan 23, 2014

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@esc24 This seems unnecessarily strict ... shouldn't this just be,

delta is None and surface_air_pressure is None

Copy link
Member Author

@esc24 esc24 Jan 27, 2014

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I've left this unchanged. My thinking is that if delta and sigma are both None you cannot construct the derived coord, regardless of the surface_air_pressure, as sigma acts as a scaling which is zero if undefined (as stated in cf spec).

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 ' \
Copy link
Member

@bjlittle bjlittle Jan 23, 2014

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@esc24 Why the need to check sigma here? I'd argue that we can relax this and only check delta and surface_air_pressure here.

'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):
Expand All @@ -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):
Expand Down Expand Up @@ -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)

Expand All @@ -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)

Expand All @@ -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)
4 changes: 2 additions & 2 deletions lib/iris/cube.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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])

Expand Down
45 changes: 34 additions & 11 deletions lib/iris/etc/pp_save_rules.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Copy link
Member

@bjlittle bjlittle Jan 23, 2014

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Nice! But ... how can you guarantee that:

  • the cube contains a model_level_number coordinate
  • all the members of the aux factory are non None

Copy link
Member

@bjlittle bjlittle Jan 23, 2014

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Erk! 💥

>>> fname = '/data/local/dataZoo/PP/COLPEX/small_colpex_theta_p_alt.pp'
>>> cube = iris.load_cube(fname, 'air_potential_temperature')
>>> print cube
air_potential_temperature / (K)     (time: 6; model_level_number: 10; grid_latitude: 83; grid_longitude: 83)
     Dimension coordinates:
          time                           x                      -                  -                   -
          model_level_number             -                      x                  -                   -
          grid_latitude                  -                      -                  x                   -
          grid_longitude                 -                      -                  -                   x
     Auxiliary coordinates:
          forecast_period                x                      -                  -                   -
          level_height                   -                      x                  -                   -
          sigma                          -                      x                  -                   -
          surface_altitude               -                      -                  x                   x
     Derived coordinates:
          altitude                       -                      x                  x                   x
     Scalar coordinates:
          forecast_reference_time: 2009-09-09 22:00:00
     Attributes:
          STASH: m01s00i004
          source: Data from Met Office Unified Model 7.04
>>> cube.remove_coord('model_level_number')
>>> print cube
air_potential_temperature / (K)     (time: 6; -- : 10; grid_latitude: 83; grid_longitude: 83)
     Dimension coordinates:
          time                           x       -                  -                   -
          grid_latitude                  -       -                  x                   -
          grid_longitude                 -       -                  -                   x
     Auxiliary coordinates:
          forecast_period                x       -                  -                   -
          level_height                   -       x                  -                   -
          sigma                          -       x                  -                   -
          surface_altitude               -       -                  x                   x
     Derived coordinates:
          altitude                       -       x                  x                   x
     Scalar coordinates:
          forecast_reference_time: 2009-09-09 22:00:00
     Attributes:
          STASH: m01s00i004
          source: Data from Met Office Unified Model 7.04
>>> iris.save(cube, 'wibble.pp')
Failed to get value ('NoneType' object has no attribute 'points') to execute:     pp.lblev = scalar_coord(cm, 'model_level_number').points[0]
Failed to get value ('NoneType' object has no attribute 'points') to execute:     pp.lblev = scalar_coord(cm, 'model_level_number').points[0]
Failed to get value ('NoneType' object has no attribute 'points') to execute:     pp.lblev = scalar_coord(cm, 'model_level_number').points[0]
Failed to get value ('NoneType' object has no attribute 'points') to execute:     pp.lblev = scalar_coord(cm, 'model_level_number').points[0]
Failed to get value ('NoneType' object has no attribute 'points') to execute:     pp.lblev = scalar_coord(cm, 'model_level_number').points[0]
Failed to get value ('NoneType' object has no attribute 'points') to execute:     pp.lblev = scalar_coord(cm, 'model_level_number').points[0]
...

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)
Copy link
Member

@bjlittle bjlittle Jan 23, 2014

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

As above ...

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.
Copy link
Member

@bjlittle bjlittle Jan 23, 2014

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Oh my ... gulp

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
Expand Down
10 changes: 5 additions & 5 deletions lib/iris/fileformats/grib/load_rules.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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.
#
Expand All @@ -15,8 +15,8 @@
# You should have received a copy of the GNU Lesser General Public License
# along with Iris. If not, see <http://www.gnu.org/licenses/>.

# 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
Copy link
Member

@bjlittle bjlittle Jan 23, 2014

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@esc24 Really? Since when ? 😖

I'm all in favour for folding segments of iris-code-generators back into iris ... but I don't think this is the time to do it. Not unless you know something I don't?


import warnings

Expand Down Expand Up @@ -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))
Expand All @@ -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)
36 changes: 34 additions & 2 deletions lib/iris/fileformats/pp_rules.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -15,8 +15,8 @@
# You should have received a copy of the GNU Lesser General Public License
# along with Iris. If not, see <http://www.gnu.org/licenses/>.

# 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

Expand Down Expand Up @@ -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))
Expand Down Expand Up @@ -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)
Loading