FIX: Change PlateCarree Vector Handling

lib/cartopy/crs.py

@@ -479,14 +479,34 @@ def transform_vectors(self, src_proj, x, y, u, v):
             raise ValueError('x, y, u and v arrays must be the same shape')
         if x.ndim not in (1, 2):
             raise ValueError('x, y, u and v must be 1 or 2 dimensional')
-        # Transform the coordinates to the target projection.
-        proj_xyz = self.transform_points(src_proj, x, y)
-        target_x, target_y = proj_xyz[..., 0], proj_xyz[..., 1]
-        # Rotate the input vectors to the projection.
-        #
-        # 1: Find the magnitude and direction of the input vectors.
-        vector_magnitudes = (u**2 + v**2)**0.5
+
+        # 1a: Find the magnitude of the input vectors
+        vector_magnitudes = np.sqrt(u**2 + v**2)
+
+        # Handle special case of PlateCarree projection
+        # This assumes we were given lon/lat coordinates and
+        # east/north vector components (NOT PlateCarree vector components)
+        # When projecting the components we need to scale the vector
+        # component for the transformation to maintain the proper angles
+        if isinstance(src_proj, PlateCarree):
+            # First handle points close to the pole by cutting off the
+            # latitude a small distance away
+            pole_cutoff = 90 - 1e-2
+            pole_points = np.abs(y) >= pole_cutoff
+            if np.any(pole_points):
+                # Warn for points close to the poles, as these are
+                # unreliable for vector directions under transformation
+                warnings.warn('Vector transforms near the pole '
+                              'may not have been transformed correctly')
+                # Move the pole points away from the pole ever so slightly
+                y[pole_points] = np.sign(y[pole_points]) * pole_cutoff
+
+            scale_factor = np.cos(y / 180 * np.pi)
+            v = v * scale_factor
+
+        # 1b: Find the direction of the input vectors
         vector_angles = np.arctan2(v, u)
+
         # 2: Find a point in the direction of the original vector that is
         #    a small distance away from the base point of the vector (near
         #    the poles the point may have to be in the opposite direction
@@ -495,6 +515,7 @@ def transform_vectors(self, src_proj, x, y, u, v):
         delta = (src_proj.x_limits[1] - src_proj.x_limits[0]) / factor
         x_perturbations = delta * np.cos(vector_angles)
         y_perturbations = delta * np.sin(vector_angles)
+
         # 3: Handle points that are invalid. These come from picking a new
         #    point that is outside the domain of the CRS. The first step is
         #    to apply the native transform to the input coordinates to make
@@ -535,21 +556,45 @@ def transform_vectors(self, src_proj, x, y, u, v):
             source_x + x_perturbations < src_proj.x_limits[0] - eps,
             source_x + x_perturbations > src_proj.x_limits[1] + eps)
         if problem_points.any():
-            warnings.warn('Some vectors at source domain corners '
+            warnings.warn('Some vectors at the source domain corners '
                           'may not have been transformed correctly')
+
         # 4: Transform this set of points to the projection coordinates and
         #    find the angle between the base point and the perturbed point
         #    in the projection coordinates (reversing the direction at any
         #    points where the original was reversed in step 3).
+        proj_xyz = self.transform_points(src_proj, x, y)
+        target_x, target_y = proj_xyz[..., 0], proj_xyz[..., 1]
+
         proj_xyz = self.transform_points(src_proj,
                                          source_x + x_perturbations,
                                          source_y + y_perturbations)
         target_x_perturbed = proj_xyz[..., 0]
         target_y_perturbed = proj_xyz[..., 1]
+
+        x_delta = target_x_perturbed - target_x
+        if isinstance(self, PlateCarree):
+            # Warn for points close to the poles, as these are
+            # unreliable for vector directions under transformation
+            pole_cutoff = 90 - 1e-2
+            pole_points = np.abs(target_y) >= pole_cutoff
+            if np.any(pole_points):
+                warnings.warn('Vector transforms near the pole '
+                              'may not have been transformed correctly')
+                # Move the pole points away from the pole ever so slightly
+                target_y[pole_points] = (np.sign(target_y[pole_points]) *
+                                         pole_cutoff)
+            # Account for scaling to the North/South components in the
+            # opposite manner as the source coordinate handling
+            # (i.e. inverse scale factor), with target_y being
+            # the latitudes we are going to
+            scale_factor = np.cos(target_y / 180 * np.pi)
+            x_delta = x_delta * scale_factor
         projected_angles = np.arctan2(target_y_perturbed - target_y,
-                                      target_x_perturbed - target_x)
+                                      x_delta)
         if reversed_vectors.any():
             projected_angles[reversed_vectors] += np.pi
+
         # 5: Form the projected vector components, preserving the magnitude
         #    of the original vectors.
         projected_u = vector_magnitudes * np.cos(projected_angles)

lib/cartopy/tests/test_vector_transform.py

@@ -6,6 +6,7 @@
 
 import numpy as np
 from numpy.testing import assert_array_almost_equal, assert_array_equal
+import pytest
 
 import cartopy.crs as ccrs
 import cartopy.vector_transform as vec_trans
@@ -144,11 +145,11 @@ def test_with_transform(self):
                                     [7.5, 7.5, 7.5, 7.5, 7.5],
                                     [10., 10., 10., 10., 10]])
         expected_u_grid = np.array([[np.nan, np.nan, np.nan, np.nan, np.nan],
-                                    [np.nan, 2.3838, 3.5025, 2.6152, np.nan],
-                                    [2, 3.0043, 4, 2.9022, 2]])
+                                    [np.nan, 2.3893, 3.5097, 2.6194, np.nan],
+                                    [2, 3.0005, 4, 2.8977, 2]])
         expected_v_grid = np.array([[np.nan, np.nan, np.nan, np.nan, np.nan],
-                                    [np.nan, 2.6527, 2.1904, 2.4192, np.nan],
-                                    [5.5, 4.6483, 4, 4.47, 5.5]])
+                                    [np.nan, 2.6486, 2.1878, 2.4138, np.nan],
+                                    [5.5, 4.6497, 4, 4.4702, 5.5]])
 
         x_grid, y_grid, u_grid, v_grid = vec_trans.vector_scalar_to_grid(
             src_crs, target_crs, (5, 3), x_nps, y_nps, u_nps, v_nps)
@@ -222,3 +223,92 @@ def test_with_scalar_field_non_ndarray_data(self):
         assert_array_almost_equal(u_grid, expected_u_grid)
         assert_array_almost_equal(v_grid, expected_v_grid)
         assert_array_almost_equal(s_grid, expected_s_grid)
+
+
+class TestTransformVectors:
+
+    def test_transform(self):
+        # Test some simple vectors to make sure they are transformed
+        # correctly.
+        rlons = np.array([-90., 0, 90., 180.])
+        rlats = np.array([0., 0., 0., 0.])
+        src_proj = ccrs.PlateCarree()
+        target_proj = ccrs.Stereographic(central_latitude=90,
+                                         central_longitude=0)
+        # transform grid eastward vectors
+        ut, vt = target_proj.transform_vectors(src_proj,
+                                               rlons,
+                                               rlats,
+                                               np.ones([4]),
+                                               np.zeros([4]))
+        assert_array_almost_equal(ut, np.array([0, 1, 0, -1]), decimal=2)
+        assert_array_almost_equal(vt, np.array([-1, 0, 1, 0]), decimal=2)
+        # transform grid northward vectors
+        ut, vt = target_proj.transform_vectors(src_proj,
+                                               rlons,
+                                               rlats,
+                                               np.zeros([4]),
+                                               np.ones([4]))
+        assert_array_almost_equal(ut, np.array([1, 0, -1, 0]), decimal=2)
+        assert_array_almost_equal(vt, np.array([0, 1, 0, -1]), decimal=2)
+        # transform grid north-eastward vectors
+        ut, vt = target_proj.transform_vectors(src_proj,
+                                               rlons,
+                                               rlats,
+                                               np.ones([4]),
+                                               np.ones([4]))
+        assert_array_almost_equal(ut, np.array([1, 1, -1, -1]), decimal=2)
+        assert_array_almost_equal(vt, np.array([-1, 1, 1, -1]), decimal=2)
+
+    def test_transform_and_inverse(self):
+        # Check a full circle transform back to the native projection.
+        x = np.arange(-60, 42.5, 2.5)
+        y = np.arange(30, 72.5, 2.5)
+        x2d, y2d = np.meshgrid(x, y)
+        u = np.cos(np.deg2rad(y2d))
+        v = np.cos(2. * np.deg2rad(x2d))
+        src_proj = ccrs.PlateCarree()
+        target_proj = ccrs.Stereographic(central_latitude=90,
+                                         central_longitude=0)
+        proj_xyz = target_proj.transform_points(src_proj, x2d, y2d)
+        xt, yt = proj_xyz[..., 0], proj_xyz[..., 1]
+        ut, vt = target_proj.transform_vectors(src_proj, x2d, y2d, u, v)
+        utt, vtt = src_proj.transform_vectors(target_proj, xt, yt, ut, vt)
+        assert_array_almost_equal(u, utt, decimal=4)
+        assert_array_almost_equal(v, vtt, decimal=4)
+
+    def test_invalid_input_domain(self):
+        # If an input coordinate is outside the input projection domain
+        # we should be able to handle it correctly.
+        rlon = np.array([270.])
+        rlat = np.array([0.])
+        u = np.array([1.])
+        v = np.array([0.])
+        src_proj = ccrs.PlateCarree()
+        target_proj = ccrs.Stereographic(central_latitude=90,
+                                         central_longitude=0)
+        ut, vt = target_proj.transform_vectors(src_proj, rlon, rlat, u, v)
+        assert_array_almost_equal(ut, np.array([0]), decimal=2)
+        assert_array_almost_equal(vt, np.array([-1]), decimal=2)
+
+    @pytest.mark.parametrize('x, y, u, v, expected', [
+        pytest.param(180, 0, 1, 0, [-1, 0], id='x non-corner'),
+        pytest.param(0, 90, 0, 1, [0, 1], id='y non-corner'),
+        pytest.param(180, 90, 1, 1, [-1, -1], id='xy corner'),
+        pytest.param(180, 90, 1, -1, [-1, 1], id='x corner'),
+        pytest.param(180, 90, -1, 1, [1, -1], id='y corner')])
+    def test_transform_on_border(self, x, y, u, v, expected):
+        src_proj = ccrs.PlateCarree()
+        target_proj = ccrs.Stereographic(central_latitude=90,
+                                         central_longitude=0)
+
+        x, y = np.array([x]), np.array([y])
+        u, v = np.array([u]), np.array([v])
+        # We expect to warn the user only when y is on the pole
+        if y[0] == 90:
+            with pytest.warns(UserWarning):
+                ut, vt = target_proj.transform_vectors(src_proj, x, y, u, v)
+        else:
+            ut, vt = target_proj.transform_vectors(src_proj, x, y, u, v)
+
+        assert_array_almost_equal([ut[0], vt[0]], expected, decimal=2)