Trouble with cross_section and 1D lat/lon

[dopplershift]

I'm not quite sure why this doesn't work:

qt = xr.open_dataset('ta.all.ERAINT_MAM2018_region.nc').metpy.parse_cf().squeeze()
qh = xr.open_dataset('q.all.ERAINT_MAM2018_region.nc').metpy.parse_cf().squeeze()
start1 = (-4, 22)
end2 = (4, 42)
crosst = cross_section(qt, start1, end2)
crossq = cross_section(qh, start1, end2)
rh = mpcalc.relative_humidity_from_specific_humidity(qt['level'], qt['ta'], qh['q'])

I'm getting:

Found valid latitude/longitude coordinates, assuming latitude_longitude for projection grid_mapping variable
Found valid latitude/longitude coordinates, assuming latitude_longitude for projection grid_mapping variable
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
~/repos/metpy/src/metpy/interpolate/slices.py in cross_section(data, start, end, steps, interp_type)
    166             crs_data = data.metpy.pyproj_crs
--> 167             x = data.metpy.x
    168         except AttributeError:

~/repos/metpy/src/metpy/xarray.py in x(self)
    430         """Return the x coordinate."""
--> 431         return self._axis('x')
    432 

~/repos/metpy/src/metpy/xarray.py in _axis(self, axis)
    384             else:
--> 385                 raise AttributeError(axis + ' attribute is not available.')
    386         else:

AttributeError: x attribute is not available.

During handling of the above exception, another exception occurred:

ValueError                                Traceback (most recent call last)
/var/folders/r8/bhrhn4m56lgdckj491rm6_w40000gp/T/ipykernel_74405/515362755.py in <module>
      3 start1 = (-4, 22)
      4 end2 = (4, 42)
----> 5 crosst = cross_section(qt, start1, end2)
      6 crossq = cross_section(qh, start1, end2)
      7 rh = mpcalc.relative_humidity_from_specific_humidity(qt['level'], qt['ta'], qh['q'])

~/repos/metpy/src/metpy/interpolate/slices.py in cross_section(data, start, end, steps, interp_type)
    154     if isinstance(data, xr.Dataset):
    155         # Recursively apply to dataset
--> 156         return data.map(cross_section, True, (start, end), steps=steps,
    157                         interp_type=interp_type)
    158     elif data.ndim == 0:

~/miniconda3/envs/test-cartopy/lib/python3.9/site-packages/xarray/core/dataset.py in map(self, func, keep_attrs, args, **kwargs)
   5106         if keep_attrs is None:
   5107             keep_attrs = _get_keep_attrs(default=False)
-> 5108         variables = {
   5109             k: maybe_wrap_array(v, func(v, *args, **kwargs))
   5110             for k, v in self.data_vars.items()

~/miniconda3/envs/test-cartopy/lib/python3.9/site-packages/xarray/core/dataset.py in <dictcomp>(.0)
   5107             keep_attrs = _get_keep_attrs(default=False)
   5108         variables = {
-> 5109             k: maybe_wrap_array(v, func(v, *args, **kwargs))
   5110             for k, v in self.data_vars.items()
   5111         }

~/repos/metpy/src/metpy/interpolate/slices.py in cross_section(data, start, end, steps, interp_type)
    167             x = data.metpy.x
    168         except AttributeError:
--> 169             raise ValueError('Data missing required coordinate information. Verify that '
    170                              'your data have been parsed by MetPy with proper x and y '
    171                              'dimension coordinates and added crs coordinate of the '

ValueError: Data missing required coordinate information. Verify that your data have been parsed by MetPy with proper x and y dimension coordinates and added crs coordinate of the correct projection for each variable.

Here are the data files.

Thoughts @jthielen? (cc: @kimutaijoy)

[jthielen]

The issue I see with these data are the time_bnds variables, which are not three-dimensional (or higher) arrays with horizontal dimension coordinates, and so you can't take a cross section through them (a requirement stated under the data parameter in the docstring for cross_section).

Luckily, there are a couple easy fixes! I'd consider time_bnds to be more of an auxiliary coordinate than a data variable, so I'd change this to

qt = xr.open_dataset('ta.all.ERAINT_MAM2018_region.nc').metpy.parse_cf().squeeze().set_coords('time_bnds')
qh = xr.open_dataset('q.all.ERAINT_MAM2018_region.nc').metpy.parse_cf().squeeze().set_coords('time_bnds')
start1 = (-4, 22)
end2 = (4, 42)
crosst = cross_section(qt, start1, end2)
crossq = cross_section(qh, start1, end2)

Alternatively, you could just drop time_bnds or only take the cross section of the DataArray of interest from each Dataset.

[dopplershift]

Would it make sense when processing a Dataset to ignore variables that don't work when we have some that will? Warning of course in that case (which could be eliminated by the workarounds you mentioned).

[jthielen]

I could see it either way.

In most cases, yes, we would just have to either manually merge instead of using map to handle the "doesn't work" case or be okay with cross_section being a no-op on a DataArray with missing coordinates (although that latter option could cause more confusion). Though, if somehow a dataset has a variable with partial horizontal coordinate data (i.e., only y or x but not both), we won't be able to get a merged result back out and will have to fail entirely anyways.

[dopplershift]

Oh, this is coming from map. That makes sense. Is there any way we could avoid mapping over DataArrays that don't make sense (e.g. are 2D)?

[jthielen]

By "avoid" do you mean drop them from the result? That should be pretty easy. Or do you mean ignoring and passing them through to the final merge? That gets into the possible issues mentioned above.

[dopplershift]

Either from the final merge or up front to avoid passing to map. (I may also not be fully understanding all the plumbing.) I'm just thinking if the result of trying something was always going to be non-sensical, and we know it, might was automatically not do it.