@log_to_metadata
def ccd_process(ccd, oscan=None, trim=None, error=False, master_bias=None,
                dark_frame=None, master_flat=None, bad_pixel_mask=None,
                gain=None, readnoise=None, oscan_median=True, oscan_model=None,
                min_value=None, dark_exposure=None, data_exposure=None,
                exposure_key=None, exposure_unit=None,
                dark_scale=False, gain_corrected=True):
    """Perform basic processing on ccd data.

    The following steps can be included:

    * overscan correction (:func:`subtract_overscan`)
    * trimming of the image (:func:`trim_image`)
    * create deviation frame (:func:`create_deviation`)
    * gain correction (:func:`gain_correct`)
    * add a mask to the data
    * subtraction of master bias (:func:`subtract_bias`)
    * subtraction of a dark frame (:func:`subtract_dark`)
    * correction of flat field (:func:`flat_correct`)

    The task returns a processed `~astropy.nddata.CCDData` object.

    Parameters
    ----------
    ccd : `~astropy.nddata.CCDData`
        Frame to be reduced.

    oscan : `~astropy.nddata.CCDData`, str or None, optional
        For no overscan correction, set to None. Otherwise provide a region
        of ccd from which the overscan is extracted, using the FITS
        conventions for index order and index start, or a
        slice from ccd that contains the overscan.
        Default is ``None``.

    trim : str or None, optional
        For no trim correction, set to None. Otherwise provide a region
        of ccd from which the image should be trimmed, using the FITS
        conventions for index order and index start.
        Default is ``None``.

    error : bool, optional
        If True, create an uncertainty array for ccd.
        Default is ``False``.

    master_bias : `~astropy.nddata.CCDData` or None, optional
        A master bias frame to be subtracted from ccd. The unit of the
        master bias frame should match the unit of the image **after
        gain correction** if ``gain_corrected`` is True.
        Default is ``None``.

    dark_frame : `~astropy.nddata.CCDData` or None, optional
        A dark frame to be subtracted from the ccd. The unit of the
        master dark frame should match the unit of the image **after
        gain correction** if ``gain_corrected`` is True.
        Default is ``None``.

    master_flat : `~astropy.nddata.CCDData` or None, optional
        A master flat frame to be divided into ccd. The unit of the
        master flat frame should match the unit of the image **after
        gain correction** if ``gain_corrected`` is True.
        Default is ``None``.

    bad_pixel_mask : `numpy.ndarray` or None, optional
        A bad pixel mask for the data. The bad pixel mask should be in given
        such that bad pixels have a value of 1 and good pixels a value of 0.
        Default is ``None``.

    gain : `~astropy.units.Quantity` or None, optional
        Gain value to multiple the image by to convert to electrons.
        Default is ``None``.

    readnoise : `~astropy.units.Quantity` or None, optional
        Read noise for the observations. The read noise should be in
        electrons.
        Default is ``None``.

    oscan_median : bool, optional
        If true, takes the median of each line. Otherwise, uses the mean.
        Default is ``True``.

    oscan_model : `~astropy.modeling.Model` or None, optional
        Model to fit to the data. If None, returns the values calculated
        by the median or the mean.
        Default is ``None``.

    min_value : float or None, optional
        Minimum value for flat field. The value can either be None and no
        minimum value is applied to the flat or specified by a float which
        will replace all values in the flat by the min_value.
        Default is ``None``.

    dark_exposure : `~astropy.units.Quantity` or None, optional
        Exposure time of the dark image; if specified, must also provided
        ``data_exposure``.
        Default is ``None``.

    data_exposure : `~astropy.units.Quantity` or None, optional
        Exposure time of the science image; if specified, must also provided
        ``dark_exposure``.
        Default is ``None``.

    exposure_key : `~ccdproc.Keyword`, str or None, optional
        Name of key in image metadata that contains exposure time.
        Default is ``None``.

    exposure_unit : `~astropy.units.Unit` or None, optional
        Unit of the exposure time if the value in the meta data does not
        include a unit.
        Default is ``None``.

    dark_scale : bool, optional
        If True, scale the dark frame by the exposure times.
        Default is ``False``.

    gain_corrected : bool, optional
        If True, the ``master_bias``, ``master_flat``, and ``dark_frame``
        have already been gain corrected.  Default is ``True``.

    Returns
    -------
    occd : `~astropy.nddata.CCDData`
        Reduded ccd.

    Examples
    --------
    1. To overscan, trim and gain correct a data set::

        >>> import numpy as np
        >>> from astropy import units as u
        >>> from astropy.nddata import CCDData
        >>> from ccdproc import ccd_process
        >>> ccd = CCDData(np.ones([100, 100]), unit=u.adu)
        >>> nccd = ccd_process(ccd, oscan='[1:10,1:100]',
        ...                    trim='[10:100, 1:100]', error=False,
        ...                    gain=2.0*u.electron/u.adu)
    """
    # make a copy of the object
    nccd = ccd.copy()

    # apply the overscan correction
    if isinstance(oscan, CCDData):
        nccd = subtract_overscan(nccd, overscan=oscan,
                                 median=oscan_median,
                                 model=oscan_model)
    elif isinstance(oscan, str):
        nccd = subtract_overscan(nccd, fits_section=oscan,
                                 median=oscan_median,
                                 model=oscan_model)
    elif oscan is None:
        pass
    else:
        raise TypeError('oscan is not None, a string, or CCDData object.')

    # apply the trim correction
    if isinstance(trim, str):
        nccd = trim_image(nccd, fits_section=trim)
    elif trim is None:
        pass
    else:
        raise TypeError('trim is not None or a string.')

    # create the error frame
    if error and gain is not None and readnoise is not None:
        nccd = create_deviation(nccd, gain=gain, readnoise=readnoise)
    elif error and (gain is None or readnoise is None):
        raise ValueError(
            'gain and readnoise must be specified to create error frame.')

    # apply the bad pixel mask
    if isinstance(bad_pixel_mask, np.ndarray):
        nccd.mask = bad_pixel_mask
    elif bad_pixel_mask is None:
        pass
    else:
        raise TypeError('bad_pixel_mask is not None or numpy.ndarray.')

    # apply the gain correction
    if not (gain is None or isinstance(gain, Quantity)):
        raise TypeError('gain is not None or astropy.units.Quantity.')

    if gain is not None and gain_corrected:
        nccd = gain_correct(nccd, gain)

    # subtracting the master bias
    if isinstance(master_bias, CCDData):
        nccd = subtract_bias(nccd, master_bias)
    elif master_bias is None:
        pass
    else:
        raise TypeError(
            'master_bias is not None or a CCDData object.')

    # subtract the dark frame
    if isinstance(dark_frame, CCDData):
        nccd = subtract_dark(nccd, dark_frame, dark_exposure=dark_exposure,
                             data_exposure=data_exposure,
                             exposure_time=exposure_key,
                             exposure_unit=exposure_unit,
                             scale=dark_scale)
    elif dark_frame is None:
        pass
    else:
        raise TypeError(
            'dark_frame is not None or a CCDData object.')

    # test dividing the master flat
    if isinstance(master_flat, CCDData):
        nccd = flat_correct(nccd, master_flat, min_value=min_value)
    elif master_flat is None:
        pass
    else:
        raise TypeError(
            'master_flat is not None or a CCDData object.')

    # apply the gain correction only at the end if gain_corrected is False
    if gain is not None and not gain_corrected:
        nccd = gain_correct(nccd, gain)

    return nccd def cosmicray_median(ccd, error_image=None, thresh=5, mbox=11, gbox=0,
                     rbox=0):
    """
    Identify cosmic rays through median technique. The median technique
    identifies cosmic rays by identifying pixels by subtracting a median image
    from the initial data array.

    Parameters
    ----------
    ccd : `~astropy.nddata.CCDData`, `numpy.ndarray` or `numpy.ma.MaskedArray`
        Data to have cosmic ray cleaned.

    thresh : float, optional
        Threshold for detecting cosmic rays.
        Default is ``5``.

    error_image : `numpy.ndarray`, float or None, optional
        Error level. If None, the task will use the standard
        deviation of the data. If an ndarray, it should have the same shape
        as data.
        Default is ``None``.

    mbox : int, optional
        Median box for detecting cosmic rays.
        Default is ``11``.

    gbox : int, optional
        Box size to grow cosmic rays. If zero, no growing will be done.
        Default is ``0``.

    rbox : int, optional
        Median box for calculating replacement values. If zero, no pixels will
        be replaced.
        Default is ``0``.

    Notes
    -----
    Similar implementation to crmedian in iraf.imred.crutil.crmedian.

    Returns
    -------
    nccd : `~astropy.nddata.CCDData` or `numpy.ndarray`
        An object of the same type as ccd is returned. If it is a
        `~astropy.nddata.CCDData`, the mask attribute will also be updated with
        areas identified with cosmic rays masked.

    nccd : `numpy.ndarray`
        If an `numpy.ndarray` is provided as ccd, a boolean ndarray with the
        cosmic rays identified will also be returned.

    Examples
    --------
    1) Given an numpy.ndarray object, the syntax for running
       cosmicray_median would be:

       >>> newdata, mask = cosmicray_median(data, error_image=error,
       ...                                  thresh=5, mbox=11,
       ...                                  rbox=11, gbox=5)   # doctest: +SKIP

       where error is an array that is the same shape as data but
       includes the pixel error. This would return a data array, newdata,
       with the bad pixels replaced by the local median from a box of 11
       pixels; and it would return a mask indicating the bad pixels.

    2) Given an `~astropy.nddata.CCDData` object with an uncertainty frame, the syntax
       for running cosmicray_median would be:

       >>> newccd = cosmicray_median(ccd, thresh=5, mbox=11,
       ...                           rbox=11, gbox=5)   # doctest: +SKIP

       The newccd object will have bad pixels in its data array replace and the
       mask of the object will be created if it did not previously exist or be
       updated with the detected cosmic rays.
    """
    if isinstance(ccd, np.ndarray):
        data = ccd

        if error_image is None:
            error_image = data.std()
        else:
            if not isinstance(error_image, (float, np.ndarray)):
                raise TypeError('error_image is not a float or ndarray.')

        # create the median image
        marr = ndimage.median_filter(data, size=(mbox, mbox))

        # Only look at the data array
        if isinstance(data, np.ma.MaskedArray):
            data = data.data

        # Find the residual image
        rarr = (data - marr) / error_image

        # identify all sources
        crarr = (rarr > thresh)

        # grow the pixels
        if gbox > 0:
            crarr = ndimage.maximum_filter(crarr, gbox)

        # replace bad pixels in the image
        ndata = data.copy()
        if rbox > 0:
            data = np.ma.masked_array(data, (crarr == 1))
            mdata = ndimage.median_filter(data, rbox)
            ndata[crarr == 1] = mdata[crarr == 1]

        return ndata, crarr
    elif isinstance(ccd, CCDData):

        # set up the error image
        if error_image is None and ccd.uncertainty is not None:
            error_image = ccd.uncertainty.array
        if ccd.data.shape != error_image.shape:
            raise ValueError('error_image is not the same shape as data.')

        data, crarr = cosmicray_median(ccd.data, error_image=error_image,
                                       thresh=thresh, mbox=mbox, gbox=gbox,
                                       rbox=rbox)

        # create the new ccd data object
        nccd = ccd.copy()
        nccd.data = data
        if nccd.mask is None:
            nccd.mask = crarr
        else:
            nccd.mask = nccd.mask + crarr
        return nccd

    else:
        raise TypeError('ccd is not an numpy.ndarray or a CCDData object.') @deprecated_renamed_argument('pssl', None, '2.3.0',
                             arg_in_kwargs=True,
                             alternative='The pssl keyword will be removed in '
                                'ccdproc 3.0. Use inbkg instead to have '
                                'astroscrappy temporarily remove the background '
                                'during processing.')
def cosmicray_lacosmic(ccd, sigclip=4.5, sigfrac=0.3,
                       objlim=5.0, gain=1.0, readnoise=6.5,
                       satlevel=65535.0, pssl=0.0, niter=4,
                       sepmed=True, cleantype='meanmask', fsmode='median',
                       psfmodel='gauss', psffwhm=2.5, psfsize=7,
                       psfk=None, psfbeta=4.765, verbose=False,
                       gain_apply=True,
                       inbkg=None, invar=None):
    r"""
    Identify cosmic rays through the L.A. Cosmic technique. The L.A. Cosmic
    technique identifies cosmic rays by identifying pixels based on a variation
    of the Laplacian edge detection. The algorithm is an implementation of the
    code describe in van Dokkum (2001) [1]_ as implemented by McCully (2014)
    [2]_. If you use this algorithm, please cite these two works.


    Parameters
    ----------
    ccd : `~astropy.nddata.CCDData` or `numpy.ndarray`
        Data to have cosmic ray cleaned.

    gain_apply : bool, optional
        If ``True``, **return gain-corrected data**, with correct units,
        otherwise do not gain-correct the data. Default is ``True`` to
        preserve backwards compatibility.

    sigclip : float, optional
        Laplacian-to-noise limit for cosmic ray detection. Lower values will
        flag more pixels as cosmic rays. Default: 4.5.

    sigfrac : float, optional
        Fractional detection limit for neighboring pixels. For cosmic ray
        neighbor pixels, a Laplacian-to-noise detection limit of
        sigfrac * sigclip will be used. Default: 0.3.

    objlim : float, optional
        Minimum contrast between Laplacian image and the fine structure image.
        Increase this value if cores of bright stars are flagged as cosmic
        rays. Default: 5.0.

    inbkg : float numpy array, optional
        A pre-determined background image, to be subtracted from ``indat``
        before running the main detection algorithm.
        This is used primarily with spectroscopic data, to remove
        sky lines and the cross-section of an object continuum during
        iteration, "protecting" them from spurious rejection (see the above
        paper). This background is not removed from the final, cleaned output
        (``cleanarr``). This should be in units of "counts", the same units of indat.
        This inbkg should be free from cosmic rays. When estimating the cosmic-ray
        free noise of the image, we will treat ``inbkg`` as a constant Poisson
        contribution to the variance.

    invar : float numpy array, optional
        A pre-determined estimate of the data variance (ie. noise squared) in
        each pixel, generated by previous processing of ``indat``. If provided,
        this is used in place of an internal noise model based on ``indat``,
        ``gain`` and ``readnoise``. This still gets median filtered and cleaned
        internally, to estimate what the noise in each pixel *would* be in the
        absence of cosmic rays. This should be in units of "counts" squared.

    pssl : float, optional
        Previously subtracted sky level in ADU. We always need to work in
        electrons for cosmic ray detection, so we need to know the sky level
        that has been subtracted so we can add it back in. Default: 0.0.

    gain : float or `~astropy.units.Quantity`, optional
        Gain of the image (electrons / ADU). We always need to work in
        electrons for cosmic ray detection. Default: 1.0

    readnoise : float, optional
        Read noise of the image (electrons). Used to generate the noise model
        of the image. Default: 6.5.

    satlevel : float, optional
        Saturation level of the image (electrons). This value is used to
        detect saturated stars and pixels at or above this level are added to
        the mask. Default: 65535.0.

    niter : int, optional
        Number of iterations of the LA Cosmic algorithm to perform. Default: 4.

    sepmed : bool, optional
        Use the separable median filter instead of the full median filter.
        The separable median is not identical to the full median filter, but
        they are approximately the same, the separable median filter is
        significantly faster, and still detects cosmic rays well. Note, this is
        a performance feature, and not part of the original L.A. Cosmic.
        Default: True

    cleantype : str, optional
        Set which clean algorithm is used:

        - ``"median"``: An unmasked 5x5 median filter.
        - ``"medmask"``: A masked 5x5 median filter.
        - ``"meanmask"``: A masked 5x5 mean filter.
        - ``"idw"``: A masked 5x5 inverse distance weighted interpolation.

        Default: ``"meanmask"``.

    fsmode : str, optional
        Method to build the fine structure image:

        - ``"median"``: Use the median filter in the standard LA Cosmic \
          algorithm.
        - ``"convolve"``: Convolve the image with the psf kernel to calculate \
          the fine structure image.

        Default: ``"median"``.

    psfmodel : str, optional
        Model to use to generate the psf kernel if fsmode == 'convolve' and
        psfk is None. The current choices are Gaussian and Moffat profiles:

        - ``"gauss"`` and ``"moffat"`` produce circular PSF kernels.
        - The ``"gaussx"`` and ``"gaussy"`` produce Gaussian kernels in the x \
          and y directions respectively.

        Default: ``"gauss"``.

    psffwhm : float, optional
        Full Width Half Maximum of the PSF to use to generate the kernel.
        Default: 2.5.

    psfsize : int, optional
        Size of the kernel to calculate. Returned kernel will have size
        psfsize x psfsize. psfsize should be odd. Default: 7.

    psfk : `numpy.ndarray` (with float dtype) or None, optional
        PSF kernel array to use for the fine structure image if
        ``fsmode == 'convolve'``. If None and ``fsmode == 'convolve'``, we
        calculate the psf kernel using ``psfmodel``. Default: None.

    psfbeta : float, optional
        Moffat beta parameter. Only used if ``fsmode=='convolve'`` and
        ``psfmodel=='moffat'``. Default: 4.765.

    verbose : bool, optional
        Print to the screen or not. Default: False.

    Notes
    -----
    Implementation of the cosmic ray identification L.A.Cosmic:
    http://www.astro.yale.edu/dokkum/lacosmic/

    Returns
    -------
    nccd : `~astropy.nddata.CCDData` or `numpy.ndarray`
        An object of the same type as ccd is returned. If it is a
        `~astropy.nddata.CCDData`, the mask attribute will also be updated with
        areas identified with cosmic rays masked. **By default, the image is
        multiplied by the gain.** You can control this behavior with the
        ``gain_apply`` argument.

    crmask : `numpy.ndarray`
        If an `numpy.ndarray` is provided as ccd, a boolean ndarray with the
        cosmic rays identified will also be returned.

    References
    ----------
    .. [1] van Dokkum, P; 2001, "Cosmic-Ray Rejection by Laplacian Edge
       Detection". The Publications of the Astronomical Society of the
       Pacific, Volume 113, Issue 789, pp. 1420-1427.
       doi: 10.1086/323894

    .. [2] McCully, C., 2014, "Astro-SCRAPPY",
       https://github.com/astropy/astroscrappy

    Examples
    --------
    1) Given an numpy.ndarray object, the syntax for running
       cosmicrar_lacosmic would be:

       >>> newdata, mask = cosmicray_lacosmic(data, sigclip=5)  #doctest: +SKIP

       where the error is an array that is the same shape as data but
       includes the pixel error. This would return a data array, newdata,
       with the bad pixels replaced by the local median from a box of 11
       pixels; and it would return a mask indicating the bad pixels.

    2) Given an `~astropy.nddata.CCDData` object with an uncertainty frame, the syntax
       for running cosmicrar_lacosmic would be:

       >>> newccd = cosmicray_lacosmic(ccd, sigclip=5)   # doctest: +SKIP

       The newccd object will have bad pixels in its data array replace and the
       mask of the object will be created if it did not previously exist or be
       updated with the detected cosmic rays.
    """
    from astroscrappy import detect_cosmics
    from astroscrappy import __version__ as asy_version

    # If we didn't get a quantity, put them in, with unit specified by the
    # documentation above.
    if not isinstance(gain, u.Quantity):
        # Gain will change the value, so use the proper units
        gain = gain * u.electron / u.adu

    # Set the units of readnoise to electrons, as specified in the
    # documentation, if no unit is present.
    if not isinstance(readnoise, u.Quantity):
        readnoise = readnoise * u.electron

    # Handle transition from old astroscrappy interface to new
    old_astroscrappy_interface = (pkgversion.parse(asy_version) <
                                  pkgversion.parse('1.1.0'))

    # Use this dictionary to define which keyword arguments are actually
    # passed to astroscrappy.
    asy_background_kwargs = {}

    # This is for handling the transition in astroscrappy versions
    data_offset = 0

    # Handle setting up the keyword arguments for both interfaces
    if old_astroscrappy_interface: # pragma: no cover
        new_args = dict(inbkg=inbkg, invar=invar)
        bad_args = []
        for k, v in new_args.items():
            if v is not None:
                bad_args.append(k)

        if bad_args:
            s = 's' if len(bad_args) > 1 else ''
            bads = ', '.join(bad_args)
            raise TypeError(f'The argument{s} {bads} only valid for astroscrappy '
                            '1.1.0 or higher.')

        if pssl != 0:
            asy_background_kwargs = dict(pssl=pssl)

    else:
        if pssl != 0:
            if (inbkg is not None) or (invar is not None):
                raise ValueError('Cannot set both pssl and inbkg')

            # The old version of astroscrappy added the bkg back in
            # if pssl was provided. The new one does not, so set an offset
            # here that we later add in then take out.
            data_offset = pssl

        asy_background_kwargs = dict(inbkg=inbkg, invar=invar)

    if isinstance(ccd, np.ndarray):
        data = ccd

        crmask, cleanarr = detect_cosmics(
            data + data_offset, inmask=None, sigclip=sigclip,
            sigfrac=sigfrac, objlim=objlim, gain=gain.value,
            readnoise=readnoise.value, satlevel=satlevel,
            niter=niter, sepmed=sepmed, cleantype=cleantype,
            fsmode=fsmode, psfmodel=psfmodel, psffwhm=psffwhm,
            psfsize=psfsize, psfk=psfk, psfbeta=psfbeta,
            verbose=verbose,
            **asy_background_kwargs)

        cleanarr = cleanarr - data_offset
        cleanarr = _astroscrappy_gain_apply_helper(cleanarr,
                                                   gain,
                                                   gain_apply,
                                                   old_astroscrappy_interface)

        return cleanarr, crmask

    elif isinstance(ccd, CCDData):
        # Start with a check for a special case: ccd is in electron, and
        # gain and readnoise have no units. In that case we issue a warning
        # instead of raising an error to avoid crashing user's pipelines.
        if ccd.unit.is_equivalent(u.electron) and gain.value != 1.0:
            warnings.warn("Image unit is electron but gain value "
                          "is not 1.0. Data maybe end up being gain "
                          "corrected twice.")

        else:
            if ((readnoise.unit == u.electron)
                and (ccd.unit == u.electron)
                and (gain.value == 1.0)):
                gain = gain.value * u.one
            # Check unit consistency before taking the time to check for
            # cosmic rays.
            if not (gain * ccd).unit.is_equivalent(readnoise.unit):
                raise ValueError('Inconsistent units for gain ({}) '.format(gain.unit) +
                                 ' ccd ({}) and readnoise ({}).'.format(ccd.unit,
                                                                        readnoise.unit))

        crmask, cleanarr = detect_cosmics(
            ccd.data + data_offset, inmask=ccd.mask,
            sigclip=sigclip, sigfrac=sigfrac, objlim=objlim, gain=gain.value,
            readnoise=readnoise.value, satlevel=satlevel,
            niter=niter, sepmed=sepmed, cleantype=cleantype,
            fsmode=fsmode, psfmodel=psfmodel, psffwhm=psffwhm,
            psfsize=psfsize, psfk=psfk, psfbeta=psfbeta, verbose=verbose,
            **asy_background_kwargs)

        # create the new ccd data object
        nccd = ccd.copy()

        cleanarr = cleanarr - data_offset
        cleanarr = _astroscrappy_gain_apply_helper(cleanarr,
                                                   gain,
                                                   gain_apply,
                                                   old_astroscrappy_interface)

        # Fix the units if the gain is being applied.
        nccd.unit = ccd.unit * gain.unit

        nccd.data = cleanarr
        if nccd.mask is None:
            nccd.mask = crmask
        else:
            nccd.mask = nccd.mask + crmask

        return nccd

    else:
        raise TypeError('ccd is not a CCDData or ndarray object.') @log_to_metadata
def create_deviation(ccd_data, gain=None, readnoise=None, disregard_nan=False):
    """
    Create a uncertainty frame. The function will update the uncertainty
    plane which gives the standard deviation for the data. Gain is used in
    this function only to scale the data in constructing the deviation; the
    data is not scaled.

    Parameters
    ----------
    ccd_data : `~astropy.nddata.CCDData`
        Data whose deviation will be calculated.

    gain : `~astropy.units.Quantity` or None, optional
        Gain of the CCD; necessary only if ``ccd_data`` and ``readnoise``
        are not in the same units. In that case, the units of ``gain``
        should be those that convert ``ccd_data.data`` to the same units as
        ``readnoise``.
        Default is ``None``.

    readnoise : `~astropy.units.Quantity` or None, optional
        Read noise per pixel.
        Default is ``None``.

    disregard_nan: boolean
        If ``True``, any value of nan in the output array will be replaced by
        the readnoise.

    {log}

    Raises
    ------
    UnitsError
        Raised if ``readnoise`` units are not equal to product of ``gain`` and
        ``ccd_data`` units.

    Returns
    -------
    ccd : `~astropy.nddata.CCDData`
        CCDData object with uncertainty created; uncertainty is in the same
        units as the data in the parameter ``ccd_data``.

    """
    if gain is not None and not isinstance(gain, Quantity):
        raise TypeError('gain must be a astropy.units.Quantity.')

    if readnoise is None:
        raise ValueError('must provide a readnoise.')

    if not isinstance(readnoise, Quantity):
        raise TypeError('readnoise must be a astropy.units.Quantity.')

    if gain is None:
        gain = 1.0 * u.dimensionless_unscaled

    if gain.unit * ccd_data.unit != readnoise.unit:
        raise u.UnitsError("units of data, gain and readnoise do not match.")

    # Need to convert Quantity to plain number because NDData data is not
    # a Quantity. All unit checking should happen prior to this point.
    gain_value = float(gain / gain.unit)
    readnoise_value = float(readnoise / readnoise.unit)

    # remove values that might be negative or treat as nan
    data = gain_value * ccd_data.data
    mask = (data < 0)
    if disregard_nan:
        data[mask] = 0
    else:
        data[mask] = np.nan
        logging.warning('Negative values in array will be replaced with nan')

    # calculate the deviation
    var = (data + readnoise_value ** 2) ** 0.5

    # ensure uncertainty and image data have same unit
    ccd = ccd_data.copy()
    var /= gain_value
    ccd.uncertainty = StdDevUncertainty(var)
    return ccd @log_to_metadata
def flat_correct(ccd, flat, min_value=None, norm_value=None):
    """Correct the image for flat fielding.

    The flat field image is normalized by its mean or a user-supplied value
    before flat correcting.

    Parameters
    ----------
    ccd : `~astropy.nddata.CCDData`
        Data to be transformed.

    flat : `~astropy.nddata.CCDData`
        Flatfield to apply to the data.

    min_value : float or None, optional
        Minimum value for flat field. The value can either be None and no
        minimum value is applied to the flat or specified by a float which
        will replace all values in the flat by the min_value.
        Default is ``None``.

    norm_value : float or None, optional
        If not ``None``, normalize flat field by this argument rather than the
        mean of the image. This allows fixing several different flat fields to
        have the same scale. If this value is negative or 0, a ``ValueError``
        is raised. Default is ``None``.

    {log}

    Returns
    -------
    ccd : `~astropy.nddata.CCDData`
        CCDData object with flat corrected.
    """
    # Use the min_value to replace any values in the flat
    use_flat = flat
    if min_value is not None:
        flat_min = flat.copy()
        flat_min.data[flat_min.data < min_value] = min_value
        use_flat = flat_min

    # If a norm_value was input and is positive, use it to scale the flat
    if norm_value is not None and norm_value > 0:
        flat_mean_val = norm_value
    elif norm_value is not None:
        # norm_value was set to a bad value
        raise ValueError('norm_value must be greater than zero.')
    else:
        # norm_value was not set, use mean of the image.
        flat_mean_val = use_flat.data.mean()

    # Normalize the flat.
    flat_mean = flat_mean_val * use_flat.unit
    flat_normed = use_flat.divide(flat_mean)

    # divide through the flat
    flat_corrected = ccd.divide(flat_normed)

    flat_corrected.meta = ccd.meta.copy()
    return flat_corrected @log_to_metadata
def gain_correct(ccd, gain, gain_unit=None):
    """Correct the gain in the image.

    Parameters
    ----------
    ccd : `~astropy.nddata.CCDData`
      Data to have gain corrected.

    gain : `~astropy.units.Quantity` or `~ccdproc.Keyword`
      gain value for the image expressed in electrons per adu.

    gain_unit : `~astropy.units.Unit` or None, optional
        Unit for the ``gain``; used only if ``gain`` itself does not provide
        units.
        Default is ``None``.

    {log}

    Returns
    -------
    result : `~astropy.nddata.CCDData`
      CCDData object with gain corrected.
    """
    if isinstance(gain, Keyword):
        gain_value = gain.value_from(ccd.header)
    elif isinstance(gain, numbers.Number) and gain_unit is not None:
        gain_value = gain * u.Unit(gain_unit)
    else:
        gain_value = gain

    result = ccd.multiply(gain_value)
    result.meta = ccd.meta.copy()
    return result @log_to_metadata
def subtract_bias(ccd, master):
    """
    Subtract master bias from image.

    Parameters
    ----------
    ccd : `~astropy.nddata.CCDData`
        Image from which bias will be subtracted.

    master : `~astropy.nddata.CCDData`
        Master image to be subtracted from ``ccd``.

    {log}

    Returns
    -------
    result : `~astropy.nddata.CCDData`
        CCDData object with bias subtracted.
    """

    try:
        result = ccd.subtract(master)
    except ValueError as e:
        if 'operand units' in str(e):
            raise u.UnitsError("Unit '{}' of the uncalibrated image does not "
                               "match unit '{}' of the calibration "
                               "image".format(ccd.unit, master.unit))
        else:
            raise e

    result.meta = ccd.meta.copy()
    return result @log_to_metadata
def subtract_dark(ccd, master, dark_exposure=None, data_exposure=None,
                  exposure_time=None, exposure_unit=None,
                  scale=False):
    """
    Subtract dark current from an image.

    Parameters
    ----------
    ccd : `~astropy.nddata.CCDData`
        Image from which dark will be subtracted.

    master : `~astropy.nddata.CCDData`
        Dark image.

    dark_exposure : `~astropy.units.Quantity` or None, optional
        Exposure time of the dark image; if specified, must also provided
        ``data_exposure``.
        Default is ``None``.

    data_exposure : `~astropy.units.Quantity` or None, optional
        Exposure time of the science image; if specified, must also provided
        ``dark_exposure``.
        Default is ``None``.

    exposure_time : str or `~ccdproc.Keyword` or None, optional
        Name of key in image metadata that contains exposure time.
        Default is ``None``.

    exposure_unit : `~astropy.units.Unit` or None, optional
        Unit of the exposure time if the value in the meta data does not
        include a unit.
        Default is ``None``.

    scale: bool, optional
        If True, scale the dark frame by the exposure times.
        Default is ``False``.

    {log}

    Returns
    -------
    result : `~astropy.nddata.CCDData`
        Dark-subtracted image.
    """
    if ccd.shape != master.shape:
        err_str = "operands could not be subtracted with shapes {} {}".format(ccd.shape, master.shape)
        raise ValueError(err_str)

    if not (isinstance(ccd, CCDData) and isinstance(master, CCDData)):
        raise TypeError("ccd and master must both be CCDData objects.")

    if (data_exposure is not None and
            dark_exposure is not None and
            exposure_time is not None):
        raise TypeError("specify either exposure_time or "
                        "(dark_exposure and data_exposure), not both.")

    if data_exposure is None and dark_exposure is None:
        if exposure_time is None:
            raise TypeError("must specify either exposure_time or both "
                            "dark_exposure and data_exposure.")
        if isinstance(exposure_time, Keyword):
            data_exposure = exposure_time.value_from(ccd.header)
            dark_exposure = exposure_time.value_from(master.header)
        else:
            data_exposure = ccd.header[exposure_time]
            dark_exposure = master.header[exposure_time]

    if not (isinstance(dark_exposure, Quantity) and
            isinstance(data_exposure, Quantity)):
        if exposure_time:
            try:
                data_exposure *= exposure_unit
                dark_exposure *= exposure_unit
            except TypeError:
                raise TypeError("must provide unit for exposure time.")
        else:
            raise TypeError("exposure times must be astropy.units.Quantity "
                            "objects.")

    try:
        if scale:
            master_scaled = master.copy()
            # data_exposure and dark_exposure are both quantities,
            # so we can just have subtract do the scaling
            master_scaled = master_scaled.multiply(data_exposure /
                                                   dark_exposure)
            result = ccd.subtract(master_scaled)
        else:
            result = ccd.subtract(master)
    except (u.UnitsError, u.UnitConversionError, ValueError) as e:
        # Astropy LTS (v1) returns a ValueError, not a UnitsError, so catch
        # that if it appears to really be a UnitsError.
        if (isinstance(e, ValueError) and
                'operand units' not in str(e) and
                astropy.__version__.startswith('1.0')):
            raise e

        # Make the error message a little more explicit than what is returned
        # by default.
        raise u.UnitsError("Unit '{}' of the uncalibrated image does not "
                           "match unit '{}' of the calibration "
                           "image".format(ccd.unit, master.unit))

    result.meta = ccd.meta.copy()
    return result @log_to_metadata
def subtract_overscan(ccd, overscan=None, overscan_axis=1, fits_section=None,
                      median=False, model=None):
    """
    Subtract the overscan region from an image.

    Parameters
    ----------
    ccd : `~astropy.nddata.CCDData`
        Data to have overscan frame corrected.

    overscan : `~astropy.nddata.CCDData` or None, optional
        Slice from ``ccd`` that contains the overscan. Must provide either
        this argument or ``fits_section``, but not both.
        Default is ``None``.

    overscan_axis : 0, 1 or None, optional
        Axis along which overscan should combined with mean or median. Axis
        numbering follows the *python* convention for ordering, so 0 is the
        first axis and 1 is the second axis.

        If overscan_axis is explicitly set to None, the axis is set to
        the shortest dimension of the overscan section (or 1 in case
        of a square overscan).
        Default is ``1``.

    fits_section : str or None, optional
        Region of ``ccd`` from which the overscan is extracted, using the FITS
        conventions for index order and index start. See Notes and Examples
        below. Must provide either this argument or ``overscan``, but not both.
        Default is ``None``.

    median : bool, optional
        If true, takes the median of each line. Otherwise, uses the mean.
        Default is ``False``.

    model : `~astropy.modeling.Model` or None, optional
        Model to fit to the data. If None, returns the values calculated
        by the median or the mean.
        Default is ``None``.

    {log}

    Raises
    ------
    TypeError
        A TypeError is raised if either ``ccd`` or ``overscan`` are not the
        correct objects.

    Returns
    -------
    ccd : `~astropy.nddata.CCDData`
        CCDData object with overscan subtracted.

    Notes
    -----
    The format of the ``fits_section`` string follow the rules for slices that
    are consistent with the FITS standard (v3) and IRAF usage of keywords like
    TRIMSEC and BIASSEC. Its indexes are one-based, instead of the
    python-standard zero-based, and the first index is the one that increases
    most rapidly as you move through the array in memory order, opposite the
    python ordering.

    The 'fits_section' argument is provided as a convenience for those who are
    processing files that contain TRIMSEC and BIASSEC. The preferred, more
    pythonic, way of specifying the overscan is to do it by indexing the data
    array directly with the ``overscan`` argument.

    Examples
    --------
    Creating a 100x100 array containing ones just for demonstration purposes::

        >>> import numpy as np
        >>> from astropy import units as u
        >>> arr1 = CCDData(np.ones([100, 100]), unit=u.adu)

    The statement below uses all rows of columns 90 through 99 as the
    overscan::

        >>> no_scan = subtract_overscan(arr1, overscan=arr1[:, 90:100])
        >>> assert (no_scan.data == 0).all()

    This statement does the same as the above, but with a FITS-style section::

        >>> no_scan = subtract_overscan(arr1, fits_section='[91:100, :]')
        >>> assert (no_scan.data == 0).all()

    Spaces are stripped out of the ``fits_section`` string.

    """
    if not (isinstance(ccd, CCDData) or isinstance(ccd, np.ndarray)):
        raise TypeError('ccddata is not a CCDData or ndarray object.')

    if ((overscan is not None and fits_section is not None) or
            (overscan is None and fits_section is None)):
        raise TypeError('specify either overscan or fits_section, but not '
                        'both.')

    if (overscan is not None) and (not isinstance(overscan, CCDData)):
        raise TypeError('overscan is not a CCDData object.')

    if (fits_section is not None and
            not isinstance(fits_section, str)):
        raise TypeError('overscan is not a string.')

    if fits_section is not None:
        overscan = ccd[slice_from_string(fits_section, fits_convention=True)]

    if overscan_axis is None:
        overscan_axis = 0 if overscan.shape[1] > overscan.shape[0] else 1

    if median:
        oscan = np.median(overscan.data, axis=overscan_axis)
    else:
        oscan = np.mean(overscan.data, axis=overscan_axis)

    if model is not None:
        of = fitting.LinearLSQFitter()
        yarr = np.arange(len(oscan))
        oscan = of(model, yarr, oscan)
        oscan = oscan(yarr)
        if overscan_axis == 1:
            oscan = np.reshape(oscan, (oscan.size, 1))
        else:
            oscan = np.reshape(oscan, (1, oscan.size))
    else:
        if overscan_axis == 1:
            oscan = np.reshape(oscan, oscan.shape + (1,))
        else:
            oscan = np.reshape(oscan, (1,) + oscan.shape)

    subtracted = ccd.copy()

    # subtract the overscan
    subtracted.data = ccd.data - oscan
    return subtracted @log_to_metadata
def transform_image(ccd, transform_func, **kwargs):
    """Transform the image.

    Using the function specified by transform_func, the transform will
    be applied to data, uncertainty, and mask in ccd.

    Parameters
    ----------
    ccd : `~astropy.nddata.CCDData`
        Data to be transformed.

    transform_func : callable
        Function to be used to transform the data, mask and uncertainty.

    kwargs :
        Additional keyword arguments to be used by the transform_func.

    {log}

    Returns
    -------
    ccd : `~astropy.nddata.CCDData`
        A transformed CCDData object.

    Notes
    -----
    At this time, transform will be applied to the uncertainty data but it
    will only transform the data. This will not properly handle uncertainties
    that arise due to correlation between the pixels.

    These should only be geometric transformations of the images. Other
    methods should be used if the units of ccd need to be changed.

    Examples
    --------
    Given an array that is 100x100::

        >>> import numpy as np
        >>> from astropy import units as u
        >>> arr1 = CCDData(np.ones([100, 100]), unit=u.adu)

    The syntax for transforming the array using
    `scipy.ndimage.shift`::

        >>> from scipy.ndimage import shift
        >>> from ccdproc import transform_image
        >>> transformed = transform_image(arr1, shift, shift=(5.5, 8.1))
    """
    # check that it is a ccddata object
    if not isinstance(ccd, CCDData):
        raise TypeError('ccd is not a CCDData.')

    # make a copy of the object
    nccd = ccd.copy()

    # transform the image plane
    try:
        nccd.data = transform_func(nccd.data, **kwargs)
    except TypeError as exc:
        if 'is not callable' in str(exc):
            raise TypeError('transform_func is not a callable.')
        raise

    # transform the uncertainty plane if it exists
    if nccd.uncertainty is not None:
        nccd.uncertainty.array = transform_func(nccd.uncertainty.array,
                                                **kwargs)

    # transform the mask plane
    if nccd.mask is not None:
        mask = transform_func(nccd.mask, **kwargs)
        nccd.mask = mask > 0

    if nccd.wcs is not None:
        warn = 'WCS information may be incorrect as no transformation was applied to it'
        warnings.warn(warn, UserWarning)

    return nccd @log_to_metadata
def trim_image(ccd, fits_section=None):
    """
    Trim the image to the dimensions indicated.

    Parameters
    ----------
    ccd : `~astropy.nddata.CCDData`
        CCD image to be trimmed, sliced if desired.

    fits_section : str or None, optional
        Region of ``ccd`` from which the overscan is extracted; see
        `~ccdproc.subtract_overscan` for details.
        Default is ``None``.

    {log}

    Returns
    -------
    trimmed_ccd : `~astropy.nddata.CCDData`
        Trimmed image.

    Examples
    --------
    Given an array that is 100x100,

        >>> import numpy as np
        >>> from astropy import units as u
        >>> arr1 = CCDData(np.ones([100, 100]), unit=u.adu)

    the syntax for trimming this to keep all of the first index but only the
    first 90 rows of the second index is

        >>> trimmed = trim_image(arr1[:, :90])
        >>> trimmed.shape
        (100, 90)
        >>> trimmed.data[0, 0] = 2
        >>> arr1.data[0, 0]
        1.0

    This both trims *and makes a copy* of the image.

    Indexing the image directly does *not* do the same thing, quite:

        >>> not_really_trimmed = arr1[:, :90]
        >>> not_really_trimmed.data[0, 0] = 2
        >>> arr1.data[0, 0]
        2.0

    In this case, ``not_really_trimmed`` is a view of the underlying array
    ``arr1``, not a copy.
    """
    if (fits_section is not None and
            not isinstance(fits_section, str)):
        raise TypeError("fits_section must be a string.")
    trimmed = ccd.copy()
    if fits_section:
        python_slice = slice_from_string(fits_section, fits_convention=True)
        trimmed = trimmed[python_slice]

    return trimmed @log_to_metadata
def wcs_project(ccd, target_wcs, target_shape=None, order='bilinear'):
    """
    Given a CCDData image with WCS, project it onto a target WCS and
    return the reprojected data as a new CCDData image.

    Any flags, weight, or uncertainty are ignored in doing the
    reprojection.

    Parameters
    ----------
    ccd : `~astropy.nddata.CCDData`
        Data to be projected.

    target_wcs : `~astropy.wcs.WCS` object
        WCS onto which all images should be projected.

    target_shape : two element list-like or None, optional
        Shape of the output image. If omitted, defaults to the shape of the
        input image.
        Default is ``None``.

    order : str, optional
        Interpolation order for re-projection. Must be one of:

        + 'nearest-neighbor'
        + 'bilinear'
        + 'biquadratic'
        + 'bicubic'

        Default is ``'bilinear'``.

    {log}

    Returns
    -------
    ccd : `~astropy.nddata.CCDData`
        A transformed CCDData object.
    """
    from astropy.nddata.ccddata import _generate_wcs_and_update_header
    from reproject import reproject_interp

    if not (ccd.wcs.is_celestial and target_wcs.is_celestial):
        raise ValueError('one or both WCS is not celestial.')

    if target_shape is None:
        target_shape = ccd.shape

    projected_image_raw, _ = reproject_interp((ccd.data, ccd.wcs),
                                              target_wcs,
                                              shape_out=target_shape,
                                              order=order)

    reprojected_mask = None
    if ccd.mask is not None:
        reprojected_mask, _ = reproject_interp((ccd.mask, ccd.wcs),
                                               target_wcs,
                                               shape_out=target_shape,
                                               order=order)
        # Make the mask 1 if the reprojected mask pixel value is non-zero.
        # A small threshold is included to allow for some rounding in
        # reproject_interp.
        reprojected_mask = reprojected_mask > 1e-8

    # The reprojection will contain nan for any pixels for which the source
    # was outside the original image. Those should be masked also.
    output_mask = np.isnan(projected_image_raw)

    if reprojected_mask is not None:
        output_mask = output_mask | reprojected_mask

    # Need to scale counts by ratio of pixel areas
    area_ratio = (proj_plane_pixel_area(target_wcs) /
                  proj_plane_pixel_area(ccd.wcs))

    # If nothing ended up masked, don't create a mask.
    if not output_mask.any():
        output_mask = None

    # If there are any wcs keywords in the header, remove them
    hdr, _ = _generate_wcs_and_update_header(ccd.header)

    nccd = CCDData(area_ratio * projected_image_raw, wcs=target_wcs,
                   mask=output_mask,
                   header=hdr, unit=ccd.unit)

    return nccd class Keyword:
    """
    """
    def __init__(self, name, unit=None, value=None):
        self._name = name
        self._unit = unit
        self.value = value

    @property
    def name(self):
        return self._name

    @property
    def unit(self):
        return self._unit

    @property
    def value(self):
        return self._value

    @value.setter
    def value(self, value):
        if value is None:
            self._value = value
        elif isinstance(value, Quantity):
            self._unit = value.unit
            self._value = value
        elif isinstance(value, str):
            if self.unit is not None:
                raise ValueError("keyword with a unit cannot have a "
                                 "string value.")
            else:
                self._value = value
        else:
            if self.unit is None:
                raise ValueError("no unit provided. Set value with "
                                 "an astropy.units.Quantity.")
            self._value = value * self.unit

    def value_from(self, header):
        """
        Set value of keyword from FITS header.

        Parameters
        ----------
        header : `~astropy.io.fits.Header`
            FITS header containing a value for this keyword.
        """
        value_from_header = header[self.name]
        self.value = value_from_header
        return self.value def _blkavg(data, newshape):
    """
    Block average an array such that it has the new shape.

    Parameters
    ----------
    data : `numpy.ndarray` or `numpy.ma.MaskedArray`
        Data to average.

    newshape : tuple
        Tuple containing the new shape for the array.

    Returns
    -------
    output : `numpy.ndarray` or `numpy.ma.MaskedArray`
        An array with the new shape and the average of the pixels.

    Raises
    ------
    TypeError
        A type error is raised if data is not an `numpy.ndarray`.

    ValueError
        A value error is raised if the dimensions of new shape is not equal
        to data.

    Notes
    -----
    This is based on the scipy cookbook for rebinning:
    http://wiki.scipy.org/Cookbook/Rebinning
    """
    # check to see that is in a nddata type
    if not isinstance(data, np.ndarray):
        raise TypeError('data is not a ndarray object.')

    # check to see that the two arrays are going to be the same length
    if len(data.shape) != len(newshape):
        raise ValueError('newshape does not have the same dimensions as data.')

    shape = data.shape
    lenShape = len(shape)
    factor = np.asarray(shape)/np.asarray(newshape)

    evList = ['data.reshape('] + \
        ['newshape[%d],int(factor[%d]),' % (i, i) for i in range(lenShape)] + \
        [')'] + ['.mean(%d)' % (i + 1) for i in range(lenShape)]

    return eval(''.join(evList)) def ccd_data(data_size=DEFAULT_DATA_SIZE,
             data_scale=DEFAULT_DATA_SCALE,
             data_mean=DEFAULT_DATA_MEAN,
             rng_seed=DEFAULT_SEED):
    """
    Return a CCDData object with units of ADU.

    The size of the data array is 100x100 but can be changed using the marker
    @pytest.mark.data_size(N) on the test function, where N should be the
    desired dimension.

    Data values are initialized to random numbers drawn from a normal
    distribution with mean of 0 and scale 1.

    The scale can be changed with the marker @pytest.marker.scale(s) on the
    test function, where s is the desired scale.

    The mean can be changed with the marker @pytest.marker.scale(m) on the
    test function, where m is the desired mean.
    """
    size = data_size
    scale = data_scale
    mean = data_mean

    with NumpyRNGContext(rng_seed):
        data = np.random.normal(loc=mean, size=[size, size], scale=scale)

    fake_meta = {'my_key': 42, 'your_key': 'not 42'}
    ccd = CCDData(data, unit=u.adu)
    ccd.header = fake_meta
    return ccd
import numpy as np
import astropy.units as u
import astropy
import pytest
import skimage
from astropy.io import fits
from astropy.modeling import models
from astropy.units.quantity import Quantity
from astropy.wcs import WCS
from astropy.utils.exceptions import AstropyUserWarning
from astropy.nddata import StdDevUncertainty, CCDData
from numpy.testing import assert_array_equal
from ccdproc.core import (
    ccd_process, cosmicray_median, cosmicray_lacosmic, create_deviation,
    flat_correct, gain_correct, subtract_bias, subtract_dark, subtract_overscan,
    transform_image, trim_image, wcs_project, Keyword)
from ccdproc.core import _blkavg
from ccdproc.tests.pytest_fixtures import ccd_data as ccd_data_func
    from ..core import block_reduce, block_average, block_replicate

def test_subtract_overscan_fails():
    ccd_data = ccd_data_func()
    # Do we get an error if the *image* is neither CCDData nor an array?
    with pytest.raises(TypeError):
        subtract_overscan(3, np.zeros((5, 5)))
    # Do we get an error if the *overscan* is not an image or an array?
    with pytest.raises(TypeError):
        subtract_overscan(np.zeros((10, 10)), 3, median=False, model=None)
    # Do we get an error if we specify both overscan and fits_section?
    with pytest.raises(TypeError):
        subtract_overscan(ccd_data, overscan=ccd_data[0:10],
                          fits_section='[1:10]')
    # Do we raise an error if we specify neither overscan nor fits_section?
    with pytest.raises(TypeError):
        subtract_overscan(ccd_data)
    # Does a fits_section which is not a string raise an error?
    with pytest.raises(TypeError):
        subtract_overscan(ccd_data, fits_section=5)


def test_trim_image_fits_section_requires_string():
    ccd_data = ccd_data_func()
    with pytest.raises(TypeError):
        trim_image(ccd_data, fits_section=5)


@pytest.mark.parametrize('mask_data, uncertainty', [
                         (False, False),
                         (True, True)])
def test_trim_image_fits_section(mask_data, uncertainty):
    ccd_data = ccd_data_func(data_size=50)
    if mask_data:
        ccd_data.mask = np.zeros_like(ccd_data)
    if uncertainty:
        err = np.random.normal(size=ccd_data.shape)
        ccd_data.uncertainty = StdDevUncertainty(err)

    trimmed = trim_image(ccd_data, fits_section='[20:40,:]')
    # FITS reverse order, bounds are inclusive and starting index is 1-based
    assert trimmed.shape == (50, 21)
    np.testing.assert_array_equal(trimmed.data, ccd_data[:, 19:40])
    if mask_data:
        assert trimmed.shape == trimmed.mask.shape
    if uncertainty:
        assert trimmed.shape == trimmed.uncertainty.array.shape


def test_trim_image_no_section():
    ccd_data = ccd_data_func(data_size=50)
    trimmed = trim_image(ccd_data[:, 19:40])
    assert trimmed.shape == (50, 21)
    np.testing.assert_array_equal(trimmed.data, ccd_data[:, 19:40])


def test_trim_with_wcs_alters_wcs():
    ccd_data = ccd_data_func()
    # WCS construction example pulled form astropy.wcs docs
    wcs = WCS(naxis=2)
    wcs.wcs.crpix = np.array(ccd_data.shape)/2
    wcs.wcs.cdelt = np.array([-0.066667, 0.066667])
    wcs.wcs.crval = [0, -90]
    wcs.wcs.ctype = ["RA---AIR", "DEC--AIR"]
    wcs.wcs.set_pv([(2, 1, 45.0)])
    ccd_wcs = CCDData(ccd_data, wcs=wcs)
    # The trim below should subtract 10 from the 2nd element of crpix.
    # (Second element because the FITS convention for index ordering is
    #  opposite that of python)
    trimmed = trim_image(ccd_wcs[10:, :])
    assert trimmed.wcs.wcs.crpix[1] == wcs.wcs.crpix[1] - 10


def test_subtract_bias():
    ccd_data = ccd_data_func()
    data_avg = ccd_data.data.mean()
    bias_level = 5.0
    ccd_data.data = ccd_data.data + bias_level
    ccd_data.header['key'] = 'value'
    master_bias_array = np.zeros_like(ccd_data.data) + bias_level
    master_bias = CCDData(master_bias_array, unit=ccd_data.unit)
    no_bias = subtract_bias(ccd_data, master_bias, add_keyword=None)
    # Does the data we are left with have the correct average?
    np.testing.assert_almost_equal(no_bias.data.mean(), data_avg)
    # With logging turned off, metadata should not change
    assert no_bias.header == ccd_data.header
    del no_bias.header['key']
    assert 'key' in ccd_data.header
    assert no_bias.header is not ccd_data.header


def test_subtract_bias_fails():
    ccd_data = ccd_data_func(data_size=50)
    # Should fail if shapes don't match
    bias = CCDData(np.array([200, 200]), unit=u.adu)
    with pytest.raises(ValueError):
        subtract_bias(ccd_data, bias)
    # Should fail because units don't match
    bias = CCDData(np.zeros_like(ccd_data), unit=u.meter)
    with pytest.raises(u.UnitsError):
        subtract_bias(ccd_data, bias)


@pytest.mark.parametrize('exposure_keyword', [True, False])
@pytest.mark.parametrize('explicit_times', [True, False])
@pytest.mark.parametrize('scale', [True, False])
def test_subtract_dark(explicit_times, scale, exposure_keyword):
    ccd_data = ccd_data_func()
    exptime = 30.0
    exptime_key = 'exposure'
    exposure_unit = u.second
    dark_level = 1.7
    master_dark_data = np.zeros_like(ccd_data.data) + dark_level
    master_dark = CCDData(master_dark_data, unit=u.adu)
    master_dark.header[exptime_key] = 2 * exptime
    dark_exptime = master_dark.header[exptime_key]
    ccd_data.header[exptime_key] = exptime
    dark_exposure_unit = exposure_unit
    if explicit_times:
        # Test case when units of dark and data exposures are different
        dark_exposure_unit = u.minute

        dark_sub = subtract_dark(ccd_data, master_dark,