Extract normalized 2D spectrum values from SHCoeffs about shtools HOT 6 OPEN

I'll put this on my to-do list. It shouldn't be too hard, but I have a few other things to do first.

from shtools.

We do have a method that does this for the 1D case

clm.spectrum()

so I think what you are looking for is to create something like

clm.spectrum2d()

The only problem is how to format that output values. The raw coeffs are stored as two arrays (1 for negative orders, the other for positive orders) and the plot function in essence just spices these two together to make a single plot.

from shtools.

Yes i would expect the clm.spectrum2d() to give the values that are plotted in plot_spectrum2d() so:

clm = coeffs.to_array()

spectrum = np.empty((lmax + 1, 2 * lmax + 1))
mpositive = clm[0, :lmax + 1, :lmax + 1] * \
    clm[0, :lmax + 1, :lmax + 1].conj()
mpositive[~coeffs.mask[0, :lmax + 1, :lmax + 1]] = np.nan
mnegative = clm[1, :lmax + 1, :lmax + 1] * \
    clm[1, :lmax + 1, :lmax + 1].conj()
mnegative[~coeffs.mask[1, :lmax + 1, :lmax + 1]] = np.nan
spectrum[:, :lmax] = np.fliplr(mnegative)[:, :lmax]
spectrum[:, lmax:] = mpositive

and then normalized by user-settable normalization standard

from shtools.

Yeah it's not pressing, i've already mostly reimplemented it for my own use-case, but would make the codebase more consistent

from shtools.

I have a quick question on how you would like to use this feature. I don't think that it would be useful to output the array as in the above example because it would be hard to map the array indices (i, j) to (degree, order) indices. In particular, an order of m=0 would be in the middle of the array, and not the first value, and you would need to manually map the order to index j.

I see two possible ways of dealing with this.

First, we could output the array as an xarray, and then manually set the coordinate values to correspond to the actual degrees and orders. You could then access a given value by using a syntax like data.loc[3, -2] to access degree 3 and order -2.

Second, we could output the array in the same format at clm.coeffs, which is a numpy array of dimension (2, lmax+1, lmax+1). In essence we would just compute clm.coeffs**2, and then multiply this by the normalization factor.

I'm leaning towards the second option as this is the most simple and most easily related to the internal data. The only negative aspect is that it it wouldn't be easy to "plot" this (but that is what the dedicated plotting routine is for). I don't have a problem with outputting an xarray, but almost all of the shtools functions output numpy arrays, so this wouldn't be entirely consistent and might cause confusion.

from shtools.

For my use: I personally use the spectrum() technique to transform the data into plot-ready formats that i can pipe into my own matplotlib/seaborn scripts, and would expect a spectrum2d() function to give a similarly ready dataset for plotting; if i want to do analysis on the coeffs by order/degree, i would use the clm.coeffs.

I get where you're going with the consideration though, and it kind of depends on whether the expected use-case is indeed plotting or further analysis. In my eyes, a spectrum2d() function should return a 2D dataset/np array, as the dimensionality is in the function name, even though this necessitates an offset 0 order for a numpy array. However, this might be different for others.

I agree that for consistency the output should be an np array and not xarray.

Alternatively, maybe a parameter-setting that lets you choose between an np array of (2,lmax+1,lmax+1) or (lmax+1, lmax*2+1) could be a solution?

from shtools.