Skip to content

Commit a2a04f0

Browse files
pp-mopp-mo
committed
Split cube region extraction to separate bonus notebook.
1 parent 7273179 commit a2a04f0

File tree

4 files changed

+1244
-869
lines changed

4 files changed

+1244
-869
lines changed

‎notebooks/.notebook_shadow_copies/MeshCube_Extraction.md‎

Lines changed: 213 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,213 @@
1+
---
2+
jupyter:
3+
jupytext:
4+
text_representation:
5+
extension: .md
6+
format_name: markdown
7+
format_version: '1.3'
8+
jupytext_version: 1.14.4
9+
kernelspec:
10+
display_name: Python 3 (ipykernel)
11+
language: python
12+
name: python3
13+
---
14+
15+
# Extract a regional mesh-cube with Iris
16+
17+
While GeoVista provides the efficient tools for mesh region extraction, it and Iris know nothing about one another.
18+
So, to calculate a regionally-extracted _Iris cube_, GeoVista can do the hard work of determining the subset of cells required, but you must then "reconstruct" an Iris cube from that information.
19+
For now, at least, there are no ready-made tools for this (either in Iris or Geovista).
20+
21+
The process requires a few steps, which we can summarise as :
22+
1. record, on the original global PolyData, the original face indices of each of the cells
23+
1. perform extraction (by BBox or otherwise) to get a regional PolyData
24+
1. get the face-indices of the selected cells from the regional PolyData
25+
1. index the Iris cube with the selected indices, on the mesh dimension, to extract the regional parts
26+
1. construct and attach a suitable Iris mesh to represent the extracted region
27+
28+
( Note: the last step itself is not strictly necessary. It may be sufficent to have a regional data cube with a notional "mesh dimension", but which does not possess an actual Iris mesh. )
29+
30+
---
31+
32+
**First, we just repeat some of the imports + code from 'Sec_04_Region_Extraction'**
33+
34+
To get a global mesh-cube and a PolyData derived from it.
35+
36+
```python
37+
# Fetch a single cube of test data
38+
from testdata_fetching import lfric_rh_singletime_2d, switch_data
39+
switch_data(use_newer_smaller_c48_data=True)
40+
lfric_rh = lfric_rh_singletime_2d()
41+
lfric_rh
42+
```
43+
44+
```python
45+
# Convert to a PyVista.PolyData
46+
from pv_conversions import pv_from_lfric_cube
47+
pv_global_rh = pv_from_lfric_cube(lfric_rh)
48+
```
49+
50+
```python
51+
# Define a lat-lon region to extract -- EXACTLY AS BEFORE
52+
from geovista.geodesic import BBox
53+
bbox = BBox(lons=[0, 70, 70, 0], lats=[-25, -25, 45, 45])
54+
```
55+
56+
```python
57+
# 'Apply' the region to the PolyData object.
58+
pv_regional_rh = bbox.enclosed(pv_global_rh)
59+
pv_regional_rh
60+
```
61+
62+
<!-- #region -->
63+
---
64+
## Now the meshcube extraction, as a sequence of steps ...
65+
66+
67+
### Step 1 : Add an auxiliary array to our _global_ PolyData
68+
This is to record, for each PolyData cell, the original (face) index which that cell has.
69+
70+
Note : we use numpy.arange() to construct a counting sequence, and make this a new array on the PolyData object, by assigning to a index-name.
71+
<!-- #endregion -->
72+
73+
```python
74+
import numpy as np
75+
face_inds = np.arange(pv_global_rh.n_cells)
76+
# Assign PolyData[<a-string>] to create a new attached array.
77+
pv_global_rh.cell_data['original_face_indices'] = face_inds
78+
```
79+
80+
---
81+
82+
### Step 2 : Extract with your Bbox to get a regional PolyData, and show the result.
83+
This code is exactly the same as the previous time we did this.
84+
85+
```python
86+
pv_regional_rh = bbox.enclosed(pv_global_rh)
87+
pv_regional_rh
88+
```
89+
90+
You can see that the new version of the extracted (regional) data now has an ***extra*** attached data array, "`original_face_indices`", which is derived from the one we added to the global data, and which holds the face indices of the _selected_ cells.
91+
92+
---
93+
94+
**Step 3 : Fetch the indices array from the regional PolyData, by indexing with the array name.**
95+
and show the result.
96+
97+
```python
98+
# Get the remaining face indices, to use for indexing the Cube.
99+
region_indices = pv_regional_rh["original_face_indices"]
100+
region_indices
101+
```
102+
103+
This contains the original face-indices of all the cells which fall within the region, _i.e. which faces those were in the global mesh_.
104+
105+
We can now apply these via indexing, to select only those cells *from the Iris cube*.
106+
107+
**Step 4 : Apply these cells as an index to the 'mesh dimension' of the original Iris lfric-rh cube**
108+
.. and print that out.
109+
110+
```python
111+
lfric_rh_region = lfric_rh[..., region_indices]
112+
lfric_rh_region
113+
```
114+
115+
This new cube contains the mesh data within our selected region.
116+
117+
However, there is a catch here : Once indexed, our cube ***no longer has a mesh***.
118+
You can see this in the printout, which lists "Auxiliary coordinates" but no "Mesh coordinates".
119+
120+
This problem will probably be fixed in future. See [here in the Iris docs](https://scitools-iris.readthedocs.io/en/latest/further_topics/ugrid/operations.html#region-extraction) for a discussion.
121+
122+
For now, what we need to do is to re-create a mesh for the regional cube.
123+
We do that in a few further steps ...
124+
125+
---
126+
127+
**Step 5a : Get the X and Y-axis coordinates from the region cube.**
128+
Use `Cube.coords('longitude')` etc.
129+
130+
```python
131+
x_coord = lfric_rh_region.coord('longitude')
132+
y_coord = lfric_rh_region.coord('latitude')
133+
```
134+
135+
**Step 5b : Create a new `iris.experimental.ugrid.Mesh`-class object, passing the X,Y coords as arguments**
136+
137+
```python
138+
from iris.experimental.ugrid.mesh import Mesh
139+
mesh = Mesh.from_coords(x_coord, y_coord)
140+
```
141+
142+
( Step 2a : **`print()` the Mesh object**
143+
Note : `Mesh` does not provide a notebook display method.
144+
)
145+
146+
```python
147+
print(mesh)
148+
```
149+
150+
---
151+
**Step 5c : Call `Mesh.to_MeshCoords` to create a pair of `MeshCoord`s containing this mesh**
152+
Note : you must specify the keyword `location="face"` : This matches the data location of the original data -- i.e. the data cells are faces.
153+
154+
```python
155+
mesh_coords = mesh.to_MeshCoords(location="face")
156+
mesh_coords
157+
```
158+
159+
---
160+
**Step 5d : (finally!!)
161+
Use `Cube.remove_coord` and `Cube.add_aux_coord` to replace each AuxCoord with its corresponding `MeshCoord` from the previous step.** Note : for 'add_aux_coord', you also need to specify the relevant cube dimension(s) : See [`Cube.add_aux_coord` in the Iris docs](https://scitools-iris.readthedocs.io/en/latest/generated/api/iris/cube.html?highlight=add_aux_coord#iris.cube.Cube.add_aux_coord)
162+
.. and show the cube ..
163+
164+
```python
165+
lfric_rh_region.remove_coord('longitude')
166+
```
167+
168+
```python
169+
lfric_rh_region.remove_coord('latitude')
170+
```
171+
172+
```python
173+
xco, yco = mesh_coords
174+
175+
lfric_rh_region.add_aux_coord(xco, 0)
176+
lfric_rh_region.add_aux_coord(yco, 0)
177+
178+
# Result : a regional Mesh-Cube with a subset of the original faces.
179+
lfric_rh_region
180+
```
181+
182+
---
183+
184+
**Lastly, plot this to see what we got.**
185+
Use the techniques as above, converting with `pv_from_lfric_cube` and plotting.
186+
187+
188+
```python
189+
pv = pv_from_lfric_cube(lfric_rh_region)
190+
pv.plot()
191+
```
192+
193+
----
194+
195+
**Investigation:** It is useful to add some extra background information to make this more visible.
196+
197+
As a minimum you can use `plotter.add_coastlines()`.
198+
199+
Another useful one is `plotter.add_base_layer()`
200+
**Question : what does that actually do ?**
201+
202+
```python
203+
from geovista import GeoPlotter
204+
plotter = GeoPlotter()
205+
plotter.add_mesh(pv)
206+
plotter.add_coastlines()
207+
plotter.add_base_layer()
208+
plotter.show()
209+
```
210+
211+
```python
212+
213+
```

‎notebooks/.notebook_shadow_copies/Sec_05_RegionExtraction.md‎

Lines changed: 5 additions & 159 deletions
Original file line numberDiff line numberDiff line change
@@ -106,177 +106,23 @@ plotter.show()
106106
```
107107

108108
```python
109-
# Temporary : show static plot for notebook review
109+
# Additional : static plot for notebook review
110110
plotter.show(jupyter_backend='static')
111111
```
112112

113-
<!-- #region -->
114113
---
115114
## Get an Iris cube for an extracted region.
116115

117116
While GeoVista provides the efficient tools for mesh region extraction, it and Iris know nothing about one another.
118-
So, to calculate a regionally-extracted _Iris cube_, GeoVista can do the hard work of determining the subset of cells required, but you must then "reconstruct" an Iris cube from that information.
119-
For now, at least, there are no ready-made tools for this (either in Iris or Geovista).
117+
So, to calculate a regionally-extracted _Iris cube_, GeoVista can do the hard work of determining the subset of cells required, but you must then "reconstruct" an Iris cube from that information.
120118

121-
So, this task is rather involved at present, but for those interested a working code example is provided as bonus content in a [separate notebook (broken link!)](to-link-here).
119+
For now, at least, there are no ready-made tools for this (in either Iris or Geovista).
122120

121+
However, the operation is possible, and may be instructive.
122+
So, for those interested, there is an extra notebook ["MeshCube_Extraction.ipynb"](./MeshCube_Extraction.ipynb), showing how this is done.
123123

124-
## DEVELOPER-TODO : relegate all this to a separate "bonus content" notebook
125124

126125

127-
While GeoVista provides the efficient tools for mesh region extraction, it and Iris know nothing about one another.
128-
So, to calculate a regionally-extracted _Iris cube_, GeoVista can do the hard work of determining the subset of cells required, but you must then "reconstruct" an Iris cube from that information.
129-
For now, at least, there are no ready-made tools for this (either in Iris or Geovista).
130-
131-
The process requires a few steps, which we can summarise as :
132-
1. record, on the original global PolyData, the original face indices of each of the cells
133-
1. perform extraction (by BBox or otherwise) to get a regional PolyData
134-
1. get the face-indices of the selected cells from the regional PolyData
135-
1. index the Iris cube with the selected indices, on the mesh dimension, to extract the regional parts
136-
1. construct and attach a suitable Iris mesh to represent the extracted region
137-
138-
( Note: the last step itself is not strictly necessary. It may be sufficent to have a regional data cube with a notional "mesh dimension", but which does not possess an actual Iris mesh. )
139-
140-
---
141-
Let's show that operation ...
142-
143-
**Step 1 : First, add an auxiliary array to the global PolyData, recording the original (face) index of each cell.**
144-
Note : use numpy.arange() to construct a counting sequence, and assign to a named index on the PolyData object.
145-
<!-- #endregion -->
146-
147-
```python
148-
import numpy as np
149-
face_inds = np.arange(pv_global_rh.n_cells)
150-
pv_global_rh.cell_data['original_face_indices'] = face_inds
151-
```
152-
153-
---
154-
155-
**Step 2 : Extract with your Bbox to get a regional PolyData, and show the result.**
156-
This code is exactly the same as the previous time we did this.
157-
158-
```python
159-
pv_regional_rh = bbox.enclosed(pv_global_rh)
160-
pv_regional_rh
161-
```
162-
163-
You can see that the new version of the extracted (regional) data now has an ***extra*** attached data array, derived from the one we added to the global data, and which holds the selected face indices.
164-
165-
---
166-
167-
**Step 3 : Fetch the indices array from the regional PolyData, by indexing with the array name.**
168-
and show the result.
169-
170-
```python
171-
# Get the remaining face indices, to use for indexing the Cube.
172-
region_indices = pv_regional_rh["original_face_indices"]
173-
region_indices
174-
```
175-
176-
This contains the original face-indices of all the cells which fall within the region, _i.e. which faces those were in the global mesh_.
177-
178-
We can now apply these indices, to select only those cells *from the Iris cube*.
179-
180-
**Step 4 : Apply these cells as an index to the 'mesh dimension' of the original Iris lfric-rh cube**
181-
.. and print that out.
182-
183-
```python
184-
lfric_rh_region = lfric_rh[..., region_indices]
185-
lfric_rh_region
186-
```
187-
188-
This new cube contains the mesh data within our selected region.
189-
190-
However, there is a catch here : Once indexed, our cube ***no longer has a mesh***.
191-
You can see this in the printout, which lists "Auxiliary coordinates" but no "Mesh coordinates".
192-
193-
This problem will probably be fixed in future. See [here in the Iris docs](https://scitools-iris.readthedocs.io/en/latest/further_topics/ugrid/operations.html#region-extraction) for a discussion.
194-
195-
For now, what we need to do is to re-create a mesh for the regional cube.
196-
We do that in a few further steps ...
197-
198-
---
199-
200-
**Step 5a : Get the X and Y-axis coordinates from the region cube.**
201-
Use `Cube.coords('longitude')` etc.
202-
203-
```python
204-
x_coord = lfric_rh_region.coord('longitude')
205-
y_coord = lfric_rh_region.coord('latitude')
206-
```
207-
208-
**Step 5b : Create a new `iris.experimental.ugrid.Mesh`-class object, passing the X,Y coords as arguments**
209-
210-
```python
211-
from iris.experimental.ugrid.mesh import Mesh
212-
mesh = Mesh.from_coords(x_coord, y_coord)
213-
```
214-
215-
( Step 2a : **`print()` the Mesh object**
216-
Note : `Mesh` does not provide a notebook display method.
217-
)
218-
219-
```python
220-
print(mesh)
221-
```
222-
223-
---
224-
**Step 5c : Call `Mesh.to_MeshCoords` to create a pair of `MeshCoord`s containing this mesh**
225-
Note : you must specify the keyword `location="face"` : This matches the data location of the original data -- i.e. the data cells are faces.
226-
227-
```python
228-
mesh_coords = mesh.to_MeshCoords(location="face")
229-
mesh_coords
230-
```
231-
232-
---
233-
**Step 5d : (finally!!)
234-
Use `Cube.remove_coord` and `Cube.add_aux_coord` to replace each AuxCoord with its corresponding `MeshCoord` from the previous step.** Note : for 'add_aux_coord', you also need to specify the relevant cube dimension(s) : See [`Cube.add_aux_coord` in the Iris docs](https://scitools-iris.readthedocs.io/en/latest/generated/api/iris/cube.html?highlight=add_aux_coord#iris.cube.Cube.add_aux_coord)
235-
.. and show the cube ..
236-
237-
```python
238-
lfric_rh_region.remove_coord('longitude')
239-
```
240-
241-
```python
242-
lfric_rh_region.remove_coord('latitude')
243-
```
244-
245-
```python
246-
xco, yco = mesh_coords
247-
248-
lfric_rh_region.add_aux_coord(xco, 0)
249-
lfric_rh_region.add_aux_coord(yco, 0)
250-
251-
# Result : a regional Mesh-Cube with a subset of the original faces.
252-
lfric_rh_region
253-
```
254-
255-
---
256-
257-
**Lastly, plot this to see what we got.**
258-
Use the techniques as above, converting with `pv_from_lfric_cube` and plotting.
259-
260-
261-
```python
262-
pv = pv_from_lfric_cube(lfric_rh_region)
263-
pv.plot()
264-
```
265-
266-
----
267-
268-
**Investigation:** It is useful to add some extra background information to make this more visible.
269-
270-
As a minimum you can use `plotter.add_coastlines()`.
271-
272-
Another useful one is `plotter.add_base_layer()`
273-
**Question : what does that actually do ?**
274-
275-
```python
276-
plotter.add_coastlines()
277-
plotter.add_base_layer()
278-
plotter.show()
279-
```
280126

281127
```python
282128

0 commit comments

Comments
 (0)