FFT Filter

[LL-Geo]

Draft of FFT Filter

Reminders:

• [ ] Run make format and make check to make sure the code follows the style guide.
• [ ] Add tests for new features or tests that would have caught the bug that you're fixing.
• [ ] Add new public functions/methods/classes to doc/api/index.rst and the base __init__.py file for the package.
• [ ] Write detailed docstrings for all functions/classes/methods. It often helps to design better code if you write the docstrings first.
• [ ] If adding new functionality, add an example to the docstring, gallery, and/or tutorials.
• [ ] Add your full name, affiliation, and ORCID (optional) to the AUTHORS.md file (if you haven't already) in case you'd like to be listed as an author on the Zenodo archive of the next release.

[welcome[bot]]

💖 Thank you for opening your first pull request in this repository! 💖

A few things to keep in mind:

• Authorship Guidelines: Our rules for giving you credit for your contributions, including authorship on publications and Zenodo archives.
• Code of Conduct: How we expect people to interact in our projects.

⭐ No matter what, we are really grateful that you put in the effort to do this! ⭐

[LL-Geo]

There is an example of how this code works. https://nbviewer.org/github/LL-Geo/PFToolbox/blob/master/FFT_Filter.ipynb

[LL-Geo]

I'll work on the documentation (formula) & formate later... Any suggestions are welcome 😄

[RichardScottOZ]

Thanks, will see what I can do later today.

[RichardScottOZ]

Sent you a pull request with a few small updates.

[RichardScottOZ]

@santisoler work here https://github.com/fatiando/harmonica/blob/fft-derivatives/harmonica/filters/fft.py for reference

[LL-Geo]

Opps.. commit to wrong branch 🙈

[santisoler]

Opps.. commit to wrong branch see_no_evil

Ha! No big deal, this can happen to everyone! But since git is so good, you can easily solve this with cherry-pick:

1. Switch to your isostasy branch.
2. Run git cherry-pick 3e4d14104c7db3c32fb4a3a95a39e935bf1f9d3a: git will find that commit across your branches and try to include it to the current branch.
3. Check if the commit has been correctly added to isostasy with git log.
4. Push isostasy

And for this branch you can use revert:

1. Switch to your main branch.
2. Run git revert 3e4d14104c7db3c32fb4a3a95a39e935bf1f9d3a
3. And then push

[LL-Geo]

Opps.. commit to wrong branch see_no_evil

Ha! No big deal, this can happen to everyone! But since git is so good, you can easily solve this with cherry-pick:

1. Switch to your isostasy branch.
2. Run git cherry-pick 3e4d14104c7db3c32fb4a3a95a39e935bf1f9d3a: git will find that commit across your branches and try to include it to the current branch.
3. Check if the commit has been correctly added to isostasy with git log.
4. Push isostasy

And for this branch you can use revert:

1. Switch to your main branch.
2. Run git revert 3e4d14104c7db3c32fb4a3a95a39e935bf1f9d3a
3. And then push

Thanks @santisoler ! I finally remember to do this.. 😅

[LL-Geo]

I have added all filters in the new format. @santisoler Just curious how do we do the test? Do we want to do it by comparing the FFT result with the forwarding result, or do we just test the filter itself? 😆 One silly question about the code style. Do we do from_ .filters._filters import * orfrom .filters._filters import derivative_upward, derivative_easting? The list will be very long...

[LL-Geo]

We might need to make the magnetic point source work first then we could test RTP and pseudo_gravity... And how could we test the Gaussian filter?..

[RichardScottOZ]

I suppose a practical usage consideration, too - memory use compare to size of dataset?

e.g. if say pseudogravity is doing more than one thing, what will it use?

[RichardScottOZ]

@LL-Geo I noticed on the Pseudogravity

[RichardScottOZ]

transform is wrong in Y?

[santisoler]

transform is wrong in Y?

@RichardScottOZ I think there's a missing return da_filter in the pseudo_gravity_kernel function: https://github.com/LL-Geo/harmonica/blob/5416778c88a5591030c765c203aa5adc54460c7b/harmonica/filters/_filters.py#L437-L528

Let me add it and give it another try.

BTW, @RichardScottOZ, are you using this branch on production already? If so, I'll try to prioritize this PR.

[santisoler]

We might need to make the magnetic point source work first then we could test RTP and pseudo_gravity...

Got it. I'll be working on the kernels for the magnetic dipoles in Choclo. When those are up, we can used them to test these filters right away.

And how could we test the Gaussian filter?

That's trickier. I think if we test the kernel function well enough, we could be confident that the public function will apply it in the right way, since it will be using the apply_filter function we are already extensively testing. I wouldn't put a lot of thought on that, at least for now.

[RichardScottOZ]

transform is wrong in Y?

@RichardScottOZ I think there's a missing return da_filter in the pseudo_gravity_kernel function: https://github.com/LL-Geo/harmonica/blob/5416778c88a5591030c765c203aa5adc54460c7b/harmonica/filters/_filters.py#L437-L528

Let me add it and give it another try.

BTW, @RichardScottOZ, are you using this branch on production already? If so, I'll try to prioritize this PR.

Yes, been using it quite a lot - but just noticed the wrong way around a week or so ago when I fed a layer into something else - which also needs I have a 50GB pseudograv map of Australia to fix too. :)

[RichardScottOZ]

Ever since LL did the sample notebook I have been tinkering.

Let me know when you want me to test something.

[LL-Geo]

transform is wrong in Y?

@RichardScottOZ I think there's a missing return da_filter in the pseudo_gravity_kernel function: https://github.com/LL-Geo/harmonica/blob/5416778c88a5591030c765c203aa5adc54460c7b/harmonica/filters/_filters.py#L437-L528 Let me add it and give it another try. BTW, @RichardScottOZ, are you using this branch on production already? If so, I'll try to prioritize this PR.

Yes, been using it quite a lot - but just noticed the wrong way around a week or so ago when I fed a layer into something else - which also needs I have a 50GB pseudograv map of Australia to fix too. :)

Oops... Sorry about that...

Have a short break this month. I'll make some updates when I back to office tomorrow. 😄

[RichardScottOZ]

Thanks!

[RichardScottOZ]

Thanks!

so what is the difference between https://github.com/LL-Geo/PFToolbox/blob/master/FFT_Filter.ipynb and the changes here?

Don't think I see it at first or second look.

[LL-Geo]

Thanks!

so what is the difference between https://github.com/LL-Geo/PFToolbox/blob/master/FFT_Filter.ipynb and the changes here?

Don't think I see it at first or second look.

Yup. These are basically the same.

My test case didn't see the issue for the Y. Could you send me some data for this?

[RichardScottOZ]

Ok, will see what I can find.

[RichardScottOZ]

Assuming of course I did it correctly, too.

[RichardScottOZ]

gAerogmagnetico_Reducido_al_polo.zip

[RichardScottOZ]

This is a decent sized chunk of Peru - clearly flipped y.

[RichardScottOZ]

Admittedly just a default test as don't have access to the raw data details

gridpseudo = FFT_Filter(mag.squeeze()).pseudo_gravity(90,0,Im=None,Dm=None,F=50000,savefilter=False)

[LL-Geo]

gAerogmagnetico_Reducido_al_polo.zip

Yup, I see the problem now.

I think the issue come from xrft, which will return an increased order of y coordinate when doing ifft.

Find similar flip latitude in xrft example https://xrft.readthedocs.io/en/latest/DFT-iDFT_example.html

[LL-Geo]

Maybe add a warning use ascending order of y?

The result is actually not wrong (The Y and matrix are all flipped). But need to be careful that sometimes the y coordinate is changed slightly (tiny change).

[RichardScottOZ]

Yeah - not wrong in that look at the output sense - but wrong if going to do downstream modelling that has affine code type work in it. That is how I noticed it originally when something that generally worked after going through a bunch of checks didn't, so finally checked that last. Do you think would could add a kwarg that flipped it back if changed? e.g. keep original y sign as such?