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): """ 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)