Extrapolate a 2D array using Fourier Transform

I need to extrapolate a given 2D array to a larger domain, keeping the spatial frequency. This is the original field: (the data file in numpy npz format and a Jupyter notebook to plot it can be found here)

The horizontal size here is 14.864408108 (critical wavelength), and the vertical size is 14.864408108/sqrt(3)

I would like to extrapolate it to a square domain of size, let's say, ten times bigger than the original horizontal size, but fixing the number of grid points to a given number (256x256 for example).

Currently I use numpy.tile to do that:

import numpy as np
from scipy import interpolate

def enlarge(field,mult_x=5,mult_y=8):
y = np.linspace(2*mult_y*Y.min(),2*mult_y*Y.max(),int(np.round(2*mult_y*field.shape/mult_x)))
field2 = np.concatenate((field, field[::-1, :]))
field4 = np.concatenate((field2, field2[:,::-1]), axis=1)
large = np.tile(field4,(mult_y, mult_x)) # 17 just because it gives a number close to 256*10
less = large[::mult_x, ::mult_x] # skip every 10 rows/columns
#print(less.shape)
f = interpolate.interp1d(y, less,axis=0)
y_new = np.linspace(2*mult_y*Y.min(),2*mult_y*Y.max(),less.shape)
field = f(y_new)
return field

And I get this figure: Can I use Fourier transform to do that?

You can use the Fourier transform for a lot of things – it can be part of a particular step in a scaling algorithm. Still, what you're looking for are scaling algorithms.

In this particular case, your original data size is (74 × 128) – what that means in wavelengths or any physical unit doesn't matter to algorithms – and you want to scale it to (512 × 512).

So what you'd need to get 128 -> 512 would be simply inserting 3 zero rows after every image row and then low-pass filtering the result. A 128-FFT, followed by a zero-padding the 512-128 higher bins and then doing a 512-IFFT would do that – but not like you'd want it to, but with circular effects at the border; what you'd be doing mathematically would be a cyclic convolution with a 1/4 band-wide filter.

Same goes for the 74 -> 512 dimension – but here you'd need to interpolate higher first (to 37·512, to be exact, i.e. by a factor of 37·256) just to decimate down back to 512 again. And then you'd still be stuck with the circular convolution.

What you in fact want to do is apply resampling by a rational rate – the core filter for that could be applied by fast convolution, i.e. with the help of overlapping multiple FFTs, but it's not really the point of the algorithm (just like "using a screwdriver" isn't the central part in "building an airplane").

So, look for any of the existing image scaling algorithms.