Vector support in bridge.Transform.from_points

changelog/676.feature.rst

@@ -0,0 +1,2 @@
+Added support for 3D vectors in :meth:`geovista.bridge.Transform.from_points`
+ (:user:`pp-mo`)

src/geovista/bridge.py

@@ -19,6 +19,7 @@
 
 from __future__ import annotations
 
+from collections.abc import Iterable
 from pathlib import Path
 from typing import TYPE_CHECKING, TypeAlias
 import warnings
@@ -34,6 +35,7 @@
     cast_UnstructuredGrid_to_PolyData,
     nan_mask,
     to_cartesian,
+    vectors_to_cartesian,
     wrap,
 )
 from .crs import WGS84, CRSLike, to_wkt
@@ -616,6 +618,11 @@ def from_points(
         zlevel: int | ArrayLike | None = None,
         zscale: float | None = None,
         clean: bool | None = None,
+        vectors: (
+            tuple[ArrayLike, ArrayLike, ArrayLike] | tuple[ArrayLike, ArrayLike] | None
+        ) = None,
+        vectors_crs: CRSLike | None = None,
+        vectors_array_name: str | None = None,
     ) -> pv.PolyData:
         """Build a point-cloud mesh from x-values, y-values and z-levels.
 
@@ -655,6 +662,21 @@ def from_points(
             Specify whether to merge duplicate points. See
             :meth:`pyvista.PolyDataFilters.clean`. 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_crs : CRSLike, optional
+            The Coordinate Reference System of the provided `vectors`. May be anything
+            accepted by :meth:`pyproj.crs.CRS.from_user_input`. Defaults to the same
+            as 'crs'.  Note that `vectors_crs` only specifies **horizontal orientation**
+            of the vectors : their magnitudes are always spatial distance, and the Z
+            component is always radial (i.e. "upwards").
+        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"``.
 
         Returns
         -------
@@ -670,22 +692,30 @@ def from_points(
         zscale = ZLEVEL_SCALE if zscale is None else float(zscale)
 
         if crs is not None:
+            # Handle specified input CRS
             crs = pyproj.CRS.from_user_input(crs)
+            if crs == WGS84:
+                crs = None
 
-            if crs != WGS84:
-                transformed = transform_points(src_crs=crs, tgt_crs=WGS84, xs=xs, ys=ys)
-                xs, ys = transformed[:, 0], transformed[:, 1]
+        if crs is None:
+            xs_latlon, ys_latlon = xs, ys
+        else:
+            # Input CRS is not plain lat-lon : convert to that, which is our reference
+            transformed = transform_points(src_crs=crs, tgt_crs=WGS84, xs=xs, ys=ys)
+            xs_latlon, ys_latlon = transformed[:, 0], transformed[:, 1]
 
         # ensure longitudes (degrees) are in half-closed interval [-180, 180)
-        xs = wrap(xs)
+        xs_latlon = wrap(xs_latlon)
 
         # reduce any singularity points at the poles to a common longitude
-        poles = np.isclose(np.abs(ys), 90)
+        poles = np.isclose(np.abs(ys_latlon), 90)
         if np.any(poles):
-            xs[poles] = 0
+            xs_latlon[poles] = 0
 
         # convert lat/lon to cartesian xyz
-        xyz = to_cartesian(xs, ys, radius=radius, zlevel=zlevel, zscale=zscale)
+        xyz = to_cartesian(
+            xs_latlon, ys_latlon, radius=radius, zlevel=zlevel, zscale=zscale
+        )
 
         # create the point-cloud mesh
         mesh = pv.PolyData(xyz)
@@ -707,6 +737,111 @@ def from_points(
             mesh.field_data[GV_FIELD_NAME] = np.array([name])
             mesh[name] = data
 
+        if vectors is not None:
+            # TODO @pp-mo: this section is very long, consider refactor
+            if vectors_array_name is None:
+                vectors_array_name = "vectors"
+
+            if (
+                not isinstance(vectors, Iterable)  # type: ignore[redundant-expr]
+                or len(vectors) not in (2, 3)
+            ):
+                msg = "Keyword 'vectors' must be a tuple of 2 or 3 array-likes."
+                raise ValueError(msg)
+
+            in_vectors = [np.asanyarray(arr) for arr in vectors]
+
+            xx, yy = in_vectors[:2]
+            zz = in_vectors[2] if len(in_vectors) == 3 else np.zeros_like(xx)
+
+            if vectors_crs is None:
+                # Vectors crs defaults to main input CRS
+                # NOTE: vectors may have a different CRS than the input points.
+                # The only likely usage is for true-lat-lon winds with locations on a
+                # rotated grid or similar, but we probably do need to support that.
+                vectors_crs = crs
+
+            if vectors_crs is not None:
+                vectors_crs = pyproj.CRS.from_user_input(vectors_crs)
+                if vectors_crs == WGS84:
+                    # As for main crs, "None" means plain-latlon crs
+                    vectors_crs = None
+
+            if vectors_crs is None:
+                vector_xs = xs_latlon
+                vector_ys = ys_latlon
+                post_rotate_matrix = None
+            else:
+                # Supplied vector xx/yy/zz are for a "different" orientation than
+                # standard lat-lon, so will need rotating afterwards.
+                # NOTE: we can only do this for specific CRS types where we know how to
+                # find **its** pole.
+                axis_info = getattr(vectors_crs, "axis_info", [])
+                ok = (
+                    isinstance(axis_info, list)
+                    and len(axis_info) >= 2
+                    and all(hasattr(x, "name") for x in axis_info)
+                )
+                if ok:
+                    axis_names = [x.name for x in axis_info]
+                    ok = axis_names[:2] == ["Longitude", "Latitude"]
+                if not ok:
+                    msg = (
+                        "Cannot determine wind directions : Target CRS type is not "
+                        f"supported for grid orientation decoding : {vectors_crs}."
+                    )
+                    raise ValueError(msg)
+
+                # For a CRS with longitude and latitude axes, its axes determine the
+                # east- and north-wards directions of surface-oriented vectors.
+                # This means we can calculate the xyz vectors for the input points in
+                # vector coordinates, and afterwards rotate them to the display.
+                # TODO @pp-mo: support other common cases as needed, for example OSGB
+                #  or Mercator.  Sadly I think there is no general solution.
+
+                # calculate post-rotate matrix, so "vectors @= post_rotate_matrix"
+                # transform the equivalent of the x,y,z basis vectors
+                bases_x = np.array([0, 90, 0])
+                bases_y = np.array([0, 0, 90])
+                transformed_bases = transform_points(
+                    src_crs=vectors_crs, tgt_crs=WGS84, xs=bases_x, ys=bases_y
+                )
+                cartesian_bases = to_cartesian(
+                    transformed_bases[:, 0], transformed_bases[:, 1]
+                )
+                post_rotate_matrix = np.array(cartesian_bases).T
+
+                # We will also need points as lons+lats **in the vector CRS** to
+                # calculate the vectors (i.e. oriented to "its" northward + eastward).
+                vector_xs = xs
+                vector_ys = ys
+                src_crs = crs or WGS84
+                if vectors_crs != src_crs:
+                    transformed = transform_points(
+                        src_crs=src_crs or WGS84, tgt_crs=vectors_crs, xs=xs, ys=ys
+                    )
+                    vector_xs, vector_ys = transformed[:, 0], transformed[:, 1]
+
+            # 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=vector_xs,
+                lats=vector_ys,
+                vectors_uvw=(xx, yy, zz),
+            )
+            mesh_vectors = np.vstack((xx, yy, zz)).T
+
+            if post_rotate_matrix is not None:
+                # At this point, vectors are correct xyz's for the original points in
+                # 'vector_crs', but not ready for plotting since the "true" point
+                # locations are different : hence apply "post-rotation".
+                # TODO @pp-mo: replace np.dot with a '@' if a masked-array multiply bug
+                #  gets fixed : see https://github.com/numpy/numpy/issues/14992
+                mesh_vectors = np.dot(post_rotate_matrix, mesh_vectors.T).T
+
+            mesh[vectors_array_name] = mesh_vectors
+            mesh.set_active_vectors(vectors_array_name)
+
         # clean the mesh
         if clean:
             mesh.clean(inplace=True)

src/geovista/cache/__init__.py

@@ -44,7 +44,7 @@
 BASE_URL: str = "https://github.com/bjlittle/geovista-data/raw/{version}/assets/"
 """Base URL for :mod:`geovista` resources."""
 
-DATA_VERSION: str = "2025.01.5"
+DATA_VERSION: str = "2025.01.6"
 """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
@@ -74,6 +75,7 @@ tests/images/examples.test__unstructured.tri_hammer.png 55a063f42722665ae3575d70
 tests/images/examples.test__unstructured.tri.png 17706b23041a027bd52be9950a86703a59b8f1282a82790ee7686832d0bc4196
 tests/images/examples.test__vector_data.flow_directions.png 26cea271234f57b721da235b0d1b9b6d35927003621ac3f4750cd8b52e99134e
 tests/images/examples.test__vector_data.vectors.png 72f1a30dd57b8c6cc635e866aa1ffd79bdfc10dd75cbf5ea8588f74781f0c490
+tests/images/examples.test__vector_data.true_winds.png b160ff21770e81514089fd1a58ed3bee137734aa78d830f1d7e3c522267f9d93
 tests/images/examples.test__vector_data.wind_arrows.png e73f77f2cb804e7d629be36ef8b0fcef1d64d653db005d29d2a2188575734809
 tests/images/examples.test__vector_data.winds_3d.png eb08f0d00a89d3429ab4e541d1b6c791d569d6bc6e359bf98d43f2937dff6753
 tests/images/examples.test__warp.from_unstructured__fvcom.png b31dba987a706601c60df6ab9c39cf03e83ed6d07acf84d2d10a798702762f08

src/geovista/common.py

@@ -71,6 +71,7 @@
     "to_lonlat",
     "to_lonlats",
     "triangulated",
+    "vectors_to_cartesian",
     "vtk_warnings_off",
     "vtk_warnings_on",
     "wrap",
@@ -753,6 +754,90 @@ def to_cartesian(
     return np.vstack(xyz).T if stacked else np.array(xyz)
 
 
+def vectors_to_cartesian(
+    lons: ArrayLike,
+    lats: ArrayLike,
+    vectors_uvw: tuple[ArrayLike, ArrayLike, ArrayLike],
+    radius: float | None = None,
+    zlevel: float | ArrayLike | None = None,
+    zscale: float | None = None,
+) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
+    """Transform geographic-oriented vector components to cartesian ``xyz`` components.
+
+    Parameters
+    ----------
+    lons, lats : ArrayLike
+        The longitude + latitude locations of the vectors (in degrees).
+        Both shapes must be the same.
+    vectors_uvw : tuple of ArrayLike
+        The eastward, northward and upward vector components.
+        All shapes must be the same as ``lons`` and ``lats``.
+    radius : float, optional
+        The radius of the sphere. Defaults to :data:`RADIUS`.
+    zlevel : int  or :data:`~numpy.typing.ArrayLike`, default=0.0
+        The z-axis level. Used in combination with the `zscale` to offset the
+        `radius` by a proportional amount i.e., ``radius * zlevel * zscale``.
+        NOTE : non-scalar `zlevel` is not actually supported, since it is planned to
+        drop support for this from :meth:`~geovista.bridge.Transform.from_points`.
+        It will raise an error.
+    zscale : float, optional
+        The proportional multiplier for z-axis `zlevel`. Defaults to
+        :data:`ZLEVEL_SCALE`.
+
+    Returns
+    -------
+    (ndarray, ndarray, ndarray)
+        The corresponding ``xyz`` cartesian vector components.
+
+    Notes
+    -----
+    .. versionadded:: 0.5.0
+
+    """
+    radius = RADIUS if radius is None else abs(float(radius))
+    zscale = ZLEVEL_SCALE if zscale is None else float(zscale)
+    # cast as float here, as from_cartesian use-case results in float zlevels
+    # that should be dtype preserved for the purpose of precision
+    if zlevel is None:
+        zlevel = 0
+    zlevel = np.atleast_1d(zlevel).astype(float)
+    if zlevel.size > 1:
+        # TODO @pp-mo: remove multiple zlevel, when removed from "from_points".
+        msg = f"'zlevel' may not be multiple, has shape {zlevel.shape}."
+        raise ValueError(msg)
+
+    radius += radius * zlevel * zscale
+
+    if lons.shape != lats.shape:
+        msg = f"'lons' and 'lats' do not have same shape: {lons.shape} != {lats.shape}."
+        raise ValueError(msg)
+
+    if any(x.shape != lons.shape for x in vectors_uvw):
+        msg = (
+            "some of 'vectors_uvw' do not have same shape as lons : "
+            f"{(x.shape for x in vectors_uvw)} != {lons.shape}."
+        )
+        raise ValueError(msg)
+
+    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.
+    # Right now, for simplicity, we just don't bother.
+    return (radius * wx, radius * wy, radius * wz)
+
+
 def to_lonlat(
     xyz: ArrayLike,
     radians: bool | None = False,

src/geovista/examples/vector_data/flow_directions.py

@@ -0,0 +1,70 @@
+#!/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.
+
+"""
+Horizontal Flow
+---------------
+
+This example demonstrates how to display equal-sized flow direction arrows.
+
+📋 Summary
+^^^^^^^^^^
+
+The data source provides X and Y arrays containing plain longitude and
+latitude values, with eastward- and northward-going flow components,
+which is the most common case.
+
+The sample flow data is actually provided in three separate surface-oriented
+component arrays, 'U, V and W', i.e. eastward, northward and vertical components,
+but we are only using the first two.  We use the
+:meth:`geovista.bridge.Transform.from_points` method, passing the winds to the
+``vectors`` keyword, producing a mesh of scattered points with attached vectors.
+
+The arrows themselves are created from this mesh via the
+:meth:`pyvista.DataSetFilters.glyph` method.
+
+In this example we display flow arrows of a fixed length, showing direction
+only, but with a colour scale indicating magnitude.
+"""  # noqa: D205,D212,D400
+
+from __future__ import annotations
+
+import geovista as gv
+from geovista.pantry.data import lfric_winds
+
+
+def main() -> None:
+    """Demonstrate a flow direction plot with fixed-length arrows."""
+    # get sample data
+    sample = lfric_winds()
+
+    mesh = gv.Transform.from_points(
+        sample.lons,
+        sample.lats,
+        vectors=(sample.u, sample.v),
+    )
+    # Note: with "scale=False", arrow size is fixed and controlled by "factor".
+    arrows = mesh.glyph(factor=0.1, scale=False)
+
+    p = gv.GeoPlotter()
+    p.add_base_layer()
+    p.add_mesh(arrows, cmap="magma")
+    p.add_coastlines()
+    p.add_graticule()
+    p.add_axes()
+
+    p.camera.zoom(1.3)
+    selected_view = [
+        (-4.0688208659033505, -2.5462610064466777, -2.859304866708606),
+        (-0.0037798285484313965, 0.005168497562408447, -0.0031679868698120117),
+        (-0.523382090763761, -0.11174892277533728, 0.8447386372874786),
+    ]
+    p.camera_position = selected_view
+    p.show()
+
+
+if __name__ == "__main__":
+    main()