What is the right way to calculate gradient in Metpy ?

[Dearsomeone]

Hello, I have found two gradient functions in metpy, i.e. mpcalc.gradient and mpcalc.geospatial_gradient. In the past, I used the NCL function grad_latlon_cfd to calculate gradient of gridded model data, such as SST and PSL. Take the CMIP6 model variable PSL as an example, which is on a 2°x2° regular lat-lon grid, I wonder which metpy function is suitable for calculate the gradient. It seems that they can both be successfully executed, but their results are different.

import metpy.calc as mpcalc

# function 1
grad1 = mpcalc.gradient(psl)

# function 2
dx, dy = mpcalc.lat_lon_grid_deltas(ds.lon, ds.lat)
grad2 = mpcalc.geospatial_gradient(psl, dx=dx, dy=dy)

[dopplershift]

metpy.calc.gradient is just doing normal finite-difference-based approximation to the derivative, assuming the data are on a Cartesian grid. Effects of calculating on a spherical Earth are only handled by computing grid spacing using "great circle" distance between lat/lon points.

metpy.calc.geospatial_gradient explicit "corrects" plain Cartesian gradient for working in a spherical coordinate system with explicit factors like GEMPAK did. Note, you shouldn't need to pass dx/dy manually and should be able to do directly grad2 = mpcalc.geospatial_gradient(psl).

geospatial_gradient is the most correct calculation, but requires the input data to have good information/metadata regarding the coordinate reference system. I'm not sure whether grad_latlon_cfd is implemented using any of these corrections.

With that said, both calculations should produce very similar results up to 30-40 degrees latitude, with increasing differences towards the poles. If that's not the behavior you observe, we should dig in further to see what's going on and make sure there's not a bug lurking there.

[Dearsomeone]

Hello, thank you very much for your prompt reply! Here I would like to make a suggestion, I find from the source code of mpcalc.advection that：

def advection():

  return -sum(
          wind * gradient
          for wind, gradient in zip(wind_vector.values(), gradient_vector)
      )

so the formula of temperature advection is
$$advcect = -(u\frac{\partial T}{\partial x}+v\frac{\partial T}{\partial y}+w\frac{\partial T}{\partial z})$$
As a former NCL user, the documentaion of NCL function advect_variable will show the calculation formula:
$$advX = U*(dX/dlon) + V*(dX/dlat)$$
The return value sign of the two functions are opposite. So I suggest that metpy should also have an formula explanation in the function API reference.