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

Implement tilt filter #261

Merged
merged 17 commits into from Jul 15, 2016

Conversation

Projects
None yet
2 participants
@mtb-za
Contributor

mtb-za commented Mar 23, 2016

This will create a filter to determine the tilt of a potential field, as developed by Miller and Singh (1994).

I could use some suggestions on how to create a sensible test for this.

This is the first of a few potential field filters that I want to implement.

Checklist:

  • Make tests for new code
  • Create/update docstrings
  • Include relevant equations and citations in docstrings
  • Code follows PEP8 style conventions
  • Code and docs have been spellchecked
  • Include new dependencies in docs, requirements.txt, README, and .travis.yml
  • Documentation builds properly
  • All tests pass
  • Can be merged
  • Changelog entry (leave for last)
@mtb-za

This comment has been minimized.

Contributor

mtb-za commented Mar 29, 2016

@leouieda Do you know why the CI is failing with The command "pip install -q --use-mirrors coverage==3.7.1 coveralls==0.5 pep8==1.6.0 sphinx_bootstrap_theme==0.4.5" failed and exited with 2 during .

I do not think that I have touched anything that should affect that....

@leouieda

This comment has been minimized.

Member

leouieda commented Apr 1, 2016

@mtb-za thanks for taking the time to do this!

I'll take a look at the Travis log and see if I can figure out what is going on.

As for the tests, does the original paper mention any analytical solutions? Maybe for a point source or something like that. That is always a nice test case. If not, then testing for sanity is already something (like checking if values are within the allowed bounds and other conditions).

@leouieda

This comment has been minimized.

Member

leouieda commented Apr 1, 2016

@mtb-za seems that your branch is way behind master. Do you remember if you pulled in the latest changes from fatiando/fatiando before creating your tilt_filer branch?

If not, then you should update your master branch and also this branch. You can do that with the following:

  1. Backup your repository somewhere, just in case.

  2. Register this repository as a remote option. We'll call it upstream:

    git remote add upstream https://github.com/fatiando/fatiando.git
    
  3. Pull in the changes from the upstream master branch to your own:

    git checkout master
    git pull upstream master
    
  4. Now we can rebase your tilt_filer branch to be based on the updated master:

    git checkout tilt_filer
    git rebase master
    
  5. Rebasing should not give any conflicts for now. If it does, you'll need to resolve them (let me know). Rebasing will change your commits, so the ones you have now are not the ones you have on Github. You'll need to force push the new commits:

    git push -f origin tilt_filer
    

Try this out and let me know if it works.

Also, keep in mind that you should always pull changes from upstream master before creating your branches. And always make your branchs from master.

@leouieda leouieda referenced this pull request Apr 1, 2016

Closed

Minor changes #262

@mtb-za mtb-za force-pushed the mtb-za:tilt_filter branch from 15a7b44 to 39e0b31 Apr 5, 2016

@mtb-za mtb-za closed this Apr 5, 2016

@mtb-za mtb-za force-pushed the mtb-za:tilt_filter branch from 39e0b31 to 237ba1d Apr 5, 2016

@mtb-za

This comment has been minimized.

Contributor

mtb-za commented Apr 5, 2016

I have no idea how that happened. I am pretty sure that I had pulled in the upstream stuff before going ahead.

Anyway, I seem to have a clean base off master and have restarted.

@mtb-za mtb-za reopened this Apr 5, 2016

@mtb-za

This comment has been minimized.

Contributor

mtb-za commented Apr 5, 2016

It looks like it will need to be a sanity check. There is no analytical solution within the original paper.

The tilt angle overcomes this problem by dealing
with the ratio of the vertical derivative to the horizontal
gradient. Since both will be smaller for deeper sources,
the ratio will still be large when over the source, pass
through zero over, or near, the edge where the vertical
derivative is zero and the horizontal gradient is a maximum,
and will be negative outside the body where the
vertical derivative is negative. By expressing this as a
tilt angle rather than a ratio it will always be in the range
-90°< TILT< 90 °.

Maybe a relatively simple test that the tilt does return a positive value over a body, zero at the edge, and negative outside the body will work?

@leouieda

This comment has been minimized.

Member

leouieda commented Apr 5, 2016

@mtb-za git can be mysterious sometimes always. Reminds me of xkcd:

@leouieda

This comment has been minimized.

Member

leouieda commented Apr 5, 2016

Maybe a relatively simple test that the tilt does return a positive value over a body, zero at the edge, and negative outside the body will work?

That might be useful. Maybe use a single rectangular prism with data almost right on top of it. 0 height for the observations might cause problems but maybe have it a 10 meters or so. In that case the tilt should behave as advertised.

Would be good to have a sanity check that everything is within the -90, 90 degree range, even for more complex data.

@@ -262,6 +263,36 @@ def tga(x, y, data, shape, method='fd'):
return res
def tilt(x, y, data, shape):
"""
Calculates the magnetic tilt, as defined by Miller and Singh (1994):

This comment has been minimized.

@leouieda

leouieda Apr 5, 2016

Member

Please add a blank line between each section.

This comment has been minimized.

@leouieda

leouieda Apr 5, 2016

Member

Maybe not "magnetic tilt", as this can be used for other things as well? "potential field tilt angle" perhaps?

Calculates the magnetic tilt, as defined by Miller and Singh (1994):
.. math::
tilt(f) = tan^{-1}(\\frac{\\frac{dT}{dz}}{\\sqrt{\\frac{dT}{dx}^2 +
\\frac{dT}{dy}^2}})

This comment has been minimized.

@leouieda

leouieda Apr 5, 2016

Member

I think the equation code needs to be idented inside the .. math:: block.

Returns:
* tilt : 1D-array
The tilt angle of the total field.
**References**

This comment has been minimized.

@leouieda

leouieda Apr 5, 2016

Member

I had to check this because I never remember, the section headers (Parameters, Returns, References) are all normal case (no bold).

The docstring conventions need some updating. I'm looking into using the numpy style docstrings because they look a bit better. It'll be some work to convert all of them but might be worth doing.

doi:10.1016/0926-9851(94)90022-1.
"""
horiz_deriv = numpy.sqrt( derivx(x, y, data, shape) +\
derivy(x, y, data, shape)

This comment has been minimized.

@leouieda

leouieda Apr 5, 2016

Member

Missing a ) here. Also, if the statement is inside () you don't need to break the line with a \. Example:

horiz_deriv = numpy.sqrt(derivx(x, y, data, shape) +
                         derivy(x, y, data, shape))
@mtb-za

This comment has been minimized.

Contributor

mtb-za commented Apr 6, 2016

Hmmmm.... There is definitely something odd going on here.

I am getting a lot of NaN results, after making a simple example similar to the TGA cookbook. I will have to do some checking as to where those might be creeping in.

Miller, Hugh G, and Vijay Singh. 1994. "Potential Field Tilt --- a New
Concept for Location of Potential Field Sources."
Journal of Applied Geophysics 32 (2--3): 213-17.
doi:10.1016/0926-9851(94)90022-1.

This comment has been minimized.

@leouieda

leouieda Apr 6, 2016

Member

@mtb-za nice! Thanks for doing this so fast.

The x and y coordinates of the grid points
* data : 1D-array
The potential field at the grid points
* shape : tuple = (ny, nx)

This comment has been minimized.

@leouieda

leouieda Apr 6, 2016

Member

The shape is (nx, ny) for consistency with x being the North direction.

strange units and so will the potential field tilt! I strongly
recommend converting the data to SI **before** calculating the
TGA is you need the gradient in Eotvos (use one of the unit conversion
functions of :mod:`fatiando.utils`).

This comment has been minimized.

@leouieda

leouieda Apr 6, 2016

Member

Is this really a problem? The tilt is dimensionless so the units the data are in shouldn't matter much, so long as they are compatible.

doi:10.1016/0926-9851(94)90022-1.
"""
horiz_deriv = numpy.sqrt( derivx(x, y, data, shape) +
derivy(x, y, data, shape) )

This comment has been minimized.

@leouieda

leouieda Apr 6, 2016

Member

Shouldn't the derivatives be squared? Maybe this is where the NaNs are coming from?

This comment has been minimized.

@mtb-za

mtb-za Apr 7, 2016

Contributor

D'oh. That fixed it. Not sure how I missed it, given how often I have written this particular equation out.

This comment has been minimized.

@leouieda

leouieda Apr 7, 2016

Member

Don't worry, that's why we do code review. It's amazing how many bugs we miss when we're head deep in something.

@leouieda leouieda added this to the 0.5 milestone Apr 6, 2016

@leouieda

This comment has been minimized.

Member

leouieda commented Apr 6, 2016

There is a build error caused by a doctest in fatiando.inversion. I'm fixing that in #264. I'll let you know as soon as that is done so that you can git rebase this off master again.

@leouieda

This comment has been minimized.

Member

leouieda commented Apr 6, 2016

@mtb-za #264 was merged and tests are working now. When you get a chance, please do a git rebase of this branch to get the new changes.

@mtb-za mtb-za force-pushed the mtb-za:tilt_filter branch from d7621a7 to 23c7c9f Apr 7, 2016

@leouieda

This comment has been minimized.

Member

leouieda commented Apr 7, 2016

@mtb-za let me know you want me to look at this again. Nice work!

@mtb-za

This comment has been minimized.

Contributor

mtb-za commented Apr 8, 2016

I could use some help with setting up the test we discussed.

The paper does mention that there is an analytic solution, by the way, but not what it is.

@leouieda

This comment has been minimized.

Member

leouieda commented Apr 8, 2016

@mtb-za testing is hard. I haven't gotten that used to it yet as well. Specially when you don't have a good comparison.

This is the case with the prism forward modeling. It's not practical to calculate the solution by hand. So what I did was run a lot of sanity checks and compare against a pure numpy implementation. See test/test_gravmag_prism.py. A lot of those tests I wrote as consequence of a bug.

Some ideas for the tilt:

  • Should all be in a certain angle range (is it -90-90?). So you can calculate for a single prism model and check if all are in the range.
  • Should be symmetric if the prism is a cube and in the middle of the data grid. Both sides of the grid should have the same tilt right?
  • Check if the code fails when it should. The shape multiplied should be the size of the arrays. The arrays should be equal size. All these should be checked with an assert at the start of the function. You tests should pass bad input and check if an AssertionError is raised.
@leouieda

This comment has been minimized.

Member

leouieda commented Apr 13, 2016

@mtb-za just thought of this and wanted to write it down:

It would be good if the function had the option to receive the x, y, z derivatives as input. They could default to None and be calculated by derivx etc if that is the case. This would be good if the user has the derivatives somehow, maybe measured or calculated some other way (like equivalent layer).

This would also be handy when testing. We could test the function passing it analytical derivatives (calculated from the model) and check against the numerical result. It should be the same if the anomaly is error free and not truncated.

mtb-za added some commits Apr 5, 2016

Implemented basic tilt filter
This will hopefully, finally, be a working implementation for the tilt
filter.
Started testing tilt filter
Created simple cookbook recipe for tilt filter.
Corrected maths error
Blindly forgot to square the derivatives.
This seems to work now.
Removed numpy requirement in tilt cookbook recipe
Do not really need it, so have taken it out.

Also altered the inclination, which I had changed for testing (and which
was fixed in the last commit).
Removed requirement for data to be SI units
As pointed out, the tilt is dimensionless. As long as things are all in
the same system, what that system is should not matter.
Edit tilt to accept pre-calculated derivatives
Following on from
#261 (comment)

Implemented these changes.
Started implementing tests for tilt
Simple test to check that values are returned within a sane range: -90
to +90 degrees.

If a value is outside of this, something has gone wrong.
Fixed typo
Typo in calculating derivatives.
Fixed typo
Typo in calculating derivatives.

Deleted scratch test folder

This was a local test/scratch file that got pushed accidentally.
Deleting it now.

Remove test_tilt.py

Scratch file that got committed by mistake.

Implementing tilt test

Basic test for sane ranges.
Deleted scratch test folder
This was a local test/scratch file that got pushed accidentally.
Deleting it now.
Pushing latest test and cookbook files
Latest versions of testing and cookbook for tilt filters.
Additional test for tilt
Tests random points inside and outside a prism to ensure that the values
are sensible. Positive values should only be above the anomaly, negative
ones outside it.

We take a series of random points, find vertices closest to them and
then get the values of those points in the tilted data. If the values
are as expected, then the code passes.

@leouieda leouieda force-pushed the mtb-za:tilt_filter branch from ffe4398 to 6e21b2e Jul 14, 2016

@leouieda

This comment has been minimized.

Member

leouieda commented Jul 14, 2016

@mtb-za I ran a git rebase master on this branch and got a sane commit history now. This is something I've struggled with before. It's one the things that can trip you up in git. But it should be fine now.

Beware though that the tests have moved to fatiando/tests and that the tests are run with py.test. Take a look at the Makefile or the latest dev docs.

data = prism.tf(x, y, z, model, inc, dec)
tilt_angle = transform.tilt(x, y, data, shape)
mpl.figure()

This comment has been minimized.

@leouieda

leouieda Jul 14, 2016

Member

Please avoid the mpl module as a replacement for matplotlib.pyplot. This will be removed in the future. I really need to write this down somewhere or deprecate it already.

def test_potential_field_tilt_works():
"gravmag.transform tilt returns sane values"
model = [Prism(-1000, 1000, -500, 500, 0, 2000, {'magnetization': 3200})]

This comment has been minimized.

@leouieda

leouieda Jul 14, 2016

Member

Another thing that will be deprecated is passing in a scalar magnetization. This is implicitly saying that there is only induced magnetization. We'll be changing the gravmag functions to accept only vector magnetization so that this is explicit.

@leouieda

This comment has been minimized.

Member

leouieda commented Jul 15, 2016

I'll make some changes to your tests just to get things moving a bit faster. I guess you're pretty sick of looking at the tilt by now 😄

#Deal with horizontal derivatives. If we have not been given an array
#with either, we will calculate them first, then get the horizontal
#derivative.
if xderiv == None:

This comment has been minimized.

@leouieda

leouieda Jul 15, 2016

Member

Avoid comparing to None with == or !=. The proper way is to use xderiv is None or xderiv is not None.

@leouieda

This comment has been minimized.

Member

leouieda commented Jul 15, 2016

Would it be worth it returning the tilt in degrees instead of radians?

Two tests to cover 100% of tilt code
The tests check if values are not absurd (out of the 90:-90 degree
range) and if the results are the same when given analytical derivatives
of the gravity field. This covers 100% of the function code and is good
enough for such a small function.

Replace the two cookbook recipes with a gallery plot using already
reduced to the pole magnetic data. Plot the zero contour of the tilt as
well. This will help with the testing and looks good enough.
@leouieda

This comment has been minimized.

Member

leouieda commented Jul 15, 2016

@mtb-za I made slight changes to your tests. The first tests for sane results and I removed the pole reduction for simplicity. It seems that the pole reduction is adding a lot of noise to the tilt, so it's better to test without it. The second tests checks if the tilt results are the same when given analytical derivatives. I'm using gravity for this case. Both cover 100% of the function code so they are good enough.

I also transformed your cookbook recipes into a gallery plot. I tried doing the pole reduction on this but it looks terrible. The problem is with the pole reduction, not the tilt. So I'm using data at the pole for the gallery and linking it to reduce_to_pole for completeness.

Lastly, I made some PEP8 adjustments that were causing Travis to fail.

All this needs now is the changelog entry. Would you do the honors? Then we can merge this and be done with it!

Updated changelog
Added changelog entry for tilt filter (*Happy dance!*).

Also tweaked a couple typos I noticed while looking for an example of
previously added functionality.
@mtb-za

This comment has been minimized.

Contributor

mtb-za commented Jul 15, 2016

Changelog is done, even if you did most of the work. Glad to finally lay this one to rest though!

Onto the next thing: depth estimations. panics

@leouieda

This comment has been minimized.

Member

leouieda commented Jul 15, 2016

@mtb-za non-sense! You did all the preparation and investigating what the tilt should actually do. I just polished it up. Truth is, I'm also learning all this as I go. I cringe every time I look at code that's > 1 month old.

I'm looking forward to the next PR! Keep 'em coming!

@leouieda leouieda merged commit b5d3828 into fatiando:master Jul 15, 2016

2 checks passed

continuous-integration/travis-ci/pr The Travis CI build passed
Details
coverage/coveralls Coverage increased (+0.1%) to 72.622%
Details

@mtb-za mtb-za deleted the mtb-za:tilt_filter branch Jul 15, 2016

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment