Remove `wcs1d_multispec` reader.

[kelle]

Specutils should be able to handle this. Not sure if it's already been implemented in specutils or not but we should try to remove this from the AstrodbKit code base.

This is the test file:
filename = "https://s3.amazonaws.com/bdnyc/optical_spectra/U10929.fits"

def identify_wcs1d_multispec(origin, *args, **kwargs):
    """
    Identifier for WCS1D multispec
    """
    hdu = kwargs.get("hdu", 0)

    # Check if number of axes is one and dimension of WCS is greater than one
    with read_fileobj_or_hdulist(*args, **kwargs) as hdulist:
        return (
            hdulist[hdu].header.get("WCSDIM", 1) > 1
            and hdulist[hdu].header["NAXIS"] > 1
            and "WAT0_001" in hdulist[hdu].header
            and hdulist[hdu].header.get("WCSDIM", 1) == hdulist[hdu].header["NAXIS"]
            and "LINEAR" in hdulist[hdu].header.get("CTYPE1", "")
        )


@data_loader("wcs1d-multispec", identifier=identify_wcs1d_multispec, extensions=["fits"], dtype=Spectrum, priority=10)
def wcs1d_multispec_loader(file_obj, flux_unit=None, hdu=0, verbose=False, **kwargs):
    """
    Loader for multiextension spectra as wcs1d. Adapted from wcs1d_fits_loader

    Parameters
    ----------
    file_obj : str, file-like or HDUList
        FITS file name, object (provided from name by Astropy I/O Registry),
        or HDUList (as resulting from astropy.io.fits.open()).
    flux_unit : :class:`~astropy.units.Unit` or str, optional
        Units of the flux for this spectrum. If not given (or None), the unit
        will be inferred from the BUNIT keyword in the header. Note that this
        unit will attempt to convert from BUNIT if BUNIT is present.
    hdu : int
        The index of the HDU to load into this spectrum.
    verbose : bool
        Print extra info.
    **kwargs
        Extra keywords for :func:`~specutils.io.parsing_utils.read_fileobj_or_hdulist`.

    Returns
    -------
    :class:`~specutils.Spectrum`
    """

    with read_fileobj_or_hdulist(file_obj, **kwargs) as hdulist:
        header = hdulist[hdu].header
        wcs = WCS(header)

        # Load data, convert units if BUNIT and flux_unit is provided and not the same
        if "BUNIT" in header:
            data = u.Quantity(hdulist[hdu].data, unit=header["BUNIT"])
            if u.A in data.unit.bases:
                data = data * u.A / u.AA  # convert ampere to Angroms
            if flux_unit is not None:
                data = data.to(flux_unit)
        else:
            data = u.Quantity(hdulist[hdu].data, unit=flux_unit)

    if wcs.wcs.cunit[0] == "" and "WAT1_001" in header:
        # Try to extract from IRAF-style card or use Angstrom as default.
        wat_dict = dict((rec.split("=") for rec in header["WAT1_001"].split()))
        unit = wat_dict.get("units", "Angstrom")
        if hasattr(u, unit):
            wcs.wcs.cunit[0] = unit
        else:  # try with unit name stripped of excess plural 's'...
            wcs.wcs.cunit[0] = unit.rstrip("s")
        if verbose:
            print(f"Extracted spectral axis unit '{unit}' from 'WAT1_001'")
    elif wcs.wcs.cunit[0] == "":
        wcs.wcs.cunit[0] = "Angstrom"

    # Compatibility attribute for lookup_table (gwcs) WCS
    wcs.unit = tuple(wcs.wcs.cunit)

    # Identify the correct parts of the data to store
    if len(data.shape) > 1:
        flux_data = data[0]
    else:
        flux_data = data
    uncertainty = None
    if "NAXIS3" in header:
        for i in range(header["NAXIS3"]):
            if "spectrum" in header.get(f"BANDID{i+1}", ""):
                flux_data = data[i]
            if "sigma" in header.get(f"BANDID{i+1}", ""):
                uncertainty = data[i]

    # Reshape arrays if needed
    if len(flux_data) == 1 and len(flux_data.shape) > 1:
        flux_data = np.reshape(flux_data, -1)
        if uncertainty is not None:
            uncertainty = np.reshape(uncertainty, -1)

    # Convert uncertainty to StdDevUncertainty array
    if uncertainty is not None:
        uncertainty = StdDevUncertainty(uncertainty)

    # Manually generate spectral axis
    pixels = [[i] + [0] * (wcs.naxis - 1) for i in range(wcs.pixel_shape[0])]
    spectral_axis = [i[0] for i in wcs.all_pix2world(pixels, 0)] * wcs.wcs.cunit[0]

    # Store header as metadata information
    meta = {"header": header}

    return Spectrum(flux=flux_data, spectral_axis=spectral_axis, uncertainty=uncertainty, meta=meta)