Add equivalent layer for harmonic functions #78

santisoler · 2019-07-22T17:41:32Z

Add an Equivalent Layer gridder for harmonic functions, based on
verde.BaseGridder using 1/r as the basis function. By default it puts
one point source bellow each observation point at a constant relative
depth and fits the corresponding coefficient to each source. Also
a custom set of point sources can be passed. Add a gallery example
showing how to grid a small region of the magnetic anomaly data from Rio
de Janeiro.

Fixes #39

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.
Write detailed docstrings for all functions/methods.
If adding new functionality, add an example to the docstring, gallery, and/or tutorials.

Following Cooper (2000)

It will make easy to find future multiple EQL gridders by tab autocompletion.

Replace the distance_to_nearest function with verde.meadian_distance and add k_nearest as optional argument for the gridder.

This new test function is intended to be run with Numba jit disabled. A great number of masses will heavily increase the computation time if Numba jit is disabled.

santisoler · 2019-07-22T17:44:51Z

I want to use some functions from #71 when predicting interpolated values and computing distances in Cartesian coordinates. They must be incorporated after #71 is merged.
I'm still using the term "vertical" for specifying the downward coordinates. It should be changed in agreement with the term used in #71.

santisoler · 2019-07-22T17:46:08Z

It also uses some functions from the dev version of Verde. A new Verde release should be made in order to CIs to pass.

Because the vertical coordinate points downwards, the vertical component for the point sources must be increased if we want them to be bellow the observation points.

Co-Authored-By: Leonardo Uieda <leouieda@gmail.com>

leouieda · 2019-08-30T17:20:50Z

@santisoler I've been experimenting with the example and code a bit.

I isolated another anomaly that was better for visualizing. This is the result using the same parameters you were using but with a diverging colormap with 0 centered on white:

Notice that there are a lot of artifacts from the flight lines. The solution is also very sensitive to the depth_factor and k_nearest arguments.

These are the results with a larger grid spacing (400m) and a uniform depth (500m below data):

Notice that the stripes are gone. This due mainly to the larger spacing. The single depth is also easier to understand. I think the variable depth might be useful in other circumstances but probably not in the cases we're looking at right now. So it might be better to use the simpler case first and add the option of variable depth as we progress.

What do you think?

If you agree, I'll go ahead and commit my changes to the code and example.

santisoler · 2019-09-02T15:22:06Z

Hi @leouieda! I really like those results and agree to keep things simple and try more complex layouts in the future. Feel free to commit your changes!

Removed the code for variable depth based on neighbor distances. Use a different anomaly for the example and plot balancing the colormap. Use deeper sources for smoother grid and appropriate grid spacing. Simple train-test split instead of CV for quicker and easier to understand example.

stickler-ci · 2019-09-13T16:27:24Z

examples/eql_harmonic.py

+import numpy as np
+import pyproj
+import verde as vd
+import harmonica as hm


I001 isort found an import in the wrong position

stickler-ci · 2019-09-13T16:27:25Z

examples/eql_harmonic.py

+import pyproj
+import verde as vd
+import harmonica as hm
+


I004 isort found an unexpected blank line in imports

stickler-ci · 2019-09-13T16:27:25Z

examples/eql_harmonic.py

+
+# Plot gridded magnetic anomaly
+tmp = grid.magnetic_anomaly.plot.pcolormesh(ax=ax2, add_colorbar=False,
+                  vmin=-maxabs, vmax=maxabs, cmap="seismic")


E128 continuation line under-indented for visual indent

leouieda · 2019-09-13T16:28:27Z

@santisoler I pushed my modifications but it still needs some polish. I haven't checked the tests, the docstrings need some updating and polish, and the example plot could be nicer.

leouieda · 2019-09-13T16:28:41Z

I'll take a look at the prism PR tomorrow.

leouieda · 2019-09-14T11:55:36Z

@santisoler this is ready for a review. Here is what I've done:

Placed EQL examples in their own category. This one is about gridding and upward continuation.
Add Dampney 1969 citation
Remove the variable depth code to simplify this setup
Rework the testing to check interpolating from one grid to a denser grid
Example uses another anomaly, no cross-validation (to make things simpler and because we don't know if it actually works), balances the color maps, uses a larger spacing, and more text at the top explaining what EQLs do.
Class docstring lists what this EQL is good for and what it can't do. This will be useful since other layers can so other things and we can link them together this way.

leouieda · 2019-09-14T11:57:20Z

@birocoles it would be great if you took a look at this as well. Particularly the harmonica/equivalent_layer/harmonic.py and examples/eql/harmonic.py. I don't think we said anything terribly wrong but you're the real expert here 👨‍🎓.

santisoler

Nice work @leouieda! This looks awesome! Both the docstring and the gallery example are very clear and I give you one point for showing both the capabilities and the flaws of the gridder.

I would make clearer what the depth argument means: at first sight it seems that all point sources are located at a constant depth by default, but in fact they are set at a constant relative depth beneath each observation point. What do you think?

santisoler · 2019-09-16T13:39:40Z

.azure-pipelines.yml

@@ -92,6 +92,12 @@ jobs:
  - bash: echo "##vso[task.prependpath]$CONDA/bin"
    displayName: Add conda to PATH

+  # On Hosted macOS, the agent user doesn't have ownership of Miniconda's installation


Shouldn't we add these lines on a different PR?

leouieda · 2019-09-19T15:15:58Z

I give you one point for showing both the capabilities and the flaws of the gridder.

🥇 🎊

I would make clearer what the depth argument means: at first sight it seems that all point sources are located at a constant depth by default, but in fact they are set at a constant relative depth beneath each observation point. What do you think?

I completely agree. Would you mind adding this? Feel free to merge once that's done. We can always improve the docs but gotta get it out at some point ;-)

Shouldn't we add these lines on a different PR?

You caught me being lazy... 😞 I'll open the PR.

Improve docstrings to make clearer that the relative_depth argument represent the relative depth at which the point sources will be located referenced to the elevation of the data points.

santisoler · 2019-09-20T14:01:42Z

I renamed the depth argument to relative_depth and improved the docstrings to make it clearer what this argument represents. I've also add a little comment on the gallery example explaining this.
I'll merge this if CIs agree. I agree we can change things later.

santisoler added 27 commits June 28, 2019 17:13

Start experimenting with EQL

c714a85

Replace 2D for 3D

a08efbe

Change gridder for interpolator

4a0751a

Continue experimenting with EQL

1a1b3e5

Initialize data array with zeros instead of empty

01938f0

Improve some styling

d2920ef

Add docstrings for functions and methods

a7a1255

Simplify arguments for jacobian_numba

1a0fe92

Split distance computation on independent function

3db62d4

Define points at depth 3 times min distance

cd005be

Following Cooper (2000)

Change default depth to 3x dist to nearest point

17df1e7

Add depth_factor parameter

e56c8d7

Move points and depth_factor arguments to __init__

512056f

Rename point masses to sources and masses to coefs

2f798df

Rename gridder to EQLHarmonic

3d0e8a7

It will make easy to find future multiple EQL gridders by tab autocompletion.

Incorporate verde.median_distance

2138f4e

Replace the distance_to_nearest function with verde.meadian_distance and add k_nearest as optional argument for the gridder.

Import EQLHarmonic on harmonica/__init__.py

4293aa1

Use verde.base.n_1d_arrays on input coordinates

3f69ede

Change coordinates from list to tuple in predict

7aa8881

Start writing an example for EQLHarmonic

a4af096

Add cross-validation

c46f62f

Improve how grid is generated on example

3516889

Start drafting tests for EQLHarmonic

ec7834a

Add test function for jacobian_numba

6d1f2f7

Merge branch 'master' into eql

c5a009f

Add test function for EQLHarmonic with few sources

8519d78

This new test function is intended to be run with Numba jit disabled. A great number of masses will heavily increase the computation time if Numba jit is disabled.

Merge branch 'master' into eql

81bac5e

Change sign when setting default point sources

8cf45e4

Because the vertical coordinate points downwards, the vertical component for the point sources must be increased if we want them to be bellow the observation points.

Retitle API section where EQLHarmonic lives

855c469

Co-Authored-By: Leonardo Uieda <leouieda@gmail.com>

santisoler mentioned this pull request Sep 5, 2019

Implement the fast iterative EQL interpolator #95

Open

stickler-ci reviewed Sep 13, 2019

View reviewed changes

Add information to class docstring

1680a37

leouieda changed the title ~~Add Equivalent Layer interpolator for harmonic functions~~ Add equivalent layer for harmonic functions Sep 14, 2019

Separate EQL examples and add more text

63b54f8

fatiando deleted a comment from stickler-ci Sep 14, 2019

Rework tests to check the actual interpolation

9f9b9d2

Take ownership of azure conda path

89ed643

leouieda force-pushed the eql branch from 6a3d59b to 89ed643 Compare September 14, 2019 16:54

Test using custom sources as well

12f38aa

santisoler commented Sep 16, 2019

View reviewed changes

leouieda mentioned this pull request Sep 19, 2019

Take ownership of conda on Azure Mac builds #102

Merged

5 tasks

santisoler added 3 commits September 20, 2019 09:48

Merge branch 'master' into eql

5bad26e

Merge branch 'master' into eql

b1c9344

Rename depth variable to relative_depth

5129904

Improve docstrings to make clearer that the relative_depth argument represent the relative depth at which the point sources will be located referenced to the elevation of the data points.

santisoler merged commit a618988 into fatiando:master Sep 20, 2019

santisoler deleted the eql branch September 20, 2019 14:13

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Add equivalent layer for harmonic functions #78

Add equivalent layer for harmonic functions #78

santisoler commented Jul 22, 2019 •

edited

santisoler commented Jul 22, 2019

santisoler commented Jul 22, 2019

leouieda commented Aug 30, 2019

santisoler commented Sep 2, 2019

stickler-ci Sep 13, 2019

stickler-ci Sep 13, 2019

stickler-ci Sep 13, 2019

leouieda commented Sep 13, 2019

leouieda commented Sep 13, 2019

leouieda commented Sep 14, 2019

leouieda commented Sep 14, 2019

santisoler left a comment

santisoler Sep 16, 2019

leouieda commented Sep 19, 2019

santisoler commented Sep 20, 2019

Add equivalent layer for harmonic functions #78

Add equivalent layer for harmonic functions #78

Conversation

santisoler commented Jul 22, 2019 • edited

santisoler commented Jul 22, 2019

santisoler commented Jul 22, 2019

leouieda commented Aug 30, 2019

santisoler commented Sep 2, 2019

stickler-ci Sep 13, 2019

Choose a reason for hiding this comment

stickler-ci Sep 13, 2019

Choose a reason for hiding this comment

stickler-ci Sep 13, 2019

Choose a reason for hiding this comment

leouieda commented Sep 13, 2019

leouieda commented Sep 13, 2019

leouieda commented Sep 14, 2019

leouieda commented Sep 14, 2019

santisoler left a comment

Choose a reason for hiding this comment

santisoler Sep 16, 2019

Choose a reason for hiding this comment

leouieda commented Sep 19, 2019

santisoler commented Sep 20, 2019

santisoler commented Jul 22, 2019 •

edited