incorrect sign in metpy.calc.vorticity when using lambert_azimuthal_equal_area projection

[gianlisi]

Hello,

first of all, thank you for the great work. Metpy, and its integration with xarray, has been incredibly useful for my research.

I am writing about using metpy.calc.vorticity on a Lambert Azimuthal Equal Area grid around the north pole. As can be seen in the left panel of the attached picture, the sign of the computed vorticity is inverted in the upper half of the domain. The right panel shows what the vorticity should look like, as computed with xarray's differentiate.

The issue seems to originate in the azimuth dependent change of sign within the metpy.calc.lat_lon_grid_deltas function. Indeed, if I comment out line 848 and line 854 of metpy/calc/tools.py, the vorticity is returned with the correct sign over the entire domain. At the same time, I am afraid that changing the code will break the computation on other grids.

The definition of the grid I am using is

grid = { 'grid_mapping_name':'lambert_azimuthal_equal_area', 'longitude_of_projection_origin':0, 'latitude_of_projection_origin':90, 'false_easting':0, 'false_northing':0 }

Any help in understanding how to compute the correct vorticity would be greatly appreciated.

Thank you again for the great work!

Gianluca

[dopplershift]

Yikes! Can you post some sample code, especially showing your call to vorticity?

[gianlisi]

Hi dopplershift,

of course, here it is. The test dataset can be downloaded from http://svante.mit.edu/~mgl/metpy/ds.nc

import xarray as xr
import metpy.calc as mpcalc
import cartopy.crs as ccrs
import matplotlib.pyplot as plt

 
# define dataset grid projection
dsProj = {
      'grid_mapping_name':'lambert_azimuthal_equal_area',
      'longitude_of_projection_origin':0,
      'latitude_of_projection_origin':90,
      'false_easting':0,
      'false_northing':0
      }

# open dataset and add units to dataarrays
ds = xr.open_dataset('ds.nc').metpy.assign_crs(dsProj)
ds['u'] /= 100 # the original data is in cm/s rather than m/s
ds['v'] /= 100
ds.u.attrs = {'units':' m/s'}
ds.v.attrs = {'units':' m/s'}
ds.longitude.attrs['units'] = 'degrees_east'
ds.latitude.attrs['units'] = 'degrees_north'

# compute vorticity using mpcalc and xarray differentiate
vortmpcalc = mpcalc.vorticity(ds.u, ds.v)
vortxarray = ( ds.v.differentiate('x') - ds.u.differentiate('y') )

# plot both vorticities for comparison 
plt.close('all')
fig,axs = plt.subplots(1, 2, subplot_kw={'projection': ccrs.NorthPolarStereo()}, figsize=(15, 5))

vortmpcalc.plot(ax=axs[0], vmin=-1e-7, vmax=1e-7, cmap='RdBu_r')
vortxarray.plot(ax=axs[1], vmin=-1e-7, vmax=1e-7, cmap='RdBu_r')

axs[0].set_title('metpy.calc.vorticity(ds.u,ds.v)')
axs[1].set_title("ds.v.differentiate('x') - ds.u.differentiate('y')")

for ax in axs.flat:
  ax.set_extent([-180,180,70,90], crs=ccrs.PlateCarree())
  ax.coastlines()

plt.show()

[gianlisi]

Hi Ryan,

I was wondering if there was any updated regarding the sign change in the vorticity computations on the EASE grid.

Thanks!

Gianluca

[dopplershift]

We're at SciPy this week, but we have some work improving the handling of these spatial derivatives on projected coordinate systems planned for the next release (nominally end of the month).

Regardless, I've reproduced the problem so we should probably convert this to an issue.

[gianlisi]

Hi Ryan,

thanks, great to hear. Please let me know if I can be of any help.

In the meantime, I will convert this to an issue.

Thank you again for all your work, and enjoy scipy!

Gianluca