Skip to content

Commit 6be7f96

Browse files
HGWrightHGWright
authored
Merge 37c4bb5 into fb039ab 
2 parents fb039ab + 37c4bb5 commit 6be7f96

File tree

19 files changed

+1235
-242
lines changed

19 files changed

+1235
-242
lines changed

‎docs/src/conf.py‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -252,6 +252,8 @@ def _dotv(version):
252252
"scipy": ("https://docs.scipy.org/doc/scipy/", None),
253253
"pandas": ("https://pandas.pydata.org/docs/", None),
254254
"dask": ("https://docs.dask.org/en/stable/", None),
255+
"geovista": ("https://geovista.readthedocs.io/en/latest/", None),
256+
"pyvista": ("https://docs.pyvista.org/", None),
255257
}
256258

257259
# The name of the Pygments (syntax highlighting) style to use.

‎docs/src/further_topics/ugrid/images/plotting.png‎

31.6 KB
Loading

‎docs/src/further_topics/ugrid/images/plotting_basic.png‎

-41.3 KB
Binary file not shown.

‎docs/src/further_topics/ugrid/images/plotting_global.png‎

-130 KB
Binary file not shown.

‎docs/src/further_topics/ugrid/operations.rst‎

Lines changed: 79 additions & 216 deletions
Large diffs are not rendered by default.

‎lib/iris/experimental/geovista.py‎

Lines changed: 339 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,339 @@
1+
# Copyright Iris contributors
2+
#
3+
# This file is part of Iris and is released under the BSD license.
4+
# See LICENSE in the root of the repository for full licensing details.
5+
"""Experimental module for using some GeoVista operations with Iris cubes."""
6+
7+
from geovista import Transform
8+
from geovista.common import VTK_CELL_IDS, VTK_POINT_IDS
9+
10+
from iris.exceptions import CoordinateNotFoundError
11+
from iris.experimental.ugrid import Mesh
12+
13+
14+
def _get_coord(cube, axis):
15+
"""Get the axis coordinates from the cube."""
16+
try:
17+
coord = cube.coord(axis=axis, dim_coords=True)
18+
except CoordinateNotFoundError:
19+
coord = cube.coord(axis=axis)
20+
return coord
21+
22+
23+
def cube_to_polydata(cube, **kwargs):
24+
r"""Create a :class:`pyvista.PolyData` object from a :class:`~iris.cube.Cube`.
25+
26+
The resulting :class:`~pyvista.PolyData` object can be plotted using
27+
a :class:`geovista.geoplotter.GeoPlotter`.
28+
29+
Uses :class:`geovista.bridge.Transform` to parse the cube's information - one
30+
of: :meth:`~geovista.bridge.Transform.from_1d` /
31+
:meth:`~geovista.bridge.Transform.from_2d` /
32+
:meth:`~geovista.bridge.Transform.from_unstructured`.
33+
34+
Parameters
35+
----------
36+
cube : :class:`~iris.cube.Cube`
37+
The Cube containing the spatial information and data for creating the
38+
class:`~pyvista.PolyData`.
39+
40+
**kwargs : dict, optional
41+
Additional keyword arguments to be passed to the relevant
42+
:class:`~geovista.bridge.Transform` method (e.g ``zlevel``).
43+
44+
Returns
45+
-------
46+
:class:`~pyvista.PolyData`
47+
The PolyData object representing the cube's spatial information and data.
48+
49+
Raises
50+
------
51+
NotImplementedError
52+
If a :class:`~iris.cube.Cube` with too many dimensions is passed. Only
53+
the horizontal data can be represented, meaning a 2D Cube, or 1D Cube
54+
if the horizontal space is described by
55+
:class:`~iris.experimental.ugrid.MeshCoord`\ s.
56+
57+
Examples
58+
--------
59+
.. testsetup::
60+
61+
from iris import load_cube, sample_data_path
62+
from iris.experimental.ugrid import PARSE_UGRID_ON_LOAD
63+
64+
cube = load_cube(sample_data_path("air_temp.pp"))
65+
cube_w_time = load_cube(sample_data_path("A1B_north_america.nc"))
66+
with PARSE_UGRID_ON_LOAD.context():
67+
cube_mesh = load_cube(sample_data_path("mesh_C4_synthetic_float.nc"))
68+
69+
>>> from iris.experimental.geovista import cube_to_polydata
70+
71+
Converting a standard 2-dimensional :class:`~iris.cube.Cube` with
72+
1-dimensional coordinates:
73+
74+
>>> print(cube.summary(shorten=True))
75+
air_temperature / (K) (latitude: 73; longitude: 96)
76+
>>> print(cube_to_polydata(cube))
77+
PolyData (...
78+
N Cells: 7008
79+
N Points: 7178
80+
N Strips: 0
81+
X Bounds: -9.992e-01, 9.992e-01
82+
Y Bounds: -9.992e-01, 9.992e-01
83+
Z Bounds: -1.000e+00, 1.000e+00
84+
N Arrays: 4
85+
86+
Configure the conversion by passing additional keyword arguments:
87+
88+
>>> print(cube_to_polydata(cube, radius=2))
89+
PolyData (...
90+
N Cells: 7008
91+
N Points: 7178
92+
N Strips: 0
93+
X Bounds: -1.998e+00, 1.998e+00
94+
Y Bounds: -1.998e+00, 1.998e+00
95+
Z Bounds: -2.000e+00, 2.000e+00
96+
N Arrays: 4
97+
98+
Converting a :class:`~iris.cube.Cube` that has a
99+
:attr:`~iris.cube.Cube.mesh` describing its horizontal space:
100+
101+
>>> print(cube_mesh.summary(shorten=True))
102+
synthetic / (1) (-- : 96)
103+
>>> print(cube_to_polydata(cube_mesh))
104+
PolyData (...
105+
N Cells: 96
106+
N Points: 98
107+
N Strips: 0
108+
X Bounds: -1.000e+00, 1.000e+00
109+
Y Bounds: -1.000e+00, 1.000e+00
110+
Z Bounds: -1.000e+00, 1.000e+00
111+
N Arrays: 4
112+
113+
Remember to reduce the dimensionality of your :class:`~iris.cube.Cube` to
114+
just be the horizontal space:
115+
116+
>>> print(cube_w_time.summary(shorten=True))
117+
air_temperature / (K) (time: 240; latitude: 37; longitude: 49)
118+
>>> print(cube_to_polydata(cube_w_time[0, :, :]))
119+
PolyData (...
120+
N Cells: 1813
121+
N Points: 1900
122+
N Strips: 0
123+
X Bounds: -6.961e-01, 6.961e-01
124+
Y Bounds: -9.686e-01, -3.411e-01
125+
Z Bounds: 2.483e-01, 8.714e-01
126+
N Arrays: 4
127+
128+
"""
129+
if cube.mesh:
130+
if cube.ndim != 1:
131+
raise NotImplementedError("Cubes with a mesh must be one dimensional")
132+
lons, lats = cube.mesh.node_coords
133+
face_node = cube.mesh.face_node_connectivity
134+
indices = face_node.indices_by_location()
135+
136+
polydata = Transform.from_unstructured(
137+
xs=lons.points,
138+
ys=lats.points,
139+
connectivity=indices,
140+
data=cube.data,
141+
name=f"{cube.name()} / ({cube.units})",
142+
start_index=face_node.start_index,
143+
**kwargs,
144+
)
145+
# TODO: Add support for point clouds
146+
elif cube.ndim == 2:
147+
x_coord = _get_coord(cube, "X")
148+
y_coord = _get_coord(cube, "Y")
149+
transform_kwargs = dict(
150+
xs=x_coord.contiguous_bounds(),
151+
ys=y_coord.contiguous_bounds(),
152+
data=cube.data,
153+
name=f"{cube.name()} / ({cube.units})",
154+
**kwargs,
155+
)
156+
coord_system = cube.coord_system()
157+
if coord_system:
158+
transform_kwargs["crs"] = coord_system.as_cartopy_crs().proj4_init
159+
160+
if x_coord.ndim == 2 and y_coord.ndim == 2:
161+
polydata = Transform.from_2d(**transform_kwargs)
162+
163+
elif x_coord.ndim == 1 and y_coord.ndim == 1:
164+
polydata = Transform.from_1d(**transform_kwargs)
165+
166+
else:
167+
raise NotImplementedError("Only 1D and 2D coordinates are supported")
168+
else:
169+
raise NotImplementedError("Cube must have a mesh or have 2 dimensions")
170+
171+
return polydata
172+
173+
174+
def extract_unstructured_region(cube, polydata, region, **kwargs):
175+
"""Index a :class:`~iris.cube.Cube` with a :attr:`~iris.cube.Cube.mesh` to a specific region.
176+
177+
Uses :meth:`geovista.geodesic.BBox.enclosed` to identify the `cube` indices
178+
that are within the specified region (`region` being a
179+
:class:`~geovista.geodesic.BBox` class).
180+
181+
Parameters
182+
----------
183+
cube : :class:`~iris.cube.Cube`
184+
The cube to be indexed (must have a :attr:`~iris.cube.Cube.mesh`).
185+
polydata : :class:`pyvista.PolyData`
186+
A :class:`~pyvista.PolyData` representing the same horizontal space as
187+
`cube`. The region extraction is first applied to `polydata`, with the
188+
resulting indices then applied to `cube`. In many cases `polydata` can
189+
be created by applying :func:`cube_to_polydata` to `cube`.
190+
region : :class:`geovista.geodesic.BBox`
191+
A :class:`~geovista.geodesic.BBox` representing the region to be
192+
extracted.
193+
**kwargs : dict, optional
194+
Additional keyword arguments to be passed to the
195+
:meth:`geovista.geodesic.BBox.enclosed` method (e.g ``preference``).
196+
197+
Returns
198+
-------
199+
:class:`~iris.cube.Cube`
200+
The region extracted cube.
201+
202+
Raises
203+
------
204+
ValueError
205+
If `polydata` and the :attr:`~iris.cube.Cube.mesh` on `cube` do not
206+
have the same shape.
207+
208+
Examples
209+
--------
210+
.. testsetup::
211+
212+
from iris import load_cube, sample_data_path
213+
from iris.coords import AuxCoord
214+
from iris.cube import CubeList
215+
from iris.experimental.ugrid import PARSE_UGRID_ON_LOAD
216+
217+
file_path = sample_data_path("mesh_C4_synthetic_float.nc")
218+
with PARSE_UGRID_ON_LOAD.context():
219+
cube_w_mesh = load_cube(file_path)
220+
221+
level_cubes = CubeList()
222+
for height_level in range(72):
223+
height_coord = AuxCoord([height_level], standard_name="height")
224+
level_cube = cube_w_mesh.copy()
225+
level_cube.add_aux_coord(height_coord)
226+
level_cubes.append(level_cube)
227+
228+
cube_w_mesh = level_cubes.merge_cube()
229+
other_cube_w_mesh = cube_w_mesh[:20, :]
230+
231+
The parameters of :func:`extract_unstructured_region` have been designed with
232+
flexibility and reuse in mind. This is demonstrated below.
233+
234+
>>> from geovista import BBox
235+
>>> from iris.experimental.geovista import cube_to_polydata, extract_unstructured_region
236+
>>> print(cube_w_mesh.shape)
237+
(72, 96)
238+
>>> # The mesh dimension represents the horizontal space of the cube.
239+
>>> print(cube_w_mesh.shape[cube_w_mesh.mesh_dim()])
240+
96
241+
>>> cube_polydata = cube_to_polydata(cube_w_mesh[0, :])
242+
>>> extracted_cube = extract_unstructured_region(
243+
... cube=cube_w_mesh,
244+
... polydata=cube_polydata,
245+
... region=BBox(lons=[0, 70, 70, 0], lats=[-25, -25, 45, 45]),
246+
... )
247+
>>> print(extracted_cube.shape)
248+
(72, 11)
249+
250+
Now reuse the same `cube` and `polydata` to extract a different region:
251+
252+
>>> new_region = BBox(lons=[0, 35, 35, 0], lats=[-25, -25, 45, 45])
253+
>>> extracted_cube = extract_unstructured_region(
254+
... cube=cube_w_mesh,
255+
... polydata=cube_polydata,
256+
... region=new_region,
257+
... )
258+
>>> print(extracted_cube.shape)
259+
(72, 6)
260+
261+
Now apply the same region extraction to a different `cube` that has the
262+
same horizontal shape:
263+
264+
>>> print(other_cube_w_mesh.shape)
265+
(20, 96)
266+
>>> extracted_cube = extract_unstructured_region(
267+
... cube=other_cube_w_mesh,
268+
... polydata=cube_polydata,
269+
... region=new_region,
270+
... )
271+
>>> print(extracted_cube.shape)
272+
(20, 6)
273+
274+
Arbitrary keywords can be passed down to
275+
:meth:`geovista.geodesic.BBox.enclosed` (``outside`` in this example):
276+
277+
>>> extracted_cube = extract_unstructured_region(
278+
... cube=other_cube_w_mesh,
279+
... polydata=cube_polydata,
280+
... region=new_region,
281+
... outside=True,
282+
... )
283+
>>> print(extracted_cube.shape)
284+
(20, 90)
285+
286+
"""
287+
if cube.mesh:
288+
# Find what dimension the mesh is in on the cube
289+
mesh_dim = cube.mesh_dim()
290+
recreate_mesh = False
291+
292+
if cube.location == "face":
293+
polydata_length = polydata.GetNumberOfCells()
294+
indices_key = VTK_CELL_IDS
295+
recreate_mesh = True
296+
elif cube.location == "node":
297+
polydata_length = polydata.GetNumberOfPoints()
298+
indices_key = VTK_POINT_IDS
299+
else:
300+
raise NotImplementedError(
301+
f"cube.location must be `face` or `node`. Found: {cube.location}."
302+
)
303+
304+
if cube.shape[mesh_dim] != polydata_length:
305+
raise ValueError(
306+
f"The mesh on the cube and the polydata"
307+
f"must have the same shape."
308+
f" Found Mesh: {cube.shape[mesh_dim]},"
309+
f" Polydata: {polydata_length}."
310+
)
311+
312+
region_polydata = region.enclosed(polydata, **kwargs)
313+
indices = region_polydata[indices_key]
314+
if len(indices) == 0:
315+
raise IndexError("No part of `polydata` falls within `region`.")
316+
317+
my_tuple = tuple(
318+
[slice(None) if i != mesh_dim else indices for i in range(cube.ndim)]
319+
)
320+
321+
region_cube = cube[my_tuple]
322+
323+
if recreate_mesh:
324+
coords_on_mesh_dim = region_cube.coords(dimensions=mesh_dim)
325+
new_mesh = Mesh.from_coords(
326+
*[c for c in coords_on_mesh_dim if c.has_bounds()]
327+
)
328+
329+
new_mesh_coords = new_mesh.to_MeshCoords(cube.location)
330+
331+
for coord in new_mesh_coords:
332+
region_cube.remove_coord(coord.name())
333+
region_cube.add_aux_coord(coord, mesh_dim)
334+
335+
# TODO: Support unstructured point based data without a mesh
336+
else:
337+
raise ValueError("Cube must have a mesh")
338+
339+
return region_cube

‎lib/iris/experimental/ugrid/mesh.py‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -863,6 +863,7 @@ def check_shape(array_name):
863863

864864
#####
865865
# TODO: remove axis assignment once Mesh supports arbitrary coords.
866+
# TODO: consider filtering coords as the first action in this method.
866867
axes_present = [guess_coord_axis(coord) for coord in coords]
867868
axes_required = ("X", "Y")
868869
if all([req in axes_present for req in axes_required]):

‎lib/iris/tests/integration/experimental/geovista/__init__.py‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,5 @@
1+
# Copyright Iris contributors
2+
#
3+
# This file is part of Iris and is released under the BSD license.
4+
# See LICENSE in the root of the repository for full licensing details.
5+
"""Integration tests for the :mod:`iris.experimental.geovista` module."""

0 commit comments

Comments
 (0)