Assign CF Compliant Topology

This example assigns the missing CF Compliant topology information to allow for subsetting with the same workflow as the fvcom example, but with a netcdf file that does not contain the necessary metadata.

Open In Colab

# Install required libraries
# !pip install -q xarray_subset_grid@git+https://github.com/asascience-open/xarray-subset-grid.git
# !pip install -q s3fs cftime xarray cf-xarray fsspec dask h5netcdf

import cf_xarray #noqa
import fsspec
import xarray as xr

import xarray_subset_grid

# Open the dataset from NODD s3 bucket directly
fs = fsspec.filesystem("s3", anon=True)
ds = xr.open_dataset(
    fs.open(
        "s3://noaa-nos-ofs-pds/sscofs/netcdf/202403/sscofs.t21z.20240320.fields.f064.nc"
    ),
    chunks={},
    drop_variables=["Itime2"],
)
ds

<xarray.Dataset> Size: 199MB
Dimensions:             (nele: 433410, node: 239734, siglay: 10, siglev: 11,
                         three: 3, time: 1, maxnode: 11, maxelem: 9, four: 4)
Coordinates:
    lon                 (node) float32 959kB dask.array<chunksize=(239734,), meta=np.ndarray>
    lat                 (node) float32 959kB dask.array<chunksize=(239734,), meta=np.ndarray>
    lonc                (nele) float32 2MB dask.array<chunksize=(433410,), meta=np.ndarray>
    latc                (nele) float32 2MB dask.array<chunksize=(433410,), meta=np.ndarray>
    siglay              (siglay, node) float32 10MB dask.array<chunksize=(10, 239734), meta=np.ndarray>
    siglev              (siglev, node) float32 11MB dask.array<chunksize=(11, 239734), meta=np.ndarray>
  * time                (time) datetime64[ns] 8B 2024-03-23T13:00:07.031249984
Dimensions without coordinates: nele, node, three, maxnode, maxelem, four
Data variables: (12/45)
    nprocs              int32 4B ...
    partition           (nele) int32 2MB dask.array<chunksize=(433410,), meta=np.ndarray>
    x                   (node) float32 959kB dask.array<chunksize=(239734,), meta=np.ndarray>
    y                   (node) float32 959kB dask.array<chunksize=(239734,), meta=np.ndarray>
    xc                  (nele) float32 2MB dask.array<chunksize=(433410,), meta=np.ndarray>
    yc                  (nele) float32 2MB dask.array<chunksize=(433410,), meta=np.ndarray>
    ...                  ...
    wet_nodes           (time, node) int32 959kB dask.array<chunksize=(1, 239734), meta=np.ndarray>
    wet_cells           (time, nele) int32 2MB dask.array<chunksize=(1, 433410), meta=np.ndarray>
    wet_nodes_prev_int  (time, node) int32 959kB dask.array<chunksize=(1, 239734), meta=np.ndarray>
    wet_cells_prev_int  (time, nele) int32 2MB dask.array<chunksize=(1, 433410), meta=np.ndarray>
    wet_cells_prev_ext  (time, nele) int32 2MB dask.array<chunksize=(1, 433410), meta=np.ndarray>
    inundation_cells    (time, nele) int32 2MB dask.array<chunksize=(1, 433410), meta=np.ndarray>
Attributes: (12/14)
    title:                       SSCOFS
    institution:                 School for Marine Science and Technology
    source:                      FVCOM_4.4.5
    history:                     model started at: 20/03/2024   21:46
    references:                  http://fvcom.smast.umassd.edu, http://codfis...
    Conventions:                 CF-1.0
    ...                          ...
    Tidal_Forcing:               Tidal Forcing Time Series Title: sscofs late...
    River_Forcing:               THERE ARE 38 RIVERS IN THIS MODEL.\nRIVER IN...
    GroundWater_Forcing:         GROUND WATER FORCING IS OFF!
    Surface_Heat_Forcing:        FVCOM variable surface heat forcing file:\nF...
    Surface_Wind_Forcing:        FVCOM variable surface Wind forcing:\nFILE N...
    Surface_PrecipEvap_Forcing:  SURFACE PRECIPITATION FORCING IS OFF

This step assigns a new variable with the mesh_toplogy cf role to add the mesh information needed for subsetting in space

ds = xarray_subset_grid.grids.ugrid.assign_ugrid_topology(ds,
                                                          face_node_connectivity='nv',
                                                          face_face_connectivity='nbe')
ds.cf['mesh_topology']

<xarray.DataArray 'mesh' ()> Size: 4B
array(0, dtype=int32)
Attributes:
    cf_role:                 mesh_topology
    topology_dimension:      2
    node_coordinates:        lon lat
    face_node_connectivity:  nv
    face_coordinates:        lonc latc
    face_face_connectivity:  nbe

The subset grid accessor now finds the ugrid topology correctly

ds.xsg.grid.name

'ugrid'