def make_mask(filename, ext, trail_coords, sublen=75, subwidth=200, order=3,
              sigma=4, pad=10, plot=False, verbose=False):
    """Create DQ mask for an image for a given satellite trail.
    This mask can be added to existing DQ data using :func:`update_dq`.

    .. note::

        Unlike :func:`detsat`, multiprocessing is not available for
        this function.

    Parameters
    ----------
    filename : str
        FITS image filename.

    ext : int, str, or tuple
        Extension for science data, as accepted by ``astropy.io.fits``.

    trail_coords : ndarray
        One of the trails returned by :func:`detsat`.
        This must be in the format of ``[[x0, y0], [x1, y1]]``.

    sublen : int, optional
        Length of strip to use as the fitting window for the trail.

    subwidth : int, optional
        Width of box to fit trail on.

    order : int, optional
        The order of the spline interpolation for image rotation.
        See :func:`skimage.transform.rotate`.

    sigma : float, optional
        Sigma of the satellite trail for detection. If points are
        a given sigma above the background in the subregion then it is
        marked as a satellite. This may need to be lowered for resolved
        trails.

    pad : int, optional
        Amount of extra padding in pixels to give the satellite mask.

    plot : bool, optional
        Plot the result.

    verbose : bool, optional
        Print extra information to the terminal, mostly for debugging.

    Returns
    -------
    mask : ndarray
        Boolean array marking the satellite trail with `True`.

    Raises
    ------
    ImportError
        Missing scipy or skimage>=0.11 packages.

    IndexError
        Invalid subarray indices.

    ValueError
        Image has no positive values, trail subarray too small, or
        trail profile not found.

    """
    if not HAS_OPDEP:
        raise ImportError('Missing scipy or skimage>=0.11 packages')

    if verbose:
        t_beg = time.time()

    fname = '{0}[{1}]'.format(filename, ext)
    image = fits.getdata(filename, ext)

    dx = image.max()
    if dx <= 0:
        raise ValueError('Image has no positive values')

    # rescale the image
    image = image / dx
    # make sure everything is at least 0
    image[image < 0] = 0

    (x0, y0), (x1, y1) = trail_coords  # p0, p1

    #  Find out how much to rotate the image
    rad = np.arctan2(y1 - y0, x1 - x0)
    newrad = (np.pi * 2) - rad
    deg = np.degrees(rad)

    if verbose:
        print('Rotation: {0}'.format(deg))

    rotate = transform.rotate(image, deg, resize=True, order=order)

    if plot and plt is not None:
        plt.ion()
        mean = np.median(image)
        stddev = image.std()
        lower = mean - stddev
        upper = mean + stddev

        fig1, ax1 = plt.subplots()
        ax1.imshow(image, vmin=lower, vmax=upper, cmap=plt.cm.gray)
        ax1.set_title(fname)

        fig2, ax2 = plt.subplots()
        ax2.imshow(rotate, vmin=lower, vmax=upper, cmap=plt.cm.gray)
        ax2.set_title('{0} rotated by {1} deg'.format(fname, deg))

        plt.draw()

    #  Will do all of this in the loop, but want to make sure there is a
    #  good point first and that there is indeed a profile to fit.
    #  get starting point
    sx, sy = _rotate_point((x0, y0), newrad, image.shape, rotate.shape)

    #  start with one subarray around p0
    dx = int(subwidth / 2)
    ix0, ix1, iy0, iy1 = _get_valid_indices(
        rotate.shape, sx - dx, sx + dx, sy - sublen, sy + sublen)
    subr = rotate[iy0:iy1, ix0:ix1]
    if len(subr) <= sublen:
        raise ValueError('Trail subarray size is {0} but expected {1} or '
                         'larger'.format(len(subr), sublen))

    # Flatten the array so we are looking along rows
    # Take median of each row, should filter out most outliers
    # This list will get appended in the loop
    medarr = np.median(subr, axis=1)
    flat = [medarr]

    # get the outliers
    # mean = biweight_location(medarr)
    mean = sigma_clipped_stats(medarr)[0]
    stddev = biweight_midvariance(medarr)

    # only flag things that are sigma from the mean
    z = np.where(medarr > (mean + (sigma * stddev)))[0]

    if plot and plt is not None:
        fig1, ax1 = plt.subplots()
        ax1.plot(medarr, 'b.')
        ax1.plot(z, medarr[z], 'r.')
        ax1.set_xlabel('Index')
        ax1.set_ylabel('Value')
        ax1.set_title('Median array in flat[0]')
        plt.draw()

    # Make sure there is something in the first pass before trying to move on
    if len(z) < 1:
        raise ValueError(
            'First look at finding a profile failed. '
            'Nothing found at {0} from background! '
            'Adjust parameters and try again.'.format(sigma))

    # get the bounds of the flagged points
    lower = z.min()
    upper = z.max()
    diff = upper - lower

    # add in a pading value to make sure all of the wings are accounted for
    lower = lower - pad
    upper = upper + pad

    # for plotting see how the profile was made (append to plot above)
    if plot and plt is not None:
        padind = np.arange(lower, upper)
        ax1.plot(padind, medarr[padind], 'yx')
        plt.draw()

    # start to create a mask
    mask = np.zeros(rotate.shape)
    lowerx, upperx, lowery, uppery = _get_valid_indices(
        mask.shape, np.floor(sx - subwidth), np.ceil(sx + subwidth),
        np.floor(sy - sublen + lower), np.ceil(sy - sublen + upper))
    mask[lowery:uppery, lowerx:upperx] = 1

    done = False
    first = True
    nextx = upperx  # np.ceil(sx + subwidth)
    centery = np.ceil(lowery + diff)  # np.ceil(sy - sublen + lower + diff)
    counter = 0

    while not done:
        # move to the right of the centerpoint first. do the same
        # as above but keep moving right until the edge is hit.
        ix0, ix1, iy0, iy1 = _get_valid_indices(
            rotate.shape, nextx - dx, nextx + dx,
            centery - sublen, centery + sublen)
        subr = rotate[iy0:iy1, ix0:ix1]

        # determines the edge, if the subr is not good, then the edge was
        # hit.
        if 0 in subr.shape:
            if verbose:
                print('Hit edge, subr shape={0}, first={1}'.format(
                    subr.shape, first))
            if first:
                first = False
                centery = sy
                nextx = sx
            else:
                done = True
            continue

        medarr = np.median(subr, axis=1)
        flat.append(medarr)

        # mean = biweight_location(medarr)
        mean = sigma_clipped_stats(medarr, sigma=sigma)[0]
        # Might give RuntimeWarning
        stddev = biweight_midvariance(medarr)
        z = np.where(medarr > (mean + (sigma * stddev)))[0]

        if len(z) < 1:
            if first:
                if verbose:
                    print('No good profile found for counter={0}. Start '
                          'moving left from starting point.'.format(counter))
                centery = sy
                nextx = sx
                first = False
            else:
                if verbose:
                    print('z={0} is less than 1, subr shape={1}, '
                          'we are done'.format(z, subr.shape))
                done = True
            continue

        # get the bounds of the flagged points
        lower = z.min()
        upper = z.max()
        diff = upper - lower

        # add in a pading value to make sure all of the wings
        # are accounted for
        lower = np.floor(lower - pad)
        upper = np.ceil(upper + pad)
        lowerx, upperx, lowery, uppery = _get_valid_indices(
            mask.shape,
            np.floor(nextx - subwidth),
            np.ceil(nextx + subwidth),
            np.floor(centery - sublen + lower),
            np.ceil(centery - sublen + upper))
        mask[lowery:uppery, lowerx:upperx] = 1

        # lower_p = (lowerx, lowery)
        upper_p = (upperx, uppery)
        # lower_t = _rotate_point(
        #     lower_p, newrad, image.shape, rotate.shape, reverse=True)
        upper_t = _rotate_point(
            upper_p, newrad, image.shape, rotate.shape, reverse=True)

        # lowy = np.floor(lower_t[1])
        highy = np.ceil(upper_t[1])
        # lowx = np.floor(lower_t[0])
        highx = np.ceil(upper_t[0])

        # Reset the next subr to be at the center of the profile
        if first:
            nextx = nextx + dx
            centery = lowery + diff  # centery - sublen + lower + diff

            if (nextx + subwidth) > rotate.shape[1]:
                if verbose:
                    print('Hit rotate edge at counter={0}'.format(counter))
                first = False
            elif (highy > image.shape[0]) or (highx > image.shape[1]):
                if verbose:
                    print('Hit image edge at counter={0}'.format(counter))
                first = False

            if not first:
                centery = sy
                nextx = sx

        # Not first, this is the pass the other way.
        else:
            nextx = nextx - dx
            centery = lowery + diff  # centery - sublen + lower + diff

            if (nextx - subwidth) < 0:
                if verbose:
                    print('Hit rotate edge at counter={0}'.format(counter))
                done = True
            elif (highy > image.shape[0]) or (highx > image.shape[1]):
                if verbose:
                    print('Hit image edge at counter={0}'.format(counter))
                done = True

        counter += 1

        # make sure it does not try to go infinetly
        if counter > 500:
            if verbose:
                print('Too many loops, exiting')
            done = True
    # End while

    rot = transform.rotate(mask, -deg, resize=True, order=1)
    ix0 = (rot.shape[1] - image.shape[1]) / 2
    iy0 = (rot.shape[0] - image.shape[0]) / 2
    lowerx, upperx, lowery, uppery = _get_valid_indices(
        rot.shape, ix0, image.shape[1] + ix0, iy0, image.shape[0] + iy0)
    mask = rot[lowery:uppery, lowerx:upperx]

    if mask.shape != image.shape:
        warnings.warn(
            'Output mask shape is {0} but input image shape is '
            '{1}'.format(mask.shape, image.shape), AstropyUserWarning)

    # Change to boolean mask
    mask = mask.astype(np.bool)

    if plot and plt is not None:
        # debugging array
        test = image.copy()
        test[mask] = 0

        mean = np.median(test)
        stddev = test.std()
        lower = mean - stddev
        upper = mean + stddev

        fig1, ax1 = plt.subplots()
        ax1.imshow(test, vmin=lower, vmax=upper, cmap=plt.cm.gray)
        ax1.set_title('Masked image')

        fig2, ax2 = plt.subplots()
        ax2.imshow(mask, cmap=plt.cm.gray)
        ax2.set_title('DQ mask')

        plt.draw()

    if verbose:
        t_end = time.time()
        print('Run time: {0} s'.format(t_end - t_beg))

    return mask
def make_mask ( filename , ext , trail_coords , sublen = 75 , subwidth = 200 , order = 3 , sigma = 4 , pad = 10 , plot = False , verbose = False ) : if not HAS_OPDEP : raise ImportError ( 'Missing scipy or skimage>=0.11 packages' ) if verbose : t_beg = time . time ( ) fname = '{0}[{1}]' . format ( filename , ext ) image = fits . getdata ( filename , ext ) dx = image . max ( ) if dx <= 0 : raise ValueError ( 'Image has no positive values' ) image = image / dx image [ image < 0 ] = 0 ( x0 , y0 ) , ( x1 , y1 ) = trail_coords rad = np . arctan2 ( y1 - y0 , x1 - x0 ) newrad = ( np . pi * 2 ) - rad deg = np . degrees ( rad ) if verbose : print ( 'Rotation: {0}' . format ( deg ) ) rotate = transform . rotate ( image , deg , resize = True , order = order ) if plot and plt is not None : plt . ion ( ) mean = np . median ( image ) stddev = image . std ( ) lower = mean - stddev upper = mean + stddev fig1 , ax1 = plt . subplots ( ) ax1 . imshow ( image , vmin = lower , vmax = upper , cmap = plt . cm . gray ) ax1 . set_title ( fname ) fig2 , ax2 = plt . subplots ( ) ax2 . imshow ( rotate , vmin = lower , vmax = upper , cmap = plt . cm . gray ) ax2 . set_title ( '{0} rotated by {1} deg' . format ( fname , deg ) ) plt . draw ( ) sx , sy = _rotate_point ( ( x0 , y0 ) , newrad , image . shape , rotate . shape ) dx = int ( subwidth / 2 ) ix0 , ix1 , iy0 , iy1 = _get_valid_indices ( rotate . shape , sx - dx , sx + dx , sy - sublen , sy + sublen ) subr = rotate [ iy0 : iy1 , ix0 : ix1 ] if len ( subr ) <= sublen : raise ValueError ( 'Trail subarray size is {0} but expected {1} or ' 'larger' . format ( len ( subr ) , sublen ) ) medarr = np . median ( subr , axis = 1 ) flat = [ medarr ] mean = sigma_clipped_stats ( medarr ) [ 0 ] stddev = biweight_midvariance ( medarr ) z = np . where ( medarr > ( mean + ( sigma * stddev ) ) ) [ 0 ] if plot and plt is not None : fig1 , ax1 = plt . subplots ( ) ax1 . plot ( medarr , 'b.' ) ax1 . plot ( z , medarr [ z ] , 'r.' ) ax1 . set_xlabel ( 'Index' ) ax1 . set_ylabel ( 'Value' ) ax1 . set_title ( 'Median array in flat[0]' ) plt . draw ( ) if len ( z ) < 1 : raise ValueError ( 'First look at finding a profile failed. ' 'Nothing found at {0} from background! ' 'Adjust parameters and try again.' . format ( sigma ) ) lower = z . min ( ) upper = z . max ( ) diff = upper - lower lower = lower - pad upper = upper + pad if plot and plt is not None : padind = np . arange ( lower , upper ) ax1 . plot ( padind , medarr [ padind ] , 'yx' ) plt . draw ( ) mask = np . zeros ( rotate . shape ) lowerx , upperx , lowery , uppery = _get_valid_indices ( mask . shape , np . floor ( sx - subwidth ) , np . ceil ( sx + subwidth ) , np . floor ( sy - sublen + lower ) , np . ceil ( sy - sublen + upper ) ) mask [ lowery : uppery , lowerx : upperx ] = 1 done = False first = True nextx = upperx centery = np . ceil ( lowery + diff ) counter = 0 while not done : ix0 , ix1 , iy0 , iy1 = _get_valid_indices ( rotate . shape , nextx - dx , nextx + dx , centery - sublen , centery + sublen ) subr = rotate [ iy0 : iy1 , ix0 : ix1 ] if 0 in subr . shape : if verbose : print ( 'Hit edge, subr shape={0}, first={1}' . format ( subr . shape , first ) ) if first : first = False centery = sy nextx = sx else : done = True continue medarr = np . median ( subr , axis = 1 ) flat . append ( medarr ) mean = sigma_clipped_stats ( medarr , sigma = sigma ) [ 0 ] stddev = biweight_midvariance ( medarr ) z = np . where ( medarr > ( mean + ( sigma * stddev ) ) ) [ 0 ] if len ( z ) < 1 : if first : if verbose : print ( 'No good profile found for counter={0}. Start ' 'moving left from starting point.' . format ( counter ) ) centery = sy nextx = sx first = False else : if verbose : print ( 'z={0} is less than 1, subr shape={1}, ' 'we are done' . format ( z , subr . shape ) ) done = True continue lower = z . min ( ) upper = z . max ( ) diff = upper - lower lower = np . floor ( lower - pad ) upper = np . ceil ( upper + pad ) lowerx , upperx , lowery , uppery = _get_valid_indices ( mask . shape , np . floor ( nextx - subwidth ) , np . ceil ( nextx + subwidth ) , np . floor ( centery - sublen + lower ) , np . ceil ( centery - sublen + upper ) ) mask [ lowery : uppery , lowerx : upperx ] = 1 upper_p = ( upperx , uppery ) upper_t = _rotate_point ( upper_p , newrad , image . shape , rotate . shape , reverse = True ) highy = np . ceil ( upper_t [ 1 ] ) highx = np . ceil ( upper_t [ 0 ] ) if first : nextx = nextx + dx centery = lowery + diff if ( nextx + subwidth ) > rotate . shape [ 1 ] : if verbose : print ( 'Hit rotate edge at counter={0}' . format ( counter ) ) first = False elif ( highy > image . shape [ 0 ] ) or ( highx > image . shape [ 1 ] ) : if verbose : print ( 'Hit image edge at counter={0}' . format ( counter ) ) first = False if not first : centery = sy nextx = sx else : nextx = nextx - dx centery = lowery + diff if ( nextx - subwidth ) < 0 : if verbose : print ( 'Hit rotate edge at counter={0}' . format ( counter ) ) done = True elif ( highy > image . shape [ 0 ] ) or ( highx > image . shape [ 1 ] ) : if verbose : print ( 'Hit image edge at counter={0}' . format ( counter ) ) done = True counter += 1 if counter > 500 : if verbose : print ( 'Too many loops, exiting' ) done = True rot = transform . rotate ( mask , - deg , resize = True , order = 1 ) ix0 = ( rot . shape [ 1 ] - image . shape [ 1 ] ) / 2 iy0 = ( rot . shape [ 0 ] - image . shape [ 0 ] ) / 2 lowerx , upperx , lowery , uppery = _get_valid_indices ( rot . shape , ix0 , image . shape [ 1 ] + ix0 , iy0 , image . shape [ 0 ] + iy0 ) mask = rot [ lowery : uppery , lowerx : upperx ] if mask . shape != image . shape : warnings . warn ( 'Output mask shape is {0} but input image shape is ' '{1}' . format ( mask . shape , image . shape ) , AstropyUserWarning ) mask = mask . astype ( np . bool ) if plot and plt is not None : test = image . copy ( ) test [ mask ] = 0 mean = np . median ( test ) stddev = test . std ( ) lower = mean - stddev upper = mean + stddev fig1 , ax1 = plt . subplots ( ) ax1 . imshow ( test , vmin = lower , vmax = upper , cmap = plt . cm . gray ) ax1 . set_title ( 'Masked image' ) fig2 , ax2 = plt . subplots ( ) ax2 . imshow ( mask , cmap = plt . cm . gray ) ax2 . set_title ( 'DQ mask' ) plt . draw ( ) if verbose : t_end = time . time ( ) print ( 'Run time: {0} s' . format ( t_end - t_beg ) ) return mask
Create DQ mask for an image for a given satellite trail.
    This mask can be added to existing DQ data using :func:`update_dq`.

    .. note::

        Unlike :func:`detsat`, multiprocessing is not available for
        this function.

    Parameters
    ----------
    filename : str
        FITS image filename.

    ext : int, str, or tuple
        Extension for science data, as accepted by ``astropy.io.fits``.

    trail_coords : ndarray
        One of the trails returned by :func:`detsat`.
        This must be in the format of ``[[x0, y0], [x1, y1]]``.

    sublen : int, optional
        Length of strip to use as the fitting window for the trail.

    subwidth : int, optional
        Width of box to fit trail on.

    order : int, optional
        The order of the spline interpolation for image rotation.
        See :func:`skimage.transform.rotate`.

    sigma : float, optional
        Sigma of the satellite trail for detection. If points are
        a given sigma above the background in the subregion then it is
        marked as a satellite. This may need to be lowered for resolved
        trails.

    pad : int, optional
        Amount of extra padding in pixels to give the satellite mask.

    plot : bool, optional
        Plot the result.

    verbose : bool, optional
        Print extra information to the terminal, mostly for debugging.

    Returns
    -------
    mask : ndarray
        Boolean array marking the satellite trail with `True`.

    Raises
    ------
    ImportError
        Missing scipy or skimage>=0.11 packages.

    IndexError
        Invalid subarray indices.

    ValueError
        Image has no positive values, trail subarray too small, or
        trail profile not found.