fatiando
harmonica.reduction_to_pole flipping the output in north-south
Using harmonica.reduction_to_pole (v 0.7.0) with input TMI data read from an ERmapper ERS grid to an xarray. Data is in the north hemisphere in a UTM projection.Plotted input using a matplotlib helper - OK. Pass the xarray to harmonica.reduction_to_pole; plot the output (rtp) with same matplotlib helper). The RTP was computed but the data is flipped north to south. Both data arrays have similar spatial attributes (pasted below). Note that the CRS is dropped in the output.
AMF_38N_TestArea_Harmonica.zip
Hi @fearaschiarrai. Thanks for opening this Issue. I was able to reproduce the bug. I'm leaving here a smaller example that reproduces it, just to make it easier to test in the future. It downloads the Lightning Creek Sill Complex magnetic grid available in Ensaio, flips the northing coordinate and computes the RTP. We can see how the RTP flips the grid as you pointed out.
import numpy as np
import xrft
import xarray as xr
import matplotlib.pyplot as plt
import ensaio
import harmonica as hm
fname = ensaio.fetch_lightning_creek_magnetic(version=1)
tf = xr.load_dataset(fname).total_field_anomaly
# Flip the northing coordinate in the data array
tf = tf[::-1, :]
print("Total field anomaly grid after flipping the northing coordinate:")
print(tf)
# Pad the grid
pad_width = {
"easting": tf.easting.size // 3,
"northing": tf.northing.size // 3,
}
tf_padded = xrft.pad(tf.drop_vars("height"), pad_width)
# Plot padded grid
scale = 0.5 * np.max(np.abs(tf_padded))
tf_padded.plot(vmin=-scale, vmax=scale, cmap="RdBu_r")
plt.gca().set_aspect("equal")
plt.title("Padded TFA grid")
plt.show()
# Compute RTP
inclination, declination = -52.98, 6.51
rtp = hm.reduction_to_pole(tf_padded, inclination=inclination, declination=declination)
# Plot RTP grid
scale = 0.5 * np.max(np.abs(rtp))
rtp.plot(vmin=-scale, vmax=scale, cmap="RdBu_r")
plt.gca().set_aspect("equal")
plt.title("RTP grid")
plt.show()
Total field anomaly grid after flipping the northing coordinate:
<xarray.DataArray 'total_field_anomaly' (northing: 370, easting: 346)> Size: 1MB
array([[ 178.79995117, 170.39995117, 160.29995117, ..., 11.39995117,
-16.00004883, -35.80004883],
[ 182.09995117, 172.59995117, 161.39995117, ..., 5.99995117,
-21.50004883, -41.00004883],
[ 182.79995117, 172.39995117, 160.79995117, ..., 0.79995117,
-24.20004883, -41.80004883],
...,
[ 37.09995117, 38.19995117, 38.59995117, ..., -103.30004883,
-102.60004883, -101.90004883],
[ 36.49995117, 37.59995117, 37.99995117, ..., -102.20004883,
-101.50004883, -100.70004883],
[ 34.99995117, 36.19995117, 36.69995117, ..., -101.10004883,
-100.40004883, -99.60004883]], shape=(370, 346))
Coordinates:
* easting (easting) float64 3kB 4.655e+05 4.656e+05 ... 4.827e+05 4.828e+05
* northing (northing) float64 3kB 7.595e+06 7.595e+06 ... 7.576e+06 7.576e+06
height (northing, easting) float64 1MB 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...
license: Creative Commons Attribution 4.0 International Licence
references: Geophysical Acquisition & Processing Section 2019. MIM Dat...
long_name: total-field magnetic anomaly
units: nT
actual_range: [-1785. 3798.]
The source of the problem is that most of the FFT filtering functions in Harmonica assume that the coordinates are in a regular grid and that the coordinates are numerically sorted. I haven't dug deep into this to know exactly at which step the flip of the grid happens, but I think it might be related to the FFT and iFTT process.
Possible solutions
A solution that I can offer you right away so you can keep working with your data is to flip the northing coordinate of your grid, in a similar way as the previous example is doing:
da_amf = da_amf[::-1, :]
You can then pad this grid and use it to compute the RTP. That should get your RTP grid in the right orientation.
On our side, I think we should either raise an error when the coordinates in the grid are not valid (non regular grid, non sorted coordinates, etc). I'll try to open a PR for this soon.
