diff --git a/changelog/676.feature.rst b/changelog/676.feature.rst new file mode 100644 index 00000000..5d962990 --- /dev/null +++ b/changelog/676.feature.rst @@ -0,0 +1,2 @@ +Added support for 3D vectors in :meth:`geovista.bridge.Transform.from_points` + (:user:`pp-mo`) diff --git a/src/geovista/bridge.py b/src/geovista/bridge.py index a623766e..874ad49d 100644 --- a/src/geovista/bridge.py +++ b/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) diff --git a/src/geovista/cache/__init__.py b/src/geovista/cache/__init__.py index 258e6abd..a1da4a46 100644 --- a/src/geovista/cache/__init__.py +++ b/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" diff --git a/src/geovista/cache/registry.txt b/src/geovista/cache/registry.txt index f855dfcf..60bc30fb 100644 --- a/src/geovista/cache/registry.txt +++ b/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 diff --git a/src/geovista/common.py b/src/geovista/common.py index a66b244c..b52651ac 100644 --- a/src/geovista/common.py +++ b/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, diff --git a/src/geovista/examples/vector_data/flow_directions.py b/src/geovista/examples/vector_data/flow_directions.py new file mode 100755 index 00000000..7d36a4e7 --- /dev/null +++ b/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() diff --git a/src/geovista/examples/vector_data/true_winds.py b/src/geovista/examples/vector_data/true_winds.py new file mode 100755 index 00000000..e2d3de6c --- /dev/null +++ b/src/geovista/examples/vector_data/true_winds.py @@ -0,0 +1,90 @@ +#!/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. + +""" +"True" Winds on a Foreign Grid +------------------------------ + +This example demonstrates how to display "true" winds, when attached to data points +given in a non-standard coordinate system. + +📋 Summary +^^^^^^^^^^ + +The example constructs some sample points on a polar stereographic grid, and assigns +synthetic vector values to them, representing wind data. These data are presented as +"true" winds, i.e. northward and eastward components, so *not* related to the sample +grid. This is a fairly typical real-world case. + +In this case, the synthetic wind is a steady 4.0 ms-1 Eastward everywhere, so the +meaning is clearly seen. + +We use the :meth:`geovista.bridge.Transform.from_points` method, passing the +winds to the ``vectors`` keyword, and specifying the different coordinate reference +systems of both the sample points (``crs``) and the vectors (``vectors_crs``). + +Please refer to other code examples, e.g. +:ref:`Wind Arrows Plot ` +for more details of the basic methods used. + +""" # noqa: D205,D212,D400 + +from __future__ import annotations + +# Use cartopy as a source for specifying Coordinate Reference Systems +# (which we can translate into PROJ via the .to_wkt() method) +import cartopy.crs as ccrs +import numpy as np + +import geovista as gv + + +def main() -> None: + """Demonstrate plotting vectors on different CRS from points.""" + # Create a mesh of individual points, adding vectors at each point. + # NOTE: this creates a mesh with 'mesh vectors' : a specific concept in PyVista. + + # Create a CRS for "ordinary latitude-longitude" + crs_latlon = ccrs.Geodetic().to_wkt() + # Create a CRS for "polar stereographic" + crs_polar = ccrs.NorthPolarStereo().to_wkt() + # Create polar grid points + xx = np.linspace(-5.0e6, 5.0e6, 20) + yy = np.linspace(-4.0e6, 2.0e6, 20) + xx, yy = np.meshgrid(xx, yy) + # Make vector component arrays matching the XX and YY arrays + xx, yy = xx.flatten(), yy.flatten() + vx, vy = 4.0 * np.ones_like(xx), np.zeros_like(yy) + + # Create a mesh of location points with attached vector information + mesh = gv.Transform.from_points( + xx, yy, vectors=(vx, vy), crs=crs_polar, vectors_crs=crs_latlon + ) + + # Create a new mesh containing arrow glyphs, from the mesh vectors. + arrows = mesh.glyph(factor=0.02) + + # Add the arrows to a Plotter with other aspects, and display + p = gv.GeoPlotter() + # Scale the base layer slightly to ensure it remains 'below' other elements. + p.add_base_layer(radius=0.99) + p.add_mesh(arrows, color="red") + p.add_coastlines() + p.add_graticule() + p.add_axes() + + # Set up a nice default view + selected_view = [ + (0.6037050543418041, -0.011033743353827528, 3.069575190259155), + (0.0, 0.0, 0.0028896447726441954), + (-0.9809865744261353, -0.019981215106321615, 0.19304427429621296), + ] + p.camera_position = selected_view + p.show() + + +if __name__ == "__main__": + main() diff --git a/src/geovista/examples/vector_data/wind_arrows.py b/src/geovista/examples/vector_data/wind_arrows.py new file mode 100755 index 00000000..7d31e0d3 --- /dev/null +++ b/src/geovista/examples/vector_data/wind_arrows.py @@ -0,0 +1,81 @@ +#!/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. + +""" +Wind Arrows Plot +---------------- + +This example demonstrates how to display horizontal winds. + +📋 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 point is a separate location in a mesh of +scattered points, and each point has an associated vector value independent of +the others. 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. + +Here we show just horizontal winds (U, V), which are usually of the most interest. +""" # noqa: D205,D212,D400 + +from __future__ import annotations + +import geovista as gv +from geovista.pantry.data import lfric_winds + + +def main() -> None: + """Demonstrate horizontal wind arrows plotting.""" + # 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: apply 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 + p = gv.GeoPlotter() + # Scale the base layer slightly to ensure it remains 'below' other elements. + p.add_base_layer(radius=0.99) + p.add_mesh(arrows, cmap="inferno") + p.add_coastlines() + p.add_graticule() + p.add_axes() + + # Set up a nice default view + p.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), + ] + p.camera_position = selected_view + p.show() + + +if __name__ == "__main__": + main() diff --git a/src/geovista/examples/vector_data/winds_3d.py b/src/geovista/examples/vector_data/winds_3d.py new file mode 100755 index 00000000..efa13243 --- /dev/null +++ b/src/geovista/examples/vector_data/winds_3d.py @@ -0,0 +1,82 @@ +#!/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 plot 3D wind vectors. + +📋 Summary +^^^^^^^^^^ + +The data source provides X and Y arrays containing plain longitude and +latitude values, which is the most common case. + +3D wind components are provided in three separate field arrays, 'U, V and W', +i.e. eastward, northward and vertical components. + +There is no connectivity provided, so each location has its own attached vector, and is +independent of the others. 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. + +Here, we display 3-dimensional wind arrows. +We have amplified the "W" components by a considerable factor, which is typical since +vertical winds are generally much smaller in magnitude and otherwise tend to be not +very visible. +""" # noqa: D205,D212,D400 + +from __future__ import annotations + +import geovista as gv +from geovista.pantry.data import lfric_winds + + +def main() -> None: + """Demonstrate 3-dimensional wind arrows plotting.""" + # 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, + # supply all three components, but with a big extra scaling on the "W" values + vectors=(sample.u, sample.v, 1500.0 * sample.w), + # offset from surface to avoid downward-pointing arrows disappearing + radius=1.1, + ) + + # Create a new mesh containing arrow glyphs, from the mesh vectors. + # NOTE: choose 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 + p = gv.GeoPlotter() + # Scale the base layer slightly to ensure it remains 'below' other elements. + p.add_base_layer(radius=0.99) + p.add_mesh(arrows, cmap="inferno") + p.add_coastlines() + p.add_graticule() + p.add_axes() + + # Set up a suitable camera view + p.camera.zoom(1.3) + selected_view = [ + (0.6917810912064826, -3.065688850990997, 0.4317999141924935), + (0.41358279170396495, 0.07362917740509836, 0.5091223320854129), + (0.8088496364623022, 0.05726400555597287, 0.5852205560833343), + ] + p.camera_position = selected_view + p.show() + + +if __name__ == "__main__": + main() diff --git a/src/geovista/pantry/data.py b/src/geovista/pantry/data.py index 557e7789..86885526 100644 --- a/src/geovista/pantry/data.py +++ b/src/geovista/pantry/data.py @@ -127,6 +127,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 | None = field(default=None) + units: str | None = field(default=None) + steps: int | None = field(default=None) + ndim: int = 1 + + def capitalise(title: str) -> str: """Capitalise each word and replace inappropriate characters. @@ -969,3 +984,38 @@ def ww3_global_tri() -> SampleUnstructuredXY: name=name, units=units, ) + + +@lru_cache(maxsize=LRU_CACHE_SIZE) +def lfric_winds() -> SampleVectorsXYUVW: + """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) diff --git a/tests/bridge/test_from_points.py b/tests/bridge/test_from_points.py index 80e9c519..3a525b09 100644 --- a/tests/bridge/test_from_points.py +++ b/tests/bridge/test_from_points.py @@ -7,6 +7,7 @@ from __future__ import annotations +import cartopy.crs as ccrs import numpy as np from pyproj import CRS, Transformer import pytest @@ -76,3 +77,156 @@ def test_scalar(): np.testing.assert_array_equal(Transform.from_points(0, [90]).points, expected) np.testing.assert_array_equal(Transform.from_points([0], 90).points, expected) np.testing.assert_array_equal(Transform.from_points([0], [90]).points, expected) + + +class TestVectors: + """Check calculations on attached vectors. + + Mostly testing against previously-obtained results, but we do have examples proving + practically correct behaviours : see src/geovista/examples/vector_data + """ + + @pytest.fixture(autouse=True) + def setup(self): + """Set useful test constants.""" + self.lats = np.array([10.0, 80.0, -70.0, 35.0]) + self.lons = np.array([5.0, 140.0, -200.0, 73]) + self.easting = 10.0e6 * np.array([1.0, 8.0, 13.0, 25.0]) + self.northing = 10.0e6 * np.array([-2.0, 5.0, 17.0, -14.0]) + self.u = np.array([10.0, 20, 15.0, 24.0]) + self.v = np.array([20.0, -17, 0.0, 33.0]) + self.w = np.array([15.0, 25, 4.0, -55.0]) + self.crs_truelatlon = WGS84 + self.crs_rotatedlatlon = ccrs.RotatedGeodetic(130.0, 65.0).to_wkt() + self.crs_northpolar = ccrs.NorthPolarStereo().to_wkt() + + def test_basic(self): + """Check basic operation on true-latlon uvw vectors.""" + mesh = Transform.from_points( + xs=self.lons, ys=self.lats, vectors=(self.u, self.v, self.w) + ) + result = mesh["vectors"].T + expected = np.array( + [ + [10.385, -29.006, -6.416, -41.658], + [10.947, -1.769, -13.627, -54.169], + [22.301, 21.668, -3.759, -4.515], + ] + ) + assert np.allclose(result, expected, atol=0.001) + + def test_nonarrays(self): + """Check basic operation with lists of numbers in place of array vectors.""" + mesh = Transform.from_points( + xs=self.lons, + ys=self.lats, + vectors=(list(self.u), list(self.v), list(self.w)), + ) + result = mesh["vectors"].T + expected = np.array( + [ + [10.385, -29.006, -6.416, -41.658], + [10.947, -1.769, -13.627, -54.169], + [22.301, 21.668, -3.759, -4.515], + ] + ) + assert np.allclose(result, expected, atol=0.001) + + def test_basic__2d_uv(self): + """Check operation with only 2 input component arrays (no W).""" + mesh = Transform.from_points( + xs=self.lons, ys=self.lats, vectors=(self.u, self.v) + ) + result = mesh["vectors"].T + expected = np.array( + [ + [-4.331, -25.681, -5.13, -28.485], + [9.659, -4.56, -14.095, -11.084], + [19.696, -2.952, 0.0, 27.032], + ] + ) + assert np.allclose(result, expected, atol=0.001) + + def test_crs(self): + """Check operation with alternate latlon-type CRS.""" + mesh = Transform.from_points( + xs=self.lons, + ys=self.lats, + vectors=(self.u, self.v), + crs=self.crs_rotatedlatlon, + ) + result = mesh["vectors"].T + expected = np.array( + [ + [-0.474, -17.651, -13.786, -32.429], + [15.592, 13.943, -5.499, 21.403], + [16.02, -13.529, -2.168, 12.461], + ] + ) + assert np.allclose(result, expected, atol=0.001) + + def test__nonlatlon_crs__fail(self): + """Check error when attempted with non-latlon CRS.""" + msg = "Cannot determine wind directions : Target CRS type is not supported.*" + with pytest.raises(ValueError, match=msg): + _ = Transform.from_points( + xs=self.easting, + ys=self.northing, + vectors=(self.u, self.v), + crs=self.crs_northpolar, + ) + + def test__nonlatloncrs__truelatlon__vectorscrs(self): + """Check ok with non-latlon CRS for points but latlon vectors.""" + mesh = Transform.from_points( + xs=self.easting, + ys=self.northing, + vectors=(self.u, self.v), + crs=self.crs_northpolar, + vectors_crs=self.crs_truelatlon, + ) + result = mesh["vectors"].T + expected = np.array( + [ + [4.722, -8.267, -9.112, -4.879], + [13.541, -24.508, -11.915, 40.408], + [17.156, -4.472, 0.0, 2.903], + ] + ) + assert np.allclose(result, expected, atol=0.001) + + def test__latlon__vectorscrs(self): + """Check operation with different specified CRS for vectors only.""" + mesh = Transform.from_points( + xs=self.lons, + ys=self.lats, + vectors=(self.u, self.v), + vectors_crs=self.crs_rotatedlatlon, + ) + result = mesh["vectors"].T + expected = np.array( + [ + [-4.066, 25.821, -8.955, -38.46], + [16.031, -2.79, -11.93, 1.87], + [15.049, 3.804, 1.578, 13.504], + ] + ) + assert np.allclose(result, expected, atol=0.001) + + def test__vectors_array_name(self): + """Check operation with alternate vectors array name.""" + mesh = Transform.from_points( + xs=self.lons, + ys=self.lats, + vectors=(self.u, self.v, self.w), + vectors_array_name="squiggle", + ) + result = mesh["squiggle"].T + expected = np.array( + [ + [10.385, -29.006, -6.416, -41.658], + [10.947, -1.769, -13.627, -54.169], + [22.301, 21.668, -3.759, -4.515], + ] + ) + assert np.allclose(result, expected, atol=0.001)