-
Notifications
You must be signed in to change notification settings - Fork 44
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Implementation of FFT geophysical filters in Harmonica #154
Comments
Hi @santisoler, thanks for your interest in
Part of it was to be compatible with prior releases. We decided to use the naming of
The documentation page (https://xrft.readthedocs.io/en/latest/) has some examples with mathematical/theoretical explanation but feel free let us know if you have ideas on how to improve it.
The
This may be related also to issue #146 ? |
Hi @roxyboy ! Thanks for the quick response.
I understand. My question is: does it worth having two different functions that are intended to do the same thing with slight differences? IMHO, I find it a little bit confusing. Nevertheless, if you find necessary to maintain both functions, I would suggest to improve their docstrings to make it very clear about their differences.
Right! My concern was about the difference between the math description on the docstring and the default behavior of the function. By default
I like the documentation! I would just add mathematical descriptions of what the
Is there any reason why it wouldn't work with
Yes! I'll bring my ideas on #146 then. Thanks again for the response! We are really looking forward to use |
Maybe @lanougue can chime to remind me of the discussion we had but I think the idea was that in order to avoid confusion between
Noted. Will do this in a new PR.
Ok, I'll add this to the documentation.
|
Hi @santisoler, "lag flag" could also work with ifft but would not have any signification. fft works together with ifft and dft works together with idft. (However, only dft/idft correctly handle the phase). |
@santisoler, |
@lanougue just jumping in here 🙂 Would it be possible to store this information in the Dataset somehow so that users don't have to pass it separately? This could be through attributes or maybe just keeping the space/time-domain coordinates in the frequency domain Dataset so that they can be restored later. |
Hi @leouieda It could be possible to store the "lag" (evaluated during dft call) but I would prefer rather not for several reasons:
As I suggested, I think it could be useful to enable passing the initial dataarray coordinates to idft in order to make idft to automattically evaluate the lag. Maybe @rabernat and @roxyboy have other pros/cons arguments. |
Thanks so much @santisoler and @leouieda for this really useful feedback. We are 💯 committed to making xrft a useful, stable dependency for downstream projects such as yours. I'll chime in with a few of my thoughts.
IMO we should strongly consider deprecating one or the other.
Our documentation basically only has examples. What is missing is the standard narrative documentation that describes the operation and theory of xrft in detail.
IMO, we should not be concerned about memory footprint for coordinates. If you are using Xarray, you already have decided you want to accept the extra memory overhead in exchange for keeping careful track of coordinates.
@lanougue: How do you know this? Have we done a user survey? Not to my knowledge. Here we have some actual users coming and requesting this feature. We can see from the harmonica workaround that this is a feature that users need. We should not be so quick to say no. As I said above, I don't think that conserving memory is a valid reason. People use xrft (and python in general) for its convenience, not its memory efficiency. In general, I think that xrft should support full "round-tripping" of data, such that
Rather than placing this burden of keeping track of coordinate information on the user, I propose we use attributes and / or extra non-dimension coordinates to provide the necessary information about lag within the dataarray returned by fft. |
PR #151 adds some of the the operation and theory, which I think we can merge.
I'll work on this in another PR. |
Or perhaps the other way around is better: |
I did not say "no" ! I only started a discussion :) |
Sure, sorry if I mischaracterized your comments @lanougue. My apologies. I appreciate everyone's viewpoint on this discussion. |
I made a PR #155 to start the deprecation cycle of |
Thanks everyone for the amazing responses. I left a few thoughs after reading all the comments in this Issue.
I totally agree on this. I think the
IMHO, if any user would like to apply the FFT to a xr.DataArray regardless of the frequency arrays or its coordinates, they can always go back to
I agree with keeping
If we are all in the same page regarding the round-tripping property of @leouieda: I think you mention something like this at some point, could you bring some light on this? It think coordinates and not-index coordinates cannot be easily removed or modified. Nevertheless which option we choose, I think we should stick with the safer one. Finally but not least, I'm willing to help opening PRs or reviewing code. Feel free to assign me any task you want me to help with! |
Hi @santisoler |
Actually this question raised up some time ago. Should xrft.fft mimic the numpy.fft (because users are used to it) or should it be the closest as possible to the mathematical Fourier transform, as you suggest. We now need to choose how xrft.fft should behave. I also raised this question in PR #155 too because it is of central interest |
This is much more clever! I totally agree, so feel free to ignore my previous proposals regarding where to place the lags. |
Hi everyone!
I'm Santiago Soler, one of the core developers of the Fatiando a Terra project: open-source tools for geophysics. We have been quite interested in the development of
xrft
, specially for managing FFT-based processes and filters that we are planning to include in Harmonica, a library for processing, forward and inverse modelling gravity and magnetic data.We are very interested in using
xrft
for managing the whole FFT and inverse FFT processes that these functions need. Our potential fields data (gravity and magnetic) is usually given as a 2D regular grids whose coordinates are projected geographical coordinates: easting and northing (Cartesian coordinates in meters). See this gallery example to get a better picture of them. In the oldfatiando
package we used to have some implementations of these methods usingnumpy.fft
directly, although we are very interested in a FFT implementation that harnessxarray
objects, since we are usingDataArray
s andDataset
s as our default structure for our 2D grids.We therefore find
xrft
the perfect solution to our problems. We really appreciate you developers for taking the time to build this great project!I personally started drafting an implementation of a directional derivatives of a regular grid using FFTs (fatiando/harmonica#238) using
xrft
under the hood for computing the FFTs and their inverse.On this process, I started getting some user experience of
xrft
and wanted to share them with you.We also found that we might have some special needs and we would like to know if you would want to include solutions to them into
xrft
. We are absolutely willing to help with the implementation and maintenance of those.Differences between
dft
andfft
We found that
xrft
hosts two different functions for handling FFT computations:dft
andfft
.We found out that
fft
is a wrapper ofdft
that passes hard-coded values to certain parameters.IMHO, it sounds slightly confusing at the beginning, and in fact I had to surf through the code to finally understand the difference between them.
Is there any strong decision for this design?
Mathematical docstring of
dft
The docstring of
dft
states that it returns the Fourier transform computed aswhere I assume that the
\overline
states the mean value of the array. Am I right?If so, the function only removes the mean value if we pass
detrend="linear"
, right?True Phase and True Amplitude
I've been struggling to understand what
true_amplitude
andtrue_phase
actually mean.I think it would help to add some mathematical background on the intended behavior of enabling and disabling them.
Coordinates shift
Most of the FFT-related functions we will implement on Harmonica can be summarized in the following workflow:
While implementing the derivatives (fatiando/harmonica#238) I've noticed that
xrft.xrft.ifft
returns aDataArray
whose coordinates are centered around zero. Although this might not be a big deal in the general case, for the Harmonica use cases it's extremely important that the output grid has the same coordinates as the original one, since they mean actual geographical locations. I tried passing a tuple to thelag
parameter with the western and southern coordinates of the original region, butifft
returns a warning saying that it ignores this parameter.So far, we manage to implement this feature on our side: https://github.com/fatiando/harmonica/blob/4e865f00f0464e24d8f1cbb8b3c825ab70aa899d/harmonica/filters/fft.py#L34-L85
Nevertheless we don't intent to maintain wrappers of the
xrft
functions in the future.Would you like to see something like this implemented in
xrft
?Plans for future padding functions
In the future we will be very interested in having functions to apply padding to our 2D grids before passing them to
fft
. We've been following the implementation ofxarray.pad
, although it currently fills withnan
s the extended coordinates which cause a lot of problems when trying to compute FFTs (see pydata/xarray#3868). Sincexarray
is not intended to be used exclusively with regular grids, extending the coordinates of aDataArray
is not trivial. But sincexrft
assumes regular grids, will you be interested in implementing specific padding functions inxrft
?Thanks again for all your work on this subject! We really appreciate it.
And sorry for my long issue. If you think it might be better to split this Issue into several ones, please let me know and I'll split it so we can chat about every topic in a more organized way.
Cheers!
cc @leouieda @aguspesce @andieie
The text was updated successfully, but these errors were encountered: