I have a real 2D Array (Nx,Ny) of some variable, let’s say absolute vorticity. I perform a 2d fft on the matrix by:

fc_2d = FFTW.r2r(arr,FFTW.R2HC)

which returns a 2d matrix of the same size of Fourier coefficients.

I want to do a 2d low pass filter by zero-ing out the scales greater than wavenumber 3 for example. I read the documentation and I’m still a little confused about the format of fc_2d. Are both dimensions in half-complex format? If so, I think I would just zero out the Fourier coefficients corresponding to wavenumbers 4 and up in both dimensions?

ind = 4
# x filter
fc_2d[ind+1:Nx-ind+1,:] .= 0
# y filter
fc_2d[:,ind+1:Nx-ind+1] .= 0

Sorry if this is really basic / does not make sense. I am new to Julia and spectral analysis. Thanks.

The R2HC format in 2d is pretty tricky to comprehend, I would recommend starting with a simple complex fft and then moving to the r2c format once you understand that in order to take advantage of the real input. (The R2HC transforms are typically slower than r2c, anyway.)

That seems a job for GMT.jl grdfft but unfortunately we don’t yet have a manual tailored for the Julia wrapper (though the grdfft is wrapped). So far, mos complete manual is still the CLI version

It seems like a gap in the API that there is not rfftfreq method for the R2HC transform. (Also that there is no convenient rfftfreq for 2D transforms.)

That kind of Manhattan metric was what @balancedflow seemed to be after, if I understood correctly their initial post. It actually be useful, for example for dealiasing in pseudo-spectral methods.

For truncating modes |\boldsymbol{k}| > k_{\text{max}}, one could do something like:

kmax² = kmax^2
for j ∈ axes(us, 2), i ∈ axes(us, 1)
k² = kx[i]^2 + ky[j]^2
if k² > kmax²
us[i, j] = 0
end
end

Or in a more generic way (also valid for dimensions other than 2):

kmax² = kmax^2
ks = (kx, ky)
for I ∈ CartesianIndices(us)
k⃗ = map(getindex, ks, Tuple(I))
k² = sum(abs2, k⃗)
if k² > kmax²
us[I] = 0
end
end

Doing that is applying a boxcar filter in the frequency domain. A boxcar filter in frequency corresponds to a sinc filtering in space domain and that will cause ringing. The way to go is to apply a cosine-taper band pass, and that is option F (or filter in GMT.jl) of grdfft does. See grdfft — GMT 6.6.0 documentation

using GMT
# Download a Geoid grid
G = grdcut("@earth_geoid_02m", region=(-40,-20,32,45));
#Filter with Boxcar (low pass wavelength > 500 km). ´f=:g´ is to infor that grid is in geogs, `N=:a` is to select best FFT dims(
Gf_box = grdfft(G, F="-/-/500000/500000", N=:a, f=:g);
# Filter with cosne taper. Pass wavelengths > 500 km and rejecting wavelengths < 100 km
Gf_ct = grdfft(G, F="-/-/500000/100000", N=:a, f=:g);

Thanks for the example code. This really helps. I was going through the documentation and trying grdfft on my data, but am having a couple of issues.

# get data
arr = nc["avo"][:,:,1] # 354x354 Matrix{Float32}
arr = mat2grid(arr); # Is it necessary to convert from a matrix to a GMT.GMTgrid here?
# do low-pass
# pass 100 km and cut out 90 km
filtered = grdfft(arr, F="-/-/100000/90000",f=:g,N=:a);

This code should pass all scales > 100 km and completely cut off scales < 90 km, and the scales between these cutoffs are slowly filtered out, correct? Though the output of this is not what I expect. When I do something like filtered = grdfft(arr, F="-/-/10000000/9000000",f=:g,N=:a); I get something that looks more physically correct. These correspond to scales 10,000 km and 9,000 km (off by a factor of 2 of what my data resolution is actually in). My data resolution is 1 km so the grid domain is 354 km x 354 km.

Am I missing something obvious here? (I guess the function doesn’t know what the resolution of my data is here, since its just a 2d array of scalars)…

Absolutely, and what you did was not enough. You need also to pass at least the referencing information (either in the x,y vectors or the the hdr vector/matrix). The mat2grid help explain those things. And that is need because the grdfft manual also explains that.

Description

grdfft will take the 2-D forward Fast Fourier Transform and perform one or more mathematical operations in the frequency domain before transforming back to the space domain. An option is provided to scale the data before writing the new values to an output file. The horizontal dimensions of the grid are assumed to be in meters. Geographical grids may be used by specifying the -fflags option that scales degrees to meters.

So, your data must be in meters or in geographical coordinates. It’s in later case that one passes the f=:g option to inform the program that it has convert from degrees to meters. But since your original post mentioned absolute vorticity I assumed that it was an oceanographic or atmospheric data. Another parameter to look at is the memory layout (row-major, column-major, top-down or botom-up). I recomend that you visualize the output of the mat2grid function to see if it makes sense.

EDIT

This makes me think that you are reading from a netCDF grid. You may save yourself the mat2grid hassle if you read the nc grid directly with gmtread or gdalread.