<commit_before>import os
import numpy as np

def sky_fibers_for_brick(survey, brickname, bands=['g','r','z'],
                         apertures_arcsec=[0.5, 0.75, 1., 1.5, 2., 3.5, 5., 7.]):
    '''
    Produces a table of possible DESI sky fiber locations for a given
    "brickname" (eg, "0001p000") read from the given LegacySurveyData object *survey*.
    '''
    import fitsio
    from astrometry.util.fits import fits_table
    from astrometry.util.util import Tan
    import photutils
    from legacypipe.utils import find_unique_pixels

    fn = survey.find_file('blobmap', brick=brickname)
    blobs = fitsio.read(fn)
    print('Blobs:', blobs.min(), blobs.max())
    header = fitsio.read_header(fn)
    wcs = Tan(header)

    goodpix = (blobs == -1)
    for band in bands:
        fn = survey.find_file('nexp', brick=brickname, band=band)
        if not os.path.exists(fn):
            # Skip
            continue
        nexp = fitsio.read(fn)
        goodpix[nexp == 0] = False

    # Cut to unique brick area... required since the blob map drops
    # blobs that are completely outside the brick's unique area, thus
    # those locations are not masked.
    brick = survey.get_brick_by_name(brickname)
    H,W = wcs.shape
    U = find_unique_pixels(wcs, W, H, None, brick.ra1, brick.ra2,
                           brick.dec1, brick.dec2)
    goodpix[U == 0] = False
    del U

    x,y,blobdist = sky_fiber_locations(goodpix)

    skyfibers = fits_table()
    skyfibers.brickid = np.zeros(len(x), np.int32) + brick.brickid
    skyfibers.brickname = np.array([brickname] * len(x))
    skyfibers.x = x.astype(np.int16)
    skyfibers.y = y.astype(np.int16)
    skyfibers.blobdist = blobdist
    skyfibers.ra,skyfibers.dec = wcs.pixelxy2radec(x+1, y+1)

    pixscale = wcs.pixel_scale()
    apertures = np.array(apertures_arcsec) / pixscale

    # Now, do aperture photometry at these points in the coadd images
    for band in bands:
        imfn = survey.find_file('image',  brick=brickname, band=band)
        ivfn = survey.find_file('invvar', brick=brickname, band=band)
        if not (os.path.exists(imfn) and os.path.exists(ivfn)):
            continue
        coimg = fitsio.read(imfn)
        coiv = fitsio.read(ivfn)

        apflux = np.zeros((len(skyfibers), len(apertures)), np.float32)
        apiv   = np.zeros((len(skyfibers), len(apertures)), np.float32)
        skyfibers.set('apflux_%s' % band, apflux)
        skyfibers.set('apflux_ivar_%s' % band, apiv)
        with np.errstate(divide='ignore', invalid='ignore'):
            imsigma = 1./np.sqrt(coiv)
            imsigma[coiv == 0] = 0
        apxy = np.vstack((skyfibers.x, skyfibers.y)).T
        for irad,rad in enumerate(apertures):
            aper = photutils.CircularAperture(apxy, rad)
            p = photutils.aperture_photometry(coimg, aper, error=imsigma)
            apflux[:,irad] = p.field('aperture_sum')
            err = p.field('aperture_sum_err')
            apiv  [:,irad] = 1. / err**2

    header = fitsio.FITSHDR()
    for i,ap in enumerate(apertures_arcsec):
        header.add_record(dict(name='AP%i' % i, value=ap, comment='Aperture radius (arcsec)'))
    skyfibers._header = header

    return skyfibers

def sky_fiber_locations(skypix, gridsize=300):
    '''Worker function for sky_fibers_for_brick--- implements the core
    trick of iteratively eroding the map of good sky locations to
    produce a distance-from-blobs map, and then return the max values
    in that map in each cell of a grid.
    '''
    # Select possible locations for sky fibers
    from scipy.ndimage.morphology import binary_dilation, binary_erosion
    from scipy.ndimage.measurements import label, find_objects, center_of_mass
    from scipy.ndimage.filters import gaussian_filter

    nerosions = np.zeros(skypix.shape, np.int16)
    nerosions += skypix
    element = np.ones((3,3), bool)
    while True:
        skypix = binary_erosion(skypix, structure=element)
        nerosions += skypix
        print('After erosion,', np.sum(skypix), 'sky pixels')
        if not np.any(skypix.ravel()):
            break

    # This is a hack to break ties in the integer 'nerosions' map.
    nerosions = gaussian_filter(nerosions.astype(np.float32), 1.0)
    peaks = (nerosions > 1)
    H,W = skypix.shape

    # find pixels that are larger than their 8 neighbors
    peaks[1:-1, 1:-1] &= (nerosions[1:-1,1:-1] >= nerosions[0:-2,1:-1])
    peaks[1:-1, 1:-1] &= (nerosions[1:-1,1:-1] >= nerosions[2:  ,1:-1])
    peaks[1:-1, 1:-1] &= (nerosions[1:-1,1:-1] >= nerosions[1:-1,0:-2])
    peaks[1:-1, 1:-1] &= (nerosions[1:-1,1:-1] >= nerosions[1:-1,2:  ])
    peaks[1:-1, 1:-1] &= (nerosions[1:-1,1:-1] >= nerosions[0:-2,0:-2])
    peaks[1:-1, 1:-1] &= (nerosions[1:-1,1:-1] >= nerosions[0:-2,2:  ])
    peaks[1:-1, 1:-1] &= (nerosions[1:-1,1:-1] >= nerosions[2:  ,0:-2])
    peaks[1:-1, 1:-1] &= (nerosions[1:-1,1:-1] >= nerosions[2:  ,2:  ])

    # Split the image into 300 x 300-pixel cells, choose the highest peak in each one
    # (note, this is ignoring the brick-to-brick margin in laying down the grid)
    sx,sy = [],[]
    xx = np.round(np.linspace(0, W, 1+np.ceil(W / gridsize))).astype(int)
    yy = np.round(np.linspace(0, H, 1+np.ceil(H / gridsize))).astype(int)
    for ylo,yhi in zip(yy, yy[1:]):
        for xlo,xhi in zip(xx, xx[1:]):
            # Find max pixel in box
            subne = nerosions[ylo:yhi, xlo:xhi]
            I = np.argmax(subne)
            # Find all pixels equal to the max and take the one closest to the center of mass.
            maxval = subne.flat[I]
            cy,cx = center_of_mass(subne == maxval)
            xg = np.arange(xhi-xlo)
            yg = np.arange(yhi-ylo)
            dd = np.exp(-((yg[:,np.newaxis] - cy)**2 + (xg[np.newaxis,:] - cx)**2))
            dd[subne != maxval] = 0
            I = np.argmax(dd.flat)
            iy,ix = np.unravel_index(I, subne.shape)
            sx.append(ix + xlo)
            sy.append(iy + ylo)

    sx = np.array(sx)
    sy = np.array(sy)
    return sx, sy, nerosions[sy,sx]

def sky_fiber_plots(survey, brickname, skyfibers, basefn, bands=['g','r','z']):
    from legacypipe.survey import get_rgb
    import fitsio
    import pylab as plt

    rgbkwargs = dict(mnmx=(-1,100.), arcsinh=1.)

    imgs = []
    for band in bands:
        fn = survey.find_file('image',  brick=brickname, band=band)
        imgs.append(fitsio.read(fn))
    rgb = get_rgb(imgs, bands, **rgbkwargs)

    ima = dict(interpolation='nearest', origin='lower')<commit_msg>adding licensing information to skyfibers code<commit_after># Licensed under a 3-clause BSD style license - see LICENSE.rst
# -*- coding: utf-8 -*-
"""
=====================
desitarget.skyfibers
=====================

Module dealing with the assignation of sky fibers at the pixel-level for target selection
"""
import os
import numpy as np

def sky_fibers_for_brick(survey, brickname, bands=['g','r','z'],
                         apertures_arcsec=[0.5, 0.75, 1., 1.5, 2., 3.5, 5., 7.]):
    '''
    Produces a table of possible DESI sky fiber locations for a given
    "brickname" (eg, "0001p000") read from the given LegacySurveyData object *survey*.
    '''
    import fitsio
    from astrometry.util.fits import fits_table
    from astrometry.util.util import Tan
    import photutils
    from legacypipe.utils import find_unique_pixels

    fn = survey.find_file('blobmap', brick=brickname)
    blobs = fitsio.read(fn)
    print('Blobs:', blobs.min(), blobs.max())
    header = fitsio.read_header(fn)
    wcs = Tan(header)

    goodpix = (blobs == -1)
    for band in bands:
        fn = survey.find_file('nexp', brick=brickname, band=band)
        if not os.path.exists(fn):
            # Skip
            continue
        nexp = fitsio.read(fn)
        goodpix[nexp == 0] = False

    # Cut to unique brick area... required since the blob map drops
    # blobs that are completely outside the brick's unique area, thus
    # those locations are not masked.
    brick = survey.get_brick_by_name(brickname)
    H,W = wcs.shape
    U = find_unique_pixels(wcs, W, H, None, brick.ra1, brick.ra2,
                           brick.dec1, brick.dec2)
    goodpix[U == 0] = False
    del U

    x,y,blobdist = sky_fiber_locations(goodpix)

    skyfibers = fits_table()
    skyfibers.brickid = np.zeros(len(x), np.int32) + brick.brickid
    skyfibers.brickname = np.array([brickname] * len(x))
    skyfibers.x = x.astype(np.int16)
    skyfibers.y = y.astype(np.int16)
    skyfibers.blobdist = blobdist
    skyfibers.ra,skyfibers.dec = wcs.pixelxy2radec(x+1, y+1)

    pixscale = wcs.pixel_scale()
    apertures = np.array(apertures_arcsec) / pixscale

    # Now, do aperture photometry at these points in the coadd images
    for band in bands:
        imfn = survey.find_file('image',  brick=brickname, band=band)
        ivfn = survey.find_file('invvar', brick=brickname, band=band)
        if not (os.path.exists(imfn) and os.path.exists(ivfn)):
            continue
        coimg = fitsio.read(imfn)
        coiv = fitsio.read(ivfn)

        apflux = np.zeros((len(skyfibers), len(apertures)), np.float32)
        apiv   = np.zeros((len(skyfibers), len(apertures)), np.float32)
        skyfibers.set('apflux_%s' % band, apflux)
        skyfibers.set('apflux_ivar_%s' % band, apiv)
        with np.errstate(divide='ignore', invalid='ignore'):
            imsigma = 1./np.sqrt(coiv)
            imsigma[coiv == 0] = 0
        apxy = np.vstack((skyfibers.x, skyfibers.y)).T
        for irad,rad in enumerate(apertures):
            aper = photutils.CircularAperture(apxy, rad)
            p = photutils.aperture_photometry(coimg, aper, error=imsigma)
            apflux[:,irad] = p.field('aperture_sum')
            err = p.field('aperture_sum_err')
            apiv  [:,irad] = 1. / err**2

    header = fitsio.FITSHDR()
    for i,ap in enumerate(apertures_arcsec):
        header.add_record(dict(name='AP%i' % i, value=ap, comment='Aperture radius (arcsec)'))
    skyfibers._header = header

    return skyfibers

def sky_fiber_locations(skypix, gridsize=300):
    '''Worker function for sky_fibers_for_brick--- implements the core
    trick of iteratively eroding the map of good sky locations to
    produce a distance-from-blobs map, and then return the max values
    in that map in each cell of a grid.
    '''
    # Select possible locations for sky fibers
    from scipy.ndimage.morphology import binary_dilation, binary_erosion
    from scipy.ndimage.measurements import label, find_objects, center_of_mass
    from scipy.ndimage.filters import gaussian_filter

    nerosions = np.zeros(skypix.shape, np.int16)
    nerosions += skypix
    element = np.ones((3,3), bool)
    while True:
        skypix = binary_erosion(skypix, structure=element)
        nerosions += skypix
        print('After erosion,', np.sum(skypix), 'sky pixels')
        if not np.any(skypix.ravel()):
            break

    # This is a hack to break ties in the integer 'nerosions' map.
    nerosions = gaussian_filter(nerosions.astype(np.float32), 1.0)
    peaks = (nerosions > 1)
    H,W = skypix.shape

    # find pixels that are larger than their 8 neighbors
    peaks[1:-1, 1:-1] &= (nerosions[1:-1,1:-1] >= nerosions[0:-2,1:-1])
    peaks[1:-1, 1:-1] &= (nerosions[1:-1,1:-1] >= nerosions[2:  ,1:-1])
    peaks[1:-1, 1:-1] &= (nerosions[1:-1,1:-1] >= nerosions[1:-1,0:-2])
    peaks[1:-1, 1:-1] &= (nerosions[1:-1,1:-1] >= nerosions[1:-1,2:  ])
    peaks[1:-1, 1:-1] &= (nerosions[1:-1,1:-1] >= nerosions[0:-2,0:-2])
    peaks[1:-1, 1:-1] &= (nerosions[1:-1,1:-1] >= nerosions[0:-2,2:  ])
    peaks[1:-1, 1:-1] &= (nerosions[1:-1,1:-1] >= nerosions[2:  ,0:-2])
    peaks[1:-1, 1:-1] &= (nerosions[1:-1,1:-1] >= nerosions[2:  ,2:  ])

    # Split the image into 300 x 300-pixel cells, choose the highest peak in each one
    # (note, this is ignoring the brick-to-brick margin in laying down the grid)
    sx,sy = [],[]
    xx = np.round(np.linspace(0, W, 1+np.ceil(W / gridsize))).astype(int)
    yy = np.round(np.linspace(0, H, 1+np.ceil(H / gridsize))).astype(int)
    for ylo,yhi in zip(yy, yy[1:]):
        for xlo,xhi in zip(xx, xx[1:]):
            # Find max pixel in box
            subne = nerosions[ylo:yhi, xlo:xhi]
            I = np.argmax(subne)
            # Find all pixels equal to the max and take the one closest to the center of mass.
            maxval = subne.flat[I]
            cy,cx = center_of_mass(subne == maxval)
            xg = np.arange(xhi-xlo)
            yg = np.arange(yhi-ylo)
            dd = np.exp(-((yg[:,np.newaxis] - cy)**2 + (xg[np.newaxis,:] - cx)**2))
            dd[subne != maxval] = 0
            I = np.argmax(dd.flat)
            iy,ix = np.unravel_index(I, subne.shape)
            sx.append(ix + xlo)
            sy.append(iy + ylo)

    sx = np.array(sx)
    sy = np.array(sy)
    return sx, sy, nerosions[sy,sx]

def sky_fiber_plots(survey, brickname, skyfibers, basefn, bands=['g','r','z']):
    from legacypipe.survey import get_rgb
    import pylab as plt

    rgbkwargs = dict(mnmx=(-1,100.), arcsinh=1.)

    imgs = []
    for band in bands:
        fn = survey.find_file('image',  brick=brickname, band=band)
        imgs.append(fitsio.read(fn))
    rgb = get_rgb(imgs, bands, **rgbkwargs)

    ima = dict(interpolation='nearest', origin='lower')
{'fraction_non_alphanumeric': 0.146010059304857, 'fraction_numerical': 0.025373470460175664, 'mean_word_length': 2.723309111235327, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 21, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '487056', 'n_tokens_mistral': 5100, 'n_tokens_neox': 4849, 'n_words': 1119}