Cross-Section from HRRR output

[atmosmattbrewer]

Hello,
I have been trying to use MetPy's Cross-section functionality on HRRR model output and have run into problems with handling the coordinates and dimensions. I asked a stack overflow question, Metpy HRRR Cross Section, and was directed here. I have looked into one of the other issues related to this topic, #949, but the I still haven't been able to get fix my issue.

The data file I am using is HRRR model output transformed from grib2 to netcdf using xarray and pynio.
Here is what my dataset looks like after I parse it with metpy.parse_cf()

When try to do the cross-section this is the errors I get:

So I tried parsing one variable as suggested,
RH = data.metpy.parse_cf('RH_P0_L100_GLC0')
x, y = RH.metpy.coordinates('x', 'y')
and got this error code instead:

Is there extra xarray manipulation that I am missing to get the dimensions and coordinates to be correct?

[jthielen]

@dopplershift and @jrleeman: The general problem at play here seems to be that the dimensions of the data are not present as coordinates, and having x and y dimension coordinates is required for the current implementation of cross section interpolation to work. This isn't the first time I've seen this problem come up with dimensions lacking coordinates in relation to cross sections (WRF has been notorious for it), so I'd think MetPy should probably have something to handle this. I had two ideas on how to resolve this (but neither seems ideal, so other suggestions would be great):

1. Alter the cross section interpolation to work with non-coordinate x and y dimensions so long as we have 2D lat and lon coordinates in the DataArray. However, making the indexing for this work may get messy, since DataArray.interp() only accepts arguments that are dimensions, not arbitrary coordinates.
2. Programmatically fill in the missing coordinates using latitude, longitude, and projection information. This is what I've done as a workaround when I've ran into this issue myself. However, this may not reliably work with the current parse_cf(), since the crs will come up as latitude_longitude if the grid_mapping attribute is missing but 2D lat and lon coordinates are still found, which appears to be the case with @atmosmattbrewer's file here. It also may require making unjustified assumptions about the data.

[jthielen]

Also, @atmosmattbrewer, to directly answer your last question, I don't believe there is any extra xarray manipulation that you are missing; it is simply a problem with MetPy's cross_section() not being able to handle your dataset at the moment. I've done workarounds for a similar problem with WRF output, but they required full projection and grid spacing information from the dataset, which I'm not sure yours has or not.

[sgdecker]

@jthielen : Since WRF output is a common analysis target, could you post the code you use to fill in the missing x and y coordinates? Perhaps this (or a suitably generalized variant) should be a MetPy helper function users can call when they know they have an Xarray Dataset that originated from WRF output.

[jthielen]

As requested, here is the code I used:

https://gist.github.com/jthielen/8881d32d08d625c75c906a7e8ad7583f

This issue with WRF, however, relates to something I've been wondering about for a little while now. I have some of these scripts built (and more I'd like to create) to make working with WRF files easier (a more pythonic API than what wrf-python currently has and having CF-compliance so it can play nice with MetPy). Salem came close to this already, but it looks like they are deprecating their WRF support soon in favor of wrf-python. And so, what would be the best place to contribute these efforts?

• MetPy? (Does providing this extra support for WRF open too much of a can of worms?)
• wrf-python? (Seems like the location where effort is converging, but this might be too far afield from their current approach)
• A new package?

[jthielen]

Also, tagging @tjwixtrom, since we talked at AMS about this and he has done some similar work.

[atmosmattbrewer]

Has anyone had any success in using the metpy cross_section for HRRR or GFS grib data yet? I am now coming back to this problem and still have not had any success in making a cross section. Still pretty much running into the same issues as mentioned above and haven't found any examples any where online.

[jthielen]

@atmosmattbrewer Now with MetPy's assign_crs and assign_y_x helpers (as of version 1.0, which is still unreleased, and since RC1 is somewhat incomplete, I'd recommend using current master branch), everything should be ready to go for the option 2 referenced above (#1004 (comment)).

Would you have a sample grib2 dataset that you'd be willing to share, along with its grid projection parameters? I can take a look and put together a quick demo based off of that.

[atmosmattbrewer]

@jthielen Right now just for testing purposes I have been using the daily native level HRRR grib2 output (https://nomads.ncep.noaa.gov/pub/data/nccf/com/hrrr/prod/hrrr.20200805/conus/)

The projection code I use is:
lon0=-97.5 lat0=38.5 lat1=38.5 proj = ccrs.LambertConformal(central_longitude=lon0, central_latitude=lat0, standard_parallels=(lat0,lat1))
Got this from the namelist.

[jthielen]

@atmosmattbrewer Here's the demo: https://nbviewer.jupyter.org/urls/dl.dropbox.com/s/ymmj7nrznbqdicp/hrrr_cross_section_demo.ipynb.

[dopplershift]

I'd be interesting in including something like this in the docs if it takes care of it (unless the data file is unreasonable which case maybe python-training.)

We should also look at whatever we can to improve the docs to cross_section to convey this.

[jthielen]

This particular example with HRRR is a 519 MB grib2 file, so probably not best to include in the docs or maybe even python-training. I have some post-processed FV3-LAM output on a coarser (25km) version of the same grid that could probably function the same way for such a demo, and it is only 13 MB.

I'll also be sure to keep assign_y_x / cross_section stuff this in mind with the "MetPy 1.0 and xarray" doc updates...as soon as I can even get the xarray compat stuff done...

[dopplershift]

Probably should be some combination of at least:

• Documenting all metadata requirements for cross_section
• Putting this in narrative documentation somewhere that's easily discoverable (maybe xarray tutorial though I have concerns that it's actually used)
• Adding an example, maybe using live data, to Python Training.

[fmaussion]

I just discovered this issue - thanks for tagging the salem issue @jthielen .

Yes, I do believe that something better than wrf-python and separate from MetPy should exist, I personally don't have the bandwith for this.

I still believe that wrf-python's API was not well designed because they seem to just have translated their NCL code to python. Now they have a lot of legacy to deal with, and on top of that the library's development does not look very active. I seems that a strong political decision at NCAR or Unidata is needed here to put more effort into python for WRF tools.

(sorry for the out-of-topic rant)

[jthielen]

I just discovered this issue - thanks for tagging the salem issue @jthielen .

Yes, I do believe that something better than wrf-python and separate from MetPy should exist, I personally don't have the bandwith for this.

I still believe that wrf-python's API was not well designed because they seem to just have translated their NCL code to python. Now they have a lot of legacy to deal with, and on top of that the library's development does not look very active. I seems that a strong political decision at NCAR or Unidata is needed here to put more effort into python for WRF tools.

(sorry for the out-of-topic rant)

Thanks for the reply! After I gave a presentation at the AMS annual meeting back in January, there was some informal support from the NCAR GeoCAT team for an overhaul of wrf-python to be more pythonic and integrate with the pangeo ecosystem, but effort would have to come from the community in the immediate future, since the GeoCAT team's limited resources were prioritized elsewhere. I had hoped to spearhead that effort, but other research priorities (and the pandemic) got in the way since then. It is still something I'd like to see happen, but I unfortunately can't prioritize it myself either right now. Though, I'd be glad to discuss ideas with anyone else who wants to move it forward!

[tjwixtrom]

I have an initial stab at this which can be found here. The basic idea is to wrap xarray's open_dataset by opening the input file, fixing metadata, and providing some basic diagnostic computations. Arrays are held in Dask, allowing for out-of-memory and delayed operations, so the code seems to be relatively efficient. Computations are currently mainly either vectorized adaptation of RIP/wrf-python or a modification of MetPy. I'm currently waiting for a few PRs to be merged in MetPy but the idea is to depend on MetPy for general computations and handle the model specific routines within the wrf-ens-tools code. We are currently very early in the development for this so input is definitely welcomed.

[fmaussion]

Great talk! salem is still attracting some users disappointed with wrf-python btw. I would much rather have them go to a better supported package...

@tjwixtrom great! There are some ideas about that in salem as well: https://github.com/fmaussion/salem/blob/697762b3b8daf195119889959535220726b12248/salem/sio.py#L951 - one of the fundamental ideas was to unstagger and add diagnostic variables on the go by mimicking netcdf variables for xarray. Maybe this is also what you are doing? In all cases I wish your efforts to take off, it would be great to help the community to converge.

[mvoss-wx]

I've tried to follow this thread, but I'm stuck. I asked Unidata support for help and they said they were able to create a cross section with my data, which perhaps that tells you about my coding skills. Anyway, as university faculty, I would like to use this great tool in my mesoscale meteorology course if anyone can point me in the right direction.

I'm trying to plot a hrrr model cross section. I can use the sample metpy code
to plot the NARR example, and I can modify this to plot GFS data. But HRRR is non CF compliant and when I try to assign crs manually, it seems work, but then gives me a strange projection. I think I'm missing something obvious here.

My code
sample ouput file from this code.
My hrrrr data file

Since this is the HRRR model and my data set is the entire domain, it's clear the data lat/lon is shifted north and east, so it seems like the assign crs is not working as intended.

Any advice would be appreciated.

[dcamron]

Just want to update here in case it helps out anyone else! Thanks @mvoss-wx for reaching out via eSupport as well.

The main issue was in making sure the projection definition lined up with the provided coordinates. In this case, using the operational HRRR namelist to define the appropriate grid_mapping parameters (as you've done!), you also have to make account for the offset in x and y given ref_x and ref_y. The grid definition here relies on the data origin being in a certin place. Once you do this, the rest of that example falls into place. Here's an example using that provided file. Note that units, offsets, projection definitions, etc. may differ depending on where you obtain/process your data.

[mvoss-wx]

@dcamron , Thanks so much for checking this out. I tested and this works perfectly for the hrrr isobaric files.
Thanks for adding the helpful comments in the code as well. In retrospect, I see I needed to dig one layer deeper on a few of those features to actually understand what is going on.
Thanks again, this is much appreciated!

[meihan33]

Hi, dcamron, mvoss-wx, and others,

I am trying to make a cross section from a HRRR pressure level file I got from NCEP AWS. I figured out to use cfgrib engine to directly read the grib2 file by xarray. However, I don't get the data structure like what's in your example, where you have nc data file with y, x, plevel as dimensions and dimension coordinates. May I ask how to get this type of nc data for HRRR and with these dimension coordinates?

I just have limited knowledge with xarray or metpy. Please see the screenshot, which shows the data structure I got:

Thank you for any suggestions you may have.

[dopplershift]

@meihan33 You can still use the approach taken in the example notebook that @jthielen posted above. I think this will work:

cf_attrs = {
    'grid_mapping_name': 'lambert_conformal_conic',
    'standard_parallel': 38.5,
    'longitude_of_central_meridian': -97.5,
    'latitude_of_projection_origin': 38.5
}
ds = ds.metpy.assign_crs(cf_attrs).metpy.assign_y_x()

[meihan33]

@dopplershift Thank you very much for your suggestion. I added one more key and value pair in the dictionary:
'earth_radius': 6371229.0
It works now.

There are only some warnings left for each variable, it shows:
/miniconda3/lib/python3.8/site-packages/metpy/xarray.py:355: UserWarning: More than one time coordinate present for variable "gh".
warnings.warn('More than one ' + axis + ' coordinate present for variable'

I could probably find a way to drop these extra coordinate or axis. But the information may be useful actually, so I would just keep them by ignoring a page of warning messages.

Thank you for any further suggestions you may have.