Skip to content

Conversation

@mmcpartland
Copy link

Adding MIT Lincoln Laboratory code to produce Turbulence Index based on Ellrod 1997. This class is intended for estimates above 1500 meters. The principal class, TurbulenceIndexAbove1500m_USAF, is based on Ellrod and a USAF algorithm description.
Note that produced values are typically small on the order of 1e-7 and are in units of 1/second^2 (i.e., s-2).

Data for acceptance tests have been provided via PR metoppv/improver_test_data#105

Please advise of any necessary changes.
-Michael

Testing:

  • Ran tests and they passed OK
  • Added new tests for the new feature(s)

CLA

  • If a new developer, signed up to CLA

mmcpartland and others added 9 commits September 2, 2025 14:57
…ured the contents of the produced data are of a type float32 per the IMPROVER standard.
Note, the code passes all (both) test but the environmental variables that had to be set for the test to function include:
IMPROVER_IGNORE_CHECKSUMS set to false
IMPROVER_ACC_TEST_DIR set to blah..blah/usaf_improver/improver_test_data/data
….py, which creates the inputs, output, and SHA256 checksums for the turbulence acceptance checks, modified the file SHA256SUMS to contain the correlated checksums.

Corrected two file names in the test code.
The acceptance test (2) each passed.
Added back the changes to SHA file. These changes were made because of unexplained merge conflicts with the MET OFFICE master branch. This was the only feasible resolution.
@MoseleyS MoseleyS self-assigned this Oct 8, 2025
Copy link
Member

@MoseleyS MoseleyS left a comment

Choose a reason for hiding this comment

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

Thanks @mmcpartland. This PR meets most of our definition-of-done. I've made some suggestions that I hope you will see the benefit of, including pointing to some existing IMPROVER code that may help simplify this method.

Comment on lines +229 to +230
Taken from Ellrod and Knapp 1997: An objective clear-air turbulence forecasting technique: verification and
operational use. Wea. Forecasting, 7, 150-165.
Copy link
Member

Choose a reason for hiding this comment

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

Thanks for including the scientific reference. To make it more discoverable through the read-the-docs pages, please move this to a References section in the class doc-string.

Turbulence Index values are typically small on the order of 1e-7 and are in units of 1/second^2 (i.e., s-2).
"""

def _verify_compatible(self, cube1:Cube, cube2:Cube, ignore: Optional[Union[str, list[str]]]=[]) ->None:
Copy link
Member

Choose a reason for hiding this comment

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

Some of this method can be replaced with from improver.utilities.cube_checker import spatial_coords_match

I don't think the units check is needed as you are also converting units to the ones you need elsewhere.


# Grab data and specify consistent units for comparisons and later maths.
# Grab U components
u_wind_constraint = iris.Constraint(cube_func=lambda cube: cube.name().startswith("UWindComponent"))
Copy link
Member

Choose a reason for hiding this comment

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

The IMPROVER code base is written to use CF-conventions, we would expect to see the name "grid_eastward_wind" here (and "grid_northward_wind" below). This would then be compatible with the outputs of the ResolveWindComponents class.

v_winds = cubes.extract(v_wind_constraint)
self._convert_cubelist_units(v_winds, 'm s-1')

geopot_constraint = iris.Constraint(cube_func=lambda cube: cube.name().startswith("GeopotentialHeightAt"))
Copy link
Member

Choose a reason for hiding this comment

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

This CF name is "geopotential_height" or "geopotential_height_at_*"

# Grab data and specify consistent units for comparisons and later maths.
# Grab U components
u_wind_constraint = iris.Constraint(cube_func=lambda cube: cube.name().startswith("UWindComponent"))
u_winds = cubes.extract(u_wind_constraint)
Copy link
Member

Choose a reason for hiding this comment

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

It should be possible to combine the separated pressure levels into one Cube with a pressure coordinate. This makes the checks elsewhere simpler. Firstly, two identical pressure levels are not permitted, secondly, the pressure coordinates on each cube can be compared directly.

Suggested change
u_winds = cubes.extract(u_wind_constraint)
u_winds = cubes.extract(u_wind_constraint).merge_cube()

However, all other IMPROVER code expects pressure to already be a coordinate with decreasing data points, so I think we should be expecting data of that standard throughout this class.

# x = turbulence_index * 1e7 ; plt.close(); p = plt.imshow(x, vmax=15); plt.colorbar(p)

pressure_as_mb = u_wind_high_press.coord("pressure").cell(0).point
name_str = f"TurbulenceIndexAbove1500mAt{int(pressure_as_mb)}mb"
Copy link
Member

Choose a reason for hiding this comment

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

The output name should also follow cf-conventions, although this item doesn't appear there.

Suggested change
name_str = f"TurbulenceIndexAbove1500mAt{int(pressure_as_mb)}mb"
name_str = "atmosphere_turbulence"

The pressure item should be on a scalar pressure coordinate on the Cube. I'm not sure about the relevance of the "above 1500m" information, so I'm not sure whether the name string should be "atmosphere_turbulence_above_1500m" or if the height information should be an attribute, or a scalar height coordinate - which would be confusing alongside a pressure coordinate.

Comment on lines +101 to +104
# Coerce all pressures level references to the same units.
# Set to mb as that is used in the product name produced.
for c in cubes:
c.coord("pressure").convert_units('millibar')
Copy link
Member

Choose a reason for hiding this comment

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

Given that I have suggested later that this should not be in the name of the output cube, this may not be needed.

CubeList([u_wind_high_pres_cube, u_wind_low_pres_cube,
v_wind_high_pres_cube, v_wind_low_pres_cube,
geopotential_high_pres_cube, geopotential_low_pres_cube]) )
print(result)
Copy link
Member

Choose a reason for hiding this comment

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

We shouldn't leave a print statement in here

Suggested change
print(result)

geopotential_high_pres_cube, geopotential_low_pres_cube]) )
print(result)

assert np.allclose(result.data, expected_cube.data)
Copy link
Member

Choose a reason for hiding this comment

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

How many significant figures are actually important here? We recommend setting the tolerance here so that insignificant changes that result from new versions of maths libraries don't flag up false-positive errors.

Suggested change
assert np.allclose(result.data, expected_cube.data)
assert np.allclose(result.data, expected_cube.data, rtol=1e-4)

'v_wind_high_press', 'v_wind_low_press',
'geopot_high_press', 'geopot_low_press'])

class TurbulenceIndexAbove1500m_USAF(PostProcessingPlugin):
Copy link
Member

Choose a reason for hiding this comment

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

I recommend removing USAF from the class name. The bit that is relevant is that this is Ellrod's method, so if your thinking is that more turbulence classes may come in the future and you want to disambiguate now, you could call this TurbulenceIndexEllrod.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

3 participants