harmonica icon indicating copy to clipboard operation
harmonica copied to clipboard

Inverted sign in upward derivative filter

Open leouieda opened this issue 5 months ago • 1 comments

Description of the problem:

The upward derivative returned has an inverted sign. It looks OK when inspecting plots because a lot of us are used to the sign convention of positive downwards. But since we use heights (and thus positive upward), the sign should be the opposite of what we expect. This is meant to improve consistency between small-scale and global data (which often has a "radial" component that is positive upward).

The fix should be relatively simple. Just add a minus sign to https://github.com/fatiando/harmonica/blob/main/harmonica/filters/_filters.py#L63

Full code that generated the error

import harmonica as hm
import verde as vd
import numpy as np
import matplotlib.pyplot as plt

model = [
    [35e3, 40e3, 50e3, 55e3, -6e3, -1e3],
]
magnetization = [
    hm.magnetic_angles_to_vec(20, -40, 15),
]
inc, dec = -40, 15
fe, fn, fu = hm.magnetic_angles_to_vec(1, inc, dec)
region = [0, 100e3, 0, 80e3]
coordinates = vd.grid_coordinates(region, spacing=500, extra_coords=1000)

be, bn, bu = hm.prism_magnetic(coordinates, model, np.transpose(magnetization), field="b")
anomalia = fe * be + fn * bn + fu * bu

grid = vd.make_xarray_grid(coordinates, anomalia, data_names="tfa", extra_coords_names="upward")
grid["d_up"] = hm.derivative_upward(grid.tfa)

fig, axes = plt.subplots(1, 2, layout="constrained", figsize=(10, 3))
grid.tfa.plot(ax=axes[0])
grid.d_up.plot(ax=axes[1])

Full error message

This produces the following result:

image

If we use positive upward, then the positive parts of the anomaly should be decreasing towards 0 (negative derivative) and the negative parts of the anomaly should be increasing towards 0 (positive derivative).

System information

  • Operating system: Linux
  • Python installation (Anaconda, system, ETS): Miniforge
  • Version of Python: 3.11
  • Version of this package: current main (bbaf5a1)
  • If using conda, paste the output of conda list below: Not using conda.
Output of conda list

PASTE OUTPUT OF CONDA LIST HERE

leouieda avatar Feb 01 '24 11:02 leouieda

The test didn't fail because this one https://github.com/fatiando/harmonica/blob/main/harmonica/tests/test_transformations.py#L246 compares the derivative of the potential with gz which is positive downward. The right thing to do is compare with -gz.

leouieda avatar Feb 01 '24 17:02 leouieda