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)

# Read salinity data directly in native float32 format
files_salinity = [dir_data / 'salinity_1.nc', dir_data / 'salinity_2.nc']
ds_salinity = xr.open_mfdataset(files_salinity, parallel=True)

ds_salinity.attrs = {
    "Conventions": "CF-1.8",
    "title": "Bay-Delta SCHISM model at Carquinez Strait",
    "source": "SCHISM model",
}

sxds = sx.read_schism_nc(grid, ds_salinity)

ds_salinity.attrs

{'Conventions': 'CF-1.8',
 'title': 'Bay-Delta SCHISM model at Carquinez Strait',
 'source': 'SCHISM model'}

sxds.attrs

{'Conventions': 'CF-1.8',
 'title': 'Bay-Delta SCHISM model at Carquinez Strait',
 'source': 'SCHISM model'}

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)

sxda_subset

<xarray.SxDataArray 'salinity' (time: 96, n_node: 62, n_layer: 23)> Size: 1MB
dask.array<getitem, shape=(96, 62, 23), dtype=float64, chunksize=(48, 62, 23), chunktype=numpy.ndarray>
Coordinates:
  * time                 (time) datetime64[ns] 768B 2014-12-17T16:30:00 ... 2...
    node_x               (n_node) float64 496B dask.array<chunksize=(62,), meta=np.ndarray>
    node_y               (n_node) float64 496B dask.array<chunksize=(62,), meta=np.ndarray>
    SCHISM_hgrid_node_x  (n_node) float64 496B dask.array<chunksize=(62,), meta=np.ndarray>
    SCHISM_hgrid_node_y  (n_node) float64 496B dask.array<chunksize=(62,), meta=np.ndarray>
    node_z               (n_node) float32 248B dask.array<chunksize=(62,), meta=np.ndarray>
Dimensions without coordinates: n_node, n_layer
Attributes:
    i23d:          2
    location:      node
    grid_mapping:  crs
    mesh:          SCHISM_hgrid

Calculate Volume under subset

from suxarray.grid.integrate import calc_volumes_from_projected, _calculate_volumes_uxarray, _calculate_mass

sxda_volumes = sx.SxDataArray(
    xr.DataArray(
        calc_volumes_from_projected(sxda_subset),
        dims=["time", "n_layer", "n_face"],
        coords={"time": sxda_subset.time},
        name="volume",
    ),
    sxgrid=sxda_subset.sxgrid,
 )

# Visualize total water column volume per face at first timestep
sxda_volumes.sum(dim="n_layer").isel(time=0).plot()

Calculate mass

sxda_mass = sx.SxDataArray(
    xr.DataArray(
        _calculate_mass(sxda_subset),
        dims=["time", "n_layer", "n_face"],
        coords={"time": sxda_subset.time},
        name="mass",
    ),
    sxgrid=sxda_subset.sxgrid,
)

# Visualize total water column mass per face at first timestep
sxda_mass.sum(dim="n_layer").isel(time=0).plot()

Subsetting SCHISM NetCDF data

Visualize the original data

Subset the data

Calculate Volume under subset

Compare planar vs uxarray volume

Calculate mass