Skip to content
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.

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

Full-featured Mercator coord system #3041

Merged
merged 6 commits into from
Jun 6, 2018
Merged
Show file tree
Hide file tree
Changes from 2 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
35 changes: 28 additions & 7 deletions lib/iris/coord_systems.py
Original file line number Diff line number Diff line change
Expand Up @@ -831,7 +831,9 @@ class Mercator(CoordSystem):

grid_mapping_name = "mercator"

def __init__(self, longitude_of_projection_origin=0, ellipsoid=None):
def __init__(self, longitude_of_projection_origin=0.0, ellipsoid=None,
min_latitude=-80.0, max_latitude=84.0,
true_scale_lat=0.0):
"""
Constructs a Mercator coord system.

Expand All @@ -840,17 +842,33 @@ def __init__(self, longitude_of_projection_origin=0, ellipsoid=None):
True longitude of planar origin in degrees.
* ellipsoid
:class:`GeogCS` defining the ellipsoid.
* min_latitude
the maximum southerly extent of the projection.
Defaults to -80 degrees.
Copy link
Member

Choose a reason for hiding this comment

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

Where do these come from? If cartopy, I'd discourage default at all for these kwargs - they aren't fundamental parts of the reference system AFAIKS.

Copy link
Member Author

Choose a reason for hiding this comment

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

Actually, GRIB GDT10 doesn't define min/max latitude, so we could drop these again (we only really need the true scale latitude). Otherwise yes, these come from cartopy.

* max_latitude
the maximum northerly extent of the projection.
Defaults to 84 degrees.
* true_scale_lat
Copy link
Member

Choose a reason for hiding this comment

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

Is this the name used by CF? Could you reference these in your description to speed up xreferencing please?

Copy link
Member Author

Choose a reason for hiding this comment

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

Looks like CF accepts one of two names: http://cfconventions.org/cf-conventions/cf-conventions.html#_mercator. I'm going to go with scale_factor_at_projection_origin as it's the more similar to the use-case we're targetting.

the latitude where the scale is 1. Defaults to 0 degrees.

"""

#: True longitude of planar origin in degrees.
Copy link
Member

Choose a reason for hiding this comment

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

I agree that these aren't particularly well written, but I don't want to lose the docstrings altogether.
e.g. https://scitools.org.uk/iris/docs/latest/iris/iris/coord_systems.html#iris.coord_systems.Mercator.ellipsoid

Copy link
Member Author

Choose a reason for hiding this comment

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

Sorry! I always forget about these special comments (perhaps because it's so unclear that they do even have special behaviour). In fact, I thought they were just unnecessary duplications of the docstring above...

I'll reinstate them and add another for the new kwarg I've introduced.

self.longitude_of_projection_origin = longitude_of_projection_origin
#: Ellipsoid definition.
self.ellipsoid = ellipsoid
self.min_latitude = min_latitude
self.max_latitude = max_latitude
self.true_scale_lat = true_scale_lat

def __repr__(self):
res = "Mercator(longitude_of_projection_origin={!r}, ellipsoid={!r})"
return res.format(self.longitude_of_projection_origin, self.ellipsoid)
res = ("Mercator(longitude_of_projection_origin={!r}, "
Copy link
Member

Choose a reason for hiding this comment

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

I'm completely happy with your implementation. FWIW, if I was implementing this though, I'd have done:

"Mercator(longitude_of_projection_origin={self.longitude_of_projection_origin!r}...)".format(self=self)

Copy link
Member Author

Choose a reason for hiding this comment

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

Yeah, it was line length that put me off doing that...

Copy link
Member Author

Choose a reason for hiding this comment

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

Also, I think that's actually less readable than having unnamed args in the format string.

Copy link
Member

Choose a reason for hiding this comment

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

Also, I think that's actually less readable than having unnamed args in the format string.

That horse has bolted I'm afraid... https://www.python.org/dev/peps/pep-0498/
😜

Copy link
Member Author

Choose a reason for hiding this comment

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

That horse has bolted

Yes yes I know, and it makes full sense in the case of full format strings (i.e. f''). I just don't much like this intermediate situation 😉

"ellipsoid={!r}, "
"min_latitude={!r}, "
"max_latitude={!r}, "
"true_scale_lat={!r})")
return res.format(self.longitude_of_projection_origin,
self.ellipsoid,
self.min_latitude,
self.max_latitude,
self.true_scale_lat)

def as_cartopy_crs(self):
if self.ellipsoid is not None:
Expand All @@ -860,7 +878,10 @@ def as_cartopy_crs(self):

return ccrs.Mercator(
central_longitude=self.longitude_of_projection_origin,
globe=globe)
globe=globe,
min_latitude=self.min_latitude,
max_latitude=self.max_latitude,
latitude_true_scale=self.true_scale_lat)

def as_cartopy_projection(self):
return self.as_cartopy_crs()
Expand Down
2 changes: 1 addition & 1 deletion lib/iris/tests/results/coord_systems/Mercator.xml
Original file line number Diff line number Diff line change
@@ -1,2 +1,2 @@
<?xml version="1.0" ?>
<mercator ellipsoid="GeogCS(semi_major_axis=6377563.396, semi_minor_axis=6356256.909)" longitude_of_projection_origin="90.0"/>
<mercator ellipsoid="GeogCS(semi_major_axis=6377563.396, semi_minor_axis=6356256.909)" longitude_of_projection_origin="90.0" max_latitude="84.0" min_latitude="-80.0" true_scale_lat="0.0"/>
8 changes: 4 additions & 4 deletions lib/iris/tests/results/netcdf/netcdf_merc.cml
Original file line number Diff line number Diff line change
Expand Up @@ -53,15 +53,15 @@
45.5158, 45.9993]]" shape="(192, 192)" standard_name="longitude" units="Unit('degrees')" value_type="float32" var_name="lon"/>
</coord>
<coord datadims="[1]">
<dimCoord id="a2415886" points="[-5.16102e+06, -5.10719e+06, -5.05336e+06, ...,
<dimCoord id="da580c54" points="[-5.16102e+06, -5.10719e+06, -5.05336e+06, ...,
5.01299e+06, 5.06682e+06, 5.12065e+06]" shape="(192,)" standard_name="projection_x_coordinate" units="Unit('m')" value_type="float32" var_name="x">
<mercator ellipsoid="GeogCS(6378169.0)" longitude_of_projection_origin="0.0"/>
<mercator ellipsoid="GeogCS(6378169.0)" longitude_of_projection_origin="0.0" max_latitude="84.0" min_latitude="-80.0" true_scale_lat="0.0"/>
</dimCoord>
</coord>
<coord datadims="[0]">
<dimCoord id="41466ff9" points="[5.16101e+06, 5.10718e+06, 5.05335e+06, ...,
<dimCoord id="caf0a46a" points="[5.16101e+06, 5.10718e+06, 5.05335e+06, ...,
-5.01294e+06, -5.06678e+06, -5.12061e+06]" shape="(192,)" standard_name="projection_y_coordinate" units="Unit('m')" value_type="float32" var_name="y">
<mercator ellipsoid="GeogCS(6378169.0)" longitude_of_projection_origin="0.0"/>
<mercator ellipsoid="GeogCS(6378169.0)" longitude_of_projection_origin="0.0" max_latitude="84.0" min_latitude="-80.0" true_scale_lat="0.0"/>
</dimCoord>
</coord>
<coord>
Expand Down
58 changes: 0 additions & 58 deletions lib/iris/tests/test_coordsystem.py
Original file line number Diff line number Diff line change
Expand Up @@ -33,17 +33,12 @@
from iris.coord_systems import *
Copy link
Member

Choose a reason for hiding this comment

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

Wow. These are some really old tests (and really old-school style). 😉

Copy link
Member Author

Choose a reason for hiding this comment

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

Whoa 😨





def osgb():
return TransverseMercator(latitude_of_projection_origin=49, longitude_of_central_meridian=-2,
false_easting=-400, false_northing=100,
scale_factor_at_central_meridian=0.9996012717,
ellipsoid=GeogCS(6377563.396, 6356256.909))

def merc():
return Mercator(longitude_of_projection_origin=90.0,
ellipsoid=GeogCS(6377563.396, 6356256.909))

def stereo():
return Stereographic(central_lat=-90, central_lon=-45,
Expand Down Expand Up @@ -391,58 +386,5 @@ def test_south_cutoff(self):
self.assertEqual(ccrs.cutoff, 30)


class Test_Mercator_construction(tests.IrisTest):
def test_merc(self):
tm = merc()
self.assertXMLElement(tm, ("coord_systems", "Mercator.xml"))


class Test_Mercator_repr(tests.IrisTest):
def test_merc(self):
tm = merc()
expected = "Mercator(longitude_of_projection_origin=90.0, "\
"ellipsoid=GeogCS(semi_major_axis=6377563.396, "\
"semi_minor_axis=6356256.909))"
self.assertEqual(expected, repr(tm))


class Test_Mercator_as_cartopy_crs(tests.IrisTest):
def test_as_cartopy_crs(self):
longitude_of_projection_origin = 90.0
ellipsoid = GeogCS(semi_major_axis=6377563.396,
semi_minor_axis=6356256.909)

merc_cs = Mercator(
longitude_of_projection_origin,
ellipsoid=ellipsoid)

expected = ccrs.Mercator(
central_longitude=longitude_of_projection_origin,
globe=ccrs.Globe(semimajor_axis=6377563.396,
semiminor_axis=6356256.909, ellipse=None))

res = merc_cs.as_cartopy_crs()
self.assertEqual(res, expected)


class Test_Mercator_as_cartopy_projection(tests.IrisTest):
def test_as_cartopy_projection(self):
longitude_of_projection_origin = 90.0
ellipsoid = GeogCS(semi_major_axis=6377563.396,
semi_minor_axis=6356256.909)

merc_cs = Mercator(
longitude_of_projection_origin,
ellipsoid=ellipsoid)

expected = ccrs.Mercator(
central_longitude=longitude_of_projection_origin,
globe=ccrs.Globe(semimajor_axis=6377563.396,
semiminor_axis=6356256.909, ellipse=None))

res = merc_cs.as_cartopy_projection()
self.assertEqual(res, expected)


if __name__ == "__main__":
tests.main()
116 changes: 116 additions & 0 deletions lib/iris/tests/unit/coord_systems/test_Mercator.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,116 @@
# (C) British Crown Copyright 2018, 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 the :class:`iris.coord_systems.Mercator` class."""

from __future__ import (absolute_import, division, print_function)
from six.moves import (filter, input, map, range, zip) # noqa

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

import cartopy.crs as ccrs
from iris.coord_systems import GeogCS, Mercator


class Test_Mercator__basics(tests.IrisTest):
def setUp(self):
self.tm = Mercator(longitude_of_projection_origin=90.0,
ellipsoid=GeogCS(6377563.396, 6356256.909))

def test_construction(self):
self.assertXMLElement(self.tm, ("coord_systems", "Mercator.xml"))

def test_repr(self):
expected = ("Mercator(longitude_of_projection_origin=90.0, "
"ellipsoid=GeogCS(semi_major_axis=6377563.396, "
"semi_minor_axis=6356256.909), "
"min_latitude=-80.0, max_latitude=84.0, "
"true_scale_lat=0.0)")
self.assertEqual(expected, repr(self.tm))


class Test_Mercator__as_cartopy_crs(tests.IrisTest):
def test_simple(self):
# Check that a projection set up with all the defaults is correctly
# converted to a cartopy CRS.
merc_cs = Mercator()
res = merc_cs.as_cartopy_crs()
expected = ccrs.Mercator(globe=ccrs.Globe())
self.assertEqual(res, expected)

def test_extra_kwargs(self):
# Check that a projection with non-default values is correctly
# converted to a cartopy CRS.
longitude_of_projection_origin = 90.0
min_lat = -75.0
true_scale_lat = 14.0
ellipsoid = GeogCS(semi_major_axis=6377563.396,
semi_minor_axis=6356256.909)

merc_cs = Mercator(
longitude_of_projection_origin,
ellipsoid=ellipsoid,
min_latitude=min_lat,
true_scale_lat=true_scale_lat)

expected = ccrs.Mercator(
central_longitude=longitude_of_projection_origin,
globe=ccrs.Globe(semimajor_axis=6377563.396,
semiminor_axis=6356256.909, ellipse=None),
min_latitude=min_lat,
latitude_true_scale=true_scale_lat)

res = merc_cs.as_cartopy_crs()
self.assertEqual(res, expected)


class Test_as_cartopy_projection(tests.IrisTest):
def test_simple(self):
# Check that a projection set up with all the defaults is correctly
# converted to a cartopy projection.
merc_cs = Mercator()
res = merc_cs.as_cartopy_projection()
expected = ccrs.Mercator(globe=ccrs.Globe())
self.assertEqual(res, expected)

def test_extra_kwargs(self):
longitude_of_projection_origin = 90.0
max_lat = 80.0
true_scale_lat = 14.0
ellipsoid = GeogCS(semi_major_axis=6377563.396,
semi_minor_axis=6356256.909)

merc_cs = Mercator(
longitude_of_projection_origin,
ellipsoid=ellipsoid,
max_latitude=max_lat,
true_scale_lat=true_scale_lat)

expected = ccrs.Mercator(
central_longitude=longitude_of_projection_origin,
globe=ccrs.Globe(semimajor_axis=6377563.396,
semiminor_axis=6356256.909, ellipse=None),
max_latitude=max_lat,
latitude_true_scale=true_scale_lat)

res = merc_cs.as_cartopy_projection()
self.assertEqual(res, expected)


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