Region Manifold Extraction

Region Manifold Extraction#

Introduction#

In this tutorial we use unstructured cubed-sphere sample data from the pantry to explore mesh regional extraction using a geodesic manifold.

Let’s Go!#

The Met Office is migrating to a unstructured cubed-sphere quadrilateral-cell mesh, which is a gridded cube projected onto a sphere i.e., there are 6 panels on the sphere, see Fig. 4:.

Each cubed-sphere is defined by the number of cells squared within each panel. In this tutorial we’ll use a C48 cubed-sphere, so there are 48 * 48 cells per panel, and 6 * 48 * 48 cells in a mesh overall.

First, let’s import geovista and generate some LFRic model cubed-sphere sample data using the lfric() and lfric_sst() convenience functions:

See also

For those interested in LFRic, please refer to LFRic - a modelling system fit for future computers.

import geovista as gv
from geovista.pantry.meshes import lfric, lfric_sst

As the lfric() and lfric_sst() functions both return a pyvista.PolyData we can simply use the pyvista.DataSet.plot method to render them:

c48 = lfric(resolution="c48")
c48.plot(show_edges=True)

Now create the same C48 cubed-sphere resolution mesh, but this time with cells populated with Sea Surface Temperature (SST) data. Note that the land masses are masked, as this is an oceanographic dataset. Masked cells within the mesh will have nan values:

c48_sst = lfric_sst()
c48_sst.plot(show_edges=True)

Our challenge is to extract all cubed-sphere cells within the panel covering the Americas.

We can achieve this using the geovista.geodesic.panel() function.

The panel() function is a convenience which allows you to extract cells from a mesh within 1 of the 6 named cubed-sphere panels i.e., africa, americas, antarctic, arctic, asia, and pacific.

Hint

You can use BBox directly to create generic geodesic bounding-box manifolds.

It returns a geovista.geodesic.BBox instance, which defines a geodesic bounding-box manifold for the requested spatial region.

So let’s go ahead and create the BBox for the americas cubed-sphere panel and use its mesh() method to visualise the manifold:

from geovista.geodesic import panel

americas_bbox = panel("americas")
americas_bbox.mesh.plot(color="orange")

Note that, the americas_bbox is constructed from geodesic lines i.e., great circles, and is a 3-D manifold i.e., a water-tight geometric structure. As such, we can then use it to extract points/cells from any underlying mesh.

Before doing that, first let’s render the americas_bbox and the c48_sst mesh together so that we can see their relationship to one another and convince ourselves that the americas_bbox instance is indeed covering the correct panel of the cubed-sphere:

Note

plotter is a GeoPlotter instance. It inherits all the behaviour of a pyvista.Plotter but also additional ✨ cartographic conveniences ✨ from GeoPlotterBase.

plotter = gv.GeoPlotter()

plotter.add_mesh(c48_sst, show_edges=True)
plotter.add_mesh(americas_bbox.mesh, color="orange")

plotter.view_yz()
plotter.camera.zoom(1.2)
plotter.add_axes()
plotter.show()

Seems to be right on the money! 👍

However, let’s apply some opacity to the c48_sst mesh so that we can see through the surface and view the americas_bbox from a different angle.

We’ll also use the view_poi() method to perform some cartographic camera controls on the rendered scene. Namey, move the camera to 30°E longitude.

plotter = gv.GeoPlotter()
plotter.add_mesh(c48_sst, opacity=0.3)
plotter.add_mesh(americas_bbox.mesh, color="orange")

plotter.view_poi(x=30)

plotter.camera.zoom(1.2)
plotter.add_axes()
plotter.show()

Rather than viewing the entire bounding-box, sometimes it’s more helpful to see only the boundary where the manifold intersects the surface of the mesh that it’s enclosing.

We can achieve this by using the geovista.geodesic.BBox.boundary() method:

plotter = gv.GeoPlotter()
plotter.add_mesh(c48_sst, show_edges=True)

plotter.add_mesh(americas_bbox.boundary(), color="orange", line_width=5)

plotter.view_xz()
plotter.camera.zoom(1.2)
plotter.add_axes()
plotter.show()

Okay, so let’s finally use the americas_bbox to extract cells from the c48_sst mesh by using the geovista.geodesic.BBox.enclosed() method:

region = americas_bbox.enclosed(c48_sst)

The extracted region returned by enclosed() is a pyvista.PolyData. Under the hood enclosed() uses the pyvista.DataSetFilters.select_enclosed_points() method to achieve this.

You may want to experiment with the preference kwarg of the enclosed() method to see the impact on the end result.

Anyways, let’s go ahead and see the extracted region, which should represent all the cells from the c48_sst mesh within the americas bounding-box:

plotter = gv.GeoPlotter()

plotter.add_mesh(region, show_edges=True)

plotter.view_xz()
plotter.add_axes()
plotter.show()

Hint

When you extract data from a mesh use the pyvista.Plotter.add_mesh() clim kwarg to ensure the same colorbar range.

Now bring this all together by rendering the extracted region along with a Natural Earth texture mapped base layer and coastlines:

clim = (271, 303)

plotter = gv.GeoPlotter()
plotter.add_mesh(region, show_edges=True, clim=clim)

plotter.add_base_layer(texture=gv.natural_earth_hypsometric())
plotter.add_coastlines(resolution="10m")

plotter.view_xz()
plotter.camera.zoom(1.2)
plotter.show_axes()
plotter.show()

Hint

Discover the data arrays available on a mesh with pyvista.DataSet.array_names and the active scalars with pyvista.DataSet.active_scalars_name.

As we’re not particularly interested in the masked land masses, we can easily remove them using the pyvista.DataSetFilters.threshold() method available on the region PolyData instance.

By default, the threshold() method will remove cells with nan values from the active scalars array on the mesh, which is just what we need:

sea_region = region.threshold()

Now let’s respin the render, but with the sea_region mesh and add some S.I. units to the scalar bar:

plotter = gv.GeoPlotter()

sargs = {"title": f"{sea_region.active_scalars_name} / K"}
plotter.add_mesh(sea_region, show_edges=True, clim=clim, scalar_bar_args=sargs)

plotter.add_base_layer(texture=gv.natural_earth_hypsometric())
plotter.add_coastlines(resolution="10m")
plotter.view_xz()
plotter.camera.zoom(1.2)
plotter.show_axes()
plotter.show()

Extension#

geovista offers basic support for Cylindrical and Pseudo-Cylindrical cartographic projections. We’re working on maturing this capability, and later extending to other classes of projections, such as Azimuthal and Conic.

To exercise this, let’s transform our sea_region mesh to a Robinson projection using a PROJ string to define the Coordinate Reference System (CRS):

crs = "+proj=robin lon_0=-90"

Note

If no crs is provided to GeoPlotter, it will assume geographic longitude and latitudes (WGS84), and render the mesh on a sphere.

We pass the crs as a kwarg to the GeoPlotter and render the projected scene:

plotter = gv.GeoPlotter(crs=crs)
plotter.add_mesh(sea_region, show_edges=True, clim=clim)
plotter.add_base_layer(texture=gv.natural_earth_hypsometric())
plotter.add_coastlines(resolution="10m")
plotter.view_xy()
plotter.camera.zoom(1.5)
plotter.show_axes()
plotter.show()

Note that, the base layer and coastlines are also automatically transformed to the target projection.

geovista also has an understanding of cartopy CRS’s. So let’s use cartopy to create a Mollweide CRS:

import cartopy.crs as ccrs

crs = ccrs.Mollweide(central_longitude=-90)

Before we use this CRS, let’s invert the sea_region i.e., find all cells not enclosed by the americas_bbox manifold. We can easily do this using the outside kwarg of enclosed():

outside = americas_bbox.enclosed(c48_sst, outside=True)
sea_outside = outside.threshold()

Now let’s see the final projected result:

plotter = gv.GeoPlotter(crs=crs)
plotter.add_mesh(sea_outside, show_edges=True, clim=clim)
plotter.add_base_layer(texture=gv.natural_earth_hypsometric())
plotter.add_coastlines(resolution="10m")
plotter.view_xy()
plotter.camera.zoom(1.5)
plotter.show_axes()
plotter.show()

Summary#

Key learning points of this tutorial:

using the geovista.pantry to generate meshes
creating a geodesic bounding-box manifold with panel()
introduction to the BBox class and using its boundary(), enclosed(), and mesh() methods
cartographic positioning of the camera with view_poi()
performing a cell threshold() of a mesh

Region Manifold Extraction

On this page

Region Manifold Extraction#

Introduction#

Let’s Go!#

Extension#

Summary#