Read SCHISM NetCDF binary outputs

# Import packages
from pathlib import Path
import numpy as np
import xarray as xr
import suxarray as sx

# Where the test data is located.
# The test data is a subset from a Bay-Delta SCHISM model at Carquinez Strait.
dir_data = Path('../../tests/testdata')

Read the data and create a grid object

suxarray can create a grid object from SCHISM NetCDF binary outputs, out2d_*.nc and zCoordinates_*.nc. The function to create a grid object from those files is open_grid. It accepts the list of files names for out2d and zCoordinates files.

files_out2d = [str(p) for p in [dir_data / 'out2d_1.nc', dir_data / 'out2d_2.nc']]
files_zcoords = [str(p) for p in [dir_data / 'zCoordinates_1.nc', dir_data / 'zCoordinates_2.nc']]
grid = sx.open_grid(files_out2d, files_zcoords)

Visualize the grid

The grid object can be visualized easily using plot method. The method uses holoviews in the backend by default. Note that the coordinates in the plot uses the global coordinate system, longitude and latitude. It is because uxarray uses the global coordinate system. The original projected coordinates is stored in the grid object.

grid.plot().opts(width=600, height=400)

Create a dataset

The suxarray dataset is a combination of a grid object and a xarray dataset. A dataset can be created using open_dataset function.

# We are converting the data to float64 for later calculation.
files_salinity = [dir_data / 'salinity_1.nc', dir_data / 'salinity_2.nc']
ds_salinity = xr.open_mfdataset(files_salinity, parallel=True).astype(np.float64)
sxds = sx.read_schism_nc(grid, ds_salinity)

Visualize the data

Let’s visualize the surface salinity data at the first time step from the dataset. Since the data is defined at the grid nodes, it will be visualized in a scatter plot.

sxds.salinity.isel(time=0, n_layer=-1).plot().opts(width=600, height=350)

Take the element average and visualize the salinity at elements

To visualize the data at elements, we need to take the element average first. Method topological_mean can be used to take the element average.

# Select the surface level first, then calculate the topological mean to elements (faces)
sxds_salinity_face = sxds.salinity.isel(n_layer=-1).topological_mean(destination="face")
# Visualize the result
sxds_salinity_face.isel(time=0).plot().opts(width=600, height=350)