Part 2 of LASSO-ShCu Tutorial: Quick-look 3-D Isosurface Plotting

This notebook contains an example 3-D isosurface plot using Plotly for the LASSO-ShCu wrfstat data. We follow the example from the Plotly documentation. We have separated this example out from lasso-shcu.ipynb since this plot is big and may cause issues for viewing the file contents without running the notebook oneself.

Author: William.Gustafson@pnnl.gov
Date: 13-May-2024

Objective of this notebook: Use plotting of 3-D isosurfaces as an example for working with the simulated cloud data from LASSO-ShCu wrfstat files.

# Libraries required for this tutorial...

from datetime import datetime, timedelta
import numpy as np
import xarray as xr

import plotly.graph_objects as go

First, we need to read the data into memory, similarly to how we did it in the lasso-shcu.ipynb notebook.

# Read in the wrfstat data...

# path_shcu_root = "/gpfs/wolf2/arm/atm124/world-shared/arm-summer-school-2024/lasso_tutorial/ShCu/untar/"  # on cumulus
path_shcu_root = "/data/project/ARM_Summer_School_2024_Data/lasso_tutorial/ShCu/untar"  # on Jupyter

case_date = datetime(2019, 4, 4)
sim_id = 4
hour_to_plot = 17

ds_stat = xr.open_dataset(f"{path_shcu_root}/{case_date:%Y%m%d}/sim{sim_id:04d}/raw_model/wrfstat_d01_{case_date:%Y-%m-%d_12:00:00}.nc")
ds_stat["Time"] = ds_stat["XTIME"]  # Fix the time coordinate

ERROR 1: PROJ: proj_create_from_database: Open of /opt/conda/share/proj failed

Now, we can make the plot, which in this case is a 3-D volume rendering of the cloud liquid water content. We only plot one quadrant of the data to keep the memory and compute usage down.

# Get the domain size and set the limitations on plot size...
nt, nz, ny, nx = ds_stat["CSV_LWC"].shape
nx_use = int(nx/3)
ny_use = int(ny/3)
nz_use = int(nz/4)

# Define the coordinates as a volume...
# We set this up in units of km using domain-averaged height data and
# assuming 100-m grid spacing.
Z, Y, X = np.meshgrid(ds_stat["CSP_Z"].sel(Time=f"{case_date:%Y-%m-%d} {hour_to_plot}:00").isel(bottom_top=slice(0,nz_use)).values*1e-3, 
                      np.arange(0,ny_use)*0.1,  # 0.1 is the conversion to km assuming dy=100 m
                      np.arange(0,nx_use)*0.1,  # and the same for dx
                      indexing='ij')

# Pull the data to plot...
plot_data = ds_stat["CSV_LWC"].\
    sel(Time=f"{case_date:%Y-%m-%d} {hour_to_plot}:00").\
    isel(west_east=slice(0,nx_use), south_north=slice(0,ny_use), bottom_top=slice(0,nz_use)).\
    values

# And, make the plot...
fig = go.Figure(data=go.Isosurface(
    x=X.flatten(), y=Y.flatten(), z=Z.flatten(),
    value=plot_data.flatten(),
    isomin=0.0001, isomax=0.001,
    colorscale='temps'
))
fig.show()

Note that Plotly permits the user to rotate, zoom, and pan the plot to see it from different sides. Just click and drag on the plot. The impact of the drag depends on the selected icon from the top-right of the panel.

Plotly is a great tool for doing quick-look 3-D plots. However, you will find it gets overwhelmed if using a large model domain. A more sophisticated way to do 3-D plots is with Paraview or VisIt. However, the learning curve for these programs is substantial and the setups can be a bit finicky on some architectures. The big advantage of these more sophisticated visualization tools comes in parallelization and server-client workflows. High-quality 3-D visualization quickly becomes a high-performance computing (HPC) application with large domains. It is common to have one or more of these visualization software packages installed on university and DOE HPC systems.

An intermediary level 3-D plotting tool is VAPOR. This one is a little easier to figure out and can more easily work with WRF output directly. Just be forewarned that trying to run any of these 3-D visualization programs requires a lot of memory and can easily bring even the highest-end laptops and desktops to their knees.