precipitable_water with xarray objects throws dimension error in rc2

[DanielAdriaansen]

Something changed from rc1 to rc2 that caused my code calculating precipitable water to break. Here is the traceback:

Traceback (most recent call last):
  File "moist_4panel.py", line 168, in <module>
    pw[x,y] = mpcalc.precipitable_water(ds['P'].isel(latitude=x,longitude=y),ds['TD'].isel(latitude=x,longitude=y),bottom=pbot,top=ptop).m
  File "/d1/anaconda3/envs/era5/lib/python3.8/site-packages/metpy/xarray.py", line 1174, in wrapper
    return wrapping(result, match)
  File "/d1/anaconda3/envs/era5/lib/python3.8/site-packages/metpy/xarray.py", line 1214, in _wrap_output_like_not_matching_units
    xr.DataArray(result, coords=match.coords, dims=match.dims) if output_xarray
  File "/d1/anaconda3/envs/era5/lib/python3.8/site-packages/xarray/core/dataarray.py", line 344, in __init__
    coords, dims = _infer_coords_and_dims(data.shape, coords, dims)
  File "/d1/anaconda3/envs/era5/lib/python3.8/site-packages/xarray/core/dataarray.py", line 121, in _infer_coords_and_dims
    raise ValueError(
ValueError: different number of dimensions on data and dims: 0 vs 1

I did some very light debugging inside xarray.py, and my best guess is that since precipitable water is an integrated quantity from 3D fields, there is a mis-match where the value returned is a single float but the input data have length of num_levels? I am going to try and do some differencing of rc1 vs rc2 xarray.py file and see if that might provide some further insight.

I printed the input and both ds['P'].isel(latitude=x,longitude=y) and ds['TD'].isel(latitude=x,longitude=y) have the same shape and size and bottom and top are just floats as input to precipitable_water().

[DanielAdriaansen]

I am using Xarray 0.16.1 but it looks like the info below is the same for 0.16.2 also.

Digging into the traceback further leads me to Xarray, in dataarray.py from infer_coords_and_dims():

    elif len(dims) != len(shape):
        raise ValueError(
            "different number of dimensions on data "
            "and dims: %s vs %s" % (len(shape), len(dims))
        )

Here, the shape of my precipitable water value (2.5675617212748807 millimeter) must be zero, while the number of dims is 1.

It looks like the type of object that's being passed into infer_coords_and_dims() is:
<class 'pint.quantity.build_quantity_class.<locals>.Quantity'>

dims is ('level',) and shape is empty () on the pint quantity object. Xarray is using "data.shape" to pass into infer_coords_and_dims(), and so either a single value pint quantity object incorrectly returns () instead of "1", or there is no shape attribute that applies to a pint quantity object.

[DanielAdriaansen]

Interesting that when I try to access the .shape attribute of a local pint quantity object in my code, I get this error:
AttributeError: Neither Quantity object nor its magnitude (1000.0) has attribute 'shape'

I am not sure why the Xarray code does not throw that error though, when they are attempting to pass data.shape into infer_coords_and_dims().

[DanielAdriaansen]

I assume I don't get the attribute error from Xarray, because they are doing some massaging of the data prior to calling infer_coords_and_dims() in the two preceeding lines:

            data = _check_data_shape(data, coords, dims)
            data = as_compatible_data(data)

[jthielen]

Thanks for the error report and catching this before 1.0 final!

It looks like I mistakenly set wrap_like='dewpoint' in the xarray decorator for this function in #1490. I think simply removing it so that the decorator is just @preprocess_and_wrap() (since this is function on just 1D profiles at this point, it should only return a Pint scalar) should fix it.

[DanielAdriaansen]

No problem! I verified your suggestion fixed my issue. I will work on a PR shortly.

Could you provide your thoughts on the additional error I am getting from mpcalc.isentropic_interpolation():

Traceback (most recent call last):
  File "isentropic_4panel.py", line 95, in <module>
    iso_anx = mpcalc.isentropic_interpolation(isentlevs,ds.level,ds['T'],ds['Q'],ds['R'],ds['U'],ds['V'],ds['Z'],temperature_out=True)
  File "/d1/anaconda3/envs/era5/lib/python3.8/site-packages/metpy/xarray.py", line 1460, in wrapper
    return func(*bound_args.args, **bound_args.kwargs)
  File "/d1/anaconda3/envs/era5/lib/python3.8/site-packages/metpy/xarray.py", line 1160, in wrapper
    result = func(*bound_args.args, **bound_args.kwargs)
  File "/d1/anaconda3/envs/era5/lib/python3.8/site-packages/metpy/units.py", line 233, in wrapper
    return func(*args, **kwargs)
  File "/d1/anaconda3/envs/era5/lib/python3.8/site-packages/metpy/calc/thermo.py", line 1936, in isentropic_interpolation
    others = interpolate_1d(isentlevels, pres_theta.m, *(arr[sorter] for arr in args),
  File "/d1/anaconda3/envs/era5/lib/python3.8/site-packages/metpy/calc/thermo.py", line 1936, in <genexpr>
    others = interpolate_1d(isentlevels, pres_theta.m, *(arr[sorter] for arr in args),
  File "/d1/anaconda3/envs/era5/lib/python3.8/site-packages/xarray/core/dataarray.py", line 642, in __getitem__
    return self.isel(indexers=self._item_key_to_dict(key))
  File "/d1/anaconda3/envs/era5/lib/python3.8/site-packages/xarray/core/dataarray.py", line 1039, in isel
    ds = self._to_temp_dataset()._isel_fancy(
  File "/d1/anaconda3/envs/era5/lib/python3.8/site-packages/xarray/core/dataset.py", line 2002, in _isel_fancy
    indexers_list = list(self._validate_indexers(indexers, missing_dims))
  File "/d1/anaconda3/envs/era5/lib/python3.8/site-packages/xarray/core/dataset.py", line 1835, in _validate_indexers
    raise IndexError(
IndexError: Unlabeled multi-dimensional array cannot be used for indexing: level

[jthielen]

That would great if you were able to work on a PR for that!

I'm not totally sure what is going wrong with mpcalc.isentropic_interpolation() in that case (I would probably need to see what isentlevs and ds are and experiment with it a little bit), but it appears that some xarray object is leaking into the internals of isentropic_interpolation where it shouldn't be (the DataArrays should be downcast to Pint Quantities by the preprocess_and_wrap wrapper). Would you be able to try using mpcalc.isentropic_interpolation_as_dataset() and see if that makes any difference?

[DanielAdriaansen]

mpcalc.isentropic_interpolation_as_dataset() worked just fine and was able to reproduce my results from rc1. Thanks for the suggestion!

Should I look further into the errors from mpcalc.isentropic_interpolation()? Is it possible Xarray support (i.e. for DataArray objects) was removed from that method in favor of the mpcalc.isentropic_interpolation_as_dataset()?

[jthielen]

Sorry for the delayed follow-up, but yes, it may still be good to look into the errors from isentropic_interpolation, since, while isentropic_interpolation_as_dataset is recommended since it gives a nicer data structure as output, isentropic_interpolation should still work on xarray input (just returning Quantities). isentropic_interpolation_as_dataset working when isentropic_interpolation doesn't may give a clue towards what is wrong.

[DanielAdriaansen]

No problem. I took a deeper look at isentropic_interpolation. For *args, the documentation says this:

args (array, optional) – Any additional variables will be interpolated to each isentropic level.

If args is an array, how would I structure it in the call to isentropic_interpolation? Right now, I simply list them in order like this:
iso_anx = mpcalc.isentropic_interpolation(isentlevs,ds.level,ds['T'],ds['Q'],ds['R'],ds['U'],ds['V'],ds['Z'],temperature_out=True)

However when I look at what the preprocess_and_wrap is doing, it seems like "args" is a single bound argument containing a list or array of DataArray objects that preprocess_and_wrap does not seem to unpack and convert to pint.Quantitiy objects? Am I close here? Seems like preprocess_and_wrap handles the positional arguments OK but when *args contains multiple DataArrays it's not clear to me that's being handled.

[dopplershift]

Just FYI, I'm about to put in a PR to fix precipitable_water. Looking at isentropic_interpolation right now.

[dopplershift]

PR with fixes for both incoming...

Thanks for reporting @DanielAdriaansen. This report, and your testing with the rc's, helped eliminate a couple bugs in 1.0 final!

[DanielAdriaansen]

PR with fixes for both incoming...

Thanks for reporting @DanielAdriaansen. This report, and your testing with the rc's, helped eliminate a couple bugs in 1.0 final!

No problem at all! Glad I was able to at least hone in on the potential problem area with isentropic_interpolation...I just didn't have any idea on whether I was correct or how to solve it. Really appreciate you getting these two fixes in and congrats on the release!

[dopplershift]

Thanks! Your detective work there made it much easier to find and fix.