From b68bebbf945a9f4ac455c8255e7d932040a7e769 Mon Sep 17 00:00:00 2001 From: Andrew Huang Date: Thu, 6 Jul 2023 10:22:17 -0400 Subject: [PATCH 1/5] Add VectorField project operation --- geoviews/operation/__init__.py | 2 +- geoviews/operation/projection.py | 42 ++++++++++++++++++++++++++--- geoviews/plotting/bokeh/__init__.py | 4 +-- geoviews/plotting/mpl/__init__.py | 4 +-- 4 files changed, 43 insertions(+), 9 deletions(-) diff --git a/geoviews/operation/__init__.py b/geoviews/operation/__init__.py index 2ff2bb32..23c7f153 100644 --- a/geoviews/operation/__init__.py +++ b/geoviews/operation/__init__.py @@ -6,7 +6,7 @@ from ..element import _Element from .projection import ( # noqa (API import) project_image, project_path, project_shape, project_points, - project_graph, project_quadmesh, project_geom, project) + project_graph, project_quadmesh, project_geom, project_vectorfield, project) from .resample import resample_geometry # noqa (API import) geo_ops = [contours, bivariate_kde] diff --git a/geoviews/operation/projection.py b/geoviews/operation/projection.py index d956f2f4..9fd9f14f 100644 --- a/geoviews/operation/projection.py +++ b/geoviews/operation/projection.py @@ -157,7 +157,7 @@ def _process_element(self, element): class project_points(_project_operation): - supported_types = [Points, Nodes, VectorField, HexTiles, Labels] + supported_types = [Points, Nodes, HexTiles, Labels] def _process_element(self, element): if not len(element): @@ -194,8 +194,6 @@ class project_geom(_project_operation): supported_types = [Rectangles, Segments] def _process_element(self, element): - if not len(element): - return element.clone(crs=self.p.projection) x0d, y0d, x1d, y1d = element.kdims x0, y0, x1, y1 = (element.dimension_values(i) for i in range(4)) p1 = self.p.projection.transform_points(element.crs, x0, y0) @@ -222,6 +220,42 @@ def _process_element(self, element): return element.clone(tuple(new_data[d.name] for d in element.dimensions()), crs=self.p.projection) +class project_vectorfield(_project_operation): + + supported_types = [VectorField] + + def _process_element(self, element): + if not len(element): + return element.clone(crs=self.p.projection) + + xdim, ydim, adim, mdim = element.dimensions()[:4] + xs, ys, ang, ms = (element.dimension_values(i) for i in range(4)) + coordinates = self.p.projection.transform_points(element.crs, xs, ys) + mask = np.isfinite(coordinates[:, 0]) + new_data = {k: v[mask] for k, v in element.columns().items()} + new_data[xdim.name] = coordinates[mask, 0] + new_data[ydim.name] = coordinates[mask, 1] + datatype = [element.interface.datatype]+element.datatype + us = np.sin(ang) * ms + vs = np.cos(ang) * ms + ut, vt = self.p.projection.transform_vectors(element.crs, xs, ys, us, vs) + with np.errstate(divide='ignore', invalid='ignore'): + angle = np.arctan2(vt, ut) + mag = np.sqrt(ut**2+vt**2) + + new_data[adim.name] = angle[mask] + new_data[mdim.name] = mag[mask] + + if len(new_data[xdim.name]) == 0: + self.warning('While projecting a {} element from a {} coordinate ' + 'reference system (crs) to a {} projection none of ' + 'the projected paths were contained within the bounds ' + 'specified by the projection. Ensure you have specified ' + 'the correct coordinate system for your data.'.format(type(element).__name__, type(element.crs).__name__, + type(self.p.projection).__name__)) + + return element.clone(tuple(new_data[d.name] for d in element.dimensions()), + crs=self.p.projection, datatype=datatype) class project_graph(_project_operation): @@ -445,7 +479,7 @@ class project(Operation): _operations = [project_path, project_image, project_shape, project_graph, project_quadmesh, project_points, - project_geom] + project_vectorfield, project_geom] def _process(self, element, key=None): for op in self._operations: diff --git a/geoviews/plotting/bokeh/__init__.py b/geoviews/plotting/bokeh/__init__.py index 9654275f..1da17235 100644 --- a/geoviews/plotting/bokeh/__init__.py +++ b/geoviews/plotting/bokeh/__init__.py @@ -23,7 +23,7 @@ ) from ...operation import ( project_image, project_points, project_path, project_graph, - project_quadmesh, project_geom + project_quadmesh, project_geom, project_vectorfield ) from ...tile_sources import _ATTRIBUTIONS from ...util import poly_types, line_types @@ -106,7 +106,7 @@ class GeoPointPlot(GeoPlot, PointPlot): class GeoVectorFieldPlot(GeoPlot, VectorFieldPlot): - _project_operation = project_points + _project_operation = project_vectorfield class GeoQuadMeshPlot(GeoPlot, QuadMeshPlot): diff --git a/geoviews/plotting/mpl/__init__.py b/geoviews/plotting/mpl/__init__.py index b2144c7b..b61e4681 100644 --- a/geoviews/plotting/mpl/__init__.py +++ b/geoviews/plotting/mpl/__init__.py @@ -37,7 +37,7 @@ from ...operation import ( project_points, project_path, project_graph, project_quadmesh, - project_geom + project_geom, project_vectorfield ) from .chart import WindBarbsPlot @@ -325,7 +325,7 @@ class GeoVectorFieldPlot(GeoPlot, VectorFieldPlot): apply_ranges = param.Boolean(default=True) - _project_operation = project_points + _project_operation = project_vectorfield class GeoWindBarbsPlot(GeoPlot, WindBarbsPlot): From af15bd194bd70de4c5e6b44f68e1620c6da7010c Mon Sep 17 00:00:00 2001 From: Philipp Rudiger Date: Sun, 10 Feb 2019 17:28:27 +0000 Subject: [PATCH 2/5] fix typo in angle calculation (#297) --- geoviews/operation/projection.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/geoviews/operation/projection.py b/geoviews/operation/projection.py index 9fd9f14f..2e34d842 100644 --- a/geoviews/operation/projection.py +++ b/geoviews/operation/projection.py @@ -236,8 +236,8 @@ def _process_element(self, element): new_data[xdim.name] = coordinates[mask, 0] new_data[ydim.name] = coordinates[mask, 1] datatype = [element.interface.datatype]+element.datatype - us = np.sin(ang) * ms - vs = np.cos(ang) * ms + us = np.cos(ang) * ms + vs = np.sin(ang) * ms ut, vt = self.p.projection.transform_vectors(element.crs, xs, ys, us, vs) with np.errstate(divide='ignore', invalid='ignore'): angle = np.arctan2(vt, ut) From 861753fbc08c2bf8414a1dd4178cca590b32521d Mon Sep 17 00:00:00 2001 From: Andrew Huang Date: Thu, 6 Jul 2023 12:23:57 -0400 Subject: [PATCH 3/5] Add barbs --- geoviews/operation/__init__.py | 3 ++- geoviews/operation/projection.py | 23 +++++++++++++++++------ geoviews/plotting/mpl/__init__.py | 4 ++-- 3 files changed, 21 insertions(+), 9 deletions(-) diff --git a/geoviews/operation/__init__.py b/geoviews/operation/__init__.py index 23c7f153..8a5184a3 100644 --- a/geoviews/operation/__init__.py +++ b/geoviews/operation/__init__.py @@ -6,7 +6,8 @@ from ..element import _Element from .projection import ( # noqa (API import) project_image, project_path, project_shape, project_points, - project_graph, project_quadmesh, project_geom, project_vectorfield, project) + project_graph, project_quadmesh, project_geom, + project_vectorfield, project_windbarbs, project) from .resample import resample_geometry # noqa (API import) geo_ops = [contours, bivariate_kde] diff --git a/geoviews/operation/projection.py b/geoviews/operation/projection.py index 2e34d842..05f00c30 100644 --- a/geoviews/operation/projection.py +++ b/geoviews/operation/projection.py @@ -15,7 +15,7 @@ from ..data import GeoPandasInterface from ..element import (Image, Shape, Polygons, Path, Points, Contours, RGB, Graph, Nodes, EdgePaths, QuadMesh, VectorField, - HexTiles, Labels, Rectangles, Segments) + HexTiles, Labels, Rectangles, Segments, WindBarbs) from ..util import ( project_extents, path_to_geom_dicts, polygons_to_geom_dicts, geom_dict_to_array_dict @@ -224,6 +224,9 @@ class project_vectorfield(_project_operation): supported_types = [VectorField] + def _calc_angles(self, ut, vt): + return np.pi / 2 - np.arctan2(ut, vt) + def _process_element(self, element): if not len(element): return element.clone(crs=self.p.projection) @@ -236,12 +239,12 @@ def _process_element(self, element): new_data[xdim.name] = coordinates[mask, 0] new_data[ydim.name] = coordinates[mask, 1] datatype = [element.interface.datatype]+element.datatype - us = np.cos(ang) * ms - vs = np.sin(ang) * ms + us = np.sin(ang) * -ms + vs = np.cos(ang) * -ms ut, vt = self.p.projection.transform_vectors(element.crs, xs, ys, us, vs) with np.errstate(divide='ignore', invalid='ignore'): - angle = np.arctan2(vt, ut) - mag = np.sqrt(ut**2+vt**2) + angle = self._calc_angles(ut, vt) + mag = np.hypot(ut, vt) new_data[adim.name] = angle[mask] new_data[mdim.name] = mag[mask] @@ -257,6 +260,14 @@ def _process_element(self, element): return element.clone(tuple(new_data[d.name] for d in element.dimensions()), crs=self.p.projection, datatype=datatype) +class project_windbarbs(project_vectorfield): + + supported_types = [WindBarbs] + + def _calc_angles(self, ut, vt): + return np.pi / 2 - np.arctan2(-vt, -ut) + + class project_graph(_project_operation): supported_types = [Graph] @@ -479,7 +490,7 @@ class project(Operation): _operations = [project_path, project_image, project_shape, project_graph, project_quadmesh, project_points, - project_vectorfield, project_geom] + project_vectorfield, project_windbarbs, project_geom] def _process(self, element, key=None): for op in self._operations: diff --git a/geoviews/plotting/mpl/__init__.py b/geoviews/plotting/mpl/__init__.py index b61e4681..2d2a0387 100644 --- a/geoviews/plotting/mpl/__init__.py +++ b/geoviews/plotting/mpl/__init__.py @@ -37,7 +37,7 @@ from ...operation import ( project_points, project_path, project_graph, project_quadmesh, - project_geom, project_vectorfield + project_geom, project_vectorfield, project_windbarbs ) from .chart import WindBarbsPlot @@ -335,7 +335,7 @@ class GeoWindBarbsPlot(GeoPlot, WindBarbsPlot): apply_ranges = param.Boolean(default=True) - _project_operation = project_points + _project_operation = project_windbarbs class GeometryPlot(GeoPlot): From 86327930f8190108e1a47cc62657db0a9590c26a Mon Sep 17 00:00:00 2001 From: Andrew Huang Date: Thu, 6 Jul 2023 13:33:54 -0400 Subject: [PATCH 4/5] Add tests --- geoviews/operation/projection.py | 4 ++- geoviews/tests/test_projection.py | 44 ++++++++++++++++++++++++++++++- 2 files changed, 46 insertions(+), 2 deletions(-) diff --git a/geoviews/operation/projection.py b/geoviews/operation/projection.py index 05f00c30..3c143b11 100644 --- a/geoviews/operation/projection.py +++ b/geoviews/operation/projection.py @@ -225,7 +225,8 @@ class project_vectorfield(_project_operation): supported_types = [VectorField] def _calc_angles(self, ut, vt): - return np.pi / 2 - np.arctan2(ut, vt) + # mathmetical convention + return np.arctan2(vt, ut) def _process_element(self, element): if not len(element): @@ -265,6 +266,7 @@ class project_windbarbs(project_vectorfield): supported_types = [WindBarbs] def _calc_angles(self, ut, vt): + # meteorological convention return np.pi / 2 - np.arctan2(-vt, -ut) diff --git a/geoviews/tests/test_projection.py b/geoviews/tests/test_projection.py index b6e1dbdd..c739bb43 100644 --- a/geoviews/tests/test_projection.py +++ b/geoviews/tests/test_projection.py @@ -1,7 +1,7 @@ import numpy as np import cartopy.crs as ccrs -from geoviews.element import Image +from geoviews.element import Image, VectorField, WindBarbs from geoviews.element.comparison import ComparisonTestCase from geoviews.operation import project @@ -31,3 +31,45 @@ def test_image_project_latlon_to_mercator(self): [ 0., 0., 0., 0., 0.], [ 12960., 17280., 21600., 4320., 8640.] ])) + + def test_project_vectorfield(self): + xs = np.linspace(10, 50, 2) + X, Y = np.meshgrid(xs, xs) + U, V = 5 * X, 1 * Y + A = np.arctan2(V, U) + M = np.hypot(U, V) + crs = ccrs.PlateCarree() + vectorfield = VectorField((X, Y, A, M), crs=crs) + projection = ccrs.Orthographic() + projected = project(vectorfield, projection=projection) + assert projected.crs == projection + + xs, ys, ang, ms = (vectorfield.dimension_values(i) for i in range(4)) + us = np.sin(ang) * -ms + vs = np.cos(ang) * -ms + u, v = projection.transform_vectors(crs, xs, ys, us, vs) + a, m = np.arctan2(v, u).T, np.hypot(u, v).T + + np.testing.assert_allclose(projected.dimension_values("Angle"), a.flatten()) + np.testing.assert_allclose(projected.dimension_values("Magnitude"), m.flatten()) + + def test_project_windbarbs(self): + xs = np.linspace(10, 50, 2) + X, Y = np.meshgrid(xs, xs) + U, V = 5 * X, 1 * Y + A = np.arctan2(V, U) + M = np.hypot(U, V) + crs = ccrs.PlateCarree() + windbarbs = WindBarbs((X, Y, A, M), crs=crs) + projection = ccrs.Orthographic() + projected = project(windbarbs, projection=projection) + assert projected.crs == projection + + xs, ys, ang, ms = (windbarbs.dimension_values(i) for i in range(4)) + us = np.sin(ang) * -ms + vs = np.cos(ang) * -ms + u, v = projection.transform_vectors(crs, xs, ys, us, vs) + a, m = np.pi / 2 - np.arctan2(-v, -u).T, np.hypot(u, v).T + + np.testing.assert_allclose(projected.dimension_values("Angle"), a.flatten()) + np.testing.assert_allclose(projected.dimension_values("Magnitude"), m.flatten()) From e362a9f4e9bfc14755a2b97e23053c18011235e7 Mon Sep 17 00:00:00 2001 From: Andrew Huang Date: Thu, 10 Aug 2023 09:13:09 -0400 Subject: [PATCH 5/5] Address comments --- geoviews/operation/projection.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/geoviews/operation/projection.py b/geoviews/operation/projection.py index 3c143b11..c24d91df 100644 --- a/geoviews/operation/projection.py +++ b/geoviews/operation/projection.py @@ -225,7 +225,7 @@ class project_vectorfield(_project_operation): supported_types = [VectorField] def _calc_angles(self, ut, vt): - # mathmetical convention + # mathematical convention; follows matplotlib return np.arctan2(vt, ut) def _process_element(self, element): @@ -266,7 +266,7 @@ class project_windbarbs(project_vectorfield): supported_types = [WindBarbs] def _calc_angles(self, ut, vt): - # meteorological convention + # meteorological convention; follows matplotlib return np.pi / 2 - np.arctan2(-vt, -ut)