How to apply MetPy's cross section analysis to a data set with 2-dimensional lat and lon?

Question

I have 4-dimensional data (time, depth, y, and x), but the latitude and longitude are both 2d arrays. y and x are just the indices, so just integers going from 0, 1...end etc. Very similar to the example data set provided by MetPy:

https://unidata.github.io/MetPy/latest/examples/cross_section.html

Unfortunately this is is not the most reproducible, because it's very specific to the data itself. But I am having trouble at the cross section part. I can parse the data according to metpy, but then I get an error when taking a cross section:

data = ds.metpy.parse_cf().squeeze()

print(data)

<xarray.Dataset>
Dimensions:               (y: 1347, x: 1379, deptht: 46, axis_nbounds: 2, time_counter: 12)
Coordinates:
    nav_lat               (y, x) float32 20.92 20.92 20.92 ... 68.49 68.49 68.48
    nav_lon               (y, x) float32 -78.95 -78.9 -78.85 ... -3.614 -3.546
  * deptht                (deptht) float32 3.047 9.454 ... 5.625e+03 5.875e+03
    time_centered         (time_counter) datetime64[ns] 1993-01-16T12:00:00 ....
  * time_counter          (time_counter) datetime64[ns] 1993-01-16T12:00:00 ....
Dimensions without coordinates: y, x, axis_nbounds
Data variables: (12/17)
    deptht_bounds         (deptht, axis_nbounds, y, x) float32 0.0 0.0 ... nan
    time_centered_bounds  (time_counter, axis_nbounds, y, x) datetime64[ns] 1...
    time_counter_bounds   (time_counter, axis_nbounds, y, x) datetime64[ns] 1...
    votemper              (time_counter, deptht, y, x) float32 26.63 ... nan
    vosaline              (time_counter, deptht, y, x) float32 35.88 ... nan
    sosstsst              (time_counter, y, x) float32 26.63 26.58 ... nan nan
    ...                    ...
    sohefldo              (time_counter, y, x) float32 -50.02 -45.85 ... nan nan
    somixhgt              (time_counter, y, x) float32 19.84 19.76 ... nan nan
    sowindsp              (time_counter, y, x) float32 5.591 5.48 ... nan nan
    sohefldp              (time_counter, y, x) float32 nan nan nan ... nan nan
    sowafldp              (time_counter, y, x) float32 3.94e-06 ... nan
    sobowlin              (time_counter, y, x) float32 20.04 20.04 ... nan nan
Attributes:
    name:                      1_VIKING20X.L46-KKG36107B_1d_19930101_19930704...
    description:               ocean T grid variables
    title:                     ocean T grid variables
    Conventions:               CF-1.6
    timeStamp:                 2019-Sep-11 21:02:06 GMT
    uuid:                      a091f081-2943-4c66-aa77-0917f654b802
    history:                   Thu Sep 12 22:15:40 2019: ncrcat -O -F /gfs1/w...
    NCO:                       4.4.8
    nco_openmp_thread_number:  1

Then I try the cross sectional part (just to note as well, the 2D lat and lon is labeled 'nav_lat' and 'nav_lon', instead of lat and lon):

 start = (40, -40)
 end = (50, -30)
    
 cross = cross_section(data, start, end).set_coords(('nav_lon', 'nav_lat'))
 print(cross)

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
~/miniconda3/envs/py3_std_maps/lib/python3.9/site-packages/metpy/interpolate/slices.py in cross_section(data, start, end, steps, interp_type)
    165         try:
--> 166             crs_data = data.metpy.pyproj_crs
    167             x = data.metpy.x

~/miniconda3/envs/py3_std_maps/lib/python3.9/site-packages/metpy/xarray.py in pyproj_crs(self)
    252         """Return the coordinate reference system (CRS) as a pyproj object."""
--> 253         return self.crs.to_pyproj()
    254 

~/miniconda3/envs/py3_std_maps/lib/python3.9/site-packages/metpy/xarray.py in crs(self)
    232             return self._data_array.coords['metpy_crs'].item()
--> 233         raise AttributeError('crs attribute is not available.')
    234 

AttributeError: crs attribute is not available.

During handling of the above exception, another exception occurred:

ValueError                                Traceback (most recent call last)
/tmp/ipykernel_2676704/3597335100.py in <module>
----> 1 cross = cross_section(data, start, end).set_coords(('nav_lon', 'nav_lat'))
      2 print(cross)

~/miniconda3/envs/py3_std_maps/lib/python3.9/site-packages/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/py3_std_maps/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/py3_std_maps/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         }

~/miniconda3/envs/py3_std_maps/lib/python3.9/site-packages/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.

I tried as an alternative, because perhaps parse_cf is causing issues:

data = ds.metpy.assign_latitude_longitude(force=True).squeeze()

but then I still get the same error when applying the cross section.

Any ideas on how to remedy this? Sorry again for the lack of reproducibility, but any ideas would be of great help :) Perhaps it has to do with the projection?

This is an example image of how the data looks at one time and depth instance (looking at temperature):

score 2 · Accepted Answer · answered Oct 30 '21 at 21:42

metpy.interpolate.cross_section requires that your data include both x and y dimension coordinates and the added metpy_crs coordinate (from either parse_cf or assign_crs). In this situation where these x and y dimension coordinates are missing, but you do have 2D latitude and longitude coordinates, these dimension coordinates can be calculated and added with .metpy.assign_y_x() (rather than assign_latitude_longitude which you stated you tried, which does the opposite--adding lat/lon auxillary coordinates from the y/x dimension coordinates).

And so, if your dataset has a valid CF grid mapping corresponding to your data projection, you'll have:

data = ds.metpy.parse_cf()
data = data.metpy.assign_y_x()

start = (40, -40)
end = (50, -30)
    
cross = cross_section(data, start, end).set_coords(('nav_lon', 'nav_lat'))

If it does not:

data = ds.metpy.assign_crs({
    "grid_mapping_name": "lambert_conformal_conic",
    "standard_parallel": 25.0,
    "longitude_of_central_meridian": 265.0,
    "latitude_of_projection_origin": 25.0,
})
data = data.metpy.assign_y_x()

start = (40, -40)
end = (50, -30)
    
cross = cross_section(data, start, end).set_coords(('nav_lon', 'nav_lat'))

(changing the projection attributes to those needed for your given projection)

Thank you so much for all the help! Yes... unfortunately I am using a model which has a nested grid, so I think the geometry might differ from valid grid mapping as I am still getting an error saying crs attribute is not available. I really need to look into the geometry of the model haha, but this is extremely helpful! :) — wabash, Oct 31 '21 at 09:03

How to apply MetPy's cross section analysis to a data set with 2-dimensional lat and lon?

1 Answers1