Add geopotential height to geopotential function?

[sgdecker]

I thought computing geopotential would be as simple as g*z (where z is geopotential height from model output), but apparently not:

from metpy.constants import g

# Code obtaining the Dataset of model output elided

p = 300*units('hPa')
z = nam['Geopotential_height_isobaric'].metpy.sel(vertical=p)
print(z)
print(g*z)

gives

<xarray.DataArray 'Geopotential_height_isobaric' (time2: 1, y: 129, x: 185)>
array([[[9725.858, 9726.138, ..., 9652.378, 9651.778],
        [9726.538, 9726.858, ..., 9648.658, 9647.818],
        ...,
        [8777.778, 8783.019, ..., 8703.098, 8693.258],
        [8767.658, 8772.738, ..., 8686.538, 8677.578]]], dtype=float32)
Coordinates:
    reftime    datetime64[ns] 2024-03-05T12:00:00
  * time2      (time2) datetime64[ns] 2024-03-05T18:00:00
    isobaric2  float64 3e+04
  * x          (x) float32 -4.226e+06 -4.185e+06 ... 3.21e+06 3.251e+06
  * y          (y) float32 -8.327e+05 -7.921e+05 ... 4.328e+06 4.369e+06
    metpy_crs  object Projection: lambert_conformal_conic
Attributes:
    long_name:                       Geopotential height @ Isobaric surface
    units:                           gpm
    Grib_Variable_Id:                VAR_0-3-5_L100
    Grib2_Parameter:                 [0 3 5]
    Grib2_Parameter_Discipline:      Meteorological products
    Grib2_Parameter_Category:        Mass
    Grib2_Parameter_Name:            Geopotential height
    Grib2_Level_Type:                100
    Grib2_Level_Desc:                Isobaric surface
    Grib2_Generating_Process_Type:   Forecast
    Grib2_Statistical_Process_Type:  UnknownStatType--1
    grid_mapping:                    LambertConformal_129X185-40p79N-100p33W
<xarray.DataArray 'Geopotential_height_isobaric' (time2: 1, y: 129, x: 185)>
<Quantity([[[95378.09  95380.83  95382.8   ... 94664.17  94657.49  94651.62 ]
  [95384.76  95387.9   95389.86  ... 94624.55  94621.016 94612.78 ]
  [95389.86  95393.38  95396.92  ... 94588.45  94583.36  94576.69 ]
  ...
  [86181.02  86234.37  86308.12  ... 85622.44  85526.33  85421.984]
  [86080.6   86131.99  86204.555 ... 85438.06  85348.234 85251.734]
  [85981.36  86031.17  86100.61  ... 85268.22  85185.84  85097.98 ]]], 'meter / second ** 2')>
Coordinates:
    reftime    datetime64[ns] 2024-03-05T12:00:00
  * time2      (time2) datetime64[ns] 2024-03-05T18:00:00
    isobaric2  float64 3e+04
  * x          (x) float32 -4.226e+06 -4.185e+06 ... 3.21e+06 3.251e+06
  * y          (y) float32 -8.327e+05 -7.921e+05 ... 4.328e+06 4.369e+06
    metpy_crs  object Projection: lambert_conformal_conic

where the units of g*z are wrong. It's possible I'm doing something wrong, but I'm guessing it is the "gpm" units that are confusing things?

Since it is apparently not so simple, perhaps there is value in adding a geopotential_height_to_geopotential function that returns the proper units?

I do see this comment, but as far as I can tell, this function does not yet exist.

[sgdecker]

OK, I see if I add .metpy.quantify() to the end of my definition of z, it does the right thing.
I'm not sure why I can't do .metpy.quantify() to my entire NAM dataset though. (I get ValueError: Unit expression cannot have a scaling factor.)

[sgdecker]

Yes, there is a suspicious array for sure:

print(nam['Soil_type_surface'])

<xarray.DataArray 'Soil_type_surface' (time2: 1, y: 129, x: 185)>
array([[[0., 0., ..., 0., 0.],
        [0., 0., ..., 0., 0.],
        ...,
        [0., 0., ..., 0., 0.],
        [0., 0., ..., 0., 0.]]], dtype=float32)
Coordinates:
    reftime    datetime64[ns] 2024-03-05T12:00:00
  * time2      (time2) datetime64[ns] 2024-03-05T18:00:00
  * x          (x) float32 -4.226e+06 -4.185e+06 ... 3.21e+06 3.251e+06
  * y          (y) float32 -8.327e+05 -7.921e+05 ... 4.328e+06 4.369e+06
    metpy_crs  object Projection: lambert_conformal_conic
    latitude   (y, x) float64 12.19 12.29 12.39 12.49 ... 57.49 57.39 57.29
    longitude  (y, x) float64 -133.5 -133.1 -132.8 ... -50.45 -49.92 -49.39
Attributes:
    long_name:                       Soil type @ Ground or water surface
    units:                           (Code table 4.213)
    Grib_Variable_Id:                VAR_2-3-0_L1
    Grib2_Parameter:                 [2 3 0]
    Grib2_Parameter_Discipline:      Land surface products
    Grib2_Parameter_Category:        Soil Products
    Grib2_Parameter_Name:            Soil type
    Grib2_Level_Type:                1
    Grib2_Level_Desc:                Ground or water surface
    Grib2_Generating_Process_Type:   Forecast
    Grib2_Statistical_Process_Type:  UnknownStatType--1
    grid_mapping:                    LambertConformal_129X185-40p79N-100p33W

I suspect there are others, but didn't look further. I wonder if there is a way to make quantify skip things with what I assume are "categorical" units (i.e., issue warnings for variables like this but not crash). Or is it simply a TDS bug to produce a unit like this in the first place?

[dopplershift]

THREDDS possibly could consider not putting that in the units (I'm pretty sure that's directly coming out of a GRIB table). We could be a bit more permissive when reading data we don't control.

Regardless, calling .metpy.quantify() on that dataset seems like doing needless work. You're doing work with, what, 6-10 fields? That call is going to needlessly download and parse the metadata for 50-100 variables in that NAM dataset, and opens you up to potential problems like you hit above.

[dopplershift]

I should also add, that if you're passing this data directly to MetPy, you shouldn't ever need to call .metpy.quantify() manually, so maybe just eliminating that call should make your life easier all around?

[sgdecker]

I only added the call to .metpy.quantify() so that I would get the correct units from g*z, so it seems like that is a case where it is indeed necessary to make that call. Are you saying if I were to pass the result of g*z to a MetPy function, suddenly the units would be correct? Or maybe since g*z is done before passing anything to a MetPy function, this is no longer "direct", so your statement above doesn't apply?

My thinking for trying to do .metpy.quantify() to the entire dataset was to eliminate any other unit gremlins that might be lurking, so that I wouldn't need to call .metpy.quantify() for each gremlin, but just do it once and for all. But maybe g*z is the only gremlin.

Since I agree that it would be good to hide/eliminate any manual .metpy.quantify() usages, it does seem like adding a MetPy function for this calculation is the way to go as you indicate below.

[dopplershift]

Ah. Yeah, if you're doing math, .metpy.quantify() is necessary since DataArrays don't know anything about units and will dump pretty much all metadata in creating a new DataArray for the results--and then there would be no unit information for later use (though you could add it back manually).

I'm trying to decide if the Dataset-level calls were intentionally chosen to error, or if that's just a consequence of looping over the calls.

[dopplershift]

Regarding the original question, would be happy to merge a PR adding the function.