Cross Section from GFS no x coordinate

When pulling data for GFS, the only direction coordinates that come through are for lon/lat and have no corresponding x, y coordinate values. The following code will reproduce the Attribute error AttributeError: x attribute is not available. I've attempted to create and x and y variable, but have not been able to unlock the magic. This error would also likely persist for NCEP/NCAR reanalysis data as that is a global lat/lon base. Currently, the metpy.parse_cf() does not generate a coordinate attribute since the crs is labelled longitude_latitude. Although interestingly it does not raise an error at that segment of the code, which directly precedes where this error occurs.

from datetime import datetime

import cartopy.crs as ccrs
from metpy.interpolate import cross_section
import xarray as xr

date = datetime(2011,2,1,0)

data = xr.open_dataset('http://www.ncei.noaa.gov/thredds/dodsC/gfs-g4-anl-files/'
                     '{0:%Y%m}/{0:%Y%m%d}/gfsanl_4_{0:%Y%m%d}_{0:%H}00_000.grb2'.format(date))

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

start = (37.0, -105.0)
end = (35.5, -65.0)

cross = cross_section(data, start, end)

@kgoebber

@kgoebber Would you be able to try taking the cross section along a single, known 3D field, such as

cross = cross_section(data['Temperature_isobaric'], start, end)

? Based on testing this out quick myself, it appears that the coordinates are being parsed correctly (see the axis argument on data['lon']), but the dataset contains non-coordinate variables that are not on the lat/lon grid (like data['sigma_layer_bounds']), and so, the cross section cannot be taken on the whole dataset.

@jthielen

Ah! Okay - I think I'm getting it now. Definitely need to add some further examples for Cross Section stuff. I was able to get it from both a single parameter as @jthielen suggested, and also by creating a new xarray Dataset that contains a small subset of the data available from the original file.

I think some further examples would be useful as well (I have some/will have more for the AMS poster on cross sections). Where would be a good location for those once they are ready?

Also, with this sort of problem specifically (supplying a dataset that has some variables without sufficient coordinates), should we be giving a more informative error, or is there some other expected behavior?

I'd think examples here and/or the gallery would probably be best. I think I'd also vote for a more useful error message.

Make necessary coordinate information for cross sections clearer

parent cf7d532

@jthielen

It seems this came up again in #2939. I looked at the error message changes made in this work, but @jthielen comment above here:

the dataset contains non-coordinate variables that are not on the lat/lon grid
still wasn't obvious to me from the error messages raised.

To answer this question raised above also by @jthielen:

Also, with this sort of problem specifically (supplying a dataset that has some variables without sufficient coordinates), should we be giving a more informative error, or is there some other expected behavior?

I think that the more informative error is the way to go here, maybe combined with some tweaks to the existing cross_section() example here: https://unidata.github.io/MetPy/latest/examples/cross_section.html#sphx-glr-examples-cross-section-py?

Since cross_section() requires a DataArray or Dataset, we could make it more explicit for users that it's useful to define a list of variables they would like a cross_section of rather than passing the whole Xarray object:

data = xr.open_dataset(get_test_data('narr_example.nc', False))
data = data.metpy.parse_cf().squeeze()
cross_section_vars = ['isobaric','Temperature','Specific_humidity','u_wind','v_wind']
cross_section = cross_section(data[cross_section_vars], start, end).set_coords(('lat','lon'))

So the question is whether we can do a better check here:

MetPy/metpy/interpolate/slices.py

Lines 156 to 158 in 3a14af1

    
           try: 
        
               crs_data = data.metpy.cartopy_crs 
        
               x = data.metpy.x

Maybe:

try:
  crs_data = data.metpy.cartopy_crs
  if isinstance(data,xr.core.dataarray.DataArray):
    x = data.metpy.x
  else:
    for var in data.data_vars:
      x = data[var].metpy.x
except AttributeError:
  raise 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 am happy to convert to issue and open a PR to add this if it seems useful.

Navigation Menu

Repository navigation

Cross Section from GFS no x coordinate #949

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Cross Section from GFS no x coordinate #949

Description

Activity

jthielen commented on Oct 1, 2018

kgoebber commented on Oct 2, 2018

jthielen commented on Oct 2, 2018

jrleeman commented on Oct 2, 2018

DanielAdriaansen commented on Mar 8, 2023

Metadata

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Issue actions