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/cache/__init__.py

@@ -33,7 +33,7 @@
 BASE_URL: str = "https://github.com/bjlittle/geovista-data/raw/{version}/assets/"
 """Base URL for :mod:`geovista` resources."""
 
-DATA_VERSION: str = "2024.01.2"
+DATA_VERSION: str = "2024.01.5"
 """The ``geovista-data`` repository version for :mod:`geovista` resources."""
 
 GEOVISTA_CACHEDIR: str = "GEOVISTA_CACHEDIR"

src/geovista/cache/registry.txt

@@ -17,6 +17,7 @@ pantry/data/lams/london.nc.bz2 33bdbe43308fb85ff2d4ae25ba321e30e38a1dd85ab7480a7
 pantry/data/lams/new_zealand.nc.bz2 e0eb824d172de988f684999ef0f5d5725fb744be64bf1ea9e0e694965d8455de
 pantry/data/lams/polar.nc.bz2 1be492f856e7a6bf1c00806a3f0a04b52e95dd6e5a874bd510d95717e5fafcff
 pantry/data/lams/uk.nc.bz2 53b6cee7342e8d80e838df63ed8949d14f2f6924a2229c76a0c432eda5b99658
+pantry/data/lfric_winds_sample.nc.bz2 ec56f98b4c0075cc858f4e1858e8013927ae2cf6a5dcefbed965e213de10c478
 pantry/data/qrclim.sst.ugrid.nc.bz2 1f31a867fe944cb56aa321809fbe12b1cf83ae4cddea3e162c4f93870d7d4990
 pantry/data/qrparam_shared.orog.ugrid.nc.bz2 7cdf3b7122d09751bb500c8cd4b6cf9053d76fff5fd8bb9a340cbd7b4557d30f
 pantry/data/oisst-avhrr.nc.bz2 32dd80d78783a872d70da581a866c86d9186496713d520394c3c4fd5ec0086ed

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.5.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,133 @@
+#!/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 vector data.
+
+📋 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 interpret.
+"""  # noqa: D205,D212,D400
+from __future__ import annotations
+
+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.
+# NOTE: this creates a mesh with 'mesh vectors' : a specific concept in PyVista.
+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: use an overall scaling factor to make the arrows a reasonable size.
+arrows = mesh.glyph(factor=0.02)
+
+# 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()
+
+
+# %%
+# Repeat, but now add in the 'W' vertical components.
+# To be visible, these need scaling up, since vertical wind values are typically much
+# smaller than horizontal.
+# We also apply a vertical offset ("radius"), 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 to W values
+    vectors_z_scaling=1500.0,
+    # offset from surface so avoid downward-pointing arrows disappearing
+    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.6917810912064826, -3.065688850990997, 0.4317999141924935),
+    (0.41358279170396495, 0.07362917740509836, 0.5091223320854129),
+    (0.8088496364623022, 0.05726400555597287, 0.5852205560833343),
+]
+plotter.camera_position = selected_view
+plotter.show()
+
+
+# %%
+# Finally, it sometimes makes more sense to display all arrows the same size so that
+# direction is always readable.
+# Here's an example with constant size arrows, but still colored by windspeed.
+
+mesh = gv.Transform.from_points(
+    sample.lons,
+    sample.lats,
+    vectors=(sample.u, sample.v),
+)
+# Note: with no scaling, the basic arrows size is now rather different
+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)
+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()

src/geovista/pantry/data.py

@@ -124,6 +124,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 +978,38 @@ def ww3_global_tri() -> SampleUnstructuredXY:
         name=name,
         units=units,
     )
+
+
+@lru_cache(maxsize=LRU_CACHE_SIZE)
+def lfric_winds() -> SampleStructuredXYZ:
+    """Download and cache unstructured 3D winds sample.
+
+    This data is derived from the LFRic test suite.
+
+    Returns
+    -------
+    SampleStructuredXYUV
+        data payload with XY spatial coordinates and UVW wind components.
+
+    Notes
+    -----
+    .. versionadded:: 0.5.0
+
+    """
+    fname = "lfric_winds_sample.nc"
+    processor = pooch.Decompress(method="auto", name=fname)
+    resource = CACHE.fetch(f"{PANTRY_DATA}/{fname}.bz2", processor=processor)
+    dataset = nc.Dataset(resource)
+
+    # 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)