Skip to content

Instantly share code, notes, and snippets.

@jswhit
Last active August 28, 2020 22:51
Show Gist options
  • Save jswhit/7a22af0659f8a46487d642ac86d8c328 to your computer and use it in GitHub Desktop.
Save jswhit/7a22af0659f8a46487d642ac86d8c328 to your computer and use it in GitHub Desktop.
interface to shtns spectral transform lib
import shtns
import numpy as np
class Spharmt(object):
"""
wrapper class for commonly used spectral transform operations in
atmospheric models.
"""
def __init__(self,nlons,nlats,ntrunc,rsphere=6.3712e6,gridtype='gaussian'):
self._shtns = shtns.sht(ntrunc, ntrunc, 1, \
norm=shtns.sht_orthonormal+shtns.SHT_NO_CS_PHASE)
if gridtype == 'gaussian':
#self._shtns.set_grid(nlats,nlons,shtns.sht_gauss_fly|shtns.SHT_PHI_CONTIGUOUS,1.e-10)
self._shtns.set_grid(nlats,nlons,shtns.sht_quick_init|shtns.SHT_PHI_CONTIGUOUS,1.e-10)
elif gridtype == 'regular':
self._shtns.set_grid(nlats,nlons,shtns.sht_reg_dct|shtns.SHT_PHI_CONTIGUOUS,1.e-10)
self.lats = np.arcsin(self._shtns.cos_theta)
self.lons = (2.*np.pi/nlons)*np.arange(nlons)
self.nlons = nlons
self.nlats = nlats
self.ntrunc = ntrunc
self.nlm = self._shtns.nlm
self.degree = self._shtns.l
self.lap = -self.degree*(self.degree+1.0).astype(np.complex)
self.invlap = np.zeros(self.lap.shape, self.lap.dtype)
self.invlap[1:] = 1./self.lap[1:]
self.rsphere = rsphere
self.lap = self.lap/rsphere**2
self.invlap = self.invlap*rsphere**2
def grdtospec(self,data):
"""compute spectral coefficients from gridded data"""
return self._shtns.analys(data)
def spectogrd(self,dataspec):
"""compute gridded data from spectral coefficients"""
return self._shtns.synth(dataspec)
def getuv(self,vrtspec,divspec):
"""compute wind vector from spectral coeffs of vorticity and divergence"""
return self._shtns.synth((self.invlap/self.rsphere)*vrtspec, (self.invlap/self.rsphere)*divspec)
def getvrtdivspec(self,u,v):
"""compute spectral coeffs of vorticity and divergence from wind vector"""
vrtspec, divspec = self._shtns.analys(u, v)
return self.lap*self.rsphere*vrtspec, self.lap*self.rsphere*divspec
def getpsichi(self,u,v):
psispec, chispec = self._shtns.analys(u, v)
return self.rsphere*self.spectogrd(psispec), self.rsphere*self.spectogrd(chispec)
def getgrad(self,divspec):
"""compute gradient vector from spectral coeffs"""
vrtspec = np.zeros(divspec.shape, dtype=np.complex)
u,v = self._shtns.synth(vrtspec,divspec)
return u/rsphere, v/rsphere
"Spharmt.py" 50L, 2429C
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment