From 76d76833768b368a5a6b1c5b00c2c40b9ce897b0 Mon Sep 17 00:00:00 2001 From: Ryan May Date: Wed, 14 Nov 2018 21:49:19 -0700 Subject: [PATCH 1/2] Fix up the satellite projections (Fixes #1144) This refactors the Geostationary and NearSidePerspective projections so that they can individually set up their boundary, limits, etc. Doing so allows fixing the incorrect boundary being set for NearSidePerspective. The two projections can't share the math for the boundary because they are fundamentally different coordinate systems. --- lib/cartopy/crs.py | 90 ++++++++++++++++++++++++++++++++++++++++------ 1 file changed, 79 insertions(+), 11 deletions(-) diff --git a/lib/cartopy/crs.py b/lib/cartopy/crs.py index 90b254c66..a211c1194 100644 --- a/lib/cartopy/crs.py +++ b/lib/cartopy/crs.py @@ -2003,20 +2003,14 @@ def __init__(self, projection, satellite_height=35785831, proj4_params.append(('sweep', sweep_axis)) super(_Satellite, self).__init__(proj4_params, globe=globe) - # TODO: Let the globe return the semimajor axis always. - a = np.float(self.globe.semimajor_axis or WGS84_SEMIMAJOR_AXIS) - b = np.float(self.globe.semiminor_axis or a) - h = np.float(satellite_height) - max_x = h * np.arcsin(a / (a + h)) - max_y = h * np.arcsin(b / (a + h)) - + def _set_bounds(self, max_x, max_y): + false_easting = self.proj4_params['x_0'] + false_northing = self.proj4_params['y_0'] coords = _ellipse_boundary(max_x, max_y, false_easting, false_northing, 61) self._boundary = sgeom.LinearRing(coords.T) - mins = np.min(coords, axis=1) - maxs = np.max(coords, axis=1) - self._x_limits = mins[0], maxs[0] - self._y_limits = mins[1], maxs[1] + self._x_limits = -max_x + false_easting, max_x + false_easting + self._y_limits = -max_y + false_northing, max_y + false_northing self._threshold = np.diff(self._x_limits)[0] * 0.02 @property @@ -2038,12 +2032,38 @@ def y_limits(self): class Geostationary(_Satellite): """ + A view appropriate for satellites in Geostationary Earth orbit. + Perspective view looking directly down from above a point on the equator. + In this projection, the projected coordinates are scanning angles measured + from the satellite looking directly downward, multiplied by the height of + the satellite. + """ def __init__(self, central_longitude=0.0, satellite_height=35785831, false_easting=0, false_northing=0, globe=None, sweep_axis='y'): + """ + Parameters + ---------- + central_longitude: float, optional + The central longitude. Defaults to 0. + satellite_height: float, optional + The height of the satellite. Defaults to 35785831 meters + (true geostationary orbit). + false_easting: + X offset from planar origin in metres. Defaults to 0. + false_northing: + Y offset from planar origin in metres. Defaults to 0. + globe: :class:`cartopy.crs.Globe`, optional + If omitted, a default globe is created. + sweep_axis: 'x' or 'y', optional. Defaults to 'y'. + Controls which axis is scanned first, and thus which angle is + applied first. The default is appropriate for Meteosat, while + 'x' should be used for GOES. + """ + super(Geostationary, self).__init__( projection='geos', satellite_height=satellite_height, @@ -2054,15 +2074,59 @@ def __init__(self, central_longitude=0.0, satellite_height=35785831, globe=globe, sweep_axis=sweep_axis) + # TODO: Let the globe return the semimajor axis always. + a = np.float(self.globe.semimajor_axis or WGS84_SEMIMAJOR_AXIS) + b = np.float(self.globe.semiminor_axis or a) + h = np.float(satellite_height) + max_x = h * np.arcsin(a / (a + h)) + max_y = h * np.arcsin(b / (a + h)) + self._set_bounds(max_x, max_y) + class NearsidePerspective(_Satellite): """ Perspective view looking directly down from above a point on the globe. + In this projection, the projected coordinates are x and y measured from + the origin of a plane tangent to the Earth directly below the perspective + point (e.g. a satellite). + """ def __init__(self, central_longitude=0.0, central_latitude=0.0, satellite_height=35785831, false_easting=0, false_northing=0, globe=None): + """ + Parameters + ---------- + central_longitude: float, optional + The central longitude. Defaults to 0. + central_latitude: float, optional + The central latitude. Defaults to 0. + satellite_height: float, optional + The height of the satellite. Defaults to 35785831 meters + (true geostationary orbit). + false_easting: + X offset from planar origin in metres. Defaults to 0. + false_northing: + 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 "nsper" projection does not handle ' + 'elliptical globes.') + super(NearsidePerspective, self).__init__( projection='nsper', satellite_height=satellite_height, @@ -2072,6 +2136,10 @@ def __init__(self, central_longitude=0.0, central_latitude=0.0, false_northing=false_northing, globe=globe) + h = np.float(satellite_height) + max_x = a * np.sqrt(h / (2 * a + h)) + self._set_bounds(max_x, max_x) + class AlbersEqualArea(Projection): """ From 0c157fdf6627c7870e839c87dc98805700819388 Mon Sep 17 00:00:00 2001 From: Ryan May Date: Wed, 14 Nov 2018 22:11:13 -0700 Subject: [PATCH 2/2] Fix up NearsideProjection tests Since these are two different projections, they should not share as much, so I've decoupled them a bit and updated some values. --- lib/cartopy/tests/crs/test_geostationary.py | 5 +- .../tests/crs/test_nearside_perspective.py | 56 ++++++++++++------- 2 files changed, 38 insertions(+), 23 deletions(-) diff --git a/lib/cartopy/tests/crs/test_geostationary.py b/lib/cartopy/tests/crs/test_geostationary.py index c0ad34501..f61c38ffb 100644 --- a/lib/cartopy/tests/crs/test_geostationary.py +++ b/lib/cartopy/tests/crs/test_geostationary.py @@ -26,14 +26,13 @@ import cartopy.crs as ccrs -# Note: code here is now shared with the NearsidePerspective test. def check_proj4_params(name, crs, other_args): expected = other_args | {'proj={}'.format(name), 'units=m', 'no_defs'} pro4_params = set(crs.proj4_init.lstrip('+').split(' +')) assert expected == pro4_params -class GeostationaryTestsMixin(object): +class TestGeostationary(object): test_class = ccrs.Geostationary expected_proj_name = 'geos' @@ -84,8 +83,6 @@ def test_eastings(self): 10434177.81588539, 5309177.81588539), decimal=4) - -class TestGeostationary(GeostationaryTestsMixin, object): def test_sweep(self): geos = ccrs.Geostationary(sweep_axis='x') other_args = {'ellps=WGS84', 'h=35785831', 'lat_0=0.0', 'lon_0=0.0', diff --git a/lib/cartopy/tests/crs/test_nearside_perspective.py b/lib/cartopy/tests/crs/test_nearside_perspective.py index 843238f5f..37b293f5d 100644 --- a/lib/cartopy/tests/crs/test_nearside_perspective.py +++ b/lib/cartopy/tests/crs/test_nearside_perspective.py @@ -23,28 +23,46 @@ from numpy.testing import assert_almost_equal -from cartopy.tests.crs.test_geostationary import (GeostationaryTestsMixin, - check_proj4_params) +from cartopy.tests.crs.test_geostationary import check_proj4_params -from cartopy.crs import NearsidePerspective +import cartopy.crs as ccrs -class TestEquatorialDefault(GeostationaryTestsMixin, object): - # Check that it behaves just like Geostationary, in the absence of a - # central_latitude parameter. - test_class = NearsidePerspective - expected_proj_name = 'nsper' +def test_default(): + geos = ccrs.NearsidePerspective() + other_args = {'a=6378137.0', 'h=35785831', 'lat_0=0.0', 'lon_0=0.0', + 'x_0=0', 'y_0=0'} + check_proj4_params('nsper', geos, other_args) -class TestOwnSpecifics(object): - def test_central_latitude(self): - # Check the effect of the added 'central_latitude' key. - geos = NearsidePerspective(central_latitude=53.7) - other_args = {'ellps=WGS84', 'h=35785831', 'lat_0=53.7', 'lon_0=0.0', - 'x_0=0', 'y_0=0'} - check_proj4_params('nsper', geos, other_args) + assert_almost_equal(geos.boundary.bounds, + (-5476336.098, -5476336.098, + 5476336.098, 5476336.098), + decimal=3) - assert_almost_equal(geos.boundary.bounds, - (-5434177.81588539, -5434177.81588539, - 5434177.81588539, 5434177.81588539), - decimal=4) + +def test_offset(): + geos = ccrs.NearsidePerspective(false_easting=5000000, + false_northing=-123000,) + other_args = {'a=6378137.0', 'h=35785831', 'lat_0=0.0', 'lon_0=0.0', + 'x_0=5000000', 'y_0=-123000'} + + check_proj4_params('nsper', geos, other_args) + + assert_almost_equal(geos.boundary.bounds, + (-476336.098, -5599336.098, + 10476336.098, 5353336.098), + decimal=4) + + +def test_central_latitude(): + # Check the effect of the added 'central_latitude' key. + geos = ccrs.NearsidePerspective(central_latitude=53.7) + other_args = {'a=6378137.0', 'h=35785831', 'lat_0=53.7', 'lon_0=0.0', + 'x_0=0', 'y_0=0'} + check_proj4_params('nsper', geos, other_args) + + assert_almost_equal(geos.boundary.bounds, + (-5476336.098, -5476336.098, + 5476336.098, 5476336.098), + decimal=3)