pythongradientpython-xarraymetpylaplacian

The laplacian of metpy can't convert the units of temperature advection

i received recently a work from university to calculate all terms of Sutcliffe equation, and now i trying to do this in python. But, for the Temperature Advection term, i need to do the laplacian of temperature advection, and this error is blocking me

raise DimensionalityError(
pint.errors.DimensionalityError: Cannot convert from 'kelvin / second ** 3' ([temperature] / [time] ** 3) to 'kelvin / meter ** 2 / second' ([temperature] / [length] ** 2 / [time])

My code:

import matplotlib.pyplot as plt
import metpy as mp
import metpy.calc as mpcalc
import metpy.constants as mpconstants
from metpy.units import units
import numpy as np 
import xarray as xr 


#Importando os dados
data = xr.open_dataset('Datain\data1_sutcliff.nc').metpy.parse_cf()

#extraindo os eixos
lat = data['latitude'][:]
lon = data['longitude'][:]
level = data['level'][:]

#calculando pontos de grades
dx,dy = mpcalc.lat_lon_grid_deltas(lon,lat)

#Colocando constantes
g = 9.8 * units('m/s^2')
f = mpcalc.coriolis_parameter(np.deg2rad(lat))*units('1/s')
R = mpconstants.Rd

#Extraindo as variaveis
#Calculando a advecção no dataset
advData = mpcalc.advection(data['t'],data['u'],data['v'])

#Extraindo dos níveis de interresse  e considerando em 500 a espessura da camada 

adv500 = advData[:,8,:,:]

lap500 = mpcalc.laplacian(adv500)

For this job, i use the NetCDF ERA5 dataset pressure level

Now, i tried to do two gradients to get laplacian, but the new issue appers for the ndims

grad  = mpcalc.gradient(adv500)
lapla = mpcalc.gradient(grad)

AttributeError: 'tuple' object has no attribute 'ndim'

Solution

  • The problem here is that laplacian will try to calculate the sum of second derivatives along ALL axes of the passed in array. In your case this would include time, which is incorrect. Thankfully the unit framework disallowed this, but the error is hard to understand. I've opened an issue to evaluate having MetPy exclude the time dimension by default (since it could determine this from the dimensionality on the xarray coordinates).

    In the meanwhile, you can get the behavior you want by manually specifying the horizontal spatial dims:

    lap500 = mpcalc.laplacian(adv500, axes=('longitude', 'latitude'))

    or just

    lap500 = mpcalc.laplacian(adv500, axes=('x', 'y'))

    The reason the two gradient calculations is failing is that gradient returns a tuple of first derivatives, one array for each dimension. If you wanted to approach it this way, you'd need to do:

    ddt, ddy, ddx = mpcalc.gradient(adv500)
lapla = mpcalc.gradient(ddy)[1] + mpcalc.gradient(ddx)[2]