From e5eee73d3f179ffd1ea382394e8c4a6d9de220cd Mon Sep 17 00:00:00 2001 From: Bill Little Date: Mon, 18 Sep 2023 14:37:50 +0100 Subject: [PATCH] support cartographic add_points (#442) * decorate dunder inits with return types * add_points * fix doc-strings * minor refactor * refactor earthquakes example * add_points doc-string and traps * fix mocker patch pattern * test coverage --- src/geovista/examples/earthquakes.py | 7 +- src/geovista/geoplotter.py | 173 +++++++++++++++++++++++- tests/geometry/test_coastlines.py | 28 ++-- tests/geoplotter/test_add_points.py | 150 ++++++++++++++++++++ tests/transform/test_transform_point.py | 6 +- 5 files changed, 340 insertions(+), 24 deletions(-) create mode 100644 tests/geoplotter/test_add_points.py diff --git a/src/geovista/examples/earthquakes.py b/src/geovista/examples/earthquakes.py index 3fde1208..689ab8f5 100755 --- a/src/geovista/examples/earthquakes.py +++ b/src/geovista/examples/earthquakes.py @@ -11,7 +11,6 @@ from warnings import warn import geovista as gv -from geovista.common import to_cartesian from geovista.pantry import usgs_earthquakes import geovista.theme # noqa: F401 @@ -56,14 +55,12 @@ def main() -> None: warn(wmsg, stacklevel=2) return - # convert coordinate to cartesian - points = to_cartesian(lons=sample.lons, lats=sample.lats) - # plot the mesh plotter = gv.GeoPlotter() sargs = {"title": "Magnitude", "shadow": True} plotter.add_points( - points, + xs=sample.lons, + ys=sample.lats, cmap="fire_r", render_points_as_spheres=True, scalars=sample.data, diff --git a/src/geovista/geoplotter.py b/src/geovista/geoplotter.py index 96997327..545175a4 100644 --- a/src/geovista/geoplotter.py +++ b/src/geovista/geoplotter.py @@ -13,9 +13,12 @@ from typing import Any from warnings import warn +import numpy.typing as npt from pyproj import CRS import pyvista as pv +import pyvista.core.utilities.helpers as helpers +from .bridge import Transform from .common import ( GV_FIELD_ZSCALE, GV_REMESH_POINT_IDS, @@ -32,8 +35,10 @@ from .core import add_texture_coords, resize, slice_mesh from .crs import ( WGS84, + CRSLike, from_wkt, get_central_meridian, + has_wkt, projected, set_central_meridian, to_wkt, @@ -51,6 +56,9 @@ __all__ = ["GeoPlotter"] +#: The valid 'style' options for adding points. +ADD_POINTS_STYLE: list[str, ...] = ["points", "points_gaussian"] + #: Proportional multiplier for z-axis levels/offsets of base-layer mesh. BASE_ZLEVEL_SCALE: int = 1.0e-3 @@ -436,7 +444,7 @@ def add_graticule( meridian. Defaults to :data:`geovista.gridlines.LONGITUDE_START`. lon_stop : float, optional The last line of longitude (degrees). The graticule will include this - meridian when it is a multiple of ``lon_step``. Also see + meridian when it is a multiple of `lon_step`. Also see ``closed_interval``. Defaults to :data:`geovista.gridlines.LONGITUDE_STOP`. lon_step : float, optional The delta (degrees) between neighbouring meridians. Defaults to @@ -447,7 +455,7 @@ def add_graticule( :data:`geovista.gridlines.LATITUDE_START`. lat_stop : float, optional The last line of latitude (degrees). The graticule will include this - parallel when it is a multiple of ``lat_step``. Defaults to + parallel when it is a multiple of `lat_step`. Defaults to :data:`geovista.gridlines.LATITUDE_STOP`. lat_step : float, optional The delta (degrees) between neighbouring parallels. Defaults to @@ -509,7 +517,7 @@ def add_graticule( ) def add_mesh(self, mesh: Any, **kwargs: Any | None) -> pv.Actor: - """Add the ``mesh`` to the plotter scene. + """Add the `mesh` to the plotter scene. See :meth:`pyvista.Plotter.add_mesh`. @@ -704,7 +712,7 @@ def add_meridians( meridian. Defaults to :data:`geovista.gridlines.LONGITUDE_START`. stop : float, optional The last line of longitude (degrees). The graticule will include this - meridian when it is a multiple of ``step``. Also see ``closed_interval``. + meridian when it is a multiple of `step`. Also see ``closed_interval``. Defaults to :data:`geovista.gridlines.LONGITUDE_STOP`. step : float, optional The delta (degrees) between neighbouring meridians. Defaults to @@ -872,11 +880,11 @@ def add_parallels( ---------- start : float, optional The first line of latitude (degrees). The graticule will include this - parallel. Also see ``poles_parallel``. Defaults to + parallel. Also see `poles_parallel`. Defaults to :data:`geovista.gridlines.LATITUDE_START`. stop : float, optional The last line of latitude (degrees). The graticule will include this - parallel when it is a multiple of ``step``. Also see ``poles_parallel`. + parallel when it is a multiple of `step`. Also see `poles_parallel`. Defaults to :data:`geovista.gridlines.LATITUDE_STOP`. step : float, optional The delta (degrees) between neighbouring parallels. Defaults to @@ -965,6 +973,159 @@ def add_parallels( point_labels_args=point_labels_args, ) + def add_points( + self, + points: npt.ArrayLike | pv.PolyData | None = None, + xs: npt.ArrayLike | None = None, + ys: npt.ArrayLike | None = None, + scalars: str | npt.ArrayLike | None = None, + crs: CRSLike | None = None, + radius: float | None = None, + style: str | None = None, + zlevel: int | npt.ArrayLike | None = None, + zscale: float | None = None, + **kwargs: Any | None, + ) -> pv.Actor: + """Add points to the plotter scene. + + See :meth:`pyvista.Plotter.add_mesh`. + + Parameters + ---------- + points : ArrayLike or PolyData, optional + Array of xyz points, in canonical `crs` units, or the points of the mesh + to be rendered. + xs : ArrayLike, optional + A 1-D, 2-D or 3-D array of point-cloud x-values, in canonical `crs` units. + Must have the same shape as the `ys`. + ys : ys : ArrayLike + A 1-D, 2-D or 3-D array of point-cloud y-values, in canonical `crs` units. + Must have the same shape as the `xs`. + scalars : str or ArrayLike, optional + Values used to color the points. Either, the string name of an array that is + present on the `points` mesh or an array equal to the number of points. + Alternatively, an array of values equal to the number of points to be + rendered. If both `color` and `scalars` are ``None``, then the active + scalars on the `points` mesh are used. + crs : CRSLike, optional + The Coordinate Reference System of the provided `points`, or `xs` and `ys`. + May be anything accepted by :meth:`pyproj.CRS.from_user_input`. Defaults + to ``EPSG:4326`` i.e., ``WGS 84``. + radius : float, optional + The radius of the mesh point-cloud. Defaults to + :data:`geovista.common.RADIUS`. + style : str, optional + Visualization style of the points to be rendered. Maybe either ``points`` + or ``points_gaussian``. The ``points_gaussian`` option maybe controlled + with the ``emissive`` and ``render_points_as_spheres`` options. + zlevel : int or ArrayLike, default=0 + The z-axis level. Used in combination with the `zscale` to offset the + `radius` by a proportional amount i.e., ``radius * zlevel * zscale``. + If `zlevel` is not a scalar, then its shape must match or broadcast + with the shape of the `xs` and `ys`. + zscale : float, optional + The proportional multiplier for z-axis `zlevel`. Defaults to + :data:`geovista.common.ZLEVEL_SCALE`. + **kwargs : dict, optional + See :meth:`pyvista.Plotter.add_mesh`. + + Returns + ------- + Actor + The rendered actor added to the plotter scene. + + Notes + ----- + .. versionadded:: 0.4.0 + + """ + if crs is not None: + # sanity check the source crs + crs = CRS.from_user_input(crs) + + if style is None: + style = "points" + + if style not in ADD_POINTS_STYLE: + options = "or ".join(f"{option!r}" for option in ADD_POINTS_STYLE) + emsg = ( + f"Invalid 'style' for 'add_points', expected {options}, " + f"got {style!r}." + ) + raise ValueError(emsg) + + if points is None and xs is None and ys is None: + emsg = ( + "Require either 'points' or both 'xs' and 'ys' to be specified, " + "got neither." + ) + raise ValueError(emsg) + + if points is not None and xs is not None and ys is not None: + emsg = ( + "Require either 'points' or both 'xs' and 'ys' to be specified, " + "got both 'points', and 'xs' and 'ys'." + ) + raise ValueError(emsg) + + if points is not None: + if xs is not None or ys is not None: + emsg = ( + "Require either 'points' or both 'xs' and 'ys' to be specified, " + "got both 'points', and 'xs' or 'ys'." + ) + raise ValueError(emsg) + + if not helpers.is_pyvista_dataset(points): + points = helpers.wrap(points) + + mesh = points + + if crs is not None: + if has_wkt(mesh): + other = from_wkt(mesh) + if other != crs: + emsg = ( + "The CRS serialized as WKT on the 'points' mesh does not " + "match the provided 'crs'." + ) + raise ValueError(emsg) + else: + # serialize the provided CRS on the points mesh as wkt + to_wkt(mesh, crs) + elif not has_wkt(mesh): + # assume CRS is WGS84 + to_wkt(mesh, WGS84) + else: + if xs is None or ys is None: + emsg = ( + "Require either 'points', or both 'xs' and 'ys' to be specified, " + "got only 'xs' or 'ys'." + ) + raise ValueError(emsg) + + if isinstance(scalars, str): + wmsg = ( + f"Ignoring the 'scalars' string name '{scalars}', as no 'points' " + "mesh was provided." + ) + warn(wmsg, stacklevel=2) + + mesh = Transform.from_points( + xs=xs, + ys=ys, + crs=crs, + radius=radius, + zlevel=zlevel, + zscale=zscale, + ) + + # defensive kwarg pop + if "texture" in kwargs: + _ = kwargs.pop("texture") + + return self.add_mesh(mesh, style=style, scalars=scalars, **kwargs) + class GeoPlotter(GeoPlotterBase, pv.Plotter): """A geospatial aware plotter. diff --git a/tests/geometry/test_coastlines.py b/tests/geometry/test_coastlines.py index b5c7ffa8..cebb4374 100644 --- a/tests/geometry/test_coastlines.py +++ b/tests/geometry/test_coastlines.py @@ -8,11 +8,13 @@ def test_defaults(mocker): """Test expected defaults are honoured.""" mesh = mocker.sentinel.mesh - fetch = mocker.patch("geovista.geometry.fetch_coastlines", return_value=mesh) - resize = mocker.patch("geovista.geometry.resize") + _ = mocker.patch("geovista.geometry.fetch_coastlines", return_value=mesh) + _ = mocker.patch("geovista.geometry.resize") result = coastlines() assert result == mesh - fetch.assert_called_once_with(resolution=COASTLINES_RESOLUTION) + from geovista.geometry import fetch_coastlines, resize + + fetch_coastlines.assert_called_once_with(resolution=COASTLINES_RESOLUTION) resize.assert_called_once_with( mesh, radius=None, zlevel=1, zscale=None, inplace=True ) @@ -25,23 +27,27 @@ def test_resize_kwarg_pass_thru(mocker): radius = mocker.sentinel.radius zscale = mocker.sentinel.zscale zlevel = mocker.sentinel.zlevel - fetch = mocker.patch("geovista.geometry.fetch_coastlines", return_value=mesh) - resize = mocker.patch("geovista.geometry.resize") + _ = mocker.patch("geovista.geometry.fetch_coastlines", return_value=mesh) + _ = mocker.patch("geovista.geometry.resize") kwargs = {"radius": radius, "zscale": zscale, "zlevel": zlevel} result = coastlines(resolution=resolution, **kwargs) assert result == mesh - fetch.assert_called_once_with(resolution=resolution) + from geovista.geometry import fetch_coastlines, resize + + fetch_coastlines.assert_called_once_with(resolution=resolution) resize.assert_called_once_with(mesh, **kwargs, inplace=True) def test_fetch_exception(mocker): """Test coastlines are loaded on cache fetch failure.""" mesh = mocker.sentinel.mesh - fetch = mocker.patch("geovista.geometry.fetch_coastlines", side_effect=ValueError) - load = mocker.patch("geovista.geometry.load_coastlines", return_value=mesh) - resize = mocker.patch("geovista.geometry.resize") + _ = mocker.patch("geovista.geometry.fetch_coastlines", side_effect=ValueError) + _ = mocker.patch("geovista.geometry.load_coastlines", return_value=mesh) + _ = mocker.patch("geovista.geometry.resize") result = coastlines() assert result == mesh - assert fetch.call_count == 1 + from geovista.geometry import fetch_coastlines, load_coastlines, resize + + assert fetch_coastlines.call_count == 1 assert resize.call_count == 1 - load.assert_called_once_with(resolution=COASTLINES_RESOLUTION) + load_coastlines.assert_called_once_with(resolution=COASTLINES_RESOLUTION) diff --git a/tests/geoplotter/test_add_points.py b/tests/geoplotter/test_add_points.py new file mode 100644 index 00000000..85fd8823 --- /dev/null +++ b/tests/geoplotter/test_add_points.py @@ -0,0 +1,150 @@ +"""Unit-tests for :meth:`geovista.geoplotter.GeoPlotter.add_points`.""" +from __future__ import annotations + +from pyproj import CRS +from pyproj.exceptions import CRSError +import pytest + +from geovista.crs import WGS84, from_wkt, has_wkt, to_wkt +from geovista.geoplotter import GeoPlotter + + +def test_style_fail(): + """Test trap of invalid style request.""" + plotter = GeoPlotter() + emsg = "Invalid 'style' for 'add_points'" + with pytest.raises(ValueError, match=emsg): + _ = plotter.add_points(style="dolce and gabbana") + + +def test_no_points_fail(): + """Test trap of no spatial points provided.""" + plotter = GeoPlotter() + emsg = "got neither" + with pytest.raises(ValueError, match=emsg): + _ = plotter.add_points() + + +def test_duplicate_points_fail(mocker): + """Test trap of duplicate spatial points provided.""" + plotter = GeoPlotter() + points, xs, ys = mocker.sentinel.points, mocker.sentinel.xs, mocker.sentinel.ys + emsg = "got both 'points', and 'xs' and 'ys'" + with pytest.raises(ValueError, match=emsg): + _ = plotter.add_points(points=points, xs=xs, ys=ys) + + +@pytest.mark.parametrize("xs, ys", [(True, None), (None, True)]) +def test_over_specified_points_fail(mocker, xs, ys): + """Test trap of points and xs or ys provided.""" + plotter = GeoPlotter() + points = mocker.sentinel.points + emsg = "got both 'points', and 'xs' or 'ys'" + with pytest.raises(ValueError, match=emsg): + _ = plotter.add_points(points=points, xs=xs, ys=ys) + + +def test_points_mesh_crs_fail(sphere): + """Test trap of crs mismatch with points mesh serialized crs as wkt.""" + plotter = GeoPlotter() + emsg = "The CRS serialized as WKT on the 'points' mesh does not match" + to_wkt(sphere, WGS84) + with pytest.raises(ValueError, match=emsg): + _ = plotter.add_points(points=sphere, crs="+proj=eqc") + + +@pytest.mark.parametrize("crs", [None, CRS.from_user_input("+proj=eqc")]) +def test_points_mesh_with_no_crs(mocker, sphere, crs): + """Test points mesh with no serialized crs.""" + assert not has_wkt(sphere) + plotter = GeoPlotter() + actor = mocker.sentinel.actor + _ = mocker.patch("geovista.geoplotter.GeoPlotterBase.add_mesh", return_value=actor) + _ = plotter.add_points(sphere, crs=crs) + plotter.add_mesh.assert_called_once_with(sphere, style="points", scalars=None) + assert has_wkt(sphere) + expected = WGS84 if crs is None else crs + assert from_wkt(sphere) == expected + + +def test_crs_invalid_fail(): + """Test trap of crs validity sanity check.""" + plotter = GeoPlotter() + emsg = "Invalid projection" + with pytest.raises(CRSError, match=emsg): + _ = plotter.add_points(crs="invalid") + + +def test_points_mesh_with_scalars(mocker, sphere): + """Test points mesh with scalars pass-through.""" + plotter = GeoPlotter() + _ = mocker.patch("geovista.geoplotter.GeoPlotterBase.add_mesh") + scalars = mocker.sentinel.scalars + _ = plotter.add_points(sphere, scalars=scalars) + plotter.add_mesh.assert_called_once_with(sphere, style="points", scalars=scalars) + + +def test_texture_kwarg_pop(mocker, sphere): + """Test defensive texture kwarg purge.""" + plotter = GeoPlotter() + _ = mocker.patch("geovista.geoplotter.GeoPlotterBase.add_mesh") + texture = mocker.sentinel.texture + _ = plotter.add_points(sphere, texture=texture) + plotter.add_mesh.assert_called_once_with(sphere, style="points", scalars=None) + + +@pytest.mark.parametrize("xs, ys", [(None, True), (True, None)]) +def test_xs_ys_incomplete_fail(xs, ys): + """Test trap of missing spatial points.""" + plotter = GeoPlotter() + emsg = "got only 'xs' or 'ys'" + with pytest.raises(ValueError, match=emsg): + _ = plotter.add_points(xs=xs, ys=ys) + + +def test_xs_ys(mocker): + """Test xs/ys points.""" + plotter = GeoPlotter() + mesh = mocker.sentinel.mesh + _ = mocker.patch("geovista.Transform.from_points", return_value=mesh) + actor = mocker.sentinel.actor + _ = mocker.patch("geovista.geoplotter.GeoPlotterBase.add_mesh", return_value=actor) + xs, ys = mocker.sentinel.xs, mocker.sentinel.ys + crs, scalars = WGS84, mocker.sentinel.scalars + radius, zlevel, zscale = ( + mocker.sentinel.radius, + mocker.sentinel.zlevel, + mocker.sentinel.zscale, + ) + result = plotter.add_points( + xs=xs, + ys=ys, + scalars=scalars, + crs=crs, + radius=radius, + zlevel=zlevel, + zscale=zscale, + ) + assert result == actor + from geovista.bridge import Transform + + Transform.from_points.assert_called_once_with( + xs=xs, ys=ys, crs=crs, radius=radius, zlevel=zlevel, zscale=zscale + ) + plotter.add_mesh.assert_called_once_with(mesh, style="points", scalars=scalars) + + +def test_xy_ys_scalars_name_warnings(mocker): + """Test xy/ys points invalid scalars array name warning is raised.""" + plotter = GeoPlotter() + mesh = mocker.sentinel.mesh + _ = mocker.patch("geovista.Transform.from_points", return_value=mesh) + _ = mocker.patch("geovista.geoplotter.GeoPlotterBase.add_mesh") + xs, ys = mocker.sentinel.xs, mocker.sentinel.ys + scalars = "invalid" + wmsg = ( + f"Ignoring the 'scalars' string name '{scalars}', " + "as no 'points' mesh was provided" + ) + with pytest.warns(UserWarning, match=wmsg): + _ = plotter.add_points(xs=xs, ys=ys, scalars=scalars) diff --git a/tests/transform/test_transform_point.py b/tests/transform/test_transform_point.py index c966e32e..1533bec7 100644 --- a/tests/transform/test_transform_point.py +++ b/tests/transform/test_transform_point.py @@ -18,11 +18,13 @@ def test_shape_fail(mocker, bad): tgt_crs = mocker.sentinel.tgt_crs x, y, z = mocker.sentinel.x, mocker.sentinel.y, mocker.sentinel.z trap = mocker.sentinel.trap - patcher = mocker.patch("geovista.transform.transform_points", return_value=bad) + _ = mocker.patch("geovista.transform.transform_points", return_value=bad) emsg = "Cannot transform point, got unexpected shape" with pytest.raises(AssertionError, match=emsg): _ = transform_point(src_crs=src_crs, tgt_crs=tgt_crs, x=x, y=y, z=z, trap=trap) - patcher.assert_called_once_with( + from geovista.transform import transform_points + + transform_points.assert_called_once_with( src_crs=src_crs, tgt_crs=tgt_crs, xs=x,