I want to invert the Laplacian on a 2d fixed grid in python, which was possible to do in ncl
using the ilapsf
function. I don't really want to try and do this from scratch using a convolution and thought this would be the bread and butter of pyngl, metview or metpy, but it seems ilapsf
didn't make the transfer to python and metpy
only does the easier Laplacian calculation, rather than solving the inverse Laplacian. I also didn't find anything in scipy
.
The old code in ncl
was as simple as this:
chi = ilapsF_Wrap(nei,0) # inverse laplace
gradchi = grad_latlon_cfd(chi, chi&latitude, chi&longitude, True, False) # calculate the gradient
I managed using the dv2uv
function in cdo, but it is the fudge of the century, involved conversion to a spectral field and then tricking it by turning the flux convergence into a divergence field, setting up a dummy zero vorticity field and then turning u and v into energy fluxes... (plus the meridional flux does weird things at the poles, which the ncl solution can handle). You would think that there would be a more generic inverse Laplacian available in python (in fact I'm sure there is, I'm just not looking in the right place).
from cdo import Cdo
cdo=Cdo()
# convert flux divergence into GG and the spectral G:
cdo.gp2sp(input="-remapcon,T511grid -selvar,nei nei.nc",output="neisp.nc")
# fudge CDO into thinking this is pure divergence (code=155)
cdo.chname("nei,sd",input="neisp.nc",output="tmp_vd.nc")
# make up a zero field for vorticity
cdo.chname("nei,svo",input="-gec,1e36 neisp.nc",output="tmp_vo.nc")
# the result is energy flux, not winds, so we will need to change the output names (and eventually metadata):
cdo.dv2uv(input="-merge tmp_vo.nc tmp_vd.nc",output="tmp_fluxes.nc")
cdo.chname("v,vflux",input="-chname,u,uflux tmp_fluxes.nc",output="fluxes.nc")