Implement NCL's spectra functionality: specx_anal, specx_ci and specxy_anal

[phillips-ad]

NCL's spectra functions are some of the most highly used by my research group (CAS/NCAR). I am not familiar with other python packages, but I would guess that a number of them already have these spectra capabilities. It would be great, following what was done with the new eof routines, if the GeoCAT group could port this capability to geocat-comp.

NCL links specx_anal specxy_anal specx_ci

[andersy005]

I'm not familiar with NCL's spectra functionality, but I recommend to anyone looking into porting this in geocat-comp to look into the xrft package: https://xrft.readthedocs.io/en/latest/index.html

[adamrher]

related issue: https://github.com/NCAR/geocat-comp/issues/43

[erogluorhan]

Thanks all for your contribution!

We will soon start working on this.

[ccalvosa]

Hello @erogluorhan Has there been an update on the implementation in Python of specx_anal? Thank you!

[erogluorhan]

Hi @ccalvosa

Thanks for asking! We have assigned these functions to one of our teammates, @pilotchute . He has been investigating and could give further insights.

[pilotchute]

@ccalvosa This is on my schedule to start within the next ten workdays. I hope to have it included in the 2021.06.0 or 2021.07.0 release of geocat.-comp

[pilotchute]

Specx<?> notes:

input:

• x: uniformly spaced time series data. (example use is monthly data from climate simulations)
• iopt: series mean, trend, or nothing, removed from data prior to primary calculation
• jave: window width for deniell smoothing. (example user always used 0.1)
• pet: percent of series data to be tapered. (not used by example user)
• scl: scaling factor. (user experience not in notes)

output:

• frq: frequency axis values for result. (uniformly spaced at resolution of algorithm)
• specx: spectral intensities for each value in frq. (I think this is the same as an FFT but the source algorithm mentioned in the docs is blackmon-tukey, not cooley-tukey. This may be due to the implementor replacing cooley with their name, when implementing tukey's FFT alg. Or it could be different. At this point the algorithm is unknown, but does not appear to return complex numbers, so it's not a conventional FFT)
• nspe: len(x)/2+1 (the number of result points below the nyquist frequency)
• sinfo: vector of length 10, containing the following:
  • 1: 'equivalent degrees of freedom' (unknown meaning, unused by our example user)(might be the number of spectral peaks above some calculated value)
  • 2: lag1 autocorrelation of the detrended and tapered input data.(used by our example user to determine if the time series data was irregular enough to be experimental results)
  • 3: same as 2, reserved for use with second axis of 2d datasets
  • 4: unused
  • 5: frequency spacing between result values. (the spectral resolution of the result)
  • 6: 'band width', = sfinfo(1)sfinfo(5)0.5. (this seems to be dependent on the degrees of freedom, which is not yet understood)
  • 7: multiplication factor for 5% limit (not used by example user, meaning unknown)
  • 8: multiplication factor for 95% limit (not used by example user, meaning unknown)
  • 9: 'degrees of freedom for series via spectrum shape blackmon-tukey'. (not used by example user, meaning unknown, this may be what earlier documentation references to blakmon-tukey referred to, which could mean that the primary calculation may be cooley-tukey's fft)
  • 10: unused.

End of notes

[pilotchute]

Additional info from Adam Phillips:

Here is a direct link to the CVDP-LE output files I was speaking of this morning: http://webext.cgd.ucar.edu/Multi-Case/CVDP-LE_repository/CMIP6_Historical_1900-2014/CMIP6_Historical_1900-2014.cvdp_data.tar Each model run (and ensemble mean) have their own .nc file. Choose any single model simulation file, and then grab the following arrays: nino34 (=raw nino3.4 timeseries. Each individual model timeseries is plotted here.) nino34_spectra (= spectra from specx_anal (0,:), then from specx_ci Markov red noise curve (1,:), 95% confidence bounds for spectrum (2,:) and 99% confidence bounds for spectrum. Each individual spectra is plotted here.)

[dcherian]

blackmon-tukey

This is a typo: it's blackman-tukey https://ocw.mit.edu/courses/earth-atmospheric-and-planetary-sciences/12-864-inference-from-data-and-models-spring-2005/lecture-notes/tsamsfmt_1_12.pdf

I'm amazed anyone uses it in this millenium :P!

[pilotchute]

blackmon-tukey

This is a typo: it's blackman-tukey https://ocw.mit.edu/courses/earth-atmospheric-and-planetary-sciences/12-864-inference-from-data-and-models-spring-2005/lecture-notes/tsamsfmt_1_12.pdf

I'm amazed anyone uses it in this millenium :P!

Thanks, I had assumed that the fortran's spelling was correct since I had found this refence to a a 1959 paper: https://link.springer.com/chapter/10.1007%2F978-1-4757-2514-8_8

[pilotchute]

• [ ] I'm thinking we replace the blackman-tukey calculation for the specx output with amplitude result of a standard FFT.
  • this is in line with modernizing the functional components, while providing the same utility.
  • unlikely that blackman-tukey has any advantage over a calculation of intensity done via a vector length result of the complex values returned by numpy's FFT.
• [ ] Does anyone have input on the return values that we don't have an established understanding of?
  • sfinfo[1,7,8,9]
• [ ] Does anyone have suggestions about keeping of removing the detrending and tapering that are not required when doing an FFT, that were required for the BT alg?

[ccalvosa]

Hi @pilotchute, Is there any update of these functions? Thx!

[pilotchute]

Hi @pilotchute, Is there any update of these functions? Thx!

I've been out on medical leave, I'm resuming work part time. Though I think that our analysis of the specx_<> function showed that they were developed prior to current signal processing techniques, and could (should?) be replaced by standard sp library functions.

My personal investigations (and interrogatory discussions with one of our scientists who uses the function) suggested that numpy's FFT served all the qualitative needs, and was more precise for any quantitative needs.

If you have a specific case of the functions output you would like to see in geocat, let's talk about it, and I would be happy to implement something from first principles in python to match your needed functionality.

[pilotchute]

@ccalvosa could I get an update from you? Is there something specific that specx_anal does that isn't currently available from a public source?

[ccalvosa]

@ccalvosa could I get an update from you? Is there something specific that specx_anal does that isn't currently available from a public source?

Hi @pilotchute, Last weeks we tried to analyze the datasets with numpy and scipy functions and finally we achieved it. Thank you for everything!

[pilotchute]

reducing priority as functionality exists in numpy/scipy. Will not close as a single call function within geoat-comp to do the signal processing would be a good enhancement.

[pilotchute]

@ccalvosa would you be willing to show me what your solution was? It would help us add an enhancement to make your process of recreating specx_anal into a single function call.

[Happy-fish-yu]

Hi, @pilotchute Is there any update of these functions? I‘m looking forward to using these functions in the future :)

[pilotchute]

@Happy-fish-yu The specx_<> functions are old (from the 1970s) and similar (and inferior) to the DFT and FFT functions developed after them.

If you would like to tell me more about your use case, I may be able to either implement a function that simplifies the use the FFT (or DFT for irregularly spaced data), but for now we are not working on this family of functions.

[rejianup]

Hii,@pilotchute, As you mentioned, other components of specx_ can be calculated with scipy or NumPy. However, calculating the theoretical Markov spectrum and lower and upper confidence curves using the lag1 autocorrelation mentioned "specx_ci' is not possible with any function in python. confidence curves are important in explaining the significance of a scientific study.

[jukent]

Using Scipy.signal module I was able to perform spec_1 and 2 from https://www.ncl.ucar.edu/Applications/spec.shtml

There are some minor differences in the results that are affected by the tapering method, smoothing method, window used during the spectral analysis, etc. The NCL code is not always descriptive about what equation it uses for which part. I tried to demonstrate how to do it as close as I can tell to NCL, but will also show other scipy parameter options in my GeoCAT.applications notebook.

Ultimately the art and science of spectral analysis has changed enough since NCL was written that I don't think I need to achieve bit-for-bit precision. I demonstrate the options and the scientist needs to know what options they'll want for their specific application.

I'll clean up and push the code tomorrow.