Python looping over fits datacube to fit spectra

Question

I have a fits datacube in Python where I need to fit the spectra of each pixel in the cube. My code for fitting just one spectrum works great, but for some reason I can't figure out how to translate that into fitting each pixel (100 pixels total for a 10by10).

ngc2617 = fits.open('NGC2617.fits')[1]
cube2617 = SpectralCube.read(ngc2617)
data2617 = ngc2617.data
nz,ny,nx = np.shape(data2617)

obs_2617 = cube2617.spectral_axis
z_2617 = 0.014326 #simbad
rest_2617 = np.asarray(obs_2617 / (1+z_2617))

flux_2617 = data2617[:,172,142]
flux_2617[np.isnan(flux_2617)]=0

halpha_chans = np.asarray(np.where((rest_2617 > 6365) & (rest_2617 < 6715)))
halphawaves = np.asarray(rest_2617[halpha_chans]).T
halphaflux = np.asarray(flux_2617[halpha_chans]).T
cont_chans = np.asarray(np.where((rest_2617 > 5500) & (rest_2617 < 5510)))
cont_2617 = np.asarray(flux_2617[cont_chans]).T


x_halpha = halphawaves.flatten()
y_halpha = halphaflux.flatten()
c = 2.998 * 10**5 #km/s

plt.figure(figsize=[10,8])
plt.subplot(2,1,1)
plt.plot(x_halpha,y_halpha, '--', label = 'data')
def Gauss(x, A, x0, fwhm):
    sigma = (fwhm * x0) / (2.355 * c)
    return A * np.exp(-(x - x0) ** 2 / (2 * sigma ** 2))

def func(x, A0, A1, x0, fwhm, fwhm1, cont):
    param1 = cont
    param2 = Gauss(x, A0, x0 - 14.75, fwhm) + Gauss(x, A1, x0, fwhm) + Gauss(x, A0*3, x0 + 20.66, fwhm) + Gauss(x, A1, x0, fwhm1)
    return param1 + param2

cont_init = np.average(cont_2617)
A0_init = np.max(y_halpha[0:150])-cont_init
A1_init = np.max(y_halpha)-cont_init
x0_init = x_halpha[np.where(y_halpha==np.max(y_halpha))][0]
fwhm_init = (x_halpha[np.where(y_halpha==np.max(y_halpha))][0] - x_halpha[np.where((np.min(abs(np.max(y_halpha)/2 - y_halpha))))][0])*2
fwhm1_init = (x_halpha[np.where(y_halpha==np.max(y_halpha))][0] - x_halpha[np.where((np.min(abs(np.max(y_halpha)/2 - y_halpha[0:150]))))][0])*2
parameters, covariance = curve_fit(func, x_halpha, y_halpha, p0=[A0_init, A1_init, x0_init, fwhm_init, fwhm1_init, cont_init])

fit_y = func(x_halpha, *parameters)
fit1_NII = Gauss(x_halpha, A0_init, parameters[2]-14.75, parameters[3]) + parameters[5]
fit2_NII = Gauss(x_halpha, A0_init, parameters[2]+20.66, parameters[3]) + parameters[5]
fit_Halph = Gauss(x_halpha, A1_init, parameters[2], parameters[3]) + parameters[5]
plt.plot(x_halpha, fit_y, '-', label='full fit')
plt.plot(x_halpha, fit1_NII, '-', label='NII 6549.86 fit')
plt.plot(x_halpha, fit2_NII, '-', label='NII 6585.27 fit')
plt.plot(x_halpha, fit_Halph, '-', label='H\u03B1 fit')
plt.title('NGC2617 H\u03B1 & NII Spectra of Central Point with Gaussian Fit', fontsize = 15)
plt.xlabel('wavelength [$\AA$]', fontsize = 12)
plt.ylabel(r'flux [10$^{-20}$ $erg * s^{-1}$ * $cm^{-2}$ * $\AA^{-1}$]', fontsize = 12)
plt.grid(True)
plt.legend()

plt.subplot(2,1,2)
plt.plot(np.linspace(0,0,6800), color = 'r', label = '0 point')
plt.scatter(x_halpha, y_halpha-fit_y, marker = 'o', label='residuals')
plt.xlim(halphawaves[0],halphawaves[-1])
plt.xlabel('wavelength [$\AA$]', fontsize = 12)
plt.ylabel('residuals', fontsize = 12)
plt.grid(True)
plt.legend()
plt.tight_layout()
plt.show()

this is my code for the single spectrum...can anyone help me create a loop that will do this but for a full cube? Thanks!