Region Manifold Extraction

Region Manifold Extraction#

Tutorial ~30 minute read

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 LFRic Model - For further information refer to LFRic - a modelling system fit for future computers. cubed-sphere sample data using the lfric() and lfric_sst() convenience functions:

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)

Attention

Our goal 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.

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

Hint

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

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:

p = gv.GeoPlotter()

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

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

Seems to be right on the money! 👍

Note

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

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.

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

p.view_poi(x=30)

p.camera.zoom(1.2)
p.add_axes()
p.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:

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

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

p.view_xz()
p.camera.zoom(1.2)
p.add_axes()
p.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() You may want to experiment with the preference kwarg of the enclosed() method to see the impact on the end result. is a pyvista.PolyData. Under the hood enclosed() uses the pyvista.DataSetFilters.select_enclosed_points() method to achieve this.

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:

p = gv.GeoPlotter()

p.add_mesh(region, show_edges=True)

p.view_xz()
p.add_axes()
p.show()

Now bring this all together by rendering the extracted region When you extract data from a mesh use the add_mesh() clim kwarg to ensure the same colorbar range. Compare this render of the region mesh with the previous render without clim. Notice the difference? along with a Natural Earth texture mapped base layer and coastlines:

clim = (271, 303)

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

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

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

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 Discover the data arrays available on a mesh with array_names, and the active scalars with active_scalars_name. 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:

p = gv.GeoPlotter()

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

p.add_base_layer(texture=gv.natural_earth_hypsometric())
p.add_coastlines(resolution="10m")
p.view_xz()
p.camera.zoom(1.2)
p.show_axes()
p.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"

We pass the crs as a kwarg to the GeoPlotter If no crs is provided to GeoPlotter, it will assume geographic longitude and latitudes (WGS84), and render the mesh on a sphere. and render the projected scene:

p = gv.GeoPlotter(crs=crs)
p.add_mesh(sea_region, show_edges=True, clim=clim)
p.add_base_layer(texture=gv.natural_earth_hypsometric())
p.add_coastlines(resolution="10m")
p.view_xy()
p.camera.zoom(1.5)
p.show_axes()
p.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:

p = gv.GeoPlotter(crs=crs)
p.add_mesh(sea_outside, show_edges=True, clim=clim)
p.add_base_layer(texture=gv.natural_earth_hypsometric())
p.add_coastlines(resolution="10m")
p.view_xy()
p.camera.zoom(1.5)
p.show_axes()
p.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