SHTOOLS
SHExpandLSQ for complex values - combine separate coefficient arrays
Hi,
I want to apply SHExpandLSQ to a complex valued array. I have understood that SHExpandLSQ only takes real input, so I'm doing separate fits on the real and imaginary parts.
However I assumed that I could just combine the returned coefficient arrays and get the same results when for expand, but this doesn't work:
hilm_real, chi2_real = pyshtools.expand.SHExpandLSQ(np.real(reim_data), inlat, inlon, sphorder)
hilm_imag, chi2_imag = pyshtools.expand.SHExpandLSQ(np.imag(reim_data), inlat, inlon, sphorder)
hlm_real = pyshtools.SHCoeffs.from_array(hilm_real)
reconstructedReal = hlm_real.expand(lat=inlat, lon=inlon)
deviation_real = np.real(reim_data) - reconstructedReal
print('deviation_real for real: min = ' + str(np.min(np.abs(deviation_real))) + ', max = ' + str(np.max(np.abs(deviation_real))) + ', mean = ' + str(np.mean(np.abs(deviation_real))) + ', chi = ' + str(np.sum(np.square(deviation_real))))
hlm_imag = pyshtools.SHCoeffs.from_array(hilm_imag)
reconstructedimag = hlm_imag.expand(lat=inlat, lon=inlon)
deviation_imag = np.imag(reim_data) - reconstructedimag
print('deviation_imag for imag: min = ' + str(np.min(np.abs(deviation_imag))) + ', max = ' + str(np.max(np.abs(deviation_imag))) + ', mean = ' + str(np.mean(np.abs(deviation_imag))) + ', chi = ' + str(np.sum(np.square(deviation_imag))))
hilm_complex = hilm_real + 1j*hilm_imag
hlm_complex = pyshtools.SHCoeffs.from_array(hilm_complex)
reconstructedComplex = hlm_complex.expand(lat=inlat, lon=inlon)
deviation_to_separate = reconstructedComplex - (reconstructedReal + 1j*reconstructedimag)
print('deviation_to_separate: min = ' + str(np.min(np.abs(deviation_to_separate))) + ', max = ' + str(np.max(np.abs(deviation_to_separate))) + ', mean = ' + str(np.mean(np.abs(deviation_to_separate))) + ', chi = ' + str(np.sum(np.square(deviation_to_separate))))
which gives
deviation_real for real: min = 1.8123470792164031e-07, max = 0.02399276818309759, mean = 0.0015879964062471847, chi = 0.13405475009741488
deviation_imag for imag: min = 1.4086760584397506e-08, max = 0.025954209358143007, mean = 0.001817595279754483,
deviation_to_separate: min = 0.0, max = 0.2844782293230475, mean = 0.07088919414456854, chi = (2.6888420851680928+32.486839201184736j)
I expected that deviation_to_separate is zero (or at least very, very close to).
What is the right way to combine these two sets of coefficients into a single set of coefficients for the complex function?
...nevermind, figured it out myself!
For future reference and everyone else wondering:
Direct addition of the raw coefficient arrays does not work, because converting the coefficient arrays from real to complex is a nontrivial step.
You need to make use of the builtin capabilities of the SHCoeffs class:
hlm_complex = hlm_real.convert(kind='complex') + 1j*hlm_imag.convert(kind='complex')
Perhaps it would be useful to add a section about this in the documentation?
I'm not quite sure where to document this. Do you have any ideas?
- The best solution to this problem would be add a new fortran function that deals with complex data. I am all for someone contributing this! I rarely work with complex data, so this is not high on priority list.
- A hacky solution would be to add a new python function that wraps the intermediate steps you used. It's not ideal, but at least it would work.
For info, I've made some improvements to the least squares inversion routines: #463.
I didn't yet implement complex coefficients, but this should be easier to do now than before (the changes would go in the new function shlsq). I would welcome a PR on this as I don't have time to look into this now.