Vector support in bridge.Transform.from_points

src/geovista/bridge.py

@@ -29,6 +29,7 @@
     ZLEVEL_SCALE,
     nan_mask,
     to_cartesian,
+    vectors_to_cartesian,
     wrap,
 )
 from .crs import WGS84, CRSLike, to_wkt
@@ -545,6 +546,9 @@ def from_points(
         zlevel: int | ArrayLike | None = None,
         zscale: float | None = None,
         clean: bool | None = None,
+        vectors: ArrayLike | tuple(ArrayLike) | None = None,
+        vectors_array_name: str | None = None,
+        vectors_z_scaling: float | None = None,
     ) -> pv.PolyData:
         """Build a point-cloud mesh from x-values, y-values and z-levels.
 
@@ -582,6 +586,18 @@ def from_points(
         clean : bool, optional
             Specify whether to merge duplicate points. Defaults to
             :data:`BRIDGE_CLEAN`.
+        vectors : tuple(ArrayLike), optional
+            If present, a tuple of 2 or 3 arrays of the same shape as `xs` and `ys`.
+            These give eastward, northward and (optionally) vertical vectors, which are
+            converted to an [N, 3] array of 3-D vectors attached to the result as a
+            points array ``mesh["vectors"]``.  This can be used to generate glyphs
+            (such as arrows) and streamlines.
+        vectors_array_name : str, optional
+            Specifies an alternate name for the points array to store the vectors.
+            Also set as the active vectors name.   Defaults to "vectors".
+        vectors_z_scaling : float, optional
+            scaling factor to apply to vertical vectors (i.e. relative to the eastward
+            and northward components).  Defaults to 1.0
 
         Returns
         -------
@@ -636,6 +652,42 @@ def from_points(
             mesh.field_data[GV_FIELD_NAME] = np.array([name])
             mesh[name] = data
 
+        if vectors is not None:
+            if vectors_array_name is None:
+                vectors_array_name = "vectors"
+
+            if not isinstance(vectors, tuple) or not all(
+                isinstance(arr, np.ndarray) for arr in vectors
+            ):
+                msg = 'Keyword "vectors" must be a tuple of array-like.'
+                raise ValueError(msg)
+
+            n_vecs = len(vectors)
+            if n_vecs not in (2, 3):
+                msg = (
+                    'Keyword "vectors" must be a tuple of 2 or 3 arrays : '
+                    f"got {n_vecs}."
+                )
+                raise ValueError(msg)
+
+            vectors = [np.asanyarray(vecdata) for vecdata in vectors]
+            xx, yy = vectors[:2]
+            zz = vectors[2] if n_vecs > 2 else np.zeros_like(xx)
+
+            if vectors_z_scaling is not None:
+                zz = zz * vectors_z_scaling
+
+            # TODO @pp-mo: should we pass flattened arrays here, and reshape
+            # as-per the inputs
+            #  (and xyz)?  not clear if multidimensional input is used or needed
+            xx, yy, zz = vectors_to_cartesian(
+                lons_lats=(xs, ys),
+                vectors_uvw=(xx, yy, zz),
+            )
+            vectors = np.vstack((xx, yy, zz)).T
+            mesh[vectors_array_name] = vectors
+            mesh.set_active_vectors(vectors_array_name)
+
         # clean the mesh
         if clean:
             mesh.clean(inplace=True)

src/geovista/common.py

@@ -750,6 +750,51 @@ def to_cartesian(
     return np.vstack(xyz).T if stacked else np.array(xyz)
 
 
+def vectors_to_cartesian(
+    lons_lats: (ArrayLike, ArrayLike),
+    vectors_uvw: (ArrayLike, ArrayLike, ArrayLike),
+) -> (np.ndarray, np.ndarray, np.ndarray):
+    """Convert geographic-oriented vectors to cartesian ``xyz`` points.
+
+    Parameters
+    ----------
+    lons_lats : pair of ArrayLike
+        The longitude + latitude locations of the vectors (in degrees).
+        Both shapes must be the same.
+    vectors_uvw : triple of ArrayLike
+        The eastward, northward and upward vector components.
+        All shapes must be the same as in ``lons_lats``.
+
+    Returns
+    -------
+    (ndarray, ndarray, ndarray)
+        The corresponding ``xyz`` cartesian vector components.
+
+    Notes
+    -----
+    .. versionadded:: 0.?.?
+
+    """
+    # TODO @pp-mo: Argument checking ???
+    lons, lats = (np.deg2rad(arr) for arr in lons_lats)
+    u, v, w = vectors_uvw
+
+    coslons = np.cos(lons)
+    sinlons = np.sin(lons)
+    coslats = np.cos(lats)
+    sinlats = np.sin(lats)
+    # N.B. the term signs are slightly unexpected here, because the viewing coord system
+    # is not quite what you may expect :  The "Y" axis goes to the right, and the "X"
+    # axis points out of the screen, towards the viewer.
+    z_factor = w * coslats - v * sinlats
+    wy = coslons * u + sinlons * z_factor
+    wx = -sinlons * u + coslons * z_factor
+    wz = v * coslats + w * sinlats
+    # NOTE: for better efficiency, we *COULD* handle the w=0 special case separately.
+
+    return wx, wy, wz
+
+
 def to_lonlat(
     xyz: ArrayLike,
     radians: bool | None = False,