Visualizing Meshes

As part of the EarthSim project, extensive support for visualizing triangular meshes was added to HoloViews and Datashader; see the websites for those projects for details of how to use this functionality. Here, we will show how to use the small utilities provided in earthsim itself, which allow you to read 3DM and mesh2d files with mesh data, and we'll then visualize the resulting structures.

In [1]:
import numpy as np
import holoviews as hv
import geoviews as gv
import cartopy.crs as ccrs

from holoviews import opts
from holoviews.operation.datashader import datashade, rasterize
import datashader as ds
from colorcet import cm_n
from earthsim.io import read_3dm_mesh, read_mesh2d

datashade.precompute = True

hv.extension('bokeh')

size = dict(width=800, height=600)
opts.defaults(
    opts.Image(**size), opts.RGB(**size), opts.VectorField(**size))

Static meshes

Static trimeshes can be loaded using the read_3dm_mesh utility, which will return its simplices and vertices. Here we will load a mesh of the Chesapeake and Delaware Bays.

Before we declare the TriMesh we will project the coordinates of the vertices to a Mercator projection for plotting, so that GeoViews does not have to re-project the coordinates each time the plot is displayed. We also declare the 'z' dimension representing the value at each vertex (water depth in this case).

In [2]:
fpath = '../data/Chesapeake_and_Delaware_Bays.3dm'
tris, verts = read_3dm_mesh(fpath)
points = gv.operation.project_points(gv.Points(verts, vdims=['z']))

Now that we have the simplices and vertices we can construct the TriMesh. Since the TriMesh is too large to display on its own in Bokeh (which provides all data to the web browser by default), we apply the datashade operation to render the data into an image before providing it to the browser. To begin with, we will visualize the underlying mesh by datashading just trimesh.edgepaths:

In [3]:
tiles = gv.WMTS('https://maps.wikimedia.org/osm-intl/{Z}/{X}/{Y}@2x.png')
trimesh = gv.TriMesh((tris, points))
chesapeake_mesh = tiles * datashade(trimesh.edgepaths)
chesapeake_mesh
Out[3]:

Here the plot may look continuous in most regions, but if you have a live, running version of this plot, you can zoom in to see that it's a wireframe with high resolution (many grid points and vertices) in the inland areas and low resolution in the open water.

Since the mesh also has associated 'z' values, we can render the value at each pixel that falls inside a triangle, interpolating from the values at the vertices when there are multiple pixels per triangle, and aggregating across vertex values (averaging in this case) when there are multiple vertices per pixel:

In [4]:
tiles = gv.WMTS('https://maps.wikimedia.org/osm-intl/{Z}/{X}/{Y}@2x.png')
trimesh = gv.TriMesh((tris, points))
chesapeake_interp = tiles * rasterize(trimesh, aggregator=ds.mean('z')).opts(colorbar=True, cmap='Viridis')
chesapeake_interp
Out[4]:

If you have a live running Python process (not a static web page), you can zoom in to the aboe plot and eventually make out the original triangles.

Meshes over time

When running environmental simulations using a tool like AdH, the mesh is also often accompanied by additional data files representing depths, velocity, and error values over time. We can load in these values and use them to visualize changes in each variable as a simulation evolves.

First, we again load a static mesh and project it from its native UTM Zone 11 coordinate system to the Mercator projection.

In [5]:
fpath = '../data/SanDiego_Mesh/SanDiego.3dm'
tris, verts = read_3dm_mesh(fpath, skiprows=2)
points = gv.operation.project_points(gv.Points(verts, vdims=['z'], crs=ccrs.UTM(11)))
trimesh = gv.TriMesh((tris, points))
san_diego = tiles * rasterize(trimesh, aggregator='mean').opts(colorbar=True, cmap='Viridis', clim=(-75, 0))
san_diego
Out[5]: