The two most general, publicly available packages for galaxy image fitting (other than GALFIT) are probably
ProFit. (Note that I am the developer of Imfit.)
There is also Lenstronomy, which is specialized for fitting gravitational lenses; this might be more relevant to your particular needs, though I know very little about it.
EDIT: I moved the last paragraph to the beginning because it's really a Python question, but I left my "manual" method answer in case it helps, too.
Forgot to add this link to Python code on Github. It helped me get started with displaying images using Python and dealing with the headers without altering them. On line 19, I changed it to "...
One thing you might try doing is to use the SkyServer Navigate interface to see if the object was imaged by SDSS. Enter the name in the "name" box in the upper left and then click on the "Resolve" button. If an image with the galaxy shows up, click on the "Object with spectra" checkbox in the "Drawing options" panel on ...
The information you need to recreate the wavelength array is in the World Coordinate System (WCS) of the header, specifically:
CRPIX1 = 1.00
CRVAL1 = 3500.0000 / central wavelength of first pixel
CDELT1 = 0.900000 / linear dispersion (Angstrom/pixel)
which lists the starting/reference pixel of the wavelength array (...
Skycat also stores the WCS information when selecting an image region from a fits file. Not sure how exactly this can be done in non-interactive mode but many features can in fact be "batched". Give it a try. Best luck!