Other grids than equi-distant cylindrical?

[milankl]

Maybe a long shot, have you ever thought about using different grids than equi-distant cylindrical (i.e. lon-lat) grids? In weather forecasting, some models also use

• reduced Gaussian grids: Fewer zonal points towards the poles to correct for the convergence of the longitudes which cause grid points to cover smaller areas towards the poles. However, I'm not sure how "fewer" is determined actually.
• octahedral grids. They are actually mathematically quite beautiful, see picture (from ECMWF)

They are constructed as follows. Get the Gaussian latitudes

julia> using FastGaussQuadrature
julia> nlat = 32
julia> latd = -asind.(FastGaussQuadrature.gausslegendre(nlat)[1])[1:nlat÷2]
16-element Vector{Float64}:
 85.7605871204438
 80.26877907225
 74.7445403686358
 69.21297616937085
 63.67863556109686
 58.14295404920328
 52.606526034345265
 47.06964205968768
 41.53246124665608
 35.99507841127159
 30.457553961152094
 24.91992862994861
 19.38223134643438
 13.84448373438486
  8.306702856518806
  2.76890300773601

Then start with 20 points for the latitude closest to the pole and increase by 4 points for every latitude towards the Equator

julia> lons = [0:360/(20+4i):360-180/(20+4i) for i in 0:length(latd)-1]

Some of the Fourier transforms get an ugly length

julia> using Primes
julia> nlons = [length(lon) for lon in lons]
julia> [factor(nlon) for nlon in nlons]
16-element Vector{Primes.Factorization{Int64}}:
 2^2 * 5
 2^3 * 3
 2^2 * 7
 2^5
 2^2 * 3^2
 2^3 * 5
 2^2 * 11
 2^4 * 3
 2^2 * 13
 2^3 * 7
 2^2 * 3 * 5
 2^6
 2^2 * 17
 2^3 * 3^2
 2^2 * 19
 2^4 * 5

but ECMWF uses FFTW successfully for that.

[jmert]

Maybe a long shot, have you ever thought about using different grids than equi-distant cylindrical (i.e. lon-lat) grids?

Yes, #68 is building up to supporting any arbitrary iso-latitude grid. I've demonstrated the ability to do ECP and HEALPix with a common transform implementation given sufficient descriptions of each grid.

[milankl]

Thanks for these efforts already, that's exciting. Would love to test the HEALPix grid in SpeedyWeather.jl (I don't think it's been used so far in weather or climate models!) Will see how much of that PR I can reuse to make this work in our application. How do you store a HEALPix map? As a 12-element vector of matrices, or differently?

[jmert]

How do you store a HEALPix map? As a 12-element vector of matrices, or differently?

HEALPix maps are just vectors with length $12N_\mathrm{side}^2$, where $N_\mathrm{side}$ is the "resolution factor" (a power of two) of a HEALPix map that basically corresponds to the number of halvings of the 12 base pixels. (There are actually two indexing schemes in the HEALPix format, so-called RING and NESTED, but for spherical harmonic transforms, the only relevant one is RING.) It's the most common full-sky mapping format used in CMB analysis, hence why it's one of the pixelization formats I've played with.

The canonical paper that describes HEALPix is Gorski et al (2004): https://arxiv.org/abs/astro-ph/0409513

[jmert]

@milankl See #74, current src/sphericalharmonics.jl, and dev docs for my new SphericalHarmonicTransforms.jl for a demonstration of implementing alternate pixelization schemes with a single version of the synthesize/analyze functions.

[milankl]

Fantastic. Many thanks, as in #73 I'll start playing around with your code, but happy to create a dependency if you happen to register SphericalHarmonicTransforms.jl

[jmert]

In my opinion, it's too early to register SphericalHarmonicTransforms.jl since there are several things I know I want to do and haven't done yet — I kind of hurriedly extracted the code I'd been working on to extract to an independent package to allow you to start exploring its features, but I think I'd rather give myself some more time to reconsider interfaces/api/etc before possibly registering generally.