Strange error for advection function

[Marston]

Hi,

I'm using a Metpy to calculate the advection field of Absolute Vorticity using the U, V and W fields passed as Xarray DataArrays as seen in the call below. Given the Xarray requirement in the doc I don't need to pass in the dx, dy, or zx. But I still got this error:

adv = mc.advection(av, u=u, v=v, w=w, x_dim=-1, y_dim=-2, vertical_dim=-3)

File ~/anaconda3/envs/contrail/lib/python3.10/site-packages/metpy/calc/tools.py:1312, in _process_deriv_args(f, axis, x, delta)
   1310     delta = np.diff(x, axis=axis)
   1311 else:
-> 1312     raise ValueError('Must specify either "x" or "delta" for value positions.')
   1314 return n, axis, delta

ValueError: Must specify either "x" or "delta" for value positions.

[dopplershift]

Can you change the call to remove the inclusion of x_dim, y_dim, and vertical_dim? Those shouldn't be needed with xarray and might be throwing off the code that's filling in the x/dx. Also, you only need to pass in w if you intend to include vertical advection (not sure if that's your intent or not).

[Marston]

Added the units per your suggestion but still getting all nans in the results:

import numpy as np
import xarray as xr
import metpy.calc as mc
import pandas as pd

def create_da(arr, time, lev, lat, lon, name, units, description):
    ts = xr.Dataset(
        {name: (("time", "lev", "lat", "lon"), arr)},
        coords={
            "time": time,
            "lev": xr.DataArray(lev, dims=('lev',), attrs={'units': 'Pa', 'axis': 'Z'}),
            "lat": xr.DataArray(lat, dims=('lat',), attrs={'units': 'degrees', 'axis': 'Y'}),
            "lon": xr.DataArray(lon, dims=('lon',), attrs={'units': 'degrees', 'axis': 'X'}),
        },
    )
    ts.attrs['units'] = units
    ts.attrs['description'] = description
    return ts.to_array()

levs = np.array([101300, 100000,  97500,  95000,  92500,  90000,  87500,  85000,
                         82500,  80000,  77500,  75000,  72500,  70000,  65000,  60000,
                         55000,  50000,  45000,  40000,  35000,  30000,  25000,  20000,
                         15000,  10000,   7000,   5000,   3000,   2000,   1000])
res = 0.25
lon = np.arange(-180, 179.75+res, res)
lat = np.arange(90, -90-res, -res)
absVort = np.random.uniform(low=-0.00265643, high=0.00315209, size=(1, 31, 721, 1440))  
vWind = np.random.uniform(low=-63.43000031, high=67.23999786, size=(1, 31, 721, 1440))
uWind = np.random.uniform(low=-48.02999878, high=115.94999695, size=(1, 31, 721, 1440))  
time = pd.date_range(start='1/1/2018', periods=1, freq='H')

    
av = create_da(absVort, time, levs, lat, lon, "absvort", "s^-1", "Absolute Vorticity")
v = create_da(vWind, time, levs, lat, lon, 'V', 'm/s', "V Wind")
u = create_da(uWind, time, levs, lat, lon, 'U', 'm/s', "U Wind")

adv = mc.advection(av, u=u, v=v)

[Marston]

<xarray.DataArray ()> <Quantity(nan, '1 / second ** 2')> <xarray.DataArray ()> <Quantity(nan, '1 / second ** 2')>

[dopplershift]

The way you're creating the DataArrays above results in an odd leading variable dimension that seems to be causing some problems (I'm not yet sure why). If I change the code above to:

av = create_da(absVort, time, levs, lat, lon, "absvort", "s^-1", "Absolute Vorticity")[0]
v = create_da(vWind, time, levs, lat, lon, 'V', 'm/s', "V Wind")[0]
u = create_da(uWind, time, levs, lat, lon, 'U', 'm/s', "U Wind")[0]

adv = mc.advection(av, u=u, v=v)

I get non-NaN results. Again, this shouldn't change the results, but for some reason it is.

Regarding the units, the expression np.array([1, 2, 3]) * units.dimensionless results in a Quantity that wraps the original data, "attaching" units to the data. You could accomplish the same thing with units.Quantity(np.array([1,2 , 3]), 'dimensionless'). The important thing from MetPy's perspective is that by turning it into a Quantity, the unit conversion succeeds. Thanks to the "fun" that is angles, 'dimensionless' and 'radians' are synonymous, so essentially we are interpreting all lat/lon coordinates as (when attached to a DataArray) as having units of radians, when our intention is to assume degrees. I'm working on a fix for that.

Thanks for stress-testing our code here! We'll get this cleaned up so everything works more easily.

[Marston]

This is good. I was not aware then to_array() was returning a list. Thanks for catching that. The code calls for the input variables to be DataArrays and it is easier to create Datasets and then convert it to the former. Maybe a suggestion for Metpy to handle both?

I think the units in Metpy is most confusing and the docs could use a little more meat on the bones. Perhaps use this example and expand on it to show how to handle units? Maybe also expanding on generating a Quantity?
Also, I find it difficult to know what units are in Metpy to match the units of the input data. Can a list/table of the available units be out in the docs?

When I add in a vertical motion W is the end product only a vertical advection of the input variable? Just want to make sure.

Metpy reminds me a lot of NCL (a great tool in many regards). I hope Metpy can grow into something like NCL.

I'm curious about the backend of Metpy that drives the calculation. Is it designed for speed? What is under the hood?
Also, one case that I looked at Metpy for and thought it fell short was this: Given a Temperature array A(lev, lat, lon) , where lev is the vertical coordinate in pressure, Pa. I would like to convert to Temperature B(lev, lat, lon) where lev is instead expressed on MSL levels. The input variables are A, geometric height GMH(lev, lat, lon) and M(lev), where M is the vector containing the new coordinates.
I also would like to control this type of interpolation, since in meteorology, not all variables are spatially contiguous. Clouds, for example would need a kind of NN interpolation technique, but Temperature could be a simple linear interpolation.
Is this something Metpy can handle, and I have missed it, or is this something in the works?

[dopplershift]

The code calls for the input variables to be DataArrays and it is easier to create Datasets and then convert it to the former.

Datasets contain one or more DataArrays (indexable by name), so fundamentally they are very different concepts. I'm curious about the "easier to create Datasets" part, though. I can convert your sample code above to create a DataArray directly:

def create_da(arr, time, lev, lat, lon, units, description):
    ts = xr.DataArray(arr,
        coords={
            "time": time,
            "lev": xr.DataArray(lev, dims=('lev',), attrs={'units': 'Pa', 'axis': 'Z'}),
            "lat": xr.DataArray(lat, dims=('lat',), attrs={'units': 'degrees', 'axis': 'Y'}),
            "lon": xr.DataArray(lon, dims=('lon',), attrs={'units': 'degrees', 'axis': 'X'}),
        },
    )
    ts.attrs['units'] = units
    ts.attrs['description'] = description
    return ts

That took minimal changes, produces identical results, and to my eye looks no more difficult than it looked before?

I think the units in Metpy is most confusing and the docs could use a little more meat on the bones. Perhaps use this example and expand on it to show how to handle units? Maybe also expanding on generating a Quantity?

Well some of the issues encountered here are actual bugs that users shouldn't be dealing with. The docs can always be improved, but from your comments I'm not sure exactly what I'd add to them to address any of the challenges you've encountered here. If you have specific additions for the xarray tutorial or units tutorial I'd be happy to see them incorporated to make things better.

Also, I find it difficult to know what units are in Metpy to match the units of the input data. Can a list/table of the available units be out in the docs?

Well the problem is that the full collection is dynamically generated from a combination of the base units and available prefixes. You can see the default set of units over in the repository for Pint, which actually provides the implementation of unit support. Also, if you've done from metpy.units import units, executing list(units) will give you a list of all the units (it's really long!)--I'm not sure how we'd actually present that usefully in our docs, and really if there was a way it would probably belong upstream in Pint's docs.

When I add in a vertical motion W is the end product only a vertical advection of the input variable? Just want to make sure.

If you "add" w, as in adv = mc.advection(av, u=u, v=v, w=w) that will be the total advection of av, including the vertical advection due to w acting on a vertical gradient of av. If you want only the vertical advection, then you'd want to only pass w, like adv = mc.advection(av, w=w).

I'm curious about the backend of Metpy that drives the calculation. Is it designed for speed? What is under the hood?

Well the beauty of using Python is you can see for yourself. Underneath what's written there we rely on the implementation of the array library. This is nominally numpy, which is written in C for speed, but could also potentially be Dask (wrapping Numpy) to facilitate distributed computing on large data. Performance is something we're going to be focusing on a bit more in the near term.

Also, one case that I looked at Metpy for and thought it fell short was this: Given a Temperature array A(lev, lat, lon) , where lev is the vertical coordinate in pressure, Pa. I would like to convert to Temperature B(lev, lat, lon) where lev is instead expressed on MSL levels.

We haven't baked that in since it hasn't been necessary for the workflows we've focused on and no one from our community has contributed an implementation. One way to accomplish it would be to:

1. Use metpy.calc.pressure_to_height_std to convert (usually 1-D) pressure to a height assuming a standard atmosphere
2. Add this height as a new coordinate on your xarray data
3. Use xarray.interp to interpolate your temperature field to a fixed set of height levels

You could go for more detail if you needed it by using instead metpy.calc.thickness_hydrostatic to get a 3D field of height. This does highlight, though, the general theme of MetPy where we want to fill in gaps of domain-specific functionality and then delegate to more general tools elsewhere in the scientific ecosystem to solve general problems. It can make things a little tricky to find sometimes, but it keeps the MetPy project limited in scope and sustainable--we'll never be one stop shopping for everything. The plus side is we get to leverage the "free work" going on elsewhere. 😆

[dopplershift]

So the original all NaN result comes from pydata/xarray#7385, coming out of our use of xr.broadcast(), but caused by the way the arrays are constructed in the sample code. When constructing a Dataset of a single variable and then calling to_array() the returned result isn't a list, but a DataArray that includes an extra variable dimension with a single coordinate that seems to correspond to the original variable. That (for reasons unknown at this point) causes broadcast to put nan in the broadcast field rather than repeating the original values.

We still need to fix in MetPy the issue with "no units lat/lon" -> "lat/lon in radians".

[Marston]

Apologies for this late reply. Thank you for support and the solution for vertical advection. Next sprint I'll be integrating Metpy into our system. I greatly appreciate your support and explanations. For the units, it is good just to give an example of how to find the right argument to Metpy for one's own variable.
With regards to converting a 3D field expressed on pressure levels to a 3D field expressed on MSL levels, I need a solution that doesn't use a standard atmosphere. But rather a solution that uese a supplied a 3D field of geometric height expressed on pressure levels.
I'll continue to explore Metpy and other tools for a working solution. In essence I would like to do for each lat and lon point and interpolation from constant pressure grid to a MSL grid and controlling the method, i.e., linear, NN, etc.

[dopplershift]

I would suggest looking at xarray.interp() for that use case.

[Marston]

Thanks for the tip.
I took a look at it, but it seems too restrictive to work. Maybe I'm missing something, but I'll look into this a bit more later.