Skip to content
Closed
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
21 changes: 15 additions & 6 deletions lib/iris/fileformats/grib/_save_rules.py
Original file line number Diff line number Diff line change
Expand Up @@ -178,7 +178,6 @@ def identification(cube, grib):


def shape_of_the_earth(cube, grib):
# assume latlon
cs = cube.coord(dimensions=[0]).coord_system

# Initially set shape_of_earth keys to missing (255 for byte, -1 for long).
Expand All @@ -200,13 +199,23 @@ def shape_of_the_earth(cube, grib):

# Spherical earth.
if ellipsoid.inverse_flattening == 0.0:
gribapi.grib_set_long(grib, "shapeOfTheEarth", 1)
gribapi.grib_set_long(grib, "scaleFactorOfRadiusOfSphericalEarth", 0)
gribapi.grib_set_long(grib, "scaledValueOfRadiusOfSphericalEarth",
ellipsoid.semi_major_axis)
# Set shapeOfTheEarth key based on value of spherical earth radius.
if int(ellipsoid.semi_major_axis) == 6367470:
gribapi.grib_set_long(grib, "shapeOfTheEarth", 0)
elif int(ellipsoid.semi_major_axis) == 6371229:
gribapi.grib_set_long(grib, "shapeOfTheEarth", 6)
else:
gribapi.grib_set_long(grib, "shapeOfTheEarth", 1)
gribapi.grib_set_long(grib,
"scaleFactorOfRadiusOfSphericalEarth", 0)
gribapi.grib_set_long(grib, "scaledValueOfRadiusOfSphericalEarth",
ellipsoid.semi_major_axis)
# Oblate spheroid earth.
else:
gribapi.grib_set_long(grib, "shapeOfTheEarth", 7)
if isinstance(cs, iris.coord_systems.OSGB):
gribapi.grib_set_long(grib, "shapeOfTheEarth", 9)
else:
gribapi.grib_set_long(grib, "shapeOfTheEarth", 7)
Copy link
Member

Choose a reason for hiding this comment

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

I know 7 was the previous default, but why not 3?

Copy link
Member

Choose a reason for hiding this comment

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

this appears to me to be due to a default in Cartopy, which expects the units to be metres, not km

However this doesn't appear to be clear in the docstrings
#1491
I will try to get clarity on this

Copy link
Member Author

Choose a reason for hiding this comment

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

I think I left 3 out deliberately because separating 3 and 7 is conceptually quite hard as they are effectively identical keys with a difference of a factor of 1000 between them. If we currently support units of m then it seems extraneous to add or change the behaviour to expect km instead...

gribapi.grib_set_long(grib, "scaleFactorOfEarthMajorAxis", 0)
gribapi.grib_set_long(grib, "scaledValueOfEarthMajorAxis",
ellipsoid.semi_major_axis)
Expand Down
2 changes: 1 addition & 1 deletion lib/iris/tests/integration/test_grib2.py
Original file line number Diff line number Diff line change
Expand Up @@ -140,7 +140,7 @@ def test_save_load(self):
'gridDefinitionTemplateNumber = 5',
'Ni = {:d}'.format(cube.shape[-1]),
'Nj = {:d}'.format(cube.shape[-2]),
'shapeOfTheEarth = 1',
'shapeOfTheEarth = 6',
'scaledValueOfRadiusOfSphericalEarth = {:d}'.format(
int(UM_DEFAULT_EARTH_RADIUS)),
'resolutionAndComponentFlags = 0',
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,6 @@ long [productionStatusOfProcessedData]: [0] != [255]
[section2Length] not found in 2nd field
long [numberOfSection] 6 out of 7 different
long [numberOfSection] 5 out of 7 different
long [shapeOfTheEarth]: [0] != [1]
scaleFactorOfRadiusOfSphericalEarth is set to missing in 1st field is not missing in 2nd field
scaledValueOfRadiusOfSphericalEarth is set to missing in 1st field is not missing in 2nd field
long [numberOfSection] 6 out of 7 different
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,6 @@ long [productionStatusOfProcessedData]: [0] != [255]
[section2Length] not found in 2nd field
long [numberOfSection] 6 out of 7 different
long [numberOfSection] 5 out of 7 different
long [shapeOfTheEarth]: [0] != [1]
scaleFactorOfRadiusOfSphericalEarth is set to missing in 1st field is not missing in 2nd field
scaledValueOfRadiusOfSphericalEarth is set to missing in 1st field is not missing in 2nd field
long [latitudeOfLastGridPoint]: [19419996] != [19419285]
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,6 @@ long [productionStatusOfProcessedData]: [0] != [255]
[section2Length] not found in 2nd field
long [numberOfSection] 6 out of 7 different
long [numberOfSection] 5 out of 7 different
long [shapeOfTheEarth]: [0] != [1]
scaleFactorOfRadiusOfSphericalEarth is set to missing in 1st field is not missing in 2nd field
scaledValueOfRadiusOfSphericalEarth is set to missing in 1st field is not missing in 2nd field
long [latitudeOfLastGridPoint]: [-89999938] != [-89999944]
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -42,25 +42,6 @@ def test__template_number(self):
grid_definition_template_0(self.test_cube, self.mock_grib)
self._check_key('gridDefinitionTemplateNumber', 0)

def test__shape_of_earth_spherical(self):
cs = GeogCS(semi_major_axis=1.23)
test_cube = self._make_test_cube(cs=cs)
grid_definition_template_0(test_cube, self.mock_grib)
self._check_key('shapeOfTheEarth', 1)
self._check_key('scaleFactorOfRadiusOfSphericalEarth', 0)
self._check_key('scaledValueOfRadiusOfSphericalEarth', 1.23)

def test__shape_of_earth_flattened(self):
cs = GeogCS(semi_major_axis=1.456,
semi_minor_axis=1.123)
test_cube = self._make_test_cube(cs=cs)
grid_definition_template_0(test_cube, self.mock_grib)
self._check_key('shapeOfTheEarth', 7)
self._check_key('scaleFactorOfEarthMajorAxis', 0)
self._check_key('scaledValueOfEarthMajorAxis', 1.456)
self._check_key('scaleFactorOfEarthMinorAxis', 0)
self._check_key('scaledValueOfEarthMinorAxis', 1.123)

def test__grid_shape(self):
test_cube = self._make_test_cube(x_points=np.arange(13),
y_points=np.arange(6))
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -52,30 +52,6 @@ def test__template_number(self):
grid_definition_template_1(self.test_cube, self.mock_grib)
self._check_key('gridDefinitionTemplateNumber', 1)

def test__shape_of_earth_spherical(self):
ellipsoid = GeogCS(1.23)
cs = RotatedGeogCS(grid_north_pole_latitude=90.0,
grid_north_pole_longitude=0.0,
ellipsoid=ellipsoid)
test_cube = self._make_test_cube(cs=cs)
grid_definition_template_1(test_cube, self.mock_grib)
self._check_key('shapeOfTheEarth', 1)
self._check_key('scaleFactorOfRadiusOfSphericalEarth', 0)
self._check_key('scaledValueOfRadiusOfSphericalEarth', 1.23)

def test__shape_of_earth_flattened(self):
ellipsoid = GeogCS(semi_major_axis=1.456, semi_minor_axis=1.123)
cs = RotatedGeogCS(grid_north_pole_latitude=90.0,
grid_north_pole_longitude=0.0,
ellipsoid=ellipsoid)
test_cube = self._make_test_cube(cs=cs)
grid_definition_template_1(test_cube, self.mock_grib)
self._check_key('shapeOfTheEarth', 7)
self._check_key('scaleFactorOfEarthMajorAxis', 0)
self._check_key('scaledValueOfEarthMajorAxis', 1.456)
self._check_key('scaleFactorOfEarthMinorAxis', 0)
self._check_key('scaledValueOfEarthMinorAxis', 1.123)

def test__grid_shape(self):
test_cube = self._make_test_cube(x_points=np.arange(13),
y_points=np.arange(6))
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -61,14 +61,6 @@ def test__template_number(self):
grid_definition_template_12(self.test_cube, self.mock_grib)
self._check_key('gridDefinitionTemplateNumber', 12)

def test__shape_of_earth(self):
grid_definition_template_12(self.test_cube, self.mock_grib)
self._check_key('shapeOfTheEarth', 7)
self._check_key('scaleFactorOfEarthMajorAxis', 0)
self._check_key('scaledValueOfEarthMajorAxis', 6377563.396)
self._check_key('scaleFactorOfEarthMinorAxis', 0)
self._check_key('scaledValueOfEarthMinorAxis', 6356256.909)

def test__grid_shape(self):
test_cube = self._make_test_cube(x_points=np.arange(13),
y_points=np.arange(6))
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -56,29 +56,6 @@ def test__template_number(self):
grid_definition_template_5(self.test_cube, self.mock_grib)
self._check_key('gridDefinitionTemplateNumber', 5)

def test__shape_of_earth_spherical(self):
cs = RotatedGeogCS(grid_north_pole_latitude=90.0,
grid_north_pole_longitude=0.0,
ellipsoid=GeogCS(52431.0))
test_cube = self._make_test_cube(cs=cs)
grid_definition_template_5(test_cube, self.mock_grib)
self._check_key('shapeOfTheEarth', 1)
self._check_key('scaleFactorOfRadiusOfSphericalEarth', 0)
self._check_key('scaledValueOfRadiusOfSphericalEarth', 52431.0)

def test__shape_of_earth_flattened(self):
ellipsoid = GeogCS(semi_major_axis=1456.0, semi_minor_axis=1123.0)
cs = RotatedGeogCS(grid_north_pole_latitude=90.0,
grid_north_pole_longitude=0.0,
ellipsoid=ellipsoid)
test_cube = self._make_test_cube(cs=cs)
grid_definition_template_5(test_cube, self.mock_grib)
self._check_key('shapeOfTheEarth', 7)
self._check_key('scaleFactorOfEarthMajorAxis', 0)
self._check_key('scaledValueOfEarthMajorAxis', 1456.0)
self._check_key('scaleFactorOfEarthMinorAxis', 0)
self._check_key('scaledValueOfEarthMinorAxis', 1123.0)

def test__grid_shape(self):
test_cube = self._make_test_cube(x_points=np.arange(13),
y_points=np.arange(6))
Expand Down
Original file line number Diff line number Diff line change
@@ -0,0 +1,142 @@
# (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 <http://www.gnu.org/licenses/>.
"""
Unit tests for
:meth:`iris.fileformats.grib._save_rules.shape_of_the_earth`.

"""

from __future__ import (absolute_import, division, print_function)

# Import iris.tests first so that some things can be initialised before
# importing anything else.
import iris.tests as tests

import numpy as np

from iris.coord_systems import GeogCS, RotatedGeogCS, TransverseMercator, OSGB
from iris.exceptions import TranslationError
from iris.tests.unit.fileformats.grib.save_rules import GdtTestMixin

from iris.fileformats.grib._save_rules import shape_of_the_earth


class Test(tests.IrisTest, GdtTestMixin):
def setUp(self):
GdtTestMixin.setUp(self)

def _spherical_earth_test_common(self, radius):
self._check_key('scaleFactorOfRadiusOfSphericalEarth', 0)
self._check_key('scaledValueOfRadiusOfSphericalEarth', radius)

def _oblate_spheroid_earth_test_common(self, semi_major_axis,
semi_minor_axis):
self._check_key('scaleFactorOfEarthMajorAxis', 0)
self._check_key('scaledValueOfEarthMajorAxis', semi_major_axis)
self._check_key('scaleFactorOfEarthMinorAxis', 0)
self._check_key('scaledValueOfEarthMinorAxis', semi_minor_axis)

def test_radius_of_earth_6367470(self):
# Test setting shapeOfTheEarth = 0
radius = 6367470
cs = GeogCS(semi_major_axis=radius)
test_cube = self._make_test_cube(cs=cs)
shape_of_the_earth(test_cube, self.mock_grib)
self._check_key('shapeOfTheEarth', 0)
self._check_key('scaleFactorOfRadiusOfSphericalEarth', 255)
self._check_key('scaledValueOfRadiusOfSphericalEarth', -1)

def test_radius_of_earth_6371229(self):
# Test setting shapeOfTheEarth = 6
radius = 6371229
cs = GeogCS(semi_major_axis=radius)
test_cube = self._make_test_cube(cs=cs)
shape_of_the_earth(test_cube, self.mock_grib)
self._check_key('shapeOfTheEarth', 6)
self._check_key('scaleFactorOfRadiusOfSphericalEarth', 255)
self._check_key('scaledValueOfRadiusOfSphericalEarth', -1)

def test_spherical_earth(self):
# Test setting shapeOfTheEarth = 1
radius = 1.23
cs = GeogCS(semi_major_axis=radius)
test_cube = self._make_test_cube(cs=cs)
shape_of_the_earth(test_cube, self.mock_grib)
self._check_key('shapeOfTheEarth', 1)
self._spherical_earth_test_common(radius)

def test_oblate_spheroid_earth(self):
# Test setting shapeOfTheEarth = 7
semi_major_axis = 1.456
semi_minor_axis = 1.123
cs = GeogCS(semi_major_axis=semi_major_axis,
semi_minor_axis=semi_minor_axis)
test_cube = self._make_test_cube(cs=cs)
shape_of_the_earth(test_cube, self.mock_grib)
self._check_key('shapeOfTheEarth', 7)
self._oblate_spheroid_earth_test_common(semi_major_axis,
semi_minor_axis)

def test_OSGB(self):
# Test setting shapeOfTheEarth = 9
cs = OSGB()
# The following are fixed for the OSGB coord system.
semi_major_axis, semi_minor_axis = (6377563.396, 6356256.909)
test_cube = self._make_test_cube(cs=cs)
shape_of_the_earth(test_cube, self.mock_grib)
self._check_key('shapeOfTheEarth', 9)
self._oblate_spheroid_earth_test_common(semi_major_axis,
semi_minor_axis)

def test_TransverseMercator_spherical(self):
# Test setting shapeOfTheEarth = 1 with a non-GeogCS coord system.
cs = TransverseMercator(49, -2, 400000, -100000, 0.9996012717,
ellipsoid=GeogCS(6377563.396))
radius = 6377563.396
test_cube = self._make_test_cube(cs=cs)
shape_of_the_earth(test_cube, self.mock_grib)
self._check_key('shapeOfTheEarth', 1)
self._spherical_earth_test_common(radius)

def test_RotatedGeogCS_spherical(self):
# Test setting shapeOfTheEarth = 1 with a RotatedGeogCS coord system.
radius = 52431.0
cs = RotatedGeogCS(grid_north_pole_latitude=90.0,
grid_north_pole_longitude=0.0,
ellipsoid=GeogCS(radius))
test_cube = self._make_test_cube(cs=cs)
shape_of_the_earth(test_cube, self.mock_grib)
self._check_key('shapeOfTheEarth', 1)
self._spherical_earth_test_common(radius)

def test_RotatedGeogCS_oblate_spheroid(self):
# Test setting shapeOfTheEarth = 7 with a RotatedGeogCS coord system.
semi_major_axis = 1456.0
semi_minor_axis = 1123.0
ellipsoid = GeogCS(semi_major_axis=semi_major_axis,
semi_minor_axis=semi_minor_axis)
cs = RotatedGeogCS(grid_north_pole_latitude=90.0,
grid_north_pole_longitude=0.0,
ellipsoid=ellipsoid)
test_cube = self._make_test_cube(cs=cs)
shape_of_the_earth(test_cube, self.mock_grib)
self._check_key('shapeOfTheEarth', 7)
self._oblate_spheroid_earth_test_common(semi_major_axis,
semi_minor_axis)

if __name__ == "__main__":
tests.main()