# Grid transformations¶

Harmonica offers some functions to apply FFT-based (Fast Fourier Transform) and finite-differences transformations to regular grids of gravity and magnetic fields located at a constant height.

In order to apply these grid transformations, we first need a regular grid in Cartesians coordinates. Let’s download a magnetic anomaly grid over the Lightning Creek Sill Complex, Australia, readily available in `ensaio`. We can load the data file using `xarray`:

```import ensaio
import xarray as xr

fname = ensaio.fetch_lightning_creek_magnetic(version=1)
magnetic_grid
```
```<xarray.DataArray 'total_field_anomaly' (northing: 370, easting: 346)>
array([[  35.        ,   36.19995   ,   36.69995   , ..., -101.099976  ,
-100.400024  ,  -99.599976  ],
[  36.5       ,   37.599976  ,   38.        , ..., -102.20007   ,
-101.5       , -100.70007   ],
[  37.099976  ,   38.19995   ,   38.599976  , ..., -103.30005   ,
-102.599976  , -101.900024  ],
...,
[ 182.79993   ,  172.40002   ,  160.79993   , ...,    0.79992676,
-24.200073  ,  -41.80005   ],
[ 182.09998   ,  172.59998   ,  161.40002   , ...,    6.        ,
-21.5       ,  -41.        ],
[ 178.79993   ,  170.40002   ,  160.29993   , ...,   11.400024  ,
-16.        ,  -35.80005   ]], dtype=float32)
Coordinates:
* easting   (easting) float64 4.655e+05 4.656e+05 ... 4.827e+05 4.828e+05
* northing  (northing) float64 7.576e+06 7.576e+06 ... 7.595e+06 7.595e+06
height    (northing, easting) float64 500.0 500.0 500.0 ... 500.0 500.0
Attributes:
Conventions:   CF-1.8
title:         Magnetic total-field anomaly of the Lightning Creek sill c...
crs:           proj=utm zone=54 south datum=WGS84 units=m no_defs ellps=W...
source:        Interpolated from airborne magnetic line data using gradie...
references:    Geophysical Acquisition & Processing Section 2019. MIM Dat...
long_name:     total-field magnetic anomaly
units:         nT
actual_range:  [-1785.  3798.]```

And plot it:

```import matplotlib.pyplot as plt

plt.gca().set_aspect("equal")
plt.title("Magnetic anomaly grid")
plt.gca().ticklabel_format(style="sci", scilimits=(0, 0))
plt.colorbar(tmp, label="nT")
plt.show()
```

In case we have a regular grid defined in geographic coordinates (longitude, latitude) we can project them to Cartesian coordinates using the `verde.project_grid` function and a map projection like the ones available in `pyproj`.

Since all the grid transformations we are going to apply are based on FFT methods, we usually want to pad them in order their increase the accuracy. We can easily do it through the `xrft.pad` function. First we need to define how much padding we want to add along each direction. We will add one third of the width and height of the grid to each side:

```pad_width = {
"easting": magnetic_grid.easting.size // 3,
"northing": magnetic_grid.northing.size // 3,
}
```

And then we can pad it, but dropping the `height` coordinate first (this is needed by the `xrft.pad` function):

```import xrft

magnetic_grid_no_height = magnetic_grid.drop_vars("height")
```
```<xarray.DataArray 'total_field_anomaly' (northing: 616, easting: 576)>
array([[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
...,
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.]], dtype=float32)
Coordinates:
* easting   (easting) float64 4.598e+05 4.598e+05 ... 4.885e+05 4.885e+05
* northing  (northing) float64 7.57e+06 7.57e+06 ... 7.601e+06 7.601e+06
Attributes:
Conventions:   CF-1.8
title:         Magnetic total-field anomaly of the Lightning Creek sill c...
crs:           proj=utm zone=54 south datum=WGS84 units=m no_defs ellps=W...
source:        Interpolated from airborne magnetic line data using gradie...
references:    Geophysical Acquisition & Processing Section 2019. MIM Dat...
long_name:     total-field magnetic anomaly
units:         nT
actual_range:  [-1785.  3798.]```
```tmp = magnetic_grid_padded.plot(cmap="seismic", center=0, add_colorbar=False)
plt.gca().set_aspect("equal")
plt.gca().ticklabel_format(style="sci", scilimits=(0, 0))
plt.colorbar(tmp, label="nT")
plt.show()
```

Now that we have the padded grid, we can apply any grid transformation.

## Upward derivative¶

Let’s calculate the upward derivative (a.k.a. vertical derivative) of the magnetic anomaly grid using the `harmonica.derivative_upward` function:

```import harmonica as hm

deriv_upward
```
```<xarray.DataArray (northing: 616, easting: 576)>
array([[ 1.33170010e-05,  9.52432798e-06,  5.67172025e-06, ...,
2.43398647e-05,  2.07242504e-05,  1.70500748e-05],
[ 2.61432378e-05,  2.24094583e-05,  1.86184385e-05, ...,
3.70069736e-05,  3.34414963e-05,  2.98204704e-05],
[ 3.90050348e-05,  3.53287011e-05,  3.15981781e-05, ...,
4.97132608e-05,  4.61968158e-05,  4.26275223e-05],
...,
[-2.49540481e-05, -2.89307684e-05, -3.29763487e-05, ...,
-1.34304860e-05, -1.72044538e-05, -2.10453870e-05],
[-1.22311478e-05, -1.61452928e-05, -2.01251523e-05, ...,
-8.78416548e-07, -4.59825256e-06, -8.38238562e-06],
[ 5.25638946e-07, -3.32701915e-06, -7.24270765e-06, ...,
1.17116107e-05,  8.04443285e-06,  4.31601166e-06]])
Coordinates:
* northing  (northing) float64 7.57e+06 7.57e+06 ... 7.601e+06 7.601e+06
* easting   (easting) float64 4.598e+05 4.598e+05 ... 4.885e+05 4.885e+05```

This grid includes all the padding we added to the original magnetic grid, so we better unpad it using `xrft.unpad`:

```deriv_upward = xrft.unpad(deriv_upward, pad_width)
deriv_upward
```
```<xarray.DataArray (northing: 370, easting: 346)>
array([[ 0.95819797,  0.62479631,  0.65249329, ..., -1.73446106,
-1.67664073, -2.7243531 ],
[ 0.63634155,  0.21904983,  0.23107703, ..., -0.49049877,
-0.45948262, -1.68410265],
[ 0.66359221,  0.23536193,  0.2450619 , ..., -0.51034902,
-0.49225175, -1.75482593],
...,
[ 3.3946594 ,  0.9299787 ,  0.84907987, ..., -0.18739683,
-0.37947336, -1.13012159],
[ 3.28895305,  0.89679032,  0.84612464, ..., -0.15550245,
-0.36489469, -1.12153511],
[ 5.04819984,  2.91262219,  2.80733139, ...,  0.11714946,
-0.38706034, -1.26040298]])
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```

And plot it:

```tmp = deriv_upward.plot(cmap="seismic", center=0, add_colorbar=False)
plt.gca().set_aspect("equal")
plt.title("Upward derivative of the magnetic anomaly")
plt.gca().ticklabel_format(style="sci", scilimits=(0, 0))
plt.colorbar(tmp, label="nT/m")
plt.show()
```

## Horizontal derivatives¶

We can also compute horizontal derivatives over a regular grid using the `harmonica.derivative_easting` and `harmonica.derivative_northing` functions.

```deriv_easting = hm.derivative_easting(magnetic_grid)
deriv_easting
```
```<xarray.DataArray 'total_field_anomaly' (northing: 370, easting: 346)>
array([[ 0.02399902,  0.01699951,  0.00400024, ...,  0.01599976,
0.015     ,  0.01600098],
[ 0.02199951,  0.015     ,  0.00400024, ...,  0.015     ,
0.015     ,  0.01599853],
[ 0.02199951,  0.015     ,  0.00300049, ...,  0.015     ,
0.01400024,  0.01399902],
...,
[-0.20799805, -0.22      , -0.22100097, ..., -0.5530005 ,
-0.42599976, -0.35199952],
[-0.19      , -0.20699951, -0.21800049, ..., -0.6019995 ,
-0.47      , -0.39      ],
[-0.16799805, -0.185     , -0.2       , ..., -0.5969995 ,
-0.47200072, -0.39600098]], dtype=float32)
Coordinates:
* easting   (easting) float64 4.655e+05 4.656e+05 ... 4.827e+05 4.828e+05
* northing  (northing) float64 7.576e+06 7.576e+06 ... 7.595e+06 7.595e+06
height    (northing, easting) float64 500.0 500.0 500.0 ... 500.0 500.0```
```deriv_northing = hm.derivative_northing(magnetic_grid)
deriv_northing
```
```<xarray.DataArray 'total_field_anomaly' (northing: 370, easting: 346)>
array([[ 0.03      ,  0.02800049,  0.02600098, ..., -0.02200195,
-0.02199951, -0.02200195],
[ 0.02099976,  0.02      ,  0.01900024, ..., -0.02200073,
-0.02199951, -0.02300049],
[ 0.00599976,  0.005     ,  0.005     , ..., -0.02299927,
-0.02400024, -0.025     ],
...,
[-0.00199951,  0.01400024,  0.02200073, ...,  0.10300049,
0.03599976, -0.015     ],
[-0.04      , -0.02      , -0.005     , ...,  0.10600097,
0.08200073,  0.06      ],
[-0.06600098, -0.04399902, -0.02200195, ...,  0.10800049,
0.11      ,  0.10399903]], dtype=float32)
Coordinates:
* easting   (easting) float64 4.655e+05 4.656e+05 ... 4.827e+05 4.828e+05
* northing  (northing) float64 7.576e+06 7.576e+06 ... 7.595e+06 7.595e+06
height    (northing, easting) float64 500.0 500.0 500.0 ... 500.0 500.0```

And plot them:

```fig, (ax1, ax2) = plt.subplots(
nrows=1, ncols=2, sharey=True, figsize=(12, 8)
)

cbar_kwargs=dict(
)
kwargs = dict(center=0, cmap="seismic", cbar_kwargs=cbar_kwargs)

tmp = deriv_easting.plot(ax=ax1, **kwargs)
tmp = deriv_northing.plot(ax=ax2, **kwargs)

ax1.set_title("Easting derivative of the magnetic anomaly")
ax2.set_title("Northing derivative of the magnetic anomaly")
for ax in (ax1, ax2):
ax.set_aspect("equal")
ax.ticklabel_format(style="sci", scilimits=(0, 0))
plt.show()
```

By default, these two functions compute the horizontal derivatives using central finite differences methods. We can choose to use either the finite difference or the FFT-based method through the `method` argument.

For example, we can pass `method="fft"` to compute the derivatives in the frequency domain:

```deriv_easting = hm.derivative_easting(magnetic_grid_padded, method="fft")
deriv_easting
```
```<xarray.DataArray (northing: 370, easting: 346)>
array([[ 0.50048189, -0.1927814 ,  0.1357725 , ...,  0.39878718,
-0.59429492,  1.39003573],
[ 0.52115518, -0.20585515,  0.1431599 , ...,  0.40290014,
-0.601334  ,  1.40562022],
[ 0.52891843, -0.20861634,  0.14348668, ...,  0.40827736,
-0.61014921,  1.42089157],
...,
[ 2.43620524, -1.3841688 ,  0.49946329, ..., -0.3666973 ,
-0.73927341,  0.42750678],
[ 2.43711461, -1.36261579,  0.49770341, ..., -0.41819685,
-0.78503279,  0.39910907],
[ 2.40158947, -1.31498147,  0.5016298 , ..., -0.43062003,
-0.75786069,  0.32352992]])
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```
```deriv_northing = hm.derivative_northing(magnetic_grid_padded, method="fft")
deriv_northing
```
```<xarray.DataArray (northing: 370, easting: 346)>
array([[ 0.4998185 ,  0.51505215,  0.52086047, ..., -1.41297197,
-1.40373341, -1.39047911],
[-0.18344108, -0.19204698, -0.19683734, ...,  0.59467737,
0.59123353,  0.58320208],
[ 0.13266601,  0.13715393,  0.14021693, ..., -0.41085619,
-0.41006457, -0.40635707],
...,
[-0.67324956, -0.62931958, -0.58751428, ...,  0.04336819,
0.07440814,  0.09627284],
[ 1.03876533,  1.0129475 ,  0.97022432, ...,  0.19945125,
0.01275229, -0.12811101],
[-2.51808914, -2.39082358, -2.23713166, ..., -0.10561208,
0.28247267,  0.55809773]])
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```
```fig, (ax1, ax2) = plt.subplots(
nrows=1, ncols=2, sharey=True, figsize=(12, 8)
)

cbar_kwargs=dict(
)
kwargs = dict(center=0, cmap="seismic", cbar_kwargs=cbar_kwargs)

tmp = deriv_easting.plot(ax=ax1, **kwargs)
tmp = deriv_northing.plot(ax=ax2, **kwargs)

ax1.set_title("Easting derivative of the magnetic anomaly")
ax2.set_title("Northing derivative of the magnetic anomaly")
for ax in (ax1, ax2):
ax.set_aspect("equal")
ax.ticklabel_format(style="sci", scilimits=(0, 0))
plt.show()
```

Important

Horizontal derivatives through finite differences are usually more accurate and have less artifacts than their FFT-based counterpart.

## Upward continuation¶

We can also upward continue the original magnetic grid. This is, estimating the magnetic field generated by the same sources at a higher altitude. The original magnetic anomaly grid is located at 500 m above the ellipsoid, as we can see in its height coordinate. If we want to get the magnetic anomaly at 1000m above the ellipsoid, we need to upward continue it a height displacement of 500m:

```upward_continued = hm.upward_continuation(
)
```

This grid includes all the padding we added to the original magnetic grid, so we better unpad it using `xrft.unpad`:

```upward_continued = xrft.unpad(upward_continued, pad_width)
upward_continued
```
```<xarray.DataArray (northing: 370, easting: 346)>
array([[  1.5318864 ,   1.85100319,   2.13679603, ..., -33.60489015,
-31.65890937, -29.67750194],
[  1.82033329,   2.17484726,   2.49236703, ..., -35.96395426,
-33.83598794, -31.6667934 ],
[  2.07317532,   2.45927687,   2.80499305, ..., -38.27995984,
-35.97492132, -33.62307385],
...,
[ 50.44856723,  53.84378576,  57.13892691, ...,   4.05301916,
2.81272894,   1.764435  ],
[ 47.56514008,  50.69951642,  53.74614319, ...,   4.66844263,
3.44420588,   2.39520989],
[ 44.6368305 ,  47.50470955,  50.29752195, ...,   5.03756138,
3.86192233,   2.84250804]])
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```

And plot it:

```tmp = upward_continued.plot(cmap="seismic", center=0, add_colorbar=False)
plt.gca().set_aspect("equal")
plt.title("Upward continued magnetic anomaly to 1000m")
plt.gca().ticklabel_format(style="sci", scilimits=(0, 0))
plt.colorbar(tmp, label="nT")
plt.show()
```

## Reduction to the pole¶

We can also apply a reduction to the pole to any magnetic anomaly grid. This transformation consists in obtaining the magnetic anomaly of the same sources as if they were located on the North magnetic pole. We can apply it throught the `harmonica.reduction_to_pole` function.

Important

Applying reduction to the pole to low latitude regions can amplify high frequency noise.

The reduction to the pole needs information about the orientation of the geomagnetic field at the location of the survey and also the orientation of the magnetization vector of the sources.

The International Global Reference Field (IGRF) can provide us information about the inclination and declination of the geomagnetic field at the time of the survey (1990 in this case):

```inclination, declination = -52.98, 6.51
```

If we consider that the sources are magnetized in the same direction as the geomagnetic survey (hypothesis that is true in case the sources don’t have any remanenece), then we can apply the reduction to the pole passing only the `inclination` and `declination` of the geomagnetic field:

```rtp_grid = hm.reduction_to_pole(
)

# Unpad the reduced to the pole grid
rtp_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```

And plot it:

```tmp = rtp_grid.plot(cmap="seismic", center=0, add_colorbar=False)
plt.gca().set_aspect("equal")
plt.title("Magnetic anomaly reduced to the pole")
plt.gca().ticklabel_format(style="sci", scilimits=(0, 0))
plt.colorbar(tmp, label="nT")
plt.show()
```

If on the other hand we have any knowledge about the orientation of the magnetization vector of the sources, we can specify the `magnetization_inclination` and `magnetization_declination`:

```mag_inclination, mag_declination = -25, 21

tmp = rtp_grid = hm.reduction_to_pole(
inclination=inclination,
declination=declination,
magnetization_inclination=mag_inclination,
magnetization_declination=mag_declination,
)

# Unpad the reduced to the pole grid
rtp_grid
```
```<xarray.DataArray (northing: 370, easting: 346)>
array([[ -77.45213498, -114.40580183, -117.03946527, ...,  -87.74334513,
-75.65311786,   43.05422627],
[ -99.41761727, -137.48341563, -141.51463219, ...,  -31.30675001,
-20.40276016,   92.22734679],
[ -99.90958793, -139.38027341, -143.11774631, ...,  -38.52925024,
-29.73748202,   83.23546089],
...,
[  13.25373968, -170.27661708, -163.93666318, ...,  -44.43377684,
-80.78997965,  -66.80602774],
[ -28.4775795 , -187.01875283, -173.90318081, ...,  -27.82006314,
-55.31901234,  -61.09796854],
[-236.07965274, -301.62726523, -299.61036415, ...,  -28.01934017,
-18.38089699,  -22.65820682]])
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```
```tmp = rtp_grid.plot(cmap="seismic", center=0, add_colorbar=False)
plt.gca().set_aspect("equal")
plt.title("Reduced to the pole with remanence")
plt.gca().ticklabel_format(style="sci", scilimits=(0, 0))
plt.colorbar(tmp, label="nT")
plt.show()
```

## Gaussians filters¶

We can also apply Gaussians low-pass and high-pass filters to any regular grid. These two need us to select a cutoff wavelenght. The low-pass filter will remove any signal with a high spatial frequency, keeping only the signal components that have a wavelenght higher than the selected cutoff wavelength. The high-pass filter, on the other hand, removes any signal with a low spatial frequency, keeping only the components with a wavelength lower than the cutoff wavelenght. These two filters can be applied to our regular grid with the `harmonica.gaussian_lowpass` and `harmonica.gaussian_highpass`.

Let’s define a cutoff wavelength of 5 kilometers:

```cutoff_wavelength = 5e3
```

Then apply the two filters to our padded magnetic grid:

```magnetic_low_freqs = hm.gaussian_lowpass(
)
magnetic_high_freqs = hm.gaussian_highpass(
)
```

```magnetic_low_freqs = xrft.unpad(magnetic_low_freqs, pad_width)
```
```magnetic_low_freqs
```
```<xarray.DataArray (northing: 370, easting: 346)>
array([[ 4.60586414e+00,  4.72026510e+00,  4.81815067e+00, ...,
-3.68551892e+01, -3.50940481e+01, -3.33406525e+01],
[ 4.70588288e+00,  4.81992298e+00,  4.91669488e+00, ...,
-3.87398513e+01, -3.68891918e+01, -3.50465475e+01],
[ 4.78689904e+00,  4.89970003e+00,  4.99449202e+00, ...,
-4.06231159e+01, -3.86830241e+01, -3.67512054e+01],
...,
[ 5.24052381e+01,  5.51977770e+01,  5.80158127e+01, ...,
-8.46616115e-01, -1.00253896e+00, -1.15999951e+00],
[ 5.00637859e+01,  5.27307793e+01,  5.54221687e+01, ...,
1.79773740e-01, -1.44068907e-02, -2.10512248e-01],
[ 4.77049502e+01,  5.02456889e+01,  5.28097252e+01, ...,
1.08697396e+00,  8.60208050e-01,  6.31175660e-01]])
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```
```magnetic_high_freqs
```
```<xarray.DataArray (northing: 370, easting: 346)>
array([[ 30.39413586,  31.47968607,  31.8818005 , ..., -64.24478639,
-65.30597631, -66.25932308],
[ 31.79411712,  32.7800526 ,  33.08330512, ..., -63.46022197,
-64.61080822, -65.65352577],
[ 32.31307654,  33.30025114,  33.60548357, ..., -62.67693291,
-63.91695151, -65.14881903],
...,
[130.39468868, 117.20224742, 102.78411409, ...,   1.64654287,
-23.19753428, -40.64004932],
[132.03618967, 119.86919632, 105.97785568, ...,   5.82022626,
-21.48559311, -40.78948775],
[131.09497652, 120.15433554, 107.49020153, ...,  10.31305045,
-16.86020805, -36.43122449]])
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```

Let’s plot the results side by side:

```import verde as vd

fig, (ax1, ax2) = plt.subplots(
nrows=1, ncols=2, sharey=True, figsize=(12, 8)
)

maxabs = vd.maxabs(magnetic_low_freqs, magnetic_high_freqs)
kwargs = dict(cmap="seismic", vmin=-maxabs, vmax=maxabs, add_colorbar=False)

tmp = magnetic_low_freqs.plot(ax=ax1, **kwargs)
tmp = magnetic_high_freqs.plot(ax=ax2, **kwargs)

ax1.set_title("Magnetic anomaly after low-pass filter")
ax2.set_title("Magnetic anomaly after high-pass filter")
for ax in (ax1, ax2):
ax.set_aspect("equal")
ax.ticklabel_format(style="sci", scilimits=(0, 0))

plt.colorbar(
tmp,
ax=[ax1, ax2],
label="nT",
orientation="horizontal",
aspect=42,
shrink=0.8,