sz = 100
tmp = np.empty((4, sy, sz))
tmp[(0), :, :] = (1.2 - 0.8) * np.random.random_sample((sy, sz)) + 0.8
tmp[(1), :, :] = (1.2 - 0.8) * np.random.random_sample((sy, sz)) + 0.8
tmp[(2), :, :] = np.ones((sy, sz))
tmp[(3), :, :] = np.ones((sy, sz)) * np.pi / 4
param_ref = np.empty((4, sy, sz, sx))
for i in range(sx):
    param_ref[:, :, :, (i)] = tmp
x = np.linspace(0, 2 * np.pi)
factor = (1.1 - 0.9) * np.random.random_sample((sy, sz, sx)) + 0.9
data = f(x, *param_ref) * factor
func1d = lambda y, *args: optimize.curve_fit(f, xdata=x, ydata=y, *args)[0]
param = np.apply_along_axis(func1d, axis=2, arr=data)
Fitting a 3D array of data to a 1D function with numpy or scipy