Skip to content
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.

Sign up for GitHub

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

Jump to bottom

Difficulty in performing Reduction to the Pole #534

Closed
gabrahamastro opened this issue Oct 21, 2024 · 8 comments
Closed

Difficulty in performing Reduction to the Pole #534

gabrahamastro opened this issue Oct 21, 2024 · 8 comments

Comments

@gabrahamastro
Copy link

I'm trying to perform a reduction to the pole operation on some gridded magnetic field data. However, I can't work out how to remove the error of Found nan(s) on the passed grid. The grid must not have missing values before computing the Fast Fourier Transform.

To get the grid in the first place, I used the below code (using this, this and this tutorial):

data = pd.read_csv('aeromag.csv')

eql = hm.EquivalentSources(damping=0.01,
                           depth=5e3,
                          ) #worked out via a cross-validation method to estimate the parameters
eql.fit(coordinates, Mag_Anomaly) #coordinates use easting, northing (cartesian), elevation, Mag_Anomaly is the associated magnetic field anomaly

#Then I created the grid

grid_coords = vd.grid_coordinates(region=region, 
                                  spacing=4e3, 
                                  extra_coords = 300
                                  ) #region is west, east, south, north in cartesian

grid = eql.grid(grid_coords, 
                data_names=["mag_field_anomaly"],
                dim=("northings","eastings"),
                #projection=projection
                ) #projection is pyproj.Proj(init="epsg:27700")

grid = vd.convexhull_mask((DataZoom.easting,DataZoom.northing), grid=grid) #remove data from outside the area

#Following this, I added some padding to prepare the grid for RTP

pad_width = {
    "easting": grid.easting.size // 3,
    "northing": grid.northing.size // 3,
}
magnetic_grid_no_height = grid.drop_vars("upward")
magnetic_grid_padded = xrft.pad(magnetic_grid_no_height, pad_width)

#Lastly, I attempted to perform the RTP

rtp_grid = hm.reduction_to_pole(
    magnetic_grid_padded, inclination=70.252, declination=-12.08
)

However, this is where the error occurred.

I'm currently using:
Verde v1.8.1
Harmonica v0.7.0
xarray 2024.9.0

I realise that there have been some comments about Equivalent Sources (#426) not being suitable in its current form for magnetic fields, so please advise about another method of gridding the magnetic field data if this would be the best way forward.
@leouieda
Copy link
Member

@gabrahamastro could you please post the error message that you got? Without it, it's hard to know what went wrong.

@gabrahamastro
Copy link
Author

Apologies @leouieda - below is the error message:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[40], [line 4](vscode-notebook-cell:?execution_count=40&line=4)
      [1](vscode-notebook-cell:?execution_count=40&line=1) grid.dropna('northing', how='all')
      [2](vscode-notebook-cell:?execution_count=40&line=2) grid.dropna('easting', how='all')
----> [4](vscode-notebook-cell:?execution_count=40&line=4) rtp_grid = hm.reduction_to_pole(
      [5](vscode-notebook-cell:?execution_count=40&line=5)     grid, inclination=inclination, declination=declination
      [6](vscode-notebook-cell:?execution_count=40&line=6) )
      [8](vscode-notebook-cell:?execution_count=40&line=8) # Unpad the reduced to the pole grid
      [9](vscode-notebook-cell:?execution_count=40&line=9) rtp_grid = xrft.unpad(rtp_grid, pad_width)

File /usr/local/Caskroom/miniforge/base/envs/hm/lib/python3.12/site-packages/harmonica/_transformations.py:334, in reduction_to_pole(grid, inclination, declination, magnetization_inclination, magnetization_declination)
    [285](https://file+.vscode-resource.vscode-cdn.net/usr/local/Caskroom/miniforge/base/envs/hm/lib/python3.12/site-packages/harmonica/_transformations.py:285) def reduction_to_pole(
    [286](https://file+.vscode-resource.vscode-cdn.net/usr/local/Caskroom/miniforge/base/envs/hm/lib/python3.12/site-packages/harmonica/_transformations.py:286)     grid,
    [287](https://file+.vscode-resource.vscode-cdn.net/usr/local/Caskroom/miniforge/base/envs/hm/lib/python3.12/site-packages/harmonica/_transformations.py:287)     inclination,
   (...)
    [290](https://file+.vscode-resource.vscode-cdn.net/usr/local/Caskroom/miniforge/base/envs/hm/lib/python3.12/site-packages/harmonica/_transformations.py:290)     magnetization_declination=None,
    [291](https://file+.vscode-resource.vscode-cdn.net/usr/local/Caskroom/miniforge/base/envs/hm/lib/python3.12/site-packages/harmonica/_transformations.py:291) ):
    [292](https://file+.vscode-resource.vscode-cdn.net/usr/local/Caskroom/miniforge/base/envs/hm/lib/python3.12/site-packages/harmonica/_transformations.py:292)     """
    [293](https://file+.vscode-resource.vscode-cdn.net/usr/local/Caskroom/miniforge/base/envs/hm/lib/python3.12/site-packages/harmonica/_transformations.py:293)     Calculate the reduction to the pole of a magnetic field grid
    [294](https://file+.vscode-resource.vscode-cdn.net/usr/local/Caskroom/miniforge/base/envs/hm/lib/python3.12/site-packages/harmonica/_transformations.py:294) 
   (...)
    [332](https://file+.vscode-resource.vscode-cdn.net/usr/local/Caskroom/miniforge/base/envs/hm/lib/python3.12/site-packages/harmonica/_transformations.py:332)     harmonica.filters.reduction_to_pole_kernel
    [333](https://file+.vscode-resource.vscode-cdn.net/usr/local/Caskroom/miniforge/base/envs/hm/lib/python3.12/site-packages/harmonica/_transformations.py:333)     """
--> [334](https://file+.vscode-resource.vscode-cdn.net/usr/local/Caskroom/miniforge/base/envs/hm/lib/python3.12/site-packages/harmonica/_transformations.py:334)     return apply_filter(
    [335](https://file+.vscode-resource.vscode-cdn.net/usr/local/Caskroom/miniforge/base/envs/hm/lib/python3.12/site-packages/harmonica/_transformations.py:335)         grid,
    [336](https://file+.vscode-resource.vscode-cdn.net/usr/local/Caskroom/miniforge/base/envs/hm/lib/python3.12/site-packages/harmonica/_transformations.py:336)         reduction_to_pole_kernel,
    [337](https://file+.vscode-resource.vscode-cdn.net/usr/local/Caskroom/miniforge/base/envs/hm/lib/python3.12/site-packages/harmonica/_transformations.py:337)         inclination=inclination,
    [338](https://file+.vscode-resource.vscode-cdn.net/usr/local/Caskroom/miniforge/base/envs/hm/lib/python3.12/site-packages/harmonica/_transformations.py:338)         declination=declination,
    [339](https://file+.vscode-resource.vscode-cdn.net/usr/local/Caskroom/miniforge/base/envs/hm/lib/python3.12/site-packages/harmonica/_transformations.py:339)         magnetization_inclination=magnetization_inclination,
    [340](https://file+.vscode-resource.vscode-cdn.net/usr/local/Caskroom/miniforge/base/envs/hm/lib/python3.12/site-packages/harmonica/_transformations.py:340)         magnetization_declination=magnetization_declination,
    [341](https://file+.vscode-resource.vscode-cdn.net/usr/local/Caskroom/miniforge/base/envs/hm/lib/python3.12/site-packages/harmonica/_transformations.py:341)     )

File /usr/local/Caskroom/miniforge/base/envs/hm/lib/python3.12/site-packages/harmonica/filters/_utils.py:42, in apply_filter(grid, fft_filter, **kwargs)
     [16](https://file+.vscode-resource.vscode-cdn.net/usr/local/Caskroom/miniforge/base/envs/hm/lib/python3.12/site-packages/harmonica/filters/_utils.py:16) """
     [17](https://file+.vscode-resource.vscode-cdn.net/usr/local/Caskroom/miniforge/base/envs/hm/lib/python3.12/site-packages/harmonica/filters/_utils.py:17) Apply a filter to a grid and return the transformed grid in spatial domain
     [18](https://file+.vscode-resource.vscode-cdn.net/usr/local/Caskroom/miniforge/base/envs/hm/lib/python3.12/site-packages/harmonica/filters/_utils.py:18) 
   (...)
     [39](https://file+.vscode-resource.vscode-cdn.net/usr/local/Caskroom/miniforge/base/envs/hm/lib/python3.12/site-packages/harmonica/filters/_utils.py:39)     ``grid``. Defined are in the spatial domain.
     [40](https://file+.vscode-resource.vscode-cdn.net/usr/local/Caskroom/miniforge/base/envs/hm/lib/python3.12/site-packages/harmonica/filters/_utils.py:40) """
     [41](https://file+.vscode-resource.vscode-cdn.net/usr/local/Caskroom/miniforge/base/envs/hm/lib/python3.12/site-packages/harmonica/filters/_utils.py:41) # Run sanity checks on the grid
---> [42](https://file+.vscode-resource.vscode-cdn.net/usr/local/Caskroom/miniforge/base/envs/hm/lib/python3.12/site-packages/harmonica/filters/_utils.py:42) grid_sanity_checks(grid)
     [43](https://file+.vscode-resource.vscode-cdn.net/usr/local/Caskroom/miniforge/base/envs/hm/lib/python3.12/site-packages/harmonica/filters/_utils.py:43) # Catch the dims of the grid
     [44](https://file+.vscode-resource.vscode-cdn.net/usr/local/Caskroom/miniforge/base/envs/hm/lib/python3.12/site-packages/harmonica/filters/_utils.py:44) dims = grid.dims

File /usr/local/Caskroom/miniforge/base/envs/hm/lib/python3.12/site-packages/harmonica/filters/_utils.py:87, in grid_sanity_checks(grid)
     [85](https://file+.vscode-resource.vscode-cdn.net/usr/local/Caskroom/miniforge/base/envs/hm/lib/python3.12/site-packages/harmonica/filters/_utils.py:85) # Check if the grid has nans
     [86](https://file+.vscode-resource.vscode-cdn.net/usr/local/Caskroom/miniforge/base/envs/hm/lib/python3.12/site-packages/harmonica/filters/_utils.py:86) if np.isnan(grid).any():
---> [87](https://file+.vscode-resource.vscode-cdn.net/usr/local/Caskroom/miniforge/base/envs/hm/lib/python3.12/site-packages/harmonica/filters/_utils.py:87)     raise ValueError(
     [88](https://file+.vscode-resource.vscode-cdn.net/usr/local/Caskroom/miniforge/base/envs/hm/lib/python3.12/site-packages/harmonica/filters/_utils.py:88)         "Found nan(s) on the passed grid. "
     [89](https://file+.vscode-resource.vscode-cdn.net/usr/local/Caskroom/miniforge/base/envs/hm/lib/python3.12/site-packages/harmonica/filters/_utils.py:89)         + "The grid must not have missing values before computing the "
     [90](https://file+.vscode-resource.vscode-cdn.net/usr/local/Caskroom/miniforge/base/envs/hm/lib/python3.12/site-packages/harmonica/filters/_utils.py:90)         + "Fast Fourier Transform."
     [91](https://file+.vscode-resource.vscode-cdn.net/usr/local/Caskroom/miniforge/base/envs/hm/lib/python3.12/site-packages/harmonica/filters/_utils.py:91)     )

ValueError: Found nan(s) on the passed grid. The grid must not have missing values before computing the Fast Fourier Transform.

@leouieda
Copy link
Member

@gabrahamastro thank you!

OK, the problem is that you're doing the convexhull_mask and then trying to reduce to the pole. For these FFT-based filters, you can't have missing values in the grids, which means no masks can be used. If you remove the convex hull masking step, it should work fine. If you really want to mask, I'd suggest doing it after the reduction to the pole or whatever other filter you may want to apply.

@leouieda
Copy link
Member

Also, note that in your reduction to the pole you don't specify the source magnetization direction. An induced magnetization is implied and your reduction will only be correct if this is true. If there is any source with remanent magnetization in your data then the reduction will be incorrect. Just thought it was worth point out.

@leouieda
Copy link
Member

Let me know if removing the mask works of it comes up again. I'll close this issue for now but if you're still getting an error, please feel free to reopen it.

@gabrahamastro
Copy link
Author

Hi @leouieda, thank you for the help! I've now removed convexhull_mask but have not been able to perform RTP yet, receiving the same error.

I've looked into it and, using the code which raised the error, I ran the below code:

if np.isnan(grid_coords).any():
    print("Error!!")
else:
    print("All good")

if np.isnan(grid).any():
    print("Error - something!!")
else:
    print("All good")

if np.isnan(grid.magnetic_anomaly).any():
    print("Error - magnetic_anomaly")
else:
    print("All good")

if np.isnan(grid.northing).any():
    print("Error - northings")
else:
    print("All good")

if np.isnan(grid.easting).any():
    print("Error - eastings")
else:
    print("All good")

if np.isnan(grid.upward).any():
    print("Error - upward")
else:
    print("All good")

running through all listed attributes in grid generated using eql.grid and in the object inputted into eql.grid, which I labelled grid_coords.

To these statements, I received All good except for the statement if np.isnan(grid).any() being the specific if statement used by the RTP operation performed by harmonica which raised the error.

This has confused me, since, if there are no nan values in each variable in grid, surely there are no nan values in the whole grid.

To help, below I'll write down the output of print(grid) followed by the error message produced when running hm.reduction_to_pole using grid as the input (I tested this both as is and with padding with xrft.pad, and both produce the same error - I'll also write out the code used here below).

Output of print(grid)

<xarray.Dataset> Size: 67kB
Dimensions:           (northing: 58, easting: 71)
Coordinates:
  * easting           (easting) float64 568B 1.462e+05 1.467e+05 ... 1.812e+05
  * northing          (northing) float64 464B 7.177e+05 7.182e+05 ... 7.464e+05
    upward            (northing, easting) float64 33kB 1.5e+03 ... 1.5e+03
Data variables:
    magnetic_anomaly  (northing, easting) float64 33kB -728.5 -673.9 ... 320.1
Attributes:
    metadata:  Generated by EquivalentSources(block_size=500, damping=0.01, d...

Code used to add padding (tested with and without)

pad_width = {
    "easting": grid.easting.size // 3,
    "northing": grid.northing.size // 3,
}
magnetic_grid_no_height = grid.drop_vars("upward")
magnetic_grid_padded = xrft.pad(magnetic_grid_no_height, pad_width)

Error message produced from hm.reduction_to_pole

----> [1](vscode-notebook-cell:?execution_count=60&line=1) rtp_grid = hm.reduction_to_pole(
      [2](vscode-notebook-cell:?execution_count=60&line=2)     grid, inclination=inclination, declination=declination
      [3](vscode-notebook-cell:?execution_count=60&line=3) )
      [5](vscode-notebook-cell:?execution_count=60&line=5) # Unpad the reduced to the pole grid
      [6](vscode-notebook-cell:?execution_count=60&line=6) rtp_grid = xrft.unpad(rtp_grid, pad_width)

File /usr/local/Caskroom/miniforge/base/envs/hm/lib/python3.12/site-packages/harmonica/_transformations.py:334, in reduction_to_pole(grid, inclination, declination, magnetization_inclination, magnetization_declination)
    [285](https://file+.vscode-resource.vscode-cdn.net/usr/local/Caskroom/miniforge/base/envs/hm/lib/python3.12/site-packages/harmonica/_transformations.py:285) def reduction_to_pole(
    [286](https://file+.vscode-resource.vscode-cdn.net/usr/local/Caskroom/miniforge/base/envs/hm/lib/python3.12/site-packages/harmonica/_transformations.py:286)     grid,
    [287](https://file+.vscode-resource.vscode-cdn.net/usr/local/Caskroom/miniforge/base/envs/hm/lib/python3.12/site-packages/harmonica/_transformations.py:287)     inclination,
   (...)
    [290](https://file+.vscode-resource.vscode-cdn.net/usr/local/Caskroom/miniforge/base/envs/hm/lib/python3.12/site-packages/harmonica/_transformations.py:290)     magnetization_declination=None,
    [291](https://file+.vscode-resource.vscode-cdn.net/usr/local/Caskroom/miniforge/base/envs/hm/lib/python3.12/site-packages/harmonica/_transformations.py:291) ):
    [292](https://file+.vscode-resource.vscode-cdn.net/usr/local/Caskroom/miniforge/base/envs/hm/lib/python3.12/site-packages/harmonica/_transformations.py:292)     """
    [293](https://file+.vscode-resource.vscode-cdn.net/usr/local/Caskroom/miniforge/base/envs/hm/lib/python3.12/site-packages/harmonica/_transformations.py:293)     Calculate the reduction to the pole of a magnetic field grid
    [294](https://file+.vscode-resource.vscode-cdn.net/usr/local/Caskroom/miniforge/base/envs/hm/lib/python3.12/site-packages/harmonica/_transformations.py:294) 
   (...)
    [332](https://file+.vscode-resource.vscode-cdn.net/usr/local/Caskroom/miniforge/base/envs/hm/lib/python3.12/site-packages/harmonica/_transformations.py:332)     harmonica.filters.reduction_to_pole_kernel
    [333](https://file+.vscode-resource.vscode-cdn.net/usr/local/Caskroom/miniforge/base/envs/hm/lib/python3.12/site-packages/harmonica/_transformations.py:333)     """
--> [334](https://file+.vscode-resource.vscode-cdn.net/usr/local/Caskroom/miniforge/base/envs/hm/lib/python3.12/site-packages/harmonica/_transformations.py:334)     return apply_filter(
    [335](https://file+.vscode-resource.vscode-cdn.net/usr/local/Caskroom/miniforge/base/envs/hm/lib/python3.12/site-packages/harmonica/_transformations.py:335)         grid,
    [336](https://file+.vscode-resource.vscode-cdn.net/usr/local/Caskroom/miniforge/base/envs/hm/lib/python3.12/site-packages/harmonica/_transformations.py:336)         reduction_to_pole_kernel,
    [337](https://file+.vscode-resource.vscode-cdn.net/usr/local/Caskroom/miniforge/base/envs/hm/lib/python3.12/site-packages/harmonica/_transformations.py:337)         inclination=inclination,
    [338](https://file+.vscode-resource.vscode-cdn.net/usr/local/Caskroom/miniforge/base/envs/hm/lib/python3.12/site-packages/harmonica/_transformations.py:338)         declination=declination,
    [339](https://file+.vscode-resource.vscode-cdn.net/usr/local/Caskroom/miniforge/base/envs/hm/lib/python3.12/site-packages/harmonica/_transformations.py:339)         magnetization_inclination=magnetization_inclination,
    [340](https://file+.vscode-resource.vscode-cdn.net/usr/local/Caskroom/miniforge/base/envs/hm/lib/python3.12/site-packages/harmonica/_transformations.py:340)         magnetization_declination=magnetization_declination,
    [341](https://file+.vscode-resource.vscode-cdn.net/usr/local/Caskroom/miniforge/base/envs/hm/lib/python3.12/site-packages/harmonica/_transformations.py:341)     )

File /usr/local/Caskroom/miniforge/base/envs/hm/lib/python3.12/site-packages/harmonica/filters/_utils.py:42, in apply_filter(grid, fft_filter, **kwargs)
     [16](https://file+.vscode-resource.vscode-cdn.net/usr/local/Caskroom/miniforge/base/envs/hm/lib/python3.12/site-packages/harmonica/filters/_utils.py:16) """
     [17](https://file+.vscode-resource.vscode-cdn.net/usr/local/Caskroom/miniforge/base/envs/hm/lib/python3.12/site-packages/harmonica/filters/_utils.py:17) Apply a filter to a grid and return the transformed grid in spatial domain
     [18](https://file+.vscode-resource.vscode-cdn.net/usr/local/Caskroom/miniforge/base/envs/hm/lib/python3.12/site-packages/harmonica/filters/_utils.py:18) 
   (...)
     [39](https://file+.vscode-resource.vscode-cdn.net/usr/local/Caskroom/miniforge/base/envs/hm/lib/python3.12/site-packages/harmonica/filters/_utils.py:39)     ``grid``. Defined are in the spatial domain.
     [40](https://file+.vscode-resource.vscode-cdn.net/usr/local/Caskroom/miniforge/base/envs/hm/lib/python3.12/site-packages/harmonica/filters/_utils.py:40) """
     [41](https://file+.vscode-resource.vscode-cdn.net/usr/local/Caskroom/miniforge/base/envs/hm/lib/python3.12/site-packages/harmonica/filters/_utils.py:41) # Run sanity checks on the grid
---> [42](https://file+.vscode-resource.vscode-cdn.net/usr/local/Caskroom/miniforge/base/envs/hm/lib/python3.12/site-packages/harmonica/filters/_utils.py:42) grid_sanity_checks(grid)
     [43](https://file+.vscode-resource.vscode-cdn.net/usr/local/Caskroom/miniforge/base/envs/hm/lib/python3.12/site-packages/harmonica/filters/_utils.py:43) # Catch the dims of the grid
     [44](https://file+.vscode-resource.vscode-cdn.net/usr/local/Caskroom/miniforge/base/envs/hm/lib/python3.12/site-packages/harmonica/filters/_utils.py:44) dims = grid.dims

File /usr/local/Caskroom/miniforge/base/envs/hm/lib/python3.12/site-packages/harmonica/filters/_utils.py:87, in grid_sanity_checks(grid)
     [85](https://file+.vscode-resource.vscode-cdn.net/usr/local/Caskroom/miniforge/base/envs/hm/lib/python3.12/site-packages/harmonica/filters/_utils.py:85) # Check if the grid has nans
     [86](https://file+.vscode-resource.vscode-cdn.net/usr/local/Caskroom/miniforge/base/envs/hm/lib/python3.12/site-packages/harmonica/filters/_utils.py:86) if np.isnan(grid).any():
---> [87](https://file+.vscode-resource.vscode-cdn.net/usr/local/Caskroom/miniforge/base/envs/hm/lib/python3.12/site-packages/harmonica/filters/_utils.py:87)     raise ValueError(
     [88](https://file+.vscode-resource.vscode-cdn.net/usr/local/Caskroom/miniforge/base/envs/hm/lib/python3.12/site-packages/harmonica/filters/_utils.py:88)         "Found nan(s) on the passed grid. "
     [89](https://file+.vscode-resource.vscode-cdn.net/usr/local/Caskroom/miniforge/base/envs/hm/lib/python3.12/site-packages/harmonica/filters/_utils.py:89)         + "The grid must not have missing values before computing the "
     [90](https://file+.vscode-resource.vscode-cdn.net/usr/local/Caskroom/miniforge/base/envs/hm/lib/python3.12/site-packages/harmonica/filters/_utils.py:90)         + "Fast Fourier Transform."
     [91](https://file+.vscode-resource.vscode-cdn.net/usr/local/Caskroom/miniforge/base/envs/hm/lib/python3.12/site-packages/harmonica/filters/_utils.py:91)     )

ValueError: Found nan(s) on the passed grid. The grid must not have missing values before computing the Fast Fourier Transform.

Final Remarks

I'm assuming there is a problem with the way I'm applying the equivalent layer, though I can't work out what it would be.

Let me know whether there's anything you'd like to know, and thank you again for your help!

@gabrahamastro
Copy link
Author

Added to this, I've run the following to convert it into an xarray.DataArray (as expected by hm.reduction_to_pole), before checking whether nan values exist, to which I got the response All good. This suggests to me that hm.reduction_to_pole doesn't work with xarray.Dataset objects (as produced by (hm.EquivalentSources()).grid()).

grid_to_array = grid.to_array()

if np.isnan(grid_to_array).any():
    print("Error")
else:
    print("All good")

The only problem is that the xarray.DataArray produced now has 3 dimensions, outputting the following when I run print(grid_to_array)

<xarray.DataArray (variable: 1, northing: 58, easting: 71)> Size: 33kB
array([[[-728.51206453, -673.88127469, -617.15209162, ...,
          -34.35146787,  -46.14826617,  -56.8376889 ],
        [-703.49636856, -640.17927609, -574.81974622, ...,
          -29.14635314,  -40.82317425,  -51.47099936],
        [-667.06020527, -596.51577751, -524.12277989, ...,
          -22.02620696,  -33.22700417,  -43.58516298],
        ...,
        [ 449.14425835,  396.68635615,  338.42827755, ...,
          308.58023459,  300.59462775,  290.40129421],
        [ 495.76486742,  446.99807664,  391.19060032, ...,
          325.38031269,  317.18262695,  306.59755843],
        [ 534.07701016,  489.12683834,  436.20881077, ...,
          339.55007815,  330.99121324,  320.06630346]]])
Coordinates:
  * easting   (easting) float64 568B 1.462e+05 1.467e+05 ... 1.807e+05 1.812e+05
  * northing  (northing) float64 464B 7.177e+05 7.182e+05 ... 7.464e+05
    upward    (northing, easting) float64 33kB 1.5e+03 1.5e+03 ... 1.5e+03
  * variable  (variable) object 8B 'magnetic_anomaly'
Attributes:
    metadata:  Generated by EquivalentSources(block_size=500, damping=0.01, d...

This means that it doesn't work with hm.reduction_to_pole, since it expects only 2 dimensions.

Do you think my theory about what is happening here is correct, and is there a better way to format an xarray.DataArray which would be compatible with hm.reduction_to_pole?

@gabrahamastro
Copy link
Author

Apologies for the frequent comments, but I seem to have made a breakthrough. I've managed to get an output for the command, having converted the xarray.Dataset to an xarray.DataArray before adding some padding (note that I adapted the code from this helpful tutorial).

Code

Converting grid to an xarray.DataArray, grid_DataArray.

grid_DataArray = xr.DataArray(
    data=grid['magnetic_anomaly'].values,
    coords={
        'northing': grid['northing'],
        'easting': grid['easting'],
        'upward': grid['upward']  # Assuming 'upward' is part of the dataset
    },
    dims=['northing', 'easting']
)

Adding padding

pad_width = {
    "easting": grid_DataArray.easting.size // 3,
    "northing": grid_DataArray.northing.size // 3,
}
magnetic_grid_no_height = grid_DataArray.drop_vars("upward")
magnetic_grid_padded = xrft.pad(magnetic_grid_no_height, pad_width)

Reducing to the pole

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)

Final Remarks

Having done this, I feel that this issue is now sorted, though I feel this may highlight a bug in the cross-compatibility of hm.reduction_to_pole and (hm.EquivalentSources()).grid().

Again, thank you @leouieda for your help and I'll make sure to post again if anything else arrises!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Assignees
No one assigned
Labels
None yet
Projects
None yet
Milestone
No milestone
Development

No branches or pull requests
2 participants
@leouieda @gabrahamastro