Subsetting SCHISM NetCDF data

This example shows how to subset the suxarray data.

# 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 test data
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)

# 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 original data

Let’s visualize the surface salinity data at the first time step from the dataset.

# 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)

Subset the data

Let’s subset the small blob at the top left corner of the original data set. We will use bounding_box_xy method in subset accessor. This method accepts the bounding box coordinates in the form of two tuples. The first tuple is the range of x, and the second tuple the range of y coordinates. Note that the visualization uses the longitude and latitude coordinates, but the dataset keeps the original projected x and y coordinates, and they are used by this subset method.

bounding_x = [569268., 570000.]
bounding_y = [4213114., 4213630.]
sxda_subset = sxds.salinity.subset.bounding_box_xy(bounding_x, bounding_y)
sxda_subset.isel(time=0, n_layer=-1).topological_mean(destination="face").plot().opts(width=600, height=350)