Repository navigation
Laplacian and Second Derivative #1389
Description
Problem:
Trying to reproduce the output from GEMPAK for the laplacian calculation. MetPy has a different implementation of calculating the Laplacian by using the sum of the second derivative of the values. GEMPAK takes the divergence of the gradient of the values.
If I use the gradient and divergence calculations from MetPy, I get a solution that correlates very well (0.999), but is slightly off calculation wise. If I use the Laplacian (or natively use the second_derivative function), then I get an answer that doesn't correlate nearly as well (0.73).
Interestingly, if I use the first_derivative twice to get a second derivative, I get the same solution as if I were using the gradient and divergence functions. Not surprising since both of those functions use first_derivative under the hood. But using the first_derivative twice does not give me the same solution as using second_derivative.
When I use made up 1D or 2D data, where I would know the answer to the first and second derivative, they all work out fine.
Examples below to play with.
Example 1
import metpy.calc as mpcalc
import numpy as np
y = -2*np.arange(0, 5)**2 + 5
print(mpcalc.second_derivative(y, delta=1))
print(mpcalc.first_derivative(mpcalc.first_derivative(y, delta=1), delta=1))Output:
[-4.0 -4.0 -4.0 -4.0 -4.0] dimensionless
[-4.0 -4.0 -4.0 -4.0 -4.0] dimensionless
Example 2
import metpy.calc as mpcalc
import numpy as np
x = np.arange(-5, 6)
y = np.arange(-5, 6)
xx, yy = np.meshgrid(x, y)
f = xx**2 + yy**2
print(mpcalc.second_derivative(f, delta=1))
print(mpcalc.first_derivative(mpcalc.first_derivative(f, delta=1), delta=1))Output:
[[2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0] [2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0] [2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0] [2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0] [2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0] [2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0] [2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0] [2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0] [2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0] [2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0] [2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0]] dimensionless
[[2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0] [2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0] [2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0] [2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0] [2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0] [2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0] [2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0] [2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0] [2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0] [2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0] [2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0]] dimensionless
Example 3
from datetime import datetime
import metpy.calc as mpcalc
import numpy as np
import xarray as xr
date = datetime(2020, 6, 1, 0)
gfs_data = xr.open_dataset('https://thredds.ucar.edu/thredds/dodsC/grib/NCEP/GFS/Global_onedeg_ana/'
f'GFS_Global_onedeg_ana_{date:%Y%m%d_%H}00.grib2').metpy.parse_cf()
data_var_hght = gfs_data.Geopotential_height_isobaric.metpy.sel(time=date, vertical=500 * units.hPa).metpy.quantify()
data_var_lat = gfs_data.lat
data_var_lon = gfs_data.lon
dx, dy = mpcalc.lat_lon_grid_deltas(data_var_lon, data_var_lat)
data_test = mpcalc.laplacian(data_var_hght)
dzdy, dzdx = mpcalc.gradient(data_var_hght, deltas=(dy, dx))
data_test2 = mpcalc.divergence(dzdx, dzdy, dx, dy)
print(data_test.values)
print(data_test2.m)Output (differences are not insubstantial):
[[ nan nan nan ... nan
nan nan]
[-5.50590864e-09 -5.50907902e-09 5.15984316e-09 ... -5.63441096e-09
5.03455036e-09 5.03294556e-09]
[ 2.60302066e-09 2.60458607e-09 -1.29180946e-09 ... 1.31759610e-09
-1.35130490e-09 -1.35126594e-09]
...
[-1.94440638e-09 -1.95399582e-09 -6.44458985e-10 ... -1.92992423e-09
-5.83437519e-10 -5.88291002e-10]
[-5.74953538e-09 -5.75114012e-09 -4.82535019e-10 ... -5.73834144e-09
-4.74550480e-10 -4.79325510e-10]
[ nan nan nan ... nan
nan nan]]
[[ nan nan nan ... nan
nan nan]
[-8.18597714e-09 -2.80715478e-09 1.18025786e-09 ... -1.54098051e-09
2.44806115e-09 7.82768321e-09]
[ 3.61531068e-09 1.65605756e-09 -3.12910951e-10 ... 1.55354859e-23
-6.86196855e-10 -2.02896995e-09]
...
[-2.28371231e-09 -1.62142904e-09 -9.49835163e-10 ... -1.27248125e-09
-9.25932002e-10 -2.47061318e-10]
[-7.16707052e-09 -4.51270505e-09 -3.18572965e-09 ... -3.16794873e-09
-1.84340204e-09 8.09758940e-10]
[ nan nan nan ... nan
nan nan]]
Activity
jthielen commentedon Jun 6, 2020
Thank you for digging into this (and for prompting some flashbacks to numerical analysis)!
While unfortunate, I am not surprised by this at all from a numerical analysis point-of-view, and I'm not sure I'd even label it a bug. Of course, on a mathematical level, the Laplacian (in Cartesian coordinates) is equivalently the sum of second derivatives and the divergence of the gradient. But, they naturally differ if you are using an explicit n-point higher-order derivative formulation like MetPy does for its second derivative. I'd argue that MetPy's formulation is more correct (or at least more numerically well-behaved/rigorous), and that iteratively applying first derivatives is composing in a kind of extra smoothing operation (since you're incorporating more non-local information at any given point). Also, unless I've missed something, I don't think GEMPAK even has true higher-order derivatives?
Though, all this being said, if MetPy relying on explicit second-order derivatives does not meet most users' expectations, I would also understand something being changed here.
Also, @kgoebber, your toy example data unfortunately do not have enough non-linearity to expose the differences. Here's an example that does, however, showing both
nanpropagation and value differences: https://nbviewer.jupyter.org/urls/dl.dropbox.com/s/5kt5gfkbm0kmsz7/metpy_issue1389_1.ipynbkgoebber commentedon Jun 8, 2020
Thanks for the great example @jthielen, I had a sneaking suspicion that it had something to do with the numerical differencing of non-uniform grid spacing. Yes, I would agree not a bug, I should have chosen a different initial issue tag. I don't think we have to match GEMPAK, but we'll need to make it explicit that we are doing the calculation differently and if there appears to be demand, we can always add a function to do that calculation in that way, but I would prefer to stay with the more numerically robust solution as we currently have it.
You are correct that GEMPAK does not have a true higher-order derivatives, just first derivative/gradient.
I've used the MetPy Laplacian calculation in some of my courses and I've never thought it was off or wrong, and there have only been the handful of notable exceptions between MetPy and GEMPAK calculations, so I wanted to make sure we rolled down this path to have a clear record of it. The other exceptions are wind speed (with GEMPAKs [outdated?] conversion factor) and potential temperature using the ratio of 2/7 instead of 287/1004.
kgoebber commentedon Jun 12, 2020
Okay, so here are some global maps of the differences between a native
second_derivativecalculation of the Laplacian or iterative first derivatives through use ofgradientanddivergence.This image is taking the Laplacian of the 500 hPa height field and they give largely the same answer, except that the
second_derivativecalculation is considerably noisier. This aligns well with @jthielen indicating that in a numerical analysis sense, there would be a bit of smoothing that happens when using an iterative first derivative method.Honestly, I'm a bit torn as to which way we should go. The noise in the second derivative form is not insubstantial and overwhelms the "signal" over northeast Russia (as an example). A simple test of using the 5 point smoother with one pass greatly reduces the noise and provides a similar smoothing to the iterative first derivative calculation. However, there are some prime meridian artifacts that appear (see below).
Others thoughts?
sgdecker commentedon Jul 24, 2020
See also #893
miniufo commentedon Aug 22, 2020
Not so sure if this is related to the type of the grid. Does the difference scheme take into account different type of model grids? Like five Arakawa stagged grids. For Arakawa C grid, 1st-order finite difference may have half-grid offset relative to its original location, and 2nd-order finite difference will shift back to the original location.
dopplershift commentedon Aug 22, 2020
These functions assume that the coordinate values (e.g. x) that you provide correspond directly to the scalar values, no shift. The calculated derivative estimates are valid at these locations as well--no half-grid shift.
I'm open to adding such capabilities, though I have no idea how readily that math can fit within our existing implementation.
miniufo commentedon Aug 23, 2020
I think you are right. It is not easy to generalize such calculation to various kinds of grids. Although most reanalysis data are on Arakawa A grid (no shift), many model outputs are defined on stagged Arakawa grids. Since difference schemes may depend on the underlying grid, using 1st-order derivative as a building block for higher-order derivatives may not be so good as calculating higher-order derivatives directly.
If I understand it right, 2nd-order derivative using twice 1st-order derivative involves ((q[i+2] - q[i]) - (q[i] - q[i-2])) /dx/dx/4, while it may also be calculated as ((q[i+1] - q[i]) - (q[i] - q[i-1])) /dx/dx (assuming a constant dx). The magnitude of the differences between the two ways depends on the smoothness of the q field.
So I guess the best way is to use higher-order function directly, rather than using 1st-order one recursively.
8 remaining items
jthielen commentedon Feb 9, 2022
I've been browsing through issues on the 1.3.0 milestone, and wanted to check in about how this could be resolved as a part of #893. Should we replace the implementation of
laplacianwithdivergenceof the scalar gradient,(Also, I put a reference to the side discussion about staggered grids in #1360, as I'd think any staggered grid handling in MetPy would come about through xgcm compatibility.)
dopplershift commentedon Feb 11, 2022
I think going down the "vector derivative-based
divergenceof the scalar gradient" is the right choice here since it's the truest to the physics and involves the least amount of re-building the corrections.second_derivativestill exists for the other case.sgdecker commentedon Feb 11, 2022
I concur with @dopplershift. The first option + resolving #893 should fix this issue for free.
dopplershift commentedon Dec 23, 2022
Closed by #2743, which added
geospatial_laplacian, which calculates the Laplacian as the divergence of the scalar gradient, including all of the projection/spherical corrections (if the coordinate/projection information is present).