Making Galaxy Zoo: Hubble images for GOODS full-depth imaging (snippet)

make_color_bviz_jpeg.py

def make_jpeg(gal,color_scheme='bviz',field='s',desaturate=False,show_img=True,load_from_mosaic=False,mosaics=None):

    ascales,anonlinearity,npix_final,apedestal = get_image_defaults(color_scheme=color_scheme)

    input_oid = gal['UID_MOSAIC_Z']
    id_str = 'g%s' % input_oid

    #print '----------------- %s --------------------' % id_str

    # Get the amount by which we need to resize the image
    # Consider all bands
    a_pix = max(gal['A_IMAGE_Z'], gal['A_IMAGE_I'], gal['A_IMAGE_V'], gal['A_IMAGE_B'])
    kron_r = max(gal['KRON_RADIUS_Z'], gal['KRON_RADIUS_I'], gal['KRON_RADIUS_V'], gal['KRON_RADIUS_B'])
    obj_size_pix = a_pix * kron_r

    # Scale in HST pixel size
    img_size_hstpix = 3.5 * obj_size_pix

    # Sanity checks -- don't zoom in or out too far
    img_size_hstpix = max(img_size_hstpix, 120.)
    img_size_hstpix = min(img_size_hstpix, 1000.)
    
    resize_factor = npix_final / img_size_hstpix

    # Print re-sizing parameters to screen
    print '%s -- Semi-major axis [pix]: %f\n Kron radius [pix]: %f\n Image size [pix]: %f\n Resize factor: %f' % (id_str,a_pix, kron_r*a_pix, img_size_hstpix, resize_factor)

    # Load FITS data from original GOODS mosaics

    if load_from_mosaic:
        if mosaics is None:
            osect = gal['OSECT_Z']
            # Note: this will work w/gzipped mosaic files, but is MUCH slower.
            with fits.open('%s/%s%s_v1.0_sc03_osect%i_drz.fits' % (mosaic_path,field,'b',osect)) as f:      
                img_b = np.transpose(f[0].data)
            with fits.open('%s/%s%s_v1.0_sc03_osect%i_drz.fits' % (mosaic_path,field,'v',osect)) as f:
                img_v = np.transpose(f[0].data)
            with fits.open('%s/%s%s_v1.0_sc03_osect%i_drz.fits' % (mosaic_path,field,'i',osect)) as f:
                img_i = np.transpose(f[0].data)
            with fits.open('%s/%s%s_v1.0_sc03_osect%i_drz.fits' % (mosaic_path,field,'z',osect)) as f:
                img_z = np.transpose(f[0].data)
        else:
            img_b,img_v,img_i,img_z = mosaics

        assert img_b.shape == img_v.shape == img_i.shape == img_z.shape, \
            'Mosaic images must be the same shape\n B:%s, V:%s, I:%s, Z:%s' % \
            (img_b.shape,img_v.shape,img_i.shape,img_z.shape)

        # Cut down the image to the Kron radius aperture

        xsize,ysize = img_v.shape
        xstart,ystart = int(round(gal['X_IMAGE_B'])),int(round(gal['Y_IMAGE_B']))
        halfbox_orig = np.floor(img_size_hstpix/2.).astype(int)

        halfbox = check_edge_size(xsize,ysize,xstart,ystart,halfbox_orig)

        img_b_cut = img_b[xstart-halfbox:xstart+halfbox,ystart-halfbox:ystart+halfbox]
        img_v_cut = img_v[xstart-halfbox:xstart+halfbox,ystart-halfbox:ystart+halfbox]
        img_i_cut = img_i[xstart-halfbox:xstart+halfbox,ystart-halfbox:ystart+halfbox]
        img_z_cut = img_z[xstart-halfbox:xstart+halfbox,ystart-halfbox:ystart+halfbox]

    # Load FITS data from existing cutouts

    else:
        with fits.open('%s/%s_s%s_thumb.fits' % (fits_path,id_str,'b')) as f:
            img_b_cut = f[0].data
        with fits.open('%s/%s_s%s_thumb.fits' % (fits_path,id_str,'v')) as f:
            img_v_cut = f[0].data
        with fits.open('%s/%s_s%s_thumb.fits' % (fits_path,id_str,'i')) as f:
            img_i_cut = f[0].data
        with fits.open('%s/%s_s%s_thumb.fits' % (fits_path,id_str,'z')) as f:
            img_z_cut = f[0].data

    # Check that all image sizes and shapes match

    assert img_b_cut.shape == img_v_cut.shape == img_i_cut.shape == img_z_cut.shape, \
        'Cutout images must be the same shape\n B:%s, V:%s, I:%s, Z:%s' % \
        (img_b_cut.shape,img_v_cut.shape,img_i_cut.shape,img_z_cut.shape)

    if color_scheme == 'bviz':
        rimage = img_z_cut
        gimage = img_i_cut
        bimage = np.array([img_b_cut,img_v_cut]).mean(axis=0)
    elif color_scheme == 'vz':
        rimage = img_z_cut
        gimage = np.array([img_z_cut,img_v_cut]).mean(axis=0)
        bimage = img_v_cut

    nx,ny = rimage.shape

    if load_from_mosaic:
        RGBim = np.array([np.transpose(rimage),np.transpose(gimage),np.transpose(bimage)])
    else:
        RGBim = np.array([rimage,gimage,bimage])
    
    # Scale and fit the image data
    RGBim = nw.scale_rgb(RGBim,scales=ascales)
    RGBim = nw.arcsinh_fit(RGBim,nonlinearity=anonlinearity)
    RGBim = nw.fit_to_box(RGBim)

    if desaturate:
        # optionally desaturate pixels that are dominated by a single
        # colour to avoid colourful speckled sky

        orig = deepcopy(RGBim)
        # Take the mean flux value between the three color bands
        a = RGBim.mean(axis=0)
        # mask pixels with no flux in any of the bands
        np.putmask(a, a == 0.0, 1.0)
        # create cube with each plane containing mean flux value
        acube = np.resize(a,(3,ny,nx))
        # scale each color to the mean, then divide by non-linearity factor
        bcube = (RGBim / acube) / anonlinearity
        # create a mask based on the weighted non-linear color values
        mask = np.array(bcube)
        # find the maximum weighted value per color
        w = np.max(mask,axis=0)
        # if the max value is greater than 1, replace with 1
        np.putmask(w, w > 1.0, 1.0)
        # invert mapping from 0 to 1.
        w = 1 - w
        # taper the weights with a sin function
        w = np.sin(w*np.pi/2.0)
        # multiply the original image by the normalized weights, add back the weighted mean flux, 
        # and [optionally] recolor strongest pixels
        RGBim = RGBim * w + a*(1-w) + a*(1-w)**2 * orig

    # Convert data to scaled bytes
    RGBim = (255.*RGBim).astype(int)
    RGBim = np.where(RGBim>255,255,RGBim)
    RGBim = np.where(RGBim<0,0,RGBim)

    # Add optional grey pedestal to the byte-scaled data

    RGBim += apedestal
    RGBim = np.where(RGBim>255,255,RGBim)

    R = RGBim[0,:,:]
    G = RGBim[1,:,:]
    B = RGBim[2,:,:]

    data = np.array([R.ravel(),G.ravel(),B.ravel()])
    data = np.transpose(data)
    pdata = []

    for each in data: 
        pdata.append(tuple(each))

    # Make Image in PIL format

    img = I.new('RGB',size=R.shape)
    img.putdata(pdata)

    # Rebin images to 424x424 pixels
    img_resized = img.resize((424,424),I.ANTIALIAS)

    # Reset orientation to match the native FITS (N up, E right)
    img_flipped = img_resized.transpose(I.FLIP_TOP_BOTTOM)

    img_out = img_flipped

    if show_img:
        img_out.show()

    # Save hardcopy as JPG
    out_jpg  = '%s/goods_%s_%s_%s_standard.jpg'  % (jpg_path,field,id_str,color_scheme)
    img_out.save(out_jpg,format='JPEG',quality=100)

    return img,pdata

def get_image_defaults(color_scheme='bviz'):

    if color_scheme == 'bviz':
        ascales= 1.75*np.array([45, 22, 28])
    
    if color_scheme == 'vz':
        ascales= 2.5*np.array([16.,11.,11.])
        #ascales= 2.0*np.array([15,10,8])
        #ascales= 1.5*np.array([26,22,21])

    anonlinearity= 2.5
    
    # number of pixels the final image should have (in x and y, each)
    npix_final = 424

    # Pedestal
    apedestal = 0

    return ascales,anonlinearity,npix_final,apedestal