fatiando / harmonica

Forward modeling, inversion, and processing gravity and magnetic data
https://www.fatiando.org/harmonica
BSD 3-Clause "New" or "Revised" License
208 stars 68 forks source link

Better default for depth in equivalent sources #424

Closed leouieda closed 2 months ago

leouieda commented 1 year ago

Description of the desired feature:

From chats with @indiauppal and @birocoles, I learned that Dampney (1969; the original equivalent source paper) advises using depths between 2.5 and 6 times the source spacing (assuming sources in a regular grid). For sources beneath data, like we use here, a sensible default is 4.5 times the median distance between sources. Me and @indiauppal tested this default and it works well most of the time. Tweaking is always good but this would provide a more sensible default than the fixed value we currently use.

This would break backwards compatibility. I'm not too worried since our docs say that you should always set the depth by hand. So changing the default will likely make the interpolation better for the few people who were using the classes without following our advice. But I can do the warning dance if others think it's necessary (but it will take longer to implement).

Are you willing to help implement and maintain this feature?

Yes. But very happy to let someone else do it (😉 @indiauppal)

leouieda commented 6 months ago

@santisoler would you like me to first make a PR with a warning that the default will change that can go into v0.7 and then we make the change for v0.8?

santisoler commented 4 months ago

Just to clarify, the depth should be set to 4.5 times the median distance between first neighbor sources, right?

If that's true, then we can use:

depth = np.median(verde.median_distance(coordinates, k_nearest=1))

Is this the idea?

indiauppal commented 4 months ago

Should we take the mean of the median distances? depth = np.mean(verde.median_distance(coordinates, k_nearest=1)) * 4.5

leouieda commented 4 months ago

I guess the mean would be better? Then it wouldn't favor close points. Perhaps it would be good to have a verde.mean_distance function as well.