wavenumber-frequency filtering of raw OLR data

[AndrewILWilliams]

Hi Brian!

Thanks for making this code available, it's been really helpful for me recently. I did try to code up something similar myself but couldn't understand how...

I was wondering though, do you know how to extend this code to do the 'wavenumber-frequency' filtering that they do in WK99? Specifically for their Figure 7? As I understand it, they define regions of the diagram which are representative of Kelvin wave activity (or ER, EIG etc) and then they apply this filtering to raw time-series' of twice-daily OLR data.

The idea of it makes sense, but I have no idea how one would actually implement that so just wondering if you had ideas/pointers!

Cheers, Andrew

[brianpm]

Hi @AndrewWilliams3142 -- Yes, the filtering is something I had considered adding, but didn't get around to it yet. I think the strategy is not too tricky. Basically, do the fourier decomposition, and then take the resulting spectrum and mask out all the wavenumbers and frequencies you don't want, and then do an inverse FFT from what is left. That transforms back to physical space, but with only the variance associated with the spectral region of interest. There might be some details to work through, but that's the idea. I doubt I will get to doing this soon, but maybe if I have some time I'll see what I can do.

[AndrewILWilliams]

Hi @brianpm. Thanks for the pointers! I think I've managed to get this working, I've tried to generate something akin to the NCL Kelvin-wave filtering example as a sanity check. Seems to work alright.

First I detrend the whole time-series and remove the mean, then I taper in time before doing the fourier decomposition. Then in this example I just set all frequencies outside of the [0.1, 0.5] cpd range to zero, and any wavenumbers outside of the [0,15] range to zero. Then I did np.fft.ifft2(...).real of that and presto!

I'll hopefully have some time in the new year to tidy this up. Perhaps we could put a proper package together for this?

[brianpm]

@AndrewWilliams3142 , Nice! That looks like it is working well.

If you get a chance to clean it up and want to do something with it, let me know. I'm happy to share my code in whatever way it is useful. I've heard that it might be getting used in some of the E3SM diagnostics being developed, too. (They might also be doing a similar filter as what you've done, but I haven't seen any code or plots from that.)

[AndrewILWilliams]

Cool stuff! Yes, I'll let you know once I've got a tidy version up and running. :)

I think there might be some subtleties I'm missing at the moment though. Do you have a simple explanation for what happens in your resolveWavesHayashi function? In the code it says that the frequencies/wavelengths from np.fft are not the same as those of the east/west waves, but I'm not sure I understand why that is?

[AndrewILWilliams]

For context, the reason I ask is because when I did the filtering (above), before selecting for the wavenumbers/frequencies I wanted, I first centered them using:

# `da_raw_fft` is the fft of the raw OLR anomalies
# da_raw_fft = np.fft.fft2(ffteqolr, axes=(0,2))

z_c = np.fft.fftshift(da_raw_fft, axes=(0,2)) 
centered_freq = np.fft.fftshift(da_raw_fft['frequency'])
centered_wavn = np.fft.fftshift(da_raw_fft['wavenumber'])

da_fft_centered = xr.DataArray(np.flip(z_c, axis=2), # Flip along 'wavenumber' or else things are the wrong way around
                                                      dims=['frequency', 'lat', 'wavenumber'],
                                                      coords={'frequency': centered_freq, 
                                                                     'lat': da_raw_fft['lat'], 
                                                                     'wavenumber': centered_wavn})

and then before performing the inverse fourier transform I had to 'invert' this to get the frequencies and wavenumbers back to what np.fft.ifft2 expects.

def prep_for_ifft2(da):
    """
    Undo the transformations I made to get the
    spectral power in correct wavenumber/frequency units
    for identifying eastward-/westward-moving waves
    
    The issue is that `np.fft.ifft2` expects data ordered 
    like the output of `np.fft.fft2` - however this
    ordering is not very physical...
    """
    data = np.flip(da.values, axis=2)
    
    data = np.fft.ifftshift(data)
    
    return data

This seems to work but I think that I should actually be doing something more like what is in your resolveWavesHayashi function? Sorry if this is unclear! I appreciate the help :)

[brianpm]

I believe that the complexity of resolveWavesHayashi is mainly switching things between quadrants to get the propagation convention used by Wheeler & Kiladis. I'll have to go and look more carefully to confirm, but I haven't had a chance lately. I'm hoping to work on some additional MJO diagnostics in the next few months, which should get me refreshed on this code.