gabrahamastro commented 2 weeks ago

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 commented 2 weeks ago

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

gabrahamastro commented 2 weeks ago

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 commented 1 week ago

@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 commented 1 week ago

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 commented 1 week ago

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 commented 1 week ago

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 commented 1 week ago

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 commented 1 week ago

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!

fatiando / harmonica

Difficulty in performing Reduction to the Pole #534

Output of `print(grid)`

Code used to add padding (tested with and without)

Error message produced from `hm.reduction_to_pole`

Final Remarks

Code

Converting `grid` to an `xarray.DataArray`, `grid_DataArray`.

Adding padding

Reducing to the pole

Final Remarks

fatiando / harmonica

Difficulty in performing Reduction to the Pole #534

Output of print(grid)

Code used to add padding (tested with and without)

Error message produced from hm.reduction_to_pole

Final Remarks

Code

Converting grid to an xarray.DataArray, grid_DataArray.

Adding padding

Reducing to the pole

Final Remarks

Output of `print(grid)`

Error message produced from `hm.reduction_to_pole`

Converting `grid` to an `xarray.DataArray`, `grid_DataArray`.