ENH: Azimuthal rotation for Orthographic projection (#2504)

* Added support for azimuthal rotation of Orthographic projection

Author: kenhira

lib/cartopy/crs.py

@@ -2094,9 +2094,14 @@ class Orthographic(Projection):
     _handles_ellipses = False
 
     def __init__(self, central_longitude=0.0, central_latitude=0.0,
-                 globe=None):
+                 azimuth=0.0, globe=None):
         proj4_params = [('proj', 'ortho'), ('lon_0', central_longitude),
-                        ('lat_0', central_latitude)]
+                        ('lat_0', central_latitude), ('alpha', azimuth)]
+        if pyproj.__proj_version__ < '9.5.0' and azimuth != 0.0:
+            warnings.warn(
+                'Setting azimuth is not supported with PROJ versions < 9.5.0. '
+                'Assuming azimuth=0. '
+                'Current PROJ version: %s' % pyproj.__proj_version__)
         super().__init__(proj4_params, globe=globe)
 
         # TODO: Let the globe return the semimajor axis always.

lib/cartopy/tests/crs/test_orthographic.py

@@ -9,6 +9,7 @@
 
 import numpy as np
 from numpy.testing import assert_almost_equal
+import pyproj
 import pytest
 
 import cartopy.crs as ccrs
@@ -17,7 +18,7 @@
 
 def test_default():
     ortho = ccrs.Orthographic()
-    other_args = {'a=6378137.0', 'lon_0=0.0', 'lat_0=0.0'}
+    other_args = {'a=6378137.0', 'lat_0=0.0', 'lon_0=0.0', 'alpha=0.0'}
     check_proj_params('ortho', ortho, other_args)
 
     # WGS84 radius * 0.99999
@@ -30,7 +31,7 @@ def test_default():
 def test_sphere_globe():
     globe = ccrs.Globe(semimajor_axis=1000, ellipse=None)
     ortho = ccrs.Orthographic(globe=globe)
-    other_args = {'a=1000', 'lon_0=0.0', 'lat_0=0.0'}
+    other_args = {'a=1000', 'lat_0=0.0', 'lon_0=0.0', 'alpha=0.0'}
     check_proj_params('ortho', ortho, other_args)
 
     assert_almost_equal(ortho.x_limits, [-999.99, 999.99])
@@ -44,7 +45,7 @@ def test_ellipse_globe():
         ortho = ccrs.Orthographic(globe=globe)
         assert len(w) == 1
 
-    other_args = {'ellps=WGS84', 'lon_0=0.0', 'lat_0=0.0'}
+    other_args = {'ellps=WGS84', 'lat_0=0.0', 'lon_0=0.0', 'alpha=0.0'}
     check_proj_params('ortho', ortho, other_args)
 
     # Limits are the same as default since ellipses are not supported.
@@ -60,7 +61,7 @@ def test_eccentric_globe():
         ortho = ccrs.Orthographic(globe=globe)
         assert len(w) == 1
 
-    other_args = {'a=1000', 'b=500', 'lon_0=0.0', 'lat_0=0.0'}
+    other_args = {'a=1000', 'b=500', 'lat_0=0.0', 'lon_0=0.0', 'alpha=0.0'}
     check_proj_params('ortho', ortho, other_args)
 
     # Limits are the same as spheres since ellipses are not supported.
@@ -73,7 +74,7 @@ def test_eccentric_globe():
 def test_central_params(lon, lat):
     ortho = ccrs.Orthographic(central_latitude=lat, central_longitude=lon)
     other_args = {f'lat_0={lat}', f'lon_0={lon}',
-                  'a=6378137.0'}
+                  'a=6378137.0', 'alpha=0.0'}
     check_proj_params('ortho', ortho, other_args)
 
     # WGS84 radius * 0.99999
@@ -90,7 +91,7 @@ def test_grid():
     ortho = ccrs.Orthographic(globe=globe)
     geodetic = ortho.as_geodetic()
 
-    other_args = {'a=1.0', 'b=1.0', 'lon_0=0.0', 'lat_0=0.0'}
+    other_args = {'a=1.0', 'b=1.0', 'lon_0=0.0', 'lat_0=0.0', 'alpha=0.0'}
     check_proj_params('ortho', ortho, other_args)
 
     assert_almost_equal(np.array(ortho.x_limits),
@@ -143,7 +144,7 @@ def test_sphere_transform():
                               globe=globe)
     geodetic = ortho.as_geodetic()
 
-    other_args = {'a=1.0', 'b=1.0', 'lon_0=-100.0', 'lat_0=40.0'}
+    other_args = {'a=1.0', 'b=1.0', 'lon_0=-100.0', 'lat_0=40.0', 'alpha=0.0'}
     check_proj_params('ortho', ortho, other_args)
 
     assert_almost_equal(np.array(ortho.x_limits),
@@ -156,3 +157,28 @@ def test_sphere_transform():
 
     inverse_result = geodetic.transform_point(result[0], result[1], ortho)
     assert_almost_equal(inverse_result, [-110.0, 30.0])
+
+def test_sphere_rotate():
+    globe = ccrs.Globe(semimajor_axis=1.0, semiminor_axis=1.0,
+                       ellipse=None)
+    ortho = ccrs.Orthographic(central_latitude=40.0, central_longitude=-100.0,
+                              azimuth=180.0, globe=globe)
+    geodetic = ortho.as_geodetic()
+
+    other_args = {'a=1.0', 'b=1.0', 'lon_0=-100.0', 'lat_0=40.0',
+                  'alpha=180.0'}
+    check_proj_params('ortho', ortho, other_args)
+
+    assert_almost_equal(np.array(ortho.x_limits),
+                        [-0.99999, 0.99999])
+    assert_almost_equal(np.array(ortho.y_limits),
+                        [-0.99999, 0.99999])
+
+    result = ortho.transform_point(-110.0, 30.0, geodetic)
+    if pyproj.__proj_version__ >= '9.5.0': # support for alpha (azimuthal rotation)
+        assert_almost_equal(result, np.array([ 0.1503837,  0.1651911]))
+    else:
+        assert_almost_equal(result, np.array([-0.1503837, -0.1651911]))
+
+    inverse_result = geodetic.transform_point(result[0], result[1], ortho)
+    assert_almost_equal(inverse_result, [-110.0, 30.0])