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.

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

Jump to bottom

Improve warped rectangular projections #1180

Merged
merged 6 commits into from
Nov 14, 2018
Merged

Improve warped rectangular projections #1180

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
7efcb8f
Simplify some bits of Robinson tests.
QuLogic Nov 13, 2018
aac843a
Warn about ellipses in warped rectangular projections.
QuLogic Nov 13, 2018
60065b3
Add docs for Mollweide and Robinson.
QuLogic Nov 13, 2018
17d901e
Add standard tests for Robinson projection.
QuLogic Nov 13, 2018
a33c594
Support offsets on warped rectangular projections.
QuLogic Nov 13, 2018
08982cc
Add types to Eckert docs.
QuLogic Nov 13, 2018
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
113 changes: 105 additions & 8 deletions lib/cartopy/crs.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1521,7 +1521,12 @@ def y_limits(self):


class _WarpedRectangularProjection(six.with_metaclass(ABCMeta, Projection)):
def __init__(self, proj4_params, central_longitude, globe=None):
def __init__(self, proj4_params, central_longitude,
false_easting=None, false_northing=None, globe=None):
if false_easting is not None:
proj4_params += [('x_0', false_easting)]
if false_northing is not None:
proj4_params += [('y_0', false_northing)]
super(_WarpedRectangularProjection, self).__init__(proj4_params,
globe=globe)

Expand Down Expand Up @@ -1567,15 +1572,22 @@ class _Eckert(six.with_metaclass(ABCMeta, _WarpedRectangularProjection)):

"""

def __init__(self, central_longitude=0, globe=None):
def __init__(self, central_longitude=0, false_easting=None,
false_northing=None, globe=None):
"""
Parameters
----------
central_longitude: optional
central_longitude: float, optional
The central longitude. Defaults to 0.
globe: optional
A :class:`cartopy.crs.Globe`. If omitted, a default globe is
created.
false_easting: float, optional
X offset from planar origin in metres. Defaults to 0.
false_northing: float, optional
Y offset from planar origin in metres. Defaults to 0.
globe: :class:`cartopy.crs.Globe`, optional
If omitted, a default globe is created.

.. note::
This projection does not handle elliptical globes.

"""
if globe is None:
Expand All @@ -1592,6 +1604,8 @@ def __init__(self, central_longitude=0, globe=None):
proj4_params = [('proj', self._proj_name),
('lon_0', central_longitude)]
super(_Eckert, self).__init__(proj4_params, central_longitude,
false_easting=false_easting,
false_northing=false_northing,
globe=globe)

@property
Expand Down Expand Up @@ -1676,9 +1690,51 @@ class EckertVI(_Eckert):


class Mollweide(_WarpedRectangularProjection):
def __init__(self, central_longitude=0, globe=None):
"""
A Mollweide projection.

This projection is pseudocylindrical, and equal area. Parallels are
unequally-spaced straight lines, while meridians are elliptical arcs up to
semicircles on the edges. Poles are points.

It is commonly used for world maps, or interrupted with several central
meridians.

"""

def __init__(self, central_longitude=0, globe=None,
false_easting=None, false_northing=None):
"""
Parameters
----------
central_longitude: float, optional
The central longitude. Defaults to 0.
false_easting: float, optional
X offset from planar origin in metres. Defaults to 0.
false_northing: float, optional
Y offset from planar origin in metres. Defaults to 0.
globe: :class:`cartopy.crs.Globe`, optional
If omitted, a default globe is created.

.. note::
This projection does not handle elliptical globes.

"""
if globe is None:
globe = Globe(semimajor_axis=WGS84_SEMIMAJOR_AXIS, ellipse=None)

# TODO: Let the globe return the semimajor axis always.
a = globe.semimajor_axis or WGS84_SEMIMAJOR_AXIS
b = globe.semiminor_axis or a

if b != a or globe.ellipse is not None:
warnings.warn('The proj "moll" projection does not handle '
'elliptical globes.')

proj4_params = [('proj', 'moll'), ('lon_0', central_longitude)]
super(Mollweide, self).__init__(proj4_params, central_longitude,
false_easting=false_easting,
false_northing=false_northing,
globe=globe)

@property
Expand All @@ -1687,7 +1743,35 @@ def threshold(self):


class Robinson(_WarpedRectangularProjection):
def __init__(self, central_longitude=0, globe=None):
"""
A Robinson projection.

This projection is pseudocylindrical, and a compromise that is neither
equal-area nor conformal. Parallels are unequally-spaced straight lines,
and meridians are curved lines of no particular form.

It is commonly used for "visually-appealing" world maps.

"""

def __init__(self, central_longitude=0, globe=None,
false_easting=None, false_northing=None):
"""
Parameters
----------
central_longitude: float, optional
The central longitude. Defaults to 0.
false_easting: float, optional
X offset from planar origin in metres. Defaults to 0.
false_northing: float, optional
Y offset from planar origin in metres. Defaults to 0.
globe: :class:`cartopy.crs.Globe`, optional
If omitted, a default globe is created.

.. note::
This projection does not handle elliptical globes.

"""
# Warn when using Robinson with proj 4.8 due to discontinuity at
# 40 deg N introduced by incomplete fix to issue #113 (see
# https://github.com/OSGeo/proj.4/issues/113).
Expand All @@ -1702,8 +1786,21 @@ def __init__(self, central_longitude=0, globe=None):
'projection may be unreliable and should be used '
'with caution.')

if globe is None:
globe = Globe(semimajor_axis=WGS84_SEMIMAJOR_AXIS, ellipse=None)
pelson marked this conversation as resolved.
Show resolved Hide resolved

# TODO: Let the globe return the semimajor axis always.
pelson marked this conversation as resolved.
Show resolved Hide resolved
a = globe.semimajor_axis or WGS84_SEMIMAJOR_AXIS
b = globe.semiminor_axis or a

if b != a or globe.ellipse is not None:
warnings.warn('The proj "robin" projection does not handle '
'elliptical globes.')

proj4_params = [('proj', 'robin'), ('lon_0', central_longitude)]
super(Robinson, self).__init__(proj4_params, central_longitude,
false_easting=false_easting,
false_northing=false_northing,
globe=globe)

@property
Expand Down
17 changes: 17 additions & 0 deletions lib/cartopy/tests/crs/test_eckert.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -51,6 +51,23 @@ def test_default(name, proj, lim):
assert_almost_equal(eck.y_limits, [-lim / 2, lim / 2])


@pytest.mark.parametrize('name, proj', [
pytest.param('eck1', ccrs.EckertI, id='EckertI'),
pytest.param('eck2', ccrs.EckertII, id='EckertII'),
pytest.param('eck3', ccrs.EckertIII, id='EckertIII'),
pytest.param('eck4', ccrs.EckertIV, id='EckertIV'),
pytest.param('eck5', ccrs.EckertV, id='EckertV'),
pytest.param('eck6', ccrs.EckertVI, id='EckertVI'),
])
def test_offset(name, proj):
crs = proj()
crs_offset = proj(false_easting=1234, false_northing=-4321)
other_args = {'a=6378137.0', 'lon_0=0', 'x_0=1234', 'y_0=-4321'}
check_proj4_params(name, crs_offset, other_args)
assert tuple(np.array(crs.x_limits) + 1234) == crs_offset.x_limits
assert tuple(np.array(crs.y_limits) - 4321) == crs_offset.y_limits


@pytest.mark.parametrize('name, proj, lim', [
pytest.param('eck1', ccrs.EckertI, 18460911.739778, id='EckertI'),
pytest.param('eck2', ccrs.EckertII, 18460911.739778, id='EckertII'),
Expand Down
22 changes: 9 additions & 13 deletions lib/cartopy/tests/crs/test_mollweide.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -36,7 +36,7 @@ def check_proj4_params(crs, other_args):

def test_default():
moll = ccrs.Mollweide()
other_args = {'ellps=WGS84', 'lon_0=0'}
other_args = {'a=6378137.0', 'lon_0=0'}
check_proj4_params(moll, other_args)

assert_almost_equal(np.array(moll.x_limits),
Expand All @@ -45,23 +45,19 @@ def test_default():
[-9020047.8480736, 9020047.8480736])


def test_eccentric_globe():
globe = ccrs.Globe(semimajor_axis=1000, semiminor_axis=500,
ellipse=None)
moll = ccrs.Mollweide(globe=globe)
other_args = {'a=1000', 'b=500', 'lon_0=0'}
check_proj4_params(moll, other_args)

assert_almost_equal(np.array(moll.x_limits),
[-2828.4271247, 2828.4271247])
assert_almost_equal(np.array(moll.y_limits),
[-1414.2135623731, 1414.2135623731])
def test_offset():
crs = ccrs.Mollweide()
crs_offset = ccrs.Mollweide(false_easting=1234, false_northing=-4321)
other_args = {'a=6378137.0', 'lon_0=0', 'x_0=1234', 'y_0=-4321'}
check_proj4_params(crs_offset, other_args)
assert tuple(np.array(crs.x_limits) + 1234) == crs_offset.x_limits
assert tuple(np.array(crs.y_limits) - 4321) == crs_offset.y_limits


@pytest.mark.parametrize('lon', [-10.0, 10.0])
def test_central_longitude(lon):
moll = ccrs.Mollweide(central_longitude=lon)
other_args = {'ellps=WGS84', 'lon_0={}'.format(lon)}
other_args = {'a=6378137.0', 'lon_0={}'.format(lon)}
check_proj4_params(moll, other_args)

assert_almost_equal(np.array(moll.x_limits),
Expand Down
89 changes: 69 additions & 20 deletions lib/cartopy/tests/crs/test_robinson.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -14,76 +14,125 @@
#
# You should have received a copy of the GNU Lesser General Public License
# along with cartopy. If not, see <https://www.gnu.org/licenses/>.
'''
"""
Tests for Robinson projection.

For now, mostly tests the workaround for a specific problem.
Problem report in : https://github.com/SciTools/cartopy/issues/23
Fix covered in : https://github.com/SciTools/cartopy/pull/277

'''
"""

from __future__ import (absolute_import, division, print_function)

import numpy as np
from numpy.testing import assert_array_almost_equal as assert_arr_almost_eq
from numpy.testing import assert_almost_equal, assert_array_almost_equal
import pytest

import cartopy.crs as ccrs


_NAN = float('nan')
_CRS_PC = ccrs.PlateCarree()
_CRS_ROB = ccrs.Robinson()

# Increase tolerance if using older proj releases
_TOL = -1 if ccrs.PROJ4_VERSION < (4, 9) else 7
_LIMIT_TOL = -1 # if ccrs.PROJ4_VERSION < (5, 2, 0) else 7


def check_proj_params(crs, other_args):
expected = other_args | {'proj=robin', 'no_defs'}
proj_params = set(crs.proj4_init.lstrip('+').split(' +'))
assert expected == proj_params


def test_default():
robin = ccrs.Robinson()
other_args = {'a=6378137.0', 'lon_0=0'}
check_proj_params(robin, other_args)

assert_almost_equal(robin.x_limits,
[-17005833.3305252, 17005833.3305252])
assert_almost_equal(robin.y_limits,
[-8625154.6651000, 8625154.6651000], _LIMIT_TOL)


def test_offset():
crs = ccrs.Robinson()
crs_offset = ccrs.Robinson(false_easting=1234, false_northing=-4321)
other_args = {'a=6378137.0', 'lon_0=0', 'x_0=1234', 'y_0=-4321'}
check_proj_params(crs_offset, other_args)
assert tuple(np.array(crs.x_limits) + 1234) == crs_offset.x_limits
assert tuple(np.array(crs.y_limits) - 4321) == crs_offset.y_limits


@pytest.mark.parametrize('lon', [-10.0, 10.0])
def test_central_longitude(lon):
robin = ccrs.Robinson(central_longitude=lon)
other_args = {'a=6378137.0', 'lon_0={}'.format(lon)}
check_proj_params(robin, other_args)

assert_almost_equal(robin.x_limits,
[-17005833.3305252, 17005833.3305252],
decimal=5)
assert_almost_equal(robin.y_limits,
[-8625154.6651000, 8625154.6651000], _LIMIT_TOL)


def test_transform_point():
"""
Mostly tests the workaround for a specific problem.
Problem report in: https://github.com/SciTools/cartopy/issues/23
Fix covered in: https://github.com/SciTools/cartopy/pull/277
"""

# this way has always worked
result = _CRS_ROB.transform_point(35.0, 70.0, _CRS_PC)
assert_arr_almost_eq(result, (2376187.27182751, 7275317.81573085), _TOL)
assert_array_almost_equal(result, (2376187.27182751, 7275317.81573085),
_TOL)

# this always did something, but result has altered
result = _CRS_ROB.transform_point(_NAN, 70.0, _CRS_PC)
result = _CRS_ROB.transform_point(np.nan, 70.0, _CRS_PC)
assert np.all(np.isnan(result))

# this used to crash + is now fixed
result = _CRS_ROB.transform_point(35.0, _NAN, _CRS_PC)
result = _CRS_ROB.transform_point(35.0, np.nan, _CRS_PC)
assert np.all(np.isnan(result))


def test_transform_points():
"""
Mostly tests the workaround for a specific problem.
Problem report in: https://github.com/SciTools/cartopy/issues/23
Copy link
Member

Choose a reason for hiding this comment

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

Missing a 2 at the end (issue 232, I believe)

Fix covered in: https://github.com/SciTools/cartopy/pull/277
"""

# these always worked
result = _CRS_ROB.transform_points(_CRS_PC,
np.array([35.0]),
np.array([70.0]))
assert_arr_almost_eq(result,
[[2376187.27182751, 7275317.81573085, 0]], _TOL)
assert_array_almost_equal(result,
[[2376187.27182751, 7275317.81573085, 0]], _TOL)

result = _CRS_ROB.transform_points(_CRS_PC,
np.array([35.0]),
np.array([70.0]),
np.array([0.0]))
assert_arr_almost_eq(result,
[[2376187.27182751, 7275317.81573085, 0]], _TOL)
assert_array_almost_equal(result,
[[2376187.27182751, 7275317.81573085, 0]], _TOL)

# this always did something, but result has altered
result = _CRS_ROB.transform_points(_CRS_PC,
np.array([_NAN]),
np.array([np.nan]),
np.array([70.0]))
assert np.all(np.isnan(result))

# this used to crash + is now fixed
result = _CRS_ROB.transform_points(_CRS_PC,
np.array([35.0]),
np.array([_NAN]))
np.array([np.nan]))
assert np.all(np.isnan(result))

# multipoint case
x = np.array([10.0, 21.0, 0.0, 77.7, _NAN, 0.0])
y = np.array([10.0, _NAN, 10.0, 77.7, 55.5, 0.0])
z = np.array([10.0, 0.0, 0.0, _NAN, 55.5, 0.0])
x = np.array([10.0, 21.0, 0.0, 77.7, np.nan, 0.0])
y = np.array([10.0, np.nan, 10.0, 77.7, 55.5, 0.0])
z = np.array([10.0, 0.0, 0.0, np.nan, 55.5, 0.0])
expect_result = np.array(
[[9.40422591e+05, 1.06952091e+06, 1.00000000e+01],
[11.1, 11.2, 11.3],
Expand Down
6 changes: 4 additions & 2 deletions lib/cartopy/tests/test_crs.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -175,9 +175,11 @@ def test_transform_points_xyz(self):
def test_globe(self):
Copy link
Member

Choose a reason for hiding this comment

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

Interesting that this test hasn't changed based on the fact that the globe changed. Shows that the test is of questionable value. Perhaps the first thing to fix is the name of the variables (rugby_globe -> footy_globe and footy_globe -> tennisball_globe). And then we should think about whether this test needs to remain at all...

Thoughts?

Copy link
Member Author

@QuLogic QuLogic Nov 13, 2018 •
edited
Loading

Choose a reason for hiding this comment

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

I mean, that's kind of the point; Mollweide doesn't support elliptical globes, so specifying ellipse or semiminor_axis does nothing.

This test doesn't actually check that the results are different in the end... only that they match some 'true' value, whatever that might be.

Copy link
Member Author

Choose a reason for hiding this comment

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

Oh, perhaps I was a bit unclear; that's the point of the changes in this file. As to the point of the test, I think it's just testing that the globe does something, though maybe it's become a bit redundant by the individual test_eccentric_ellipse tests in the crs directory.

Maybe we should add a test_globe for the ones that only support spheres and remove this test?

Copy link
Member

Choose a reason for hiding this comment

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

Yep. That works.

Copy link
Member

Choose a reason for hiding this comment

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

Maybe we should add a test_globe for the ones that only support spheres and remove this test?

Either that, or we keep the test but fix the variable name. Either works for me. I'll defer this though in order to move the PR forwards.

# Ensure the globe affects output.
rugby_globe = ccrs.Globe(semimajor_axis=9000000,
semiminor_axis=1000000)
semiminor_axis=9000000,
ellipse=None)
footy_globe = ccrs.Globe(semimajor_axis=1000000,
semiminor_axis=1000000)
semiminor_axis=1000000,
ellipse=None)

rugby_moll = ccrs.Mollweide(globe=rugby_globe)
footy_moll = ccrs.Mollweide(globe=footy_globe)
Expand Down