Add a minimalistic implementation of Euler Deconvolution

[dragger21]

Greetings from Indonesia.

Bringing the Euler Deconvolution back like in previous version (Python 2.7) should be interesting. It can predict/detect the depth of source anomaly in potential field (like gravity and magnetic).

[santisoler]

Hi @dragger21! Thanks for opening this issue!

Adding an Euler deconvolution implementation is in our future plans for Harmonica.
As I remember, @leouieda was having some ideas for it.

[leouieda]

Well, to start it would be great to have a minimalistic implementation: only the inversion solution, no moving windows and all that.

This should be less work to code and can then be used as a basis for more elaborate methods.

@dragger21 would be interested in implementing this? We're always happy to have more people involved in development 🙂

[dragger21]

To @leouieda and @santisoler, I'm sorry for the late reply.

I think that is very interesting to minimise implementation!
Quick question, how are we going to do that without any moving windows? Based on Reid et al. (1990) this work should include moving windows. Do you have any other literature that's pointed out that?

I actually very interested to implement this in your code. I feel very happy to be involved in the development.

[leouieda]

@dragger21 no worries, we're all busy and can relate to taking a while to respond 🙂 We'd be glad to have you involved in the development 🎉

The moving window strategy is a way to try to deal with multiple sources (in my opinion, it kind of fails to do so but that is for another discussion). But it's not intrinsic to the method and really not required if you have already isolated the target anomaly. If there is 1 anomaly in your data, it's actually worse if you do a moving window because the solutions for most windows are completely meaningless (as Alan Reid himself pointed out in several papers).

From a coding perspective, this is what I had in mind:

1. Implement a function or a class that does the inversion of Euler's homogeneity equation (no windows, only solve the equation and calculate the source position and base level). Let's call this EulerDeconv for example.
2. Implement a strategy that selects the anomalies and runs EulerDeconv on each anomaly. This could be a moving window strategy or something else entirely. But it will depend on the implementation in step 1.

So to get started, I would recommend implementing step 1 only. We can later add windows etc on top of this. It helps to break down the coding work into smaller tasks like this because:

1. Each step is quicker to implement and less code for us to review and give feedback on.
2. It's less code for you to write and more importantly, less documentation and tests (which is what most people struggle with the most when first starting out).
3. If you get busy later, someone else can continue the work and you still get credit for adding one piece of the code.

[leouieda]

What I had in mind for the EulerDeconv class is something like this:

class EulerDeconv:
    """
    ...
    """
    def __init__(self, structural_index):
        self.structural_index = structural_index
        # The estimated parameters. Start them with None
        self.location_ = None
        self.base_level_ = None

    def fit(self, coordinates, data):
        # This function runs the inversion on the input data
        # data = (field, east_deriv, north_deriv, up_deriv)
        self.location_ =  the estimated location
        self.base_level_ = the estimated base level

This class runs the inversion and stores the estimated source position and base level. It would be used like this:

# Load some data. I'm assuming that they've been projected already.
data = pandas.read_csv("some-data-file.csv")

# Calculate derivatives with equivalent sources
eqs = hm.EquivalentSources(depth=..., damping=...)
eqs.fit((data.easting, data.northing, data.height), data.total_field_anomaly)
delta = 1  # finite difference step in meters
data["east_deriv"] = (
    eqs.predict((data.easting + delta, data.northing, data.height))
    - eqs.predict((data.easting - delta, data.northing, data.height))
) / 2 * delta
data["north_deriv"] = (
    eqs.predict((data.easting, data.northing + delta, data.height))
    - eqs.predict((data.easting, data.northing - delta, data.height))
) / 2 * delta
data["up_deriv"] = (
    eqs.predict((data.easting, data.northing, data.height + delta))
    - eqs.predict((data.easting, data.northing, data.height))
) / delta

# Run Euler deconvolution. Let's say the data is a dipole anomaly of an intrusion.
euler = hm.EulerDeconv(structural_index=3)
euler.fit(
    coordinates=(data.easting, data.northing, data.height), 
    data=(data.total_field_anomaly, data.east_deriv, data.north_deriv, data.up_deriv),
)
print(euler.location_)
print(euler.base_level_)

[dragger21]

I apologize again for the long reply. Thanks for the insight, Leo & Santiago. I'm trying to seek and simplify the code as soon as possible. I'll notify you later.

[leouieda]

@dragger21 no rush! Please reach out if you need any help

[AnthonyDunk]

Has there been any progress made on adding Euler Deconvolution to Harmonica yet? I'm interested in using it to run this example, but since I have Python 3.10, I'm out of luck at the moment: https://legacy.fatiando.org/gallery/gravmag/euler_moving_window.html
Thanks.

[nikita-k0v]

bump, any progress on this?