Reduction to the pole of a magnetic anomaly grid¶

Out:

Reduced to the pole magnetic grid:
 <xarray.DataArray (northing: 370, easting: 346)>
array([[  14.15556282,   10.38429905,   10.02003492, ..., -219.81018247,
        -210.93025635, -179.38408421],
       [  -3.21243279,   -9.37121782,  -10.96584477, ..., -165.15303656,
        -158.0958179 , -133.29352054],
       [  -2.3217145 ,   -9.44938578,  -11.35228505, ..., -170.79735965,
        -165.25667548, -141.046265  ],
       ...,
       [  45.45699147,  -24.80993602,  -51.27393817, ...,  -40.42604764,
         -64.12375288,  -75.97556361],
       [  36.91047424,  -37.13717168,  -58.40649469, ...,  -34.55576242,
         -55.65612416,  -71.01718934],
       [-102.42457874, -155.67864799, -165.96649111, ...,  -36.95819436,
         -35.04014832,  -40.15060688]])
Coordinates:
  * northing  (northing) float64 7.576e+06 7.576e+06 ... 7.595e+06 7.595e+06
  * easting   (easting) float64 4.655e+05 4.656e+05 ... 4.827e+05 4.828e+05
<IPython.core.display.Image object>

import ensaio
import pygmt
import verde as vd
import xarray as xr
import xrft

import harmonica as hm

# Fetch magnetic grid over the Lightning Creek Sill Complex, Australia using
# Ensaio and load it with Xarray
fname = ensaio.fetch_lightning_creek_magnetic(version=1)
magnetic_grid = xr.load_dataarray(fname)

# Pad the grid to increase accuracy of the FFT filter
pad_width = {
    "easting": magnetic_grid.easting.size // 3,
    "northing": magnetic_grid.northing.size // 3,
}
# drop the extra height coordinate
magnetic_grid_no_height = magnetic_grid.drop_vars("height")
magnetic_grid_padded = xrft.pad(magnetic_grid_no_height, pad_width)

# Define the inclination and declination of the region by the time of the data
# acquisition (1990).
inclination, declination = -52.98, 6.51

# Apply a reduction to the pole over the magnetic anomaly grid. We will assume
# that the sources share the same inclination and declination as the
# geomagnetic field.
rtp_grid = hm.reduction_to_pole(
    magnetic_grid_padded, inclination=inclination, declination=declination
)

# Unpad the reduced to the pole grid
rtp_grid = xrft.unpad(rtp_grid, pad_width)

# Show the reduced to the pole grid
print("\nReduced to the pole magnetic grid:\n", rtp_grid)


# Plot original magnetic anomaly and the reduced to the pole
fig = pygmt.Figure()
with fig.subplot(nrows=1, ncols=2, figsize=("28c", "15c"), sharey="l"):
    # Make colormap for both plots (saturate it a little bit)
    scale = 0.5 * vd.maxabs(magnetic_grid, rtp_grid)
    pygmt.makecpt(cmap="polar+h", series=[-scale, scale], background=True)
    with fig.set_panel(panel=0):
        # Plot magnetic anomaly grid
        fig.grdimage(
            grid=magnetic_grid,
            projection="X?",
            cmap=True,
        )
        # Add colorbar
        fig.colorbar(
            frame='af+l"Magnetic anomaly [nT]"',
            position="JBC+h+o0/1c+e",
        )
    with fig.set_panel(panel=1):
        # Plot upward reduced to the pole grid
        fig.grdimage(grid=rtp_grid, projection="X?", cmap=True)
        # Add colorbar
        fig.colorbar(
            frame='af+l"Reduced to the pole [nT]"',
            position="JBC+h+o0/1c+e",
        )
fig.show()

Total running time of the script: ( 0 minutes 6.581 seconds)

Gallery generated by Sphinx-Gallery

Harmonica v0.6.0

Reduction to the pole of a magnetic anomaly grid

Reduction to the pole of a magnetic anomaly grid¶