Deltares / xugrid

Xarray and unstructured grids
https://deltares.github.io/xugrid/
MIT License
61 stars 8 forks source link

Add bounds to coordinates #146

Open Huite opened 1 year ago

Huite commented 1 year ago

The UGRID conventions describe the bounds coordinates. This is somewhat of a hassle to set now:

ds = uds.ugrid.to_dataset()
grid = uds.ugrid.grid

bounds = grid.face_node_coordinates
x_bounds = bounds[..., 0]
y_bounds = bounds[..., 1]

name_x = "mesh2d_x_bounds"
name_y = "mesh2d_y_bounds"
ds["mesh2d_node_x"].attrs["bounds"] = name_x
ds["mesh2d_node_y"].attrs["bounds"] = name_y
ds[name_x] = xr.DataArray(x_bounds, dims=(grid.face_dimension, "mesh2d_nMax_face_nodes"))
ds[name_y] = xr.DataArray(y_bounds, dims=(grid.face_dimension, "mesh2d_nMax_face_nodes"))

They could be added via a keyword argument to assign_face_coords or assign_edge_coords, e.g. bounds=True.

Huite commented 1 year ago

The bounds attr should be present in the face_x and face_y coordinate rather than the node coordinates.

The following works for visualizing in FEWS:

      grid = da_ptid.ugrid.grid

      bounds = grid.face_node_coordinates
      x_bounds = bounds[..., 0]
      y_bounds = bounds[..., 1]

      name_x = "mesh2d_face_x_bnd"
      name_y = "mesh2d_face_y_bnd"
      uda_waqgeom["mesh2d_face_x"].attrs["bounds"] = name_x
      uda_waqgeom["mesh2d_face_y"].attrs["bounds"] = name_y
      uda_waqgeom["mesh2d_face_x"].attrs["units"] = "degrees_east"
      uda_waqgeom["mesh2d_face_y"].attrs["units"] = "degrees_north"
      uda_waqgeom[name_x] = xr.DataArray(x_bounds, dims=(grid.face_dimension, "mesh2d_nMax_face_nodes"))
      uda_waqgeom[name_y] = xr.DataArray(y_bounds, dims=(grid.face_dimension, "mesh2d_nMax_face_nodes"))
sibrenloos commented 1 year ago

In addition to (in front of) the above code lines are used:

da_ptid = xu.UgridDataArray.from_structured(ptid) da_ptid = da_ptid.dropna(dim=da_ptid.ugrid.grid.face_dimension) uda_waqgeom = da_ptid.ugrid.to_dataset(optional_attributes=True)