harmonica
harmonica copied to clipboard
fatiando
LP / HP filtering of grav mag data
Hi, is there a way to filter data in terms of wavelength, so a " classical " low pass / high pass filter ? Or is this in another package (e.g. harmonica) ? So far, I didn't find any proper implementation or such a feature.
Are you willing to help implement and maintain this feature? Yes
Hi @nholz. Thanks for opening this Issue. As you might know Harmonica is still in current development stages, so we are still thinking which is the best design to implement features in the frequency domain. Meanwhile we are recommending xrft to perform the FFT and much more.
I don't know what kind of filter you want to apply, but you can find implementations of several bandpass filters using numpy and scipy. For example, the scipy-cookbook shows how to implement a Butterworth filter: https://scipy-cookbook.readthedocs.io/items/ButterworthBandpass.html
Hopes it could be useful for you.
We hope we can implement features like upward continuation in frequency domain (this was requested on the Fatiando Slack group) and other features in the near future. Nevertheless, maybe it would be better for this kind of generic filters to live inside xrft, while any specific geophysical frequency calculation could be included inside Harmonica.
Hi @nholz, as @santisoler said we've been trying to avoid including code in Harmonica that is generally applicable outside of the gravity and magnetics. So regular band-pass filters would be a great addition to xrft. It might be worth opening an issue there to ask for it and the ball rolling. Please feel free to link back to this issue as I'd love to be involved in the design and implementation of that in xrft.
@leouieda is there now a pythonic, don't need a PhD in geophysics way to approach this?
@santisoler scipy and GMT etc type explanations seem to me to be the esoteric, already need to know how to do it and what you are doing type of thing to start with.
xrft is coming along quite nicely and provides a much more intuitive way to do FFTs of xarray objects (it returns actual frequency/wavenumber and correctly ordered dimensions). So this should probably be contributed to that project instead since it's not specific to geophysics.
I would suggest starting with a Gaussian filter, which is much simpler than the Butterworth (not sure if there even is an n-dimensional version of it). This would be taking the FFT of the input array, making a Gaussian with a given center and standard deviation sampled on the same the frequency domain coordinates as the array, multiple the 2, and inverse FFT. Could be implemented in an arbitrary number of dimensions even, which will probably please the xrft developers.
Leaving this open so we can refer to it but will mark it as "wont fix" just to make it easier for us to triage later.
Sounds sensible.
Interesting idea - which did some filtering in a project earlier in the year, will check what we did.
xrft is coming along quite nicely and provides a much more intuitive way to do FFTs of xarray objects (it returns actual frequency/wavenumber and correctly ordered dimensions). So this should probably be contributed to that project instead since it's not specific to geophysics.
I would suggest starting with a Gaussian filter, which is much simpler than the Butterworth (not sure if there even is an n-dimensional version of it). This would be taking the FFT of the input array, making a Gaussian with a given center and standard deviation sampled on the same the frequency domain coordinates as the array, multiple the 2, and inverse FFT. Could be implemented in an arbitrary number of dimensions even, which will probably please the xrft developers.
-- Ok, centre based on what - not sure I get that part? Thanks!
It's been a long road since 2020.. haha. Sorry for the long waiting.. Following the release of v0.6.0 #386, we are able to run Gaussian High pass and low pass filters see example. I'll close this issue 😄
