cannot infer dimension to concatenate

[thomasarsouze]

Hi,
I'm quite interested in your package, so I wanted to give it a try.
I've got a regular ORCA1 simulation, and I just want to load some grid_U and grid_V files to begin with.

data_files
[PosixPath('/esnas/exp/ecearth/a13c/original_files/19500101/fc0/outputs/a13c_1m_20000101_20001231_grid_U.nc'),
 PosixPath('/esnas/exp/ecearth/a13c/original_files/19500101/fc0/outputs/a13c_1m_20010101_20011231_grid_U.nc'),
 PosixPath('/esnas/exp/ecearth/a13c/original_files/19500101/fc0/outputs/a13c_1m_20020101_20021231_grid_U.nc'),
 PosixPath('/esnas/exp/ecearth/a13c/original_files/19500101/fc0/outputs/a13c_1m_20030101_20031231_grid_U.nc'),
 PosixPath('/esnas/exp/ecearth/a13c/original_files/19500101/fc0/outputs/a13c_1m_20040101_20041231_grid_U.nc'),
 PosixPath('/esnas/exp/ecearth/a13c/original_files/19500101/fc0/outputs/a13c_1m_20050101_20051231_grid_U.nc'),
 PosixPath('/esnas/exp/ecearth/a13c/original_files/19500101/fc0/outputs/a13c_1m_20060101_20061231_grid_U.nc'),
 PosixPath('/esnas/exp/ecearth/a13c/original_files/19500101/fc0/outputs/a13c_1m_20070101_20071231_grid_U.nc'),
 PosixPath('/esnas/exp/ecearth/a13c/original_files/19500101/fc0/outputs/a13c_1m_20080101_20081231_grid_U.nc'),
 PosixPath('/esnas/exp/ecearth/a13c/original_files/19500101/fc0/outputs/a13c_1m_20090101_20091231_grid_U.nc'),
 PosixPath('/esnas/exp/ecearth/a13c/original_files/19500101/fc0/outputs/a13c_1m_20000101_20001231_grid_V.nc'),
 PosixPath('/esnas/exp/ecearth/a13c/original_files/19500101/fc0/outputs/a13c_1m_20010101_20011231_grid_V.nc'),
 PosixPath('/esnas/exp/ecearth/a13c/original_files/19500101/fc0/outputs/a13c_1m_20020101_20021231_grid_V.nc'),
 PosixPath('/esnas/exp/ecearth/a13c/original_files/19500101/fc0/outputs/a13c_1m_20030101_20031231_grid_V.nc'),
 PosixPath('/esnas/exp/ecearth/a13c/original_files/19500101/fc0/outputs/a13c_1m_20040101_20041231_grid_V.nc'),
 PosixPath('/esnas/exp/ecearth/a13c/original_files/19500101/fc0/outputs/a13c_1m_20050101_20051231_grid_V.nc'),
 PosixPath('/esnas/exp/ecearth/a13c/original_files/19500101/fc0/outputs/a13c_1m_20060101_20061231_grid_V.nc'),
 PosixPath('/esnas/exp/ecearth/a13c/original_files/19500101/fc0/outputs/a13c_1m_20070101_20071231_grid_V.nc'),
 PosixPath('/esnas/exp/ecearth/a13c/original_files/19500101/fc0/outputs/a13c_1m_20080101_20081231_grid_V.nc'),
 PosixPath('/esnas/exp/ecearth/a13c/original_files/19500101/fc0/outputs/a13c_1m_20090101_20091231_grid_V.nc')]

)

aux_files
[PosixPath('/esnas/releases/models/nemo/v3.6_ecearth/inidata/nemo/ORCA1L75/mesh_mask_nemo.Ec3.2_O1L75.nc')]

netcdf a13c_1m_20000101_20001231_grid_U {
dimensions:
	axis_nbounds = 2 ;
	x = 362 ;
	y = 292 ;
	depthu = 75 ;
	time_counter = UNLIMITED ; // (12 currently)
variables:
	float nav_lat(y, x) ;
	float nav_lon(y, x) ;
	float depthu(depthu) ;
	float depthu_bounds(depthu, axis_nbounds) ;
	float uo(time_counter, depthu, y, x) ;
	double time_centered(time_counter) ;
	double time_centered_bounds(time_counter, axis_nbounds) ;
	double time_counter(time_counter) ;
	double time_counter_bounds(time_counter, axis_nbounds) ;
	float u2o(time_counter, depthu, y, x) ;
	float tauuo(time_counter, y, x) ;
	float umo(time_counter, depthu, y, x) ;
}
netcdf a13c_1m_20000101_20001231_grid_V {
dimensions:
	axis_nbounds = 2 ;
	x = 362 ;
	y = 292 ;
	depthv = 75 ;
	time_counter = UNLIMITED ; // (12 currently)
variables:
	float nav_lat(y, x) ;
	float nav_lon(y, x) ;
	float depthv(depthv) ;
	float depthv_bounds(depthv, axis_nbounds) ;
	float vo(time_counter, depthv, y, x) ;
	double time_centered(time_counter) ;
	double time_centered_bounds(time_counter, axis_nbounds) ;
	double time_counter(time_counter) ;
	double time_counter_bounds(time_counter, axis_nbounds) ;
	float v2o(time_counter, depthv, y, x) ;
	float tauvo(time_counter, y, x) ;
	float vmo(time_counter, depthv, y, x) ;
}

But when I run the load_orca_dataset as in your example :

ds_xorca = load_xorca_dataset(data_files=data_files, aux_files=aux_files,
                              decode_cf=True, 
                              input_ds_chunks={"time_counter": 6, "t": 6,
                                               "z": 23, "Z": 23,
                                               "deptht": 23, "depthu": 23,
                                               "depthv": 23, "depthw": 23,
                                               "y": 511, "Y": 511,
                                               "x": 361, "X": 361},
                              target_ds_chunks={"t": 6,
                                                "z_c": 23, "z_l": 23,
                                                "y_c": 511, "y_r": 511,
                                                "x_c": 360, "x_r": 360})

(actually the chunks don't make sense, but this is not an issue here, and it's the same problem with other chunks or no chunks at all)
then it's failing :

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-29-160e480e9682> in <module>()
     10                                                 "z_c": 23, "z_l": 23,
     11                                                 "y_c": 511, "y_r": 511,
---> 12                                                 "x_c": 360, "x_r": 360})
     13 ##                              input_ds_chunks={"time_counter":12,"depthu":25,"depthv":25,"z":25},
     14 ##                              target_ds_chunks={"t":12,"z_c":25},

~/miniconda3/envs/xorca/lib/python3.6/site-packages/xorca/lib.py in load_xorca_dataset(data_files, aux_files, decode_cf, **kwargs)
    389                     **kwargs),
    390                     data_files, _data_files_chunks)),
--> 391             key=_get_first_time_step_if_any))
    392 
    393     # Add info from aux files

~/miniconda3/envs/xorca/lib/python3.6/site-packages/xarray/core/combine.py in auto_combine(datasets, concat_dim, compat, data_vars, coords)
    434         concatenated = [_auto_concat(ds, dim=dim,
    435                                      data_vars=data_vars, coords=coords)
--> 436                         for ds in grouped]
    437     else:
    438         concatenated = datasets

~/miniconda3/envs/xorca/lib/python3.6/site-packages/xarray/core/combine.py in <listcomp>(.0)
    434         concatenated = [_auto_concat(ds, dim=dim,
    435                                      data_vars=data_vars, coords=coords)
--> 436                         for ds in grouped]
    437     else:
    438         concatenated = datasets

~/miniconda3/envs/xorca/lib/python3.6/site-packages/xarray/core/combine.py in _auto_concat(datasets, dim, data_vars, coords)
    359                                  'concatenate: %s' % concat_dims)
    360             elif len(concat_dims) == 0:
--> 361                 raise ValueError('cannot infer dimension to concatenate: '
    362                                  'supply the ``concat_dim`` argument '
    363                                  'explicitly')

ValueError: cannot infer dimension to concatenate: supply the ``concat_dim`` argument explicitly

I don't see what's different in my test compared to your example...
Thanks in advance.
Thomas.

[jk-rieck]

Hi Thomas,
I had that error before. In my case the problem was that the mesh_mask file did not have the same sizes of dimensions x and y as the data_files. So I guess as a first try you could double-check, that the dimensions that all files share (except for time_counter), have the same sizes.
Cheers,
Jan

[thomasarsouze]

Thanks for your answer.
I've checked all my data_files and my mesh_mask, and everything looks ok as for the dimensions.

netcdf mesh_mask_nemo.Ec3.2_O1L75 {
dimensions:
	x = 362 ;
	y = 292 ;
	z = 75 ;
	t = UNLIMITED ; // (1 currently)
variables:
	float nav_lon(y, x) ;
	float nav_lat(y, x) ;
	float nav_lev(z) ;
	double time_counter(t) ;
	byte tmask(t, z, y, x) ;
	byte umask(t, z, y, x) ;
	byte vmask(t, z, y, x) ;
	byte fmask(t, z, y, x) ;
	byte tmaskutil(t, y, x) ;
	byte umaskutil(t, y, x) ;
	byte vmaskutil(t, y, x) ;
	byte fmaskutil(t, y, x) ;
	float glamt(t, y, x) ;
	float glamu(t, y, x) ;
	float glamv(t, y, x) ;
	float glamf(t, y, x) ;
	float gphit(t, y, x) ;
	float gphiu(t, y, x) ;
	float gphiv(t, y, x) ;
	float gphif(t, y, x) ;
	double e1t(t, y, x) ;
	double e1u(t, y, x) ;
	double e1v(t, y, x) ;
	double e1f(t, y, x) ;
	double e2t(t, y, x) ;
	double e2u(t, y, x) ;
	double e2v(t, y, x) ;
	double e2f(t, y, x) ;
	double ff(t, y, x) ;
	short mbathy(t, y, x) ;
	short misf(t, y, x) ;
	float isfdraft(t, y, x) ;
	double e3t_0(t, z, y, x) ;
	double e3u_0(t, z, y, x) ;
	double e3v_0(t, z, y, x) ;
	double e3w_0(t, z, y, x) ;
	float gdept_0(t, z, y, x) ;
	float gdepu(t, z, y, x) ;
	float gdepv(t, z, y, x) ;
	float gdepw_0(t, z, y, x) ;
	double gdept_1d(t, z) ;
	double gdepw_1d(t, z) ;
	double e3t_1d(t, z) ;
	double e3w_1d(t, z) ;
}

[willirath]

Hi Thomas, is there any chance that we can get the data files to test this ourselves? How big are they? (We could provide an upload link if necessary.)

[thomasarsouze]

One file is 1.1Go and meshmask is 400Mo.
I've reproduced the issue with 4 files (2 in U and 2 in V).
So I guess the overall thing would be ~5Go. I can also try to substract a selection of variables.
If that's fine for you, I can upload that on your server.

[willirath]

(This should have been sent yesterday already...)

We can handle the full data. Can you send me an E-Mail to wrath[at]geomar.de so that I can send you an upload link to our file drop?

[willirath]

@thomasarsouze Some progress:

If I delete all vars relying on the dim axis_nbounds, I can load the data set. This is not a satisfactory fix, because we'd like to be able to work with out-of-the-box Nemo output, but maybe it's a start to get you to play with the package?

[willirath]

With the original data in original_data/ and an empy mod_data/, this should do it (using 4 cores, change the xargs ... -P flag accordingly):

find original_data/a*nc -print0 | \
    xargs -0 -n1 -P4 -I{} bash -c \
        "ncks \
            -x -v depthv_bounds \
            -x -v depthu_bounds \
            -x -v time_centered_bounds \
            -x -v time_counter_bounds \
            {} mod_data/\`basename {}\`"

That's not true, sorry. It just worked, because I was lazy and only used U files in my tests.

[jk-rieck]

Hm, I guessed that axis_nbounds could be the problem. However, I also have some files with dim axis_nbounds and the according variables depending on it, and with these files everything works. @willirath I can give you these files as well, maybe you can figure out why the behavior is different for apparently similar files...

[willirath]

I'm not sure it's about the bounds. I'll have to check this on Wednesday, though.

[thomasarsouze]

Ok thanks for the temporary fix. I'll test that.

[willirath]

There does not seem to be a real fix (I was to quick drawing conclusions from a too simple test.).

What I found is that I could load fine with only the U or with only the V files as data files.

[thomasarsouze]

ok, sorry, I got confused with the modifications indeed...

[willirath]

But maybe the datasets created using only the U files and only the V files could point to a solution.

[willirath]

Solved... The problem was that your files use a different convention for naming the fields than we used internally. We've used vozocrtx for your uo etc. It is possible to override these names and add own variables using a keyword argument. This gist shows a solution: https://gist.github.com/willirath/301225aec5c6776a40be31e0a66e2f73

This also points to the problem that all the test data used to design xorca were from a single group (and hence covered only a subset of the conventions used in the community). I'd love to generalize this to cover more diverse conventions.

[thomasarsouze]

Thanks a lot @willirath ! Can't believe I didn't think about this when I had the orca_names.py file on my screen.
It would be great to be able to update the orca_variables dictionary, at least with CMIP variables names, but also with personal dictionaries. This would probably imply some work though if we want to skip the renaming before the calculation step (checking within a list of possible names for every variables ?) but this would be much more easy to adapt for evey group of NEMO users.
Thanks again.
Thomas.