Skip to content

Commit

Permalink
made a start -- not finished
Browse files Browse the repository at this point in the history
  • Loading branch information
ChrisBarker-NOAA committed Jun 10, 2024
1 parent fbc5c34 commit 1bccf36
Show file tree
Hide file tree
Showing 3 changed files with 179 additions and 125 deletions.
18 changes: 18 additions & 0 deletions tests/test_grids/test_ugrid.py
Original file line number Diff line number Diff line change
Expand Up @@ -285,6 +285,24 @@ def test_assign_ugrid_topology_dict():
assert mesh['face_face_connectivity'] == 'nbe'


def test_assign_ugrid_topology_start_index_one():
ds = xr.open_dataset(EXAMPLE_DATA / "SFBOFS_subset1.nc")
ds = ugrid.assign_ugrid_topology(ds, **grid_topology)
# bit fragile, but easy to write the test
assert ds['nv'].attrs['start_index'] == 1
assert ds['nbe'].attrs['start_index'] == 1


def test_assign_ugrid_topology_start_index_zero():
ds = xr.open_dataset(EXAMPLE_DATA / "small_ugrid_zero_based.nc")
ds = ugrid.assign_ugrid_topology(ds, face_node_connectivity='mesh_face_nodes')
# bit fragile, but easy to write the test
assert ds['mesh_face_nodes'].attrs['start_index'] == 0
assert ds['mesh_edge_nodes'].attrs['start_index'] == 0
assert ds['mesh_boundary_nodes'].attrs['start_index'] == 0



# NOTE: these tests are probably not complete -- but they are something.
# we really should have a complete UGRID example to test with.
def test_grid_vars():
Expand Down
151 changes: 150 additions & 1 deletion xarray_subset_grid/grids/ugrid.py
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,6 @@

from xarray_subset_grid.grid import Grid
from xarray_subset_grid.utils import (
assign_ugrid_topology,
normalize_polygon_x_coords,
ray_tracing_numpy,
)
Expand Down Expand Up @@ -168,3 +167,153 @@ def subset_polygon(
if has_face_face_connectivity:
ds_subset[mesh.face_face_connectivity][:] = face_face_new
return ds_subset

def assign_ugrid_topology(ds: xr.Dataset,
*,
face_node_connectivity: str,
face_face_connectivity: str = None,
node_coordinates: str = None,
face_coordinates: str = None,
start_index: int = None,
) -> xr.Dataset:
# Should this be "make entire dataset UGRID compliant ?"
# That would mean that the grid variables should all get metadata,
# such as "location"
# and we'd need to clean up coordinates that shouldn't be coordinates.
# ("node is one that's in the UGRID test file")
"""
Assign the UGRID topology to the dataset
Only the face_node_connectivity parameter is required.
The face_face_connectivity parameter is optional.
If the variable for face_node_connectivity is named nv, the function call should look like this:
```
ds = assign_ugrid_topology(ds, face_node_connectivity="nv")
```
This will assign a new variable to the dataset called "mesh" with the face_node_connectivity attribute. It will
also assign the node_coordinates attribute to the dataset with the lon and lat coordinate variable names, introspected
from the face_node_connectivity variables' dimensions.
You can also optionally specify any of the following mesh topology variables by passing them as keyword arguments
- node_coordinates: If not specified, the function will introspect the dataset for the longitude and latitude coordinate variable names using cf_xarray
- face_face_connectivity
- face_coordinates: If not specified, the function will introspect the dataset for the longitude and latitude coordinate variable names matching the face_node_connectivity variable
but do nothing if they are not found
Args:
ds (xr.Dataset): The dataset to assign the UGRID topology to
face_node_connectivity (str): THe variable name of the face definitions
face_face_connectivity: str = None,
node_coordinates: str = None,
face_coordinates: str = None,
(See the UGRID conventions for descriptions of these)
You can pass a dict in with all the grid topology variables:
```
grid_topology = {'node_coordinates': ('lon', 'lat'),
'face_node_connectivity': 'nv',
'node_coordinates': ('lon', 'lat'),
'face_coordinates': ('lonc', 'latc'),
}
```
"""
# All the possible attributes:
# Required:
# node_coordinates
# face_node_connectivity

# Optional:
# face_dimension
# edge_node_connectivity
# edge_dimension
# Optional attributes
# face_edge_connectivity
# face_face_connectivity
# edge_face_connectivity
# boundary_node_connectivity
# face_coordinates
# edge_coordinates

# Get the variable name for the face_node_connectivity
# face_node_connectivity = attrs.get("face_node_connectivity", None)
# if face_node_connectivity is None:
# raise ValueError("The face_node_connectivity attribute is required")
#face_face_connectivity = attrs.get("face_face_connectivity", None)

# Get the longitude and latitude coordinate variable names
# node_coords = attrs.get("node_coordinates", None)
# face_coords = attrs.get("face_coordinates", None)

node_coords = node_coordinates
face_coords = face_coordinates

if not face_coords:
try:
face_coords = ds[face_node_connectivity].cf.coordinates
face_coords = [f"{coord[0]}" for coord in face_coords.values()]
except AttributeError:
face_coords = None

if not node_coords:
try:
if face_coords:
filter = face_coords
else:
filter = []

coords = ds.cf.coordinates
node_lon = [c for c in coords["longitude"] if c not in filter][0]
node_lat = [c for c in coords["latitude"] if c not in filter][0]
node_coords = [node_lon, node_lat]
except AttributeError:
raise ValueError(
"The dataset does not have cf_compliant node coordinates longitude and latitude coordinates"
)

if not face_face_connectivity:
# face_face_connectivity var should have same dimensions as
# face_node_connectivity this is assuming that only one will match!
for var_name, var in ds.variables.items():
if var_name == face_node_connectivity:
continue
if var.dims == ds[face_node_connectivity].dims:
face_face_connectivity = var_name
break

mesh_attrs = {
"cf_role": "mesh_topology",
"topology_dimension": np.int32(2),
"node_coordinates": " ".join(node_coords),
"face_node_connectivity": face_node_connectivity,
}

if face_coords:
mesh_attrs["face_coordinates"] = " ".join(face_coords)
if face_face_connectivity:
mesh_attrs["face_face_connectivity"] = face_face_connectivity

# Assign the mesh topology to the dataset
ds = ds.assign(
mesh=((), np.int32(0), mesh_attrs),
)

# check for start_index, and set it if there.
if start_index is None:
start_index = int(ds[face_node_connectivity].min())

if not start_index in (0, 1):
raise ValueError(f"start_index must be 1 or 0, not {start_index}")

# assign the start_index to all the grid variables
for var_name in (face_node_connectivity, face_face_connectivity):
if var_name:
var = ds[var_name]
var.attrs.setdefault('start_index', start_index)

return ds
135 changes: 11 additions & 124 deletions xarray_subset_grid/utils.py
Original file line number Diff line number Diff line change
Expand Up @@ -67,135 +67,22 @@ def ray_tracing_numpy(x, y, poly):
return inside


# This should probably be defined in ugrid.py -- but importing it there
# for now.
def assign_ugrid_topology(ds: xr.Dataset,
*,
face_node_connectivity: str,
face_face_connectivity: str = None,
node_coordinates: str = None,
face_coordinates: str = None,
) -> xr.Dataset:
# Should this be "make entire dataset UGRID compliant ?"
# That would mean that the grid varaibles should all get metadata, such as "location"
# and we'd need to clean up coordinates that shouldn't be coordinates.
# ("node is one that's in the UGRID test file")
"""
Assign the UGRID topology to the dataset
Only the face_node_connectivity parameter is required.
The face_face_connectivity parameter is optional.
If the variable for face_node_connectivity is named nv, the function call should look like this:
```
ds = assign_ugrid_topology(ds, face_node_connectivity="nv")
```
This will assign a new variable to the dataset called "mesh" with the face_node_connectivity attribute. It will
also assign the node_coordinates attribute to the dataset with the lon and lat coordinate variable names, introspected
from the face_node_connectivity variables' dimensions.
You can also optionally specify any of the following mesh topology variables by passing them as keyword arguments
- node_coordinates: If not specified, the function will introspect the dataset for the longitude and latitude coordinate variable names using cf_xarray
- face_face_connectivity
- face_coordinates: If not specified, the function will introspect the dataset for the longitude and latitude coordinate variable names matching the face_node_connectivity variable
but do nothing if they are not found
Args:
ds (xr.Dataset): The dataset to assign the UGRID topology to
face_node_connectivity (str): THe variable name of the face definitions
face_face_connectivity: str = None,
node_coordinates: str = None,
face_coordinates: str = None,
(See the UGRID conventions for descriptions of these)
You can pass a dict in with all the grid topology variables:
```
grid_topology = {'node_coordinates': ('lon', 'lat'),
'face_node_connectivity': 'nv',
'node_coordinates': ('lon', 'lat'),
'face_coordinates': ('lonc', 'latc'),
}
```
"""
# Get the variable name for the face_node_connectivity
# face_node_connectivity = attrs.get("face_node_connectivity", None)
# if face_node_connectivity is None:
# raise ValueError("The face_node_connectivity attribute is required")
#face_face_connectivity = attrs.get("face_face_connectivity", None)

# Get the longitude and latitude coordinate variable names
# node_coords = attrs.get("node_coordinates", None)
# face_coords = attrs.get("face_coordinates", None)

node_coords = node_coordinates
face_coords = face_coordinates

if not face_coords:
try:
face_coords = ds[face_node_connectivity].cf.coordinates
face_coords = [f"{coord[0]}" for coord in face_coords.values()]
except AttributeError:
face_coords = None

if not node_coords:
try:
if face_coords:
filter = face_coords
else:
filter = []

coords = ds.cf.coordinates
node_lon = [c for c in coords["longitude"] if c not in filter][0]
node_lat = [c for c in coords["latitude"] if c not in filter][0]
node_coords = [node_lon, node_lat]
except AttributeError:
raise ValueError(
"The dataset does not have cf_compliant node coordinates longitude and latitude coordinates"
)

if not face_face_connectivity:
# face_face_connectivity var should have same dimensions as
# face_node_connectivity this is assuming that only one will match!
for var_name, var in ds.variables.items():
if var_name == face_node_connectivity:
continue
if var.dims == ds[face_node_connectivity].dims:
face_face_connectivity = var_name
break

mesh_attrs = {
"cf_role": "mesh_topology",
"topology_dimension": np.int32(2),
"node_coordinates": " ".join(node_coords),
"face_node_connectivity": face_node_connectivity,
}

if face_coords:
mesh_attrs["face_coordinates"] = " ".join(face_coords)
if face_face_connectivity:
mesh_attrs["face_face_connectivity"] = face_face_connectivity

# Assign the mesh topology to the dataset
ds = ds.assign(
mesh=((), np.int32(0), mesh_attrs),
)

return ds


def convert_bytes(num):
# This is defined in ugrid.py
# this placeholder for backwards compatibility for a brief period
def assign_ugrid_topology(*args, **kwargs):
warnings.warn(DeprecationWarning, "The function `assign_grid_topology` has been moved to the "
"`grids.ugrid` module. It will not be able to be called from "
"the utils `module` in the future.")
from .grids.ugrid import assign_ugrid_topology
return assign_ugrid_topology(*args, **kwargs)

def format_bytes(num):
"""
This function will convert bytes to MB.... GB... etc
from:
https://stackoverflow.com/questions/12523586/python-format-size-application-converting-b-to-kb-mb-gb-tb
"""
# stupid, but handy for demos, etc.
# not much to it, but handy for demos, etc.

step_unit = 1024 #1024 bad the size

Expand Down

0 comments on commit 1bccf36

Please sign in to comment.