Stream: python-questions

Topic: spherical harmonics from lat-lon grid

view this post on Zulip Falko Judt (Apr 26 2022 at 20:12):

Hi folks, I'm wondering if anyone has computed spherical harmonics in Python from model output. So far I have used NCL's shaeC. I see that there is the pyspharm package, and will probably try that if there isn't a GeoCAT function that mimics the old NCL shaeC.

view this post on Zulip Katie Dagon (Apr 26 2022 at 20:35):

@Alea Kootz is working on something related to this for the @geocat team after I had requested functionality to replicate NCL's gradsf function which I used to calculate lat/lon gradients on a fixed grid. Relevant Github issues are pasted below, and there may be something coming via geocat-comp:
https://github.com/NCAR/geocat-comp/issues/215
https://github.com/NCAR/geocat-comp/issues/147

view this post on Zulip Alea Kootz (Apr 26 2022 at 23:25):

@Falko Judt Yes there is a spherical harmonics decomposition and recomposition set of functions in geocat-comp. They were released last month. :)

view this post on Zulip Alea Kootz (Apr 26 2022 at 23:27):

https://geocat-comp.readthedocs.io/en/latest/user_api/generated/geocat.comp.spherical.decomposition.html#geocat.comp.spherical.decomposition

view this post on Zulip Falko Judt (Apr 27 2022 at 17:47):

Thanks @Alea Kootz and @Katie Dagon! I think this is exactly what I need.

view this post on Zulip Alea Kootz (Apr 27 2022 at 17:49):

Let me know if you have any questions.

view this post on Zulip Falko Judt (May 23 2022 at 21:34):

Hi @Alea Kootz & @geocat ! So I tried to use geocat.comp.spherical.decomposition in the same way I used sheaC in NCL. But there are some issues.

In NCL, I provided the data on a fixed lat-lon grid (this is the only way the NCL function works), geocat.comp.spherical.decomposition seems to be more general, and so one needs to create weights beforehand with geocat.comp.spherical.scale_voronoi. This is where I'm stuck...

Here's what I did:

import xarray as xr
import numpy as np

ds = xr.open_dataset('/glade/campaign/mmm/wmr/fjudt/projects/dyamond_1/3.75km/latlon/latlon.nc')
lat = np.deg2rad(ds.latitude+90.)  #---lat array needs to go from 0--pi, but in our lats go from pi/2 to -pi/2.
lon = np.deg2rad(ds.longitude)
data = ds.vorticity_200hPa

#---lat/lon need to be 2D, but in the fixed grid data file they're just 1D arrays. So make 2D arrays via broadcasting...
lats1, lats2 = xr.broadcast(lat,lon)
lons1, lons2 = xr.broadcast(lon,lat)

%%time
weights = gc.spherical.scale_voronoi(lons1, lats1)
---> get an error here (also it's slow)

%%time
gc.spherical.decomposition(data, weights, lons1, lats1)

Any idea why, and is there an example that tells us how to use geocat.comp.spherical.decomposition if you're coming from NCL's shaeC?

view this post on Zulip Alea Kootz (May 23 2022 at 21:36):

@Falko Judt Check what you are doing with the broadcast you should only need to call it once.

view this post on Zulip Alea Kootz (May 23 2022 at 21:38):

And I would suggest using numpy.meshgrid instead of xr.broadcast for your 2d lat/lon arrays.

view this post on Zulip Alea Kootz (May 23 2022 at 21:40):

see here to understand how that works from the spherical.py test setup.

Last updated: May 16 2025 at 17:14 UTC