Look into using astropy.modelling.parallel_fit_dask() #12
Description
The parallel_fit_dask() looks like a clean way to do the model-fitting.
With this function, there's no need for manipulating axes and manual fits.hdu.header items.
