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
@@ -546,6 +547,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.
 
@@ -583,6 +587,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
         -------
@@ -637,6 +653,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,

src/geovista/examples/vector_data/wind_arrows.py

@@ -0,0 +1,139 @@
+#!/usr/bin/env python3
+# Copyright (c) 2021, GeoVista Contributors.
+#
+# This file is part of GeoVista and is distributed under the 3-Clause BSD license.
+# See the LICENSE file in the package root directory for licensing details.
+
+"""
+3D Wind Arrows
+--------------
+
+This example demonstrates how to display wind vectors.
+
+📋 Summary
+^^^^^^^^^^
+
+The data source provides X and Y arrays containing plain longitude and
+latitude values, which is the most common case.  
+
+The wind information is provided in three separate field arrays, 'U, V and W'
+ -- i.e. eastward, northward and vertical components.  
+These values are coded for each location (X, Y), measured relative to the longitude,
+latitude and vertical directions at each point.
+
+There is no connectivity provided, so each location has a vector and is independent of
+the others.  Hence we use the `geovista.Transform.from_points` function, passing the
+winds to the `vectors` keyword. 
+
+Initially, we can just show the horizontal winds, as this easier to interprt
+We scale up the 'W' values, since vertical winds are typically much smaller than
+horizontal.  Coastlines and a base layer are also added for ease of viewing.
+
+"""  # noqa: D205,D212,D400
+import geovista as gv
+from geovista.pantry.data import lfric_winds
+
+# get sample data
+sample = lfric_winds()
+
+# Create a mesh of individual points, adding vectors at each point
+mesh = gv.Transform.from_points(
+    sample.lons,
+    sample.lats,
+    vectors = (sample.u, sample.v),
+)
+
+# Create a new mesh containing arrow glyphs, from the mesh vectors
+# NOTE: the 'mesh vectors' are a specific concept in PyVista
+# NOTE ALSO: the 'arrows' property is effectively a convenience for calling
+#  :meth:'~pyvista.Dataset.glyph'
+arrows = mesh.glyph(factor=0.02)  # Note the overall scaling factor
+
+# Add the arrows to a plotter with other aspects, and display
+plotter = gv.GeoPlotter()
+plotter.add_base_layer(radius=0.99)
+plotter.add_mesh(arrows, cmap="inferno")
+plotter.add_coastlines()
+plotter.add_graticule()
+plotter.add_axes()
+
+# Set up a nice default view
+plotter.camera.zoom(1.3)  # adjusts the camera view angle
+selected_view = [
+    (-4.0688208659033505, -2.5462610064466777, -2.859304866708606),
+    (-0.0037798285484313965, 0.005168497562408447, -0.0031679868698120117),
+    (-0.523382090763761, -0.11174892277533728, 0.8447386372874786)
+]
+plotter.camera_position = selected_view
+plotter.show()
+print(plotter.camera_position)
+
+
+# %%
+# Repeat, but now add in the 'W' vertical components.
+# These need scaling up, since vertical winds are typically much smaller than
+# horizontal.
+# We also use one colour, and apply a vertical offset to prevent downward-going arrows
+# from disappearing into the surface.
+
+# Create a mesh of individual points, adding vectors at each point
+mesh = gv.Transform.from_points(
+    sample.lons,
+    sample.lats,
+    # supply all three components
+    vectors = (sample.u, sample.v, sample.w),
+    # apply additional scaling and a vertical offset
+    vectors_z_scaling=1500.,
+    radius=1.1
+)
+arrows = mesh.glyph(factor=0.02)
+
+plotter = gv.GeoPlotter()
+plotter.add_base_layer()
+plotter.add_mesh(arrows, color="darkred")
+plotter.add_coastlines()
+plotter.add_graticule()
+plotter.add_axes()
+
+plotter.camera.zoom(1.3)
+selected_view = [
+    (0.9892890077409511, -2.9925812011503097, 1.008438916341214),
+    (0.456372154072792, 0.10044567821980169, 0.7120015972700701),
+    (-0.39009517660643345, 0.021012607195809205, 0.9205347486799345)
+]
+plotter.camera_position = selected_view
+plotter.show()
+print(plotter.camera_position)
+
+
+# %%
+# Finally, it sometimes makes more sense to display all arrows the same size so that
+# direction is always readable.
+# Here's an example of constant size, but still colored by windspeed.
+
+# Create a mesh of individual points, adding vectors at each point
+mesh = gv.Transform.from_points(
+    sample.lons,
+    sample.lats,
+    # supply all three components
+    vectors = (sample.u, sample.v),
+)
+# Note: the overall size scale is now different, too
+arrows = mesh.glyph(factor=0.1, scale=False)
+
+plotter = gv.GeoPlotter()
+plotter.add_base_layer()
+plotter.add_mesh(arrows, cmap="magma")
+plotter.add_coastlines()
+plotter.add_graticule()
+plotter.add_axes()
+
+plotter.camera.zoom(1.3)  # adjusts the camera view angle
+selected_view = [
+    (-4.0688208659033505, -2.5462610064466777, -2.859304866708606),
+    (-0.0037798285484313965, 0.005168497562408447, -0.0031679868698120117),
+    (-0.523382090763761, -0.11174892277533728, 0.8447386372874786)
+]
+plotter.camera_position = selected_view
+plotter.show()
+print(plotter.camera_position)

src/geovista/pantry/data.py

@@ -15,6 +15,7 @@
 from dataclasses import dataclass, field
 from enum import Enum
 from functools import lru_cache
+from pathlib import Path
 from typing import TYPE_CHECKING
 
 import lazy_loader as lazy
@@ -124,6 +125,21 @@ class SampleUnstructuredXY:
     ndim: int = 2
 
 
+@dataclass(frozen=True)
+class SampleVectorsXYUVW:
+    """Data container for vector information on unconnected points."""
+
+    lons: ArrayLike
+    lats: ArrayLike
+    u: ArrayLike
+    v: ArrayLike
+    w: ArrayLike = field(default=None)
+    name: str = field(default=None)
+    units: str = field(default=None)
+    steps: int = field(default=None)
+    ndim: int = 1
+
+
 def capitalise(title: str) -> str:
     """Capitalise each word and replacing inappropriate characters.
 
@@ -963,3 +979,50 @@ def ww3_global_tri() -> SampleUnstructuredXY:
         name=name,
         units=units,
     )
+
+
+
+@lru_cache(maxsize=LRU_CACHE_SIZE)
+def lfric_winds() -> SampleStructuredXYZ:
+    """Download and cache cloud-point sample data.
+
+    Load Met Office Unified Model (UM) ORCA2 curvilinear mesh with gradient filter.
+
+    Returns
+    -------
+    SampleStructuredXYZ
+        The gradient filtered spatial coordinates and data payload.
+
+    Notes
+    -----
+    .. versionadded:: 0.2.0
+
+    """
+    # TODO:
+    #  #1 access with pooch and cache
+    #  #2 include the data in geovista-data
+    # fname = "lfric-winds.nc"
+    # processor = pooch.Decompress(method="auto", name=fname)
+    # resource = CACHE.fetch(f"{PANTRY_DATA}/{fname}.bz2", processor=processor)
+    # dataset = nc.Dataset(resource)
+    temporary_path = (
+        Path(__file__).parent
+        / "temporary_mini_test_data"
+        / "lfric_winds_sample.nc"
+    )
+    dataset = nc.Dataset(temporary_path)
+
+    # load the lon/lat/zlevel points
+    lons = dataset.variables["Mesh2d_face_x"][:]
+    lats = dataset.variables["Mesh2d_face_y"][:]
+    # NOTE: for now, take first time and depth.
+    # This effectively 'pretends' U and V are on same levels as W -- which they aren't!
+    units = dataset.variables["u_in_w3"].units
+    u = dataset.variables["u_in_w3"][0, 0]
+    v = dataset.variables["v_in_w3"][0, 0]
+    w = dataset.variables["w_in_wth"][0, 0]
+    name = "Wind"
+
+    return SampleVectorsXYUVW(lons, lats, u, v, w, name=name, units=units)
+
+