After the last update to XROMS, the interpolation module no longer works for horizontal interpolation. But rather errors with the message:

NameError: name 'xe' is not defined

in interp.py.

It looks as though the module import for xesmf is commented out.

Hello,
I just accidentally found your code, and it seems really useful! I am a CROCO/ROMS-Agrif user, do you know if your code is compatible with its outputs? Most things are similar, but since there are some differences in the name of variables, and sigma coordinates, would it deal with it correctly?
Thank you for developing and sharing such useful tools.
Cheers
Ana Machado

I am trying to run an xroms workflow I have used successfully in the past. I upgraded recently to dask=2022.7.1 and xarray=2022.6.0 after which this is the first time I am returning to the same workflow. I am trying something like this, very basic stuff and early in the workflow before any real computations:

xrargs = dict(parallel='True',data_vars='minimal',
              coords='minimal',compat='override',drop_variables=droplist)

chunks = 'auto'

# Just reading in two files from the full list of files
ds = xroms.open_mfnetcdf(files[0:2],chunks=chunks,xrargs=xrargs)

NotImplementedError: Cannot chunk along a core dimension for a grid ufunc which has a signature which includes one of the axis positions ['inner', 'outer']

    812 ds = xr.open_mfdataset(files, chunks=chunks, **xrargsin)
    814 # modify Dataset with useful ROMS z coords and make xgcm grid operations usable.
--> 815 ds, grid = roms_dataset(ds, Vtransform=Vtransform, add_verts=add_verts, proj=proj)
    817 return ds

File ./xroms/xroms.py:194, in roms_dataset(ds, Vtransform, add_verts, proj)
    187 z_w0.attrs = {
    188     "long_name": "depth of W-points",
    189     "field": "z_w0, scalar",
    190     "units": "m",
    191 }
    193 ds.coords["z_w"] = xroms.order(z_w)
--> 194 ds.coords["z_w_u"] = grid.interp(ds.z_w, "X")
    195 ds.coords["z_w_u"].attrs = {
    196     "long_name": "depth of U-points on vertical W grid",
    197     "time": "ocean_time",
    198     "field": "z_w_u, scalar, series",
    199     "units": "m",
    200 }
    201 ds.coords["z_w_v"] = grid.interp(ds.z_w, "Y")

File /software/conda/envs/dask_2022/lib/python3.10/site-packages/xgcm/grid.py:2041, in Grid.interp(self, da, axis, **kwargs)
   1989 def interp(self, da, axis, **kwargs):
   1990     """
   1991     Interpolate neighboring points to the intermediate grid point along
   1992     this axis.
   (...)
   2039     >>> grid.interp(da, ["X", "Y"], periodic={"X": True, "Y": False})
   2040     """
-> 2041     return self._1d_grid_ufunc_dispatch("interp", da, axis, **kwargs)

File /software/conda/envs/dask_2022/lib/python3.10/site-packages/xgcm/grid.py:1836, in Grid._1d_grid_ufunc_dispatch(self, funcname, data, axis, to, keep_coords, metric_weighted, other_component, **kwargs)
   1833 else:
   1834     map_overlap = False
-> 1836 array = grid_ufunc(
   1837     self,
   1838     array,
   1839     axis=[(ax_name,)],
   1840     keep_coords=keep_coords,
   1841     dask=dask,
   1842     map_overlap=map_overlap,
   1843     other_component=other_component,
   1844     **remaining_kwargs,
   1845 )
   1847 if ax_metric_weighted:
   1848     metric = self.get_metric(array, ax_metric_weighted)

File /software/conda/envs/dask_2022/lib/python3.10/site-packages/xgcm/grid_ufunc.py:460, in GridUFunc.__call__(self, grid, axis, *args, **kwargs)
    458 map_overlap = kwargs.pop("map_overlap", self.map_overlap)
    459 pad_before_func = kwargs.pop("pad_before_func", self.pad_before_func)
--> 460 return apply_as_grid_ufunc(
    461     self.ufunc,
    462     *args,
    463     axis=axis,
    464     grid=grid,
    465     signature=self.signature,
    466     boundary_width=self.boundary_width,
    467     boundary=boundary,
    468     dask=dask,
    469     map_overlap=map_overlap,
    470     pad_before_func=pad_before_func,
    471     **kwargs,
    472 )

File /software/conda/envs/dask_2022/lib/python3.10/site-packages/xgcm/grid_ufunc.py:745, in apply_as_grid_ufunc(func, axis, grid, signature, boundary_width, boundary, fill_value, keep_coords, dask, map_overlap, pad_before_func, other_component, *args, **kwargs)
    742 # Maybe map function over chunked core dims using dask.array.map_overlap
    743 if map_overlap:
    744     # Disallow situations where shifting axis position would cause chunk size to change
--> 745     _check_if_length_would_change(sig)
    747     mapped_func = _map_func_over_core_dims(
    748         func,
    749         args,
   (...)
    753         out_dtypes,
    754     )
    755 else:

File /software/conda/envs/dask_2022/lib/python3.10/site-packages/xgcm/grid_ufunc.py:1002, in _check_if_length_would_change(signature)
    996 all_ax_positions = set(
    997     p
    998     for arg_ps in signature.in_ax_positions + signature.out_ax_positions
    999     for p in arg_ps
   1000 )
   1001 if any(pos in DISALLOWED_OVERLAP_POSITIONS for pos in all_ax_positions):
-> 1002     raise NotImplementedError(
   1003         "Cannot chunk along a core dimension for a grid ufunc which has a signature which "
   1004         f"includes one of the axis positions {DISALLOWED_OVERLAP_POSITIONS}"
   1005     )

NotImplementedError: Cannot chunk along a core dimension for a grid ufunc which has a signature which includes one of the axis positions ['inner', 'outer']

1. I have a dataset from an ocean model (ROMS).
2. I use resample to compute daily-averaged values of two variables, do some basic operations to remove the monthly climatology and compute the covariance of the anomalies.
3. When I use xroms.ddxi to get the x-derivative of the covariance, I get an error message complaining about the absence of metrics.

I have uploaed a copy of the notebook (in html format) which shows the error message here (have omitted some initial cells to create the list of filenames, etc.)

The attached notebook shows the grid attribute is present for the variable uT_atu so I am not sure why I am getting that error message.

The vertical coordinate in stratification frequency requires quite a bit of wrangling, dropping and replacing vertical coordinates. This part could be re-written to be simpler and more clear. I hope @dcherian could help with this.

The recombined resulting DataArray is on vertical w-points, as it should be, but only the inner points are defined. This will make it difficult to work well with other ROMS DataArrays that are defined on the 'outer' vertical points by default. It would be better to extrapolate values to the surface and bottom, or replace these with boundary values (i.e., zero at the bottom).

Would be nice to be able to conda install xroms from conda-forge.

Any interest in releasing on pypi?

If released on pypi we could then use https://github.com/conda/grayskull to create the conda-forge recipe...

Hello,
I'm running into an issue when trying to run roms_dataset that I believe may be related to a recent upgrade of my xarray package (version 2024.5.0) and is hopefully useful for others. Below is an example of the relevant code that I'm running.

import xarray as xr
import xroms

ds = xr.open_dataset(file_name, chunks={'ocean_time': 1})
ds = xroms.roms_dataset(ds)

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[7], line 2
      1 ds = xr.open_dataset(file_name, chunks={'ocean_time': 1})
----> 2 ds, xgrid = xroms.roms_dataset(ds)

File /global/u2/user/mambaforge/envs/mamba_env1/lib/python3.10/site-packages/xroms/xroms.py:243, in roms_dataset(ds, Vtransform, add_verts, proj, include_Z0, include_3D_metrics, include_cell_volume, include_cell_area)
    232     z_rho0.attrs = {
    233         "long_name": "depth of RHO-points",
    234         "field": "z_rho0, scalar",
    235         "units": "m",
    236     }
    237     z_w0.attrs = {
    238         "long_name": "depth of W-points",
    239         "field": "z_w0, scalar",
    240         "units": "m",
    241     }
--> 243 ds.coords["z_w"] = order(z_w)
    244 ds.coords["z_w_u"] = grid_interp(xgrid, ds["z_w"], "X")
    245 # ds.coords["z_w_u"] = xgrid.interp(ds.z_w, "X")

File /global/u2/user/mambaforge/envs/mamba_env1/lib/python3.10/site-packages/xroms/utilities.py:1572, in order(var)
   1546 """Reorder var to typical dimensional ordering.
   1547 
   1548 Parameters
   (...)
   1566 >>> xroms.order(var)
   1567 """
   1569 # this looks at the dims on var directly which tends to be more accurate
   1570 # since the DataArray can have extra coordinates included (but dropping
   1571 # can drop too many variables).
-> 1572 return var.cf.transpose(
   1573     *[
   1574         dim
   1575         for dim in ["T", "Z", "Y", "X"]
   1576         if dim in var.cf.axes and var.cf.axes[dim][0] in var.dims
   1577     ]
   1578 )

File /global/u2/user/mambaforge/envs/mamba_env1/lib/python3.10/site-packages/cf_xarray/accessor.py:693, in _getattr.<locals>.wrapper(*args, **kwargs)
    689 posargs, arguments = accessor._process_signature(
    690     func, args, kwargs, key_mappers
    691 )
    692 final_func = extra_decorator(func) if extra_decorator else func
--> 693 result = final_func(*posargs, **arguments)
    694 if wrap_classes and isinstance(result, _WRAPPED_CLASSES):
    695     result = _CFWrappedClass(result, accessor)

File /global/u2/user/mambaforge/envs/mamba_env1/lib/python3.10/site-packages/xarray/util/deprecation_helpers.py:140, in deprecate_dims.<locals>.wrapper(*args, **kwargs)
    132     emit_user_level_warning(
    133         f"The `{old_name}` argument has been renamed to `dim`, and will be removed "
    134         "in the future. This renaming is taking place throughout xarray over the "
   (...)
    137         PendingDeprecationWarning,
    138     )
    139     kwargs["dim"] = kwargs.pop(old_name)
--> 140 return func(*args, **kwargs)

File /global/u2/user/mambaforge/envs/mamba_env1/lib/python3.10/site-packages/xarray/core/dataarray.py:3058, in DataArray.transpose(self, transpose_coords, missing_dims, *dim)
   3025 """Return a new DataArray object with transposed dimensions.
   3026 
   3027 Parameters
   (...)
   3055 Dataset.transpose
   3056 """
   3057 if dim:
-> 3058     dim = tuple(infix_dims(dim, self.dims, missing_dims))
   3059 variable = self.variable.transpose(*dim)
   3060 if transpose_coords:

File /global/u2/user/mambaforge/envs/mamba_env1/lib/python3.10/site-packages/xarray/namedarray/utils.py:174, in infix_dims(dims_supplied, dims_all, missing_dims)
    172             yield d
    173 else:
--> 174     existing_dims = drop_missing_dims(dims_supplied, dims_all, missing_dims)
    175     if set(existing_dims) ^ set(dims_all):
    176         raise ValueError(
    177             f"{dims_supplied} must be a permuted list of {dims_all}, unless `...` is included"
    178         )

File /global/u2/user/mambaforge/envs/mamba_env1/lib/python3.10/site-packages/xarray/namedarray/utils.py:128, in drop_missing_dims(supplied_dims, dims, missing_dims)
    126     supplied_dims_set = {val for val in supplied_dims if val is not ...}
    127     if invalid := supplied_dims_set - set(dims):
--> 128         raise ValueError(
    129             f"Dimensions {invalid} do not exist. Expected one or more of {dims}"
    130         )
    132     return supplied_dims
    134 elif missing_dims == "warn":

ValueError: Dimensions {'X', 'T', 'Y', 'Z'} do not exist. Expected one or more of ('ocean_time', 'eta_rho', 'xi_rho', 's_w')

Hi,

A trivial detail, but slightly annoying.

There are four syntax errors in xroms.py all related to the docstrings.
For roms_dataset, open_netcdf, open_zarr the docstring starts with ''' and end with """.
For the function hgrad is opposite, starts with """ and end with '''.

Automatic running of black before commit can be discussed,
but black (and possibly other tools) will prevent this kind of errors to get into the repo.

Bjørn Ådlandsvik

Hi!

I just found your code by scanning through the PANGEO forum and I am impressed by the scope of functions that seem to be implemented. I have written a package called xcroco which I started many years ago when nothing else was available. I've had brief discussions with @jbusecke from PANGEO and came to use xgcm upon his recommendation. xcroco is what I use to load and analyse my own CROCO model output. I haven't promoted it really and not many people other than myself have used it.

So i was wondering if you could fill me in on the current developments regarding analysis of ROMS / CROCO output with xarray and xgcm. Is xroms what most people working with those outputs seem to be using in 2022? If so, I would be happy to join the effort, contribute some of my code etc. if there is interest. I'm not sure what functionality may be missing that I could help with. You could check out my xcroco example notebook to get an idea of what my package contains so far.

Hope to collaborate in the future!
I feel like the ROMS / CROCO community is in need of some standardization... a lot of people seem to still be using their own tools, written in Matlab (!) or Python.

Cheers,
Jaard

Dear All,

I started to analyze my CROCO output by using XROMS. Really I loved it, it is amazing for calculation of several physical variables in a very efficient and smart way specially for curvilinear grid system.

Recently, I tried to calculate the geostrophic current and I found that in the derived.py the surface zonal geostrophic current are calculated as
ug = (-g * dzetadxi / to_u(f, xgrid, hboundary=hboundary, hfill_value=hfill_value).

However in ROMS the xi represent the X direction and eta represent the Y direction so basically the derivative of zeta with respect to xi shall represent vg instead of ug.

It is my understanding, if some thing I am missing please correct me.

Thank you.

Anurag

I'm probably missing something, but there seems to be an issue interpolating in space the barotropic (2-d spatial., plus time) output variables. (such as zeta) in the below code I and trying to interpolate zeta to a -D set of spatial points, and throws an error.

`import xroms
import xarray as xr

lat=[-68.7712,
-69.0220,
-68.9802,
-69.5237,
-70.7780,
-72.0740,
-73.1610,
-73.2865,
-73.9972,
-74.6661]
lon=[ 43.1379,
42.3528,
41.6986,
40.9136,
40.7173,
40.3902,
40.0304,
39.3435,
38.5911,
37.5117]

url='https://tds.marine.rutgers.edu/thredds/dodsC/roms/doppio/2017_da/his/History_Best'
ds=xroms.open_netcdf(url)
dstmp=ds.isel(time=100)
test = dstmp.zeta.xroms.interpll(lon, lat,which='pairs')
`

The error is:

IndexError Traceback (most recent call last)
/tmp/ipykernel_184701/2301636929.py in
26 ds=xroms.open_netcdf(url)
27 dstmp=ds.isel(time=100)
---> 28 test = dstmp.zeta.xroms.interpll(lon, lat,which='pairs')

/cache/home/elhunter/python/packages/xroms/xroms/accessor.py in interpll(self, lons, lats, which)
1360 """
1361
-> 1362 return xroms.interpll(self.da, lons, lats, which=which)
1363
1364 def isoslice(self, iso_values, iso_array=None, axis="Z"):

/cache/home/elhunter/python/packages/xroms/xroms/interp.py in interpll(var, lons, lats, which)
102 # get z coordinates to go with interpolated output if not available
103 if zkey_varint == []:
--> 104 zkey = [coord for coord in var.coords if "z_" in coord and "0" not in coord][
105 0
106 ] # str

IndexError: list index out of range

It looks as though it is trying to find a vertical coordinate, and fails since the zeta variable does not have one. Is that correct?

Thanks,
Eli

Right now open_zarr requires a list, but it should also accept a string for a single dataset.

Hi,

First, I am not negative to pygridgen by itself. I have seen Rob giving very impressive demonstrations of the software. But I do not think it is a good idea to make xroms depend on it. I will illustrate this be giving my user experience today.

I try to run the examples, in particular the new interpolation example.

After fixing the bugs described in my previous issue. I started again.
A download of 2GB for running an example, is on the heavy side, but OK.

With the data in place, I started again. Crash, pygridgen is missing.

No problem I thought, conda is your friend, and conda install pygridgen worked nicely.

Restart kernel and try again. Now it crashed due to a bug in pygridgen:
While executing the line

  ds, grid = xroms.roms_dataset(ds, add_verts=True, proj=proj)

I get a NameError from pygridgen/grid.py saying "empty" is not defined
The offending code looks like

  1017     x = empty((Mp + 1, Lp + 1), dtype='d')

Of course, conda may have given me a bad version

pygridgen                 0.2.1            py38hc1659b7_3    conda-forge

Then I go one step further, try github and find https://github.com/rsignell-usgs/pygridgen
As this is maintained by a master, I was quite happy.
However, the last commit was in 2012, it is written in python 2 and depends on basemap
and some obscure features of matplotlib no longer existing.

I am sure you at Texas A&M have an updated version of pygridgen working nicely for you.
But I was hoping that xroms should stay away from hard-to-come-by and hard-to-install specific
C-based dependencies.

The vertices defined by Rob is a good idea that may be useful in xroms, but they do not
require the full pygridgen package.

General comment

It boils down to what you (we?) want xroms to be. A package supporting research at TA&M
(of course an important use of your own developer resources)
or a package of general use by the intersection of the ROMS and python communities. It is however
important to tell potential users the main objectives of the package.

BTW. I feel the same tension with my particle tracking code. It is tempting to tailor it for IMR research
and our own computing environment, but I try to avoid it to make the code generally useful.

Sorry for the long rant,
Bjørn Å.

New roms version has Vtransform from 1-5 and so far xroms only has 1 and 2.

Hello,

While browsing through xroms code, I noticed a possible error in the calculation of 'dz':

    ds["dz"] = grid.diff(ds.z_w, "Z")
    ds["dz"].attrs = {
        "long_name": "vertical layer thickness on vertical RHO grid",
        "time": "ocean_time",
        "field": "dz, scalar, series",
        "units": "m",
    }´

    ds["dz_w"] = grid.diff(ds.z_rho, "Z", boundary="fill")
    ds["dz_w"].attrs = {
        "long_name": "vertical layer thickness on vertical W grid",
        "time": "ocean_time",
        "field": "dz_w, scalar, series",
        "units": "m",
    }

Came across this problem when trying to load in a normal output NetCDF file. Currently I try to apply xroms to my dataset loaded in with xarray and get no error message before python kills the application. I'm running this through a Jupyter notebook, and had the same issue come up using python in a terminal window.

file = 'roms_his.nc'
ds = xarray.open_dataset(file)
ds = xroms.roms_dataset(ds)

I stepped through the roms_dataset function line by line and am first encountering an issue on line 511 which is calculating dz_w with:

ds["dz_w"] = xgrid.diff(ds.z_rho, "Z", boundary="fill")

The variable z_rho is present in the dataset and appears to be calculated normally, so I'm not sure why the diff command works above on line 503 to calculate dz and not below for dz_w.

Turns out that I missed specifying chunking originally, doing so fixes my issue

ds = xarray.open_dataset(files[0], chunks={'ocean_time': 1})
ds = xroms.roms_dataset(ds)

Hi all,

Great package! In migrating some of my workflow across from home-cooked xarray/xgcm to xroms, I have come across an issue with saving a dataset opened with xroms.open_mfnetcdf to netcdf.

When I run ds.isel(ocean_time=0).to_netcdf('filename.nc')

I get the following error to do with attributes associated with the grid object:

`TypeError: Invalid value for attr 'grid': <xgcm.Grid>
X Axis (not periodic, boundary=None):

• center xi_rho --> inner
• inner xi_u --> center
  Y Axis (not periodic, boundary=None):
• center eta_rho --> inner
• inner eta_v --> center
  Z Axis (not periodic, boundary=None):
• center s_rho --> outer
• outer s_w --> center. For serialization to netCDF files, its value must be of one of the following types: str, Number, ndarray, number, list, tuple`

Has anyone else encountered this behaviour? Unfortunately I was not able to access the files used in the example notebook, so this may be a function of my particular ROMS dataset.

Hi everyone,

I have a memory issue when plotting the Hovmoller diagram of the vertical component of surface/bottom relative vorticity (dvdx-dudz) in a large ROMS dataset, which has the dimension of {"ocean_time": 96, "s_rho": 30, "eta_v": 601, "xi_u": 677} for velocity v (@hetland ).

Here is what I did.

`
import matplotlib.pyplot as plt
import xroms
import cmocean
from glob import glob
%matplotlib inline

files = glob('./roms_*.nc')
ds = xroms.open_mfnetcdf(files)
ds, grid = xroms.roms_dataset(ds)

zeta = xroms.relative_vorticity(ds.u,ds.v,grid)
plt.pcolormesh(ds.lon_rho,ds.ocean_time,zeta.isel(s_w=-1, eta_v=0),cmap=cmocean.cm.balance)
`
In the last line, I tried to get a pcolor of zeta which should have a size of 96 by 677 (ocean_time by xi_u).
And I did this in a jupyter notebook on a local desktop which has 16GB memory. However, it got stuck at the very last command line and after some time the notebook ran out of memory. Does anyone have any idea how to solve this? Thanks.

Jinliang

Hi,

I am generally very impressed by the xisoslice function. It is a good example of what can be done
within the xarray framework. The new examples show that it is widely useful.

However, if the iso-value is found exactly in the iso-array it returns not-a-number.
I consider this a bug.

My attempt for an easy fix was to change the test for negativity to <= 0.
Then it returns a value, unfortunately in error.

There is also a more arguable bug if the iso-value is an upper or lower boundary of the
iso-array. Also here NaN is returned. This may be the expected behaviour.

Below is the test function a wrote to understand the behaviour of xisoslice.
The bug is caught by the 2nd test, test_no_interpolation.
The questionable boundary handling is exposed by the 3rd test.

------------ 8< ------------------------

import numpy as np
import xarray as xr
import xroms
import pytest

Z = np.linspace(-50, 0, 51)
A = np.sin(Z)
Zdims = ["Z"]
Zcoords = dict(Z=Z)

1D tests

def test_linear_OK():
"""Interpolation is ordinay linear interpolation"""
iso_array = xr.DataArray(Z, dims=Zdims, coords=Zcoords)
projected_array = xr.DataArray(A, dims="Z", coords=Zcoords)
iso_value = -22.3
proj_val = xroms.xisoslice(iso_array, iso_value, projected_array, "Z")
assert proj_val == projected_array.interp(Z=iso_value)

def test_no_interpolation():
"""Works correctly if the iso_value is already an element in the iso_array"""
iso_array = xr.DataArray(Z, dims="Z", coords=Zcoords)
projected_array = xr.DataArray(A, dims="Z", coords=Zcoords)
iso_value = Z[11]
proj_val = xroms.xisoslice(iso_array, iso_value, projected_array, "Z")
assert proj_val == projected_array[11]

def test_boundary():
"""Bounday values are handled correctly"""
iso_array = xr.DataArray(Z, dims=Zdims, coords=Zcoords)
projected_array = xr.DataArray(A, dims="Z", coords=Zcoords)
# Upper boundary
iso_value = Z[-1]
proj_val = xroms.xisoslice(iso_array, iso_value, projected_array, "Z")
assert proj_val == A[-1]
# Lower boundary
iso_value = Z[0]
proj_val = xroms.xisoslice(iso_array, iso_value, projected_array, "Z")
assert proj_val == A[0]

def test_value_out_of_range():
"""Should return nan if value is out of range"""
iso_array = xr.DataArray(Z, dims="Z", coords=Zcoords)
projected_array = xr.DataArray(A, dims="Z", coords=Zcoords)
iso_value = 3.14 # Too high
proj_val = xroms.xisoslice(iso_array, iso_value, projected_array, "Z")
assert np.isnan(proj_val)
iso_value = -70 # Too low
proj_val = xroms.xisoslice(iso_array, iso_value, projected_array, "Z")
assert np.isnan(proj_val)

def test_non_monotone():
"""Interpolate where unique, not-a-number if not unique"""

iso_array = xr.DataArray(-np.abs(Z + 10), dims=Zdims, coords=Zcoords)
projected_array = xr.DataArray(np.arange(51), dims="Z", coords=Zcoords)
# Unique iso-"surface", index = 40 + iso_value
# other "solution" 40 - iso_value is out of range
iso_value = -10.5
proj_val = xroms.xisoslice(iso_array, iso_value, projected_array, "Z")
assert proj_val == 40 + iso_value
# Non-unique iso-"surface", index = 40+iso_value and 40-iso_value
# Refuse to guess
iso_value = -5.5
proj_val = xroms.xisoslice(iso_array, iso_value, projected_array, "Z")
assert np.isnan(proj_val)

I use this little snippet to construct appropriate chunking dictionary for open_mfdataset

chunk = {"xi": 500, "eta": 120, "ocean_time": 20} 

chunks = {}
for sub in ["rho", "u", "v", "psi"]:
    for k, v in chunk.items():
        chunks[f"{k}_{sub}"] = v
chunks["ocean_time"] = chunk["ocean_time"]
chunks

{'xi_rho': 500,
 'eta_rho': 120,
 'xi_u': 500,
 'eta_u': 120,
 'xi_v': 500,
 'eta_v': 120,
 'xi_psi': 500,
 'eta_psi': 120,
 'ocean_time': 20}

def construct_chunks(xi=None, eta=None, s=None, ocean_time=None):
    pass

chunks= xroms.construct_chunks(xi=500, eta=120, ocean_time=20)
ds = xr.open_dataset(..., chunks=chunks)

I recently upgraded xgcm to ver 0.8.0 from 0.5.2 and xroms.open_mfnetcdf fails. If running 0.8.0 with the below code, the following error returns when horizontally interpolating the vertical coordinate. Not sure which xgcm update resulted in the error appearing.

import glob
import xroms 

path = glob.glob('/d1/shared/TXLA_ROMS/numerical_mixing/non-nest/ver1/1hr/ocean_avg_0000*.nc')
ds_avg = xroms.open_mfnetcdf(path, chunks = {'ocean_time': 1})
---------------------------------------------------------------------------
NotImplementedError                       Traceback (most recent call last)
Input In [1], in <cell line: 5>()
      2 import xroms 
      4 path = glob.glob('/d1/shared/TXLA_ROMS/numerical_mixing/non-nest/ver1/1hr/ocean_avg_0000*.nc')
----> 5 ds_avg = xroms.open_mfnetcdf(path, chunks = {'ocean_time': 1})

File ~/.local/lib/python3.9/site-packages/xroms-0.0.0-py3.9.egg/xroms/xroms.py:815, in open_mfnetcdf(files, chunks, xrargs, Vtransform, add_verts, proj)
    812 ds = xr.open_mfdataset(files, chunks=chunks, **xrargsin)
    814 # modify Dataset with useful ROMS z coords and make xgcm grid operations usable.
--> 815 ds, grid = roms_dataset(ds, Vtransform=Vtransform, add_verts=add_verts, proj=proj)
    817 return ds

File ~/.local/lib/python3.9/site-packages/xroms-0.0.0-py3.9.egg/xroms/xroms.py:194, in roms_dataset(ds, Vtransform, add_verts, proj)
    187 z_w0.attrs = {
    188     "long_name": "depth of W-points",
    189     "field": "z_w0, scalar",
    190     "units": "m",
    191 }
    193 ds.coords["z_w"] = xroms.order(z_w)
--> 194 ds.coords["z_w_u"] = grid.interp(ds.z_w, "X")
    195 ds.coords["z_w_u"].attrs = {
    196     "long_name": "depth of U-points on vertical W grid",
    197     "time": "ocean_time",
    198     "field": "z_w_u, scalar, series",
    199     "units": "m",
    200 }
    201 ds.coords["z_w_v"] = grid.interp(ds.z_w, "Y")

File ~/.conda/envs/xroms/lib/python3.9/site-packages/xgcm/grid.py:2041, in Grid.interp(self, da, axis, **kwargs)
   1989 def interp(self, da, axis, **kwargs):
   1990     """
   1991     Interpolate neighboring points to the intermediate grid point along
   1992     this axis.
   (...)
   2039     >>> grid.interp(da, ["X", "Y"], periodic={"X": True, "Y": False})
   2040     """
-> 2041     return self._1d_grid_ufunc_dispatch("interp", da, axis, **kwargs)

File ~/.conda/envs/xroms/lib/python3.9/site-packages/xgcm/grid.py:1836, in Grid._1d_grid_ufunc_dispatch(self, funcname, data, axis, to, keep_coords, metric_weighted, other_component, **kwargs)
   1833 else:
   1834     map_overlap = False
-> 1836 array = grid_ufunc(
   1837     self,
   1838     array,
   1839     axis=[(ax_name,)],
   1840     keep_coords=keep_coords,
   1841     dask=dask,
   1842     map_overlap=map_overlap,
   1843     other_component=other_component,
   1844     **remaining_kwargs,
   1845 )
   1847 if ax_metric_weighted:
   1848     metric = self.get_metric(array, ax_metric_weighted)

File ~/.conda/envs/xroms/lib/python3.9/site-packages/xgcm/grid_ufunc.py:460, in GridUFunc.__call__(self, grid, axis, *args, **kwargs)
    458 map_overlap = kwargs.pop("map_overlap", self.map_overlap)
    459 pad_before_func = kwargs.pop("pad_before_func", self.pad_before_func)
--> 460 return apply_as_grid_ufunc(
    461     self.ufunc,
    462     *args,
    463     axis=axis,
    464     grid=grid,
    465     signature=self.signature,
    466     boundary_width=self.boundary_width,
    467     boundary=boundary,
    468     dask=dask,
    469     map_overlap=map_overlap,
    470     pad_before_func=pad_before_func,
    471     **kwargs,
    472 )

File ~/.conda/envs/xroms/lib/python3.9/site-packages/xgcm/grid_ufunc.py:745, in apply_as_grid_ufunc(func, axis, grid, signature, boundary_width, boundary, fill_value, keep_coords, dask, map_overlap, pad_before_func, other_component, *args, **kwargs)
    742 # Maybe map function over chunked core dims using dask.array.map_overlap
    743 if map_overlap:
    744     # Disallow situations where shifting axis position would cause chunk size to change
--> 745     _check_if_length_would_change(sig)
    747     mapped_func = _map_func_over_core_dims(
    748         func,
    749         args,
   (...)
    753         out_dtypes,
    754     )
    755 else:

File ~/.conda/envs/xroms/lib/python3.9/site-packages/xgcm/grid_ufunc.py:1002, in _check_if_length_would_change(signature)
    996 all_ax_positions = set(
    997     p
    998     for arg_ps in signature.in_ax_positions + signature.out_ax_positions
    999     for p in arg_ps
   1000 )
   1001 if any(pos in DISALLOWED_OVERLAP_POSITIONS for pos in all_ax_positions):
-> 1002     raise NotImplementedError(
   1003         "Cannot chunk along a core dimension for a grid ufunc which has a signature which "
   1004         f"includes one of the axis positions {DISALLOWED_OVERLAP_POSITIONS}"
   1005     )

NotImplementedError: Cannot chunk along a core dimension for a grid ufunc which has a signature which includes one of the axis positions ['inner', 'outer']

xarray.Dataset
Dimensions:
tracer: 5, s_rho: 30, s_w: 31, eta_rho: 191, xi_rho: 671, xi_u: 670, eta_v: 190, ocean_time: 1080
Coordinates: (31)
Data variables: (165)
Attributes: (32)