Download this example.

Rasterizer post-processing in 3D

This example demonstrates the use of the rasterizer tool for post-processing 3D hydrodynamic outputs. The rasterizer interpolates fields defined on the unstructured mesh onto regular grids for visualization and analysis.

Several types of extractions are illustrated:

  • horizontal fields (surface, layers, bottom)

  • horizontal slices at fixed depth

  • vertical transects

  • projections to geographic (lon/lat) coordinates

The script reads the mesh and output fields, performs rasterization in both mesh and geographic coordinates, and visualizes the results using matplotlib.

import numpy as np
from matplotlib import pyplot as plt

from slim4 import mpi, post, slim, slim3d
from slim4.slim import proj_transform


rank = mpi.rank

We load the mesh and associated 3D output fields. The vertical coordinates and temperature field are read from the simulation output, and we select the last available time for visualization.

output_path = "/path/to/output/"
meshfile = "/path/to/meshfile.msh"

partitioned_mesh = slim3d.Mesh.partition_output(output_path, meshfile)
mesh = slim3d.Mesh.from_output(output_path, partitioned_mesh)

z_data = slim.DataOutput(output_path, "z", mesh)
temperature_data = slim.DataOutput(output_path, "temperature", mesh)
time = temperature_data.time[-1]

We define a small helper function to display rasterized fields on rank 0. This keeps the example compact while making the repeated plotting calls easier to read. By default, all MPI ranks have access to the full rasterized output, although the rendering/interpolation itself is performed in parallel (independently on each rank/GPU).

def show_field(field, extent, title=None, aspect=None):
    if rank != 0:
        return

    plt.figure()
    plt.imshow(field, cmap="turbo", extent=extent)
    if title is not None:
        plt.title(title)
    if aspect is not None:
        plt.gca().set_aspect(aspect)

We create a rasterizer in lon/lat coordinates and define a bounding box covering the full domain. The mesh coordinates are first projected to lon/lat, after which rasterization is performed in geographic coordinates.

If either width or height is left unspecified, the rasterizer automatically adjusts the other dimension to preserve square pixels in physical space rather than in coordinate space. This avoids distortions caused by effects such as meridian convergence (which would otherwise introduce an approximate factor 1/cos(lat) in the apparent aspect ratio).

lonlat = proj_transform(mesh.mesh2d.projection, "lonlat", mesh.mesh2d.xnodes)
lonmin = mpi.min_scalar(lonlat[0].min())
latmin = mpi.min_scalar(lonlat[1].min())
lonmax = mpi.max_scalar(lonlat[0].max())
latmax = mpi.max_scalar(lonlat[1].max())
lonlat_extent = (lonmin, lonmax, latmin, latmax)


horizontal_rasterizer = post.Rasterizer(mesh, projection="lonlat")
horizontal_rasterizer.set_bbox(lonmin, lonmax, latmin, latmax, height=1000)

Various horizontal fields are then extracted and rendered: surface values, values at a given layer, bottom values, and slices at fixed depth. If either width or height is negative, it is automatically inferred to preserve the aspect ratio of the bounding box.

# Top of the second layer (counting from 0).
z_layer_data = slim3d.DataExtractLayer(z_data, 1, topface=True)
layer_field = horizontal_rasterizer.render(z_layer_data, time).get()[0]
show_field(layer_field, lonlat_extent, title="Depth at the top of layer 1")

# Surface temperature.
surface_temperature = horizontal_rasterizer.render(temperature_data, time).get()[0]
show_field(surface_temperature, lonlat_extent, title="Surface temperature")

# Temperature slice at fixed depth.
depth_slice = horizontal_rasterizer.render(
    temperature_data,
    time,
    zcut=-50,
    z=z_data,
    render_type="z",
).get()[0]
show_field(depth_slice, lonlat_extent, title="Temperature at z = -50 m")

# Bottom elevation, using both discontinuous and continuous bottom extraction.
bottom_field = horizontal_rasterizer.render(z_data, time, render_type="bottom").get()[0]
show_field(bottom_field, lonlat_extent, title="Bottom")

continuous_bottom_field = horizontal_rasterizer.render(
    z_data,
    time,
    render_type="cgbottom",
).get()[0]
show_field(continuous_bottom_field, lonlat_extent, title="Continuous bottom")

We define a vertical transect across the domain. The rasterizer samples the 3D field along this transect and interpolates it onto a regular grid in the vertical plane, allowing visualization of vertical structure.

transect_xy = np.array(
    [
        [(4 * lonmin + lonmax) / 5, (lonmin + lonmax) / 2, (lonmin + 4 * lonmax) / 5],
        [(latmin + latmax) / 2, (latmin + latmax) / 2, (latmin + 2 * latmax) / 3],
    ]
)
zmin = -1000
zmax = 5
transect_width = 2000
transect_height = 1000

transect_rasterizer = post.Rasterizer(mesh, projection="lonlat")
transect_points, transect_length = transect_rasterizer.set_vertical_transect(
    transect_xy,
    zmin=zmin,
    zmax=zmax,
    width=transect_width,
    height=transect_height,
)

transect = transect_rasterizer.render(temperature_data, time, z=z_data).get()[0]

if rank == 0:
    show_field(surface_temperature, lonlat_extent, title="Transect position")
    plt.plot(transect_xy[0], transect_xy[1], ".-")
    plt.plot(transect_points[0], transect_points[1], "--k", alpha=0.5)

    show_field(
        transect, extent=(0, transect_length, zmin, zmax), title="Temperature transect"
    )
    plt.gca().set_aspect(
        (transect_height / (zmax - zmin)) / (transect_width / transect_length)
    )

Finally, we show the generated figures.

if rank == 0:
    plt.show()