Reprojecting, stacking and colour-combination

Reprojecting, stacking and colour-combination#

On this page we’ll reproject our pre-reduced exposures on a common pixel grid, stack them, and produce a colour image.

You can copy or write the code shown below in a script, or alternatively directly download this page as a jupyter notebook file.

To run the code, you’ll need the module dataredconfig.py, as explained here.

import dataredconfig

import numpy as np
import astropy
import astropy.visualization
from astropy import units as u

%matplotlib ipympl
import matplotlib
from matplotlib import pyplot as plt
import ccdproc

light_prered_dir = dataredconfig.work_dir / "LIGHT_PRERED"
science_files = ccdproc.ImageFileCollection(light_prered_dir, keywords=dataredconfig.ifc_header_keywords)

# Let's first get an overview of all available files:
science_files.summary

# Where to write the reprojected files:
dest_dir = dataredconfig.work_dir / "REPROJ"
dest_dir.mkdir(exist_ok=True)

Rerun the following three cells with all filters:

Note that the following is a minimalistic example. It might give you an acceptable image, but it can certainly be improved a lot. In particular, the background might have changed during the observations. You could address this by subtracting the background individually from each exposure before combining them. Also, the flux scaling has likely changed (airmass, absorption). This could be addressed by computing some scaling factors for each image (based on photometry), and passing those to the ccdproc.combine function.

# Select object and filter:
selected_object = "M 37"
selected_filter = "r"

selected_science_files = science_files.filter(object=selected_object, filter=selected_filter)
selected_science_files.summary

# Select first image as the target to project on (note: same image for *all* filters!)
target_image = ccdproc.CCDData.read(science_files.filter(object=selected_object).files[0], unit='adu')

# Looping over the images to do the reprojection (takes a while):
for ccd, filename in selected_science_files.ccds(ccd_kwargs={'unit': 'adu'}, return_fname=True):
    print(f"Reprojecting {filename}...")
    
    ccd = ccdproc.wcs_project(ccd, target_image.wcs)

    # Write to disk:
    ccd.data = ccd.data.astype('float32')
    ccd.write(dest_dir / filename, overwrite=True)

# Combine the reprojected images of that filter:

files_to_combine = ccdproc.ImageFileCollection(dest_dir).files_filtered(object=selected_object, filter=selected_filter, include_path=True)
ccd = ccdproc.combine(files_to_combine,
                        method='average', scale=None,
                        sigma_clip=True, sigma_clip_low_thresh=5, sigma_clip_high_thresh=5,
                        sigma_clip_func=np.ma.median, signma_clip_dev_func=astropy.stats.mad_std,
                        mem_limit=1e9
                    )

ccd.meta['combined'] = True
ccd.data = ccd.data.astype('float32') # Converts to float32 to save space
ccd.write(dataredconfig.work_dir / f"combi_{selected_object}_{selected_filter}.fits", overwrite=True)

At this stage you might want to:

combine images of separate filters to get a color image (see below)
create a greyscale image (png file) from a single stack (see further below)

Colour-combination#

This is the colour-combination algorithm used by SDSS. It might not be “optimal” for pretty pictures (well that’s a matter of taste), but it is scientifically interesting.

Details can be found here: https://docs.astropy.org/en/stable/visualization/rgb.html#astropy-visualization-rgb

g_ccd = ccdproc.CCDData.read(dataredconfig.work_dir / f"combi_{selected_object}_g.fits")
r_ccd = ccdproc.CCDData.read(dataredconfig.work_dir / f"combi_{selected_object}_r.fits")
i_ccd = ccdproc.CCDData.read(dataredconfig.work_dir / f"combi_{selected_object}_i.fits")

# i -> R
i_ccd.data *= 1.0
# r -> G
r_ccd.data *= 0.4
# g -> B
g_ccd.data *= 0.3

sky_levels = (np.nanmedian(i_ccd.data), np.nanmedian(r_ccd.data), np.nanmedian(g_ccd.data))

rgbimage = astropy.visualization.make_lupton_rgb(i_ccd.data, r_ccd.data, g_ccd.data,
                                                 minimum=1.0*np.array(sky_levels),
                                                 stretch=30, Q=5, 
                                                 filename=dataredconfig.work_dir/f"combi_{selected_object}.jpg")

# Note that this both returns the image as a numpy array, and writes it to jpg.

fig = plt.figure(figsize=(10, 6))
ax = fig.subplots()
ax.imshow(rgbimage , origin="lower")
fig.tight_layout()
plt.show()

Creating a greyscale png file from a stacked image#

# Read the stack we want to use:
stackccd = ccdproc.CCDData.read(dataredconfig.work_dir / f"combi_{selected_object}_{selected_filter}.fits")

# It can be helpful so see a histogram:

fig = plt.figure(figsize=(6, 4))
ax = fig.subplots()
ax.hist(stackccd.data.flatten(), bins=1000, range=[-50, 10000])
fig.tight_layout()
plt.show()

# Adjusting cuts and scale parameters is best done on an interactive figure:

fig = plt.figure(figsize=(10, 6))
ax = fig.subplots()
ax.imshow(np.log10(stackccd.data - 1800.0),
          origin="lower",
          vmin = 1.4,
          vmax = 3.5,
          #norm = matplotlib.colors.LogNorm(vmin = 1600, vmax = 4000),
          cmap = "Greys_r",
         )
fig.tight_layout()
plt.show()

# The same parameters from the figure above can now be copied in imsave to write a full resolution png:

import matplotlib.image
matplotlib.image.imsave("pretty.png", 
        np.log10(stackccd.data - 1800.0),
        origin="lower",
        vmin = 1.4,
        vmax = 3.5,
        #norm = matplotlib.colors.LogNorm(vmin = 1600, vmax = 4000),
        cmap = "Greys_r",
        )

Reprojecting, stacking and colour-combination

Contents

Reprojecting, stacking and colour-combination#

Colour-combination#

Creating a greyscale png file from a stacked image#