Calling grdfft with pygmt

Q&A PyGMT Q&A
1

I’m having trouble calling grdfft from pygmt (v0.4.0). Basically i want to loop through a directory (windows) and calculate the 2d FFT (radial average power spectrum) on a bunch of grid files and save the output to file for each grid. I use gdal to read the grids and then convert them to xarray. I use the virtualfile_from_grid to load the grid into GMT and then call grdfft with the arguments (args).

   for fn in os.listdir(griddir):
        if fn.endswith(".grd"):
            print("Doing fft on file:", fn)
            datagrid = gdal.Open(griddir + dsep + fn)
            band = datagrid.GetRasterBand(1)
            array = band.ReadAsArray()
            xgrid = xr.DataArray(array)

            with pygmt.clib.Session() as session:
                with session.virtualfile_from_grid(xgrid) as fin:
                    with pygmt.helpers.GMTTempFile() as fout:
                        args = "{} -Na -Er -V ->{}".format(fin, fout.name)
                        session.call_module('grdfft', args)
                        print(fout.read().strip())

This is the output and errror I get. Interestingly if i replace “grdfft” with “grdinfo” (with appropriate args) I get the expected output from grdinfo. So my basic pipeline seems ok, just some issue with grdfft???

grdfft [INFORMATION]: Processing input grid(s)
grdfft [INFORMATION]: Cartesian input grid
grdfft [INFORMATION]: Cartesian input grid
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "D:\xxxx\curiepointdepth.py", line 334, in runfft
    session.call_module('grdfft', args)
  File "C:\Users\xxxxx\Anaconda3\envs\cpd\lib\site-packages\pygmt\clib\session.py", line 501, in call_module
    status = c_call_module(
OSError: exception: access violation reading 0x0000000000000000
2

I think your life would be much simpler if you called grdfft directly on the grid file names and avoid all that gdal+xarrays stuff

3

ok thanks, after some experimentation I found I can in fact do as you suggest:

with pygmt.clib.Session() as session:
    fout="fft_result.txt"
    args = f"{fn} -Na -Er -G{fout}"
    session.call_module('grdfft', args)

where fn is my input grid file name.