Commit 2870ec84 authored by karpov-sv's avatar karpov-sv
Browse files

Documentation for pre-processing images, including stacking them. Also, simple...

Documentation for pre-processing images, including stacking them. Also, simple example notebook for image stacking.
parent 10401938
......@@ -52,6 +52,7 @@ Apart of Python requirements that will be installed automatically, *STDPipe* als
- [SExtractor](
- [SCAMP](
- [PSFEx](
- [SWarp](
- [Astrometry.Net](
......@@ -35,7 +35,7 @@ extensions = [
# "recommonmark",
# "sphinx_math_dollar",
# "sphinx.ext.autodoc",
# 'sphinx.ext.pngmath',
......@@ -113,7 +113,7 @@ exclude_patterns = ["_build"]
# The reST default role (used for this markup: `text`) to use for all
# documents.
# default_role = None
default_role = 'obj'
# If true, '()' will be appended to :func: etc. cross-reference text.
# add_function_parentheses = True
Pre-processing the data
The codes in *STDPipe* expect as an input the *science-ready* images, cleaned as much as possible from instrumental signatures and imaging artefacts. In practice, it means that the image should be
- bias and dark subtracted
- flat-fielded.
Also, the artefacts such as saturated stars, bleeding charges, cosmic ray hits etc have to be masked.
All these tasks are outside of *STDPipe* per se, as they are highly instrument and site specific.
On the other hand, they may usually be performed easily using standard Python/NumPy/AstroPy routines and libraries like `Astro-SCRAPPY <>`_ etc.
E.g. to pre-process the raw image using pre-computed master dark and flat, mask some common problems, and cleanup the cosmic rays you may do something like that:
.. code-block:: python
image = fits.getdata(filename).astype(np.double)
header = fits.getheader(filename)
dark = fits.getdata(darkname)
flat = fits.getdata(flatname)
image -= dark
image *= np.median(flat)/flat
saturation = header.get('SATURATE') or 50000 # Guess saturation level from FITS header
mask = np.isnan(image) # mask NaNs in the input image
mask |= image > saturation # mask saturated pixels
mask |= flat < 0.5*np.median(flat) # mask underilluminated/vignetted regions
from astropy.stats import mad_std
mask |= dark > np.median(dark) + 10.0*mad_std(dark) # mask hotter pixels
gain = header.get('GAIN') or 1.0 # Guess gain from FITS header
import astroscrappy
# mask cosmic rays using LACosmic algorithm
cmask,cimage = astroscrappy.detect_cosmics(image, mask, gain=gain, verbose=True)
mask |= cmask
Stacking the images
You may want to stack/coadd or mosaic some images before processing them. While there are dedicated large-scape packages like `Montage <>`_ that handle it *properly*, it still may be done manually with relatively little efforts using e.g. Python `reproject <>`_ package.
You may check `simple example notebook <>`_
that shows how to do it.
*STDPipe* also contains a simple wrapper for `SWarp <>`_ re-projection code, :func:`stdpipe.templates.reproject_swarp`, that is implemented to resemble the calling conventions of `reproject <>`_ package routines as much as possible - i.e. allows directly stacking the image files without loading them to memory first:
.. autofunction:: stdpipe.templates.reproject_swarp
Using STDPipe
The codes in *STDPipe* expect pre-processed, or *science-ready* images, with instrumental signatures removed or masked as much as possible.
*STDPipe* is a library of routines that operate on standard Python objects, especially the ones from Numpy (two-dimensional arrays for images) and Astropy (tables, WCS solutions, FITS files, ...). Therefore, in order to use it, you will need to import all necessary packages first. E.g. first lines of your script may look like that:
.. code-block:: python
# Matplotlib is for plotting!
import matplotlib.pyplot as plt
# NumPy for data arrays
import numpy as np
Object detection
# AstroPy WCS and FITS support
from astropy.wcs import WCS
from import fits
# Support for sky coordinates and times from AstroPy
from astropy.coordinates import SkyCoord
from astropy.time import Time
Astrometric calibration
# Disable some annoying warnings from AstroPy
import warnings
from astropy.wcs import FITSFixedWarning
warnings.simplefilter(action='ignore', category=FITSFixedWarning)
from astropy.utils.exceptions import AstropyUserWarning
warnings.simplefilter(action='ignore', category=AstropyUserWarning)
Photometric calibration
Next, you will need to import the modules from *STDPipe* itself:
Differential imaging
.. code-block:: python
# Load (most of) our sub-modules
from stdpipe import astrometry, photometry, catalogs, cutouts, templates, subtraction, plots, psf, pipeline, utils
Now you have everything imported, and may start actual data analysis!
Data processing
.. toctree::
:maxdepth: 2
Common principles
The functions included in *STDPipe* try to follow several common principles related to their behaviour and arguments. They are summarized below.
- Most of functions accept ``verbose`` argument that controls the amount of informational outputs the function produces. You may use ``verbose=True`` to see the details of what exactly the function is doing. Also, you may pass any ``print``-like function to this argument to receive the messages instead of printing - so e.g. they may be logged.
.. code-block:: python
# Define simple logging function
def print_to_file(*args, logname='/tmp/logfile', clear=False, **kwargs):
if clear and os.path.exists(logname):
print('Clearing', logname)
if len(args) or len(kwargs):
print(*args, **kwargs)
with open(logname, 'a+') as lfd:
print(file=lfd, *args, **kwargs)
# verbose = True
verbose = print_to_file
# Now use it to redirect STDPipe function output to log file
obj = photometry.get_objects_sextractor(image, mask=mask, verbose=verbose)
- Functions that produce (and then delete of course) some temporary files during its operation usually accept ``_tmpdir`` argument to manually specify the location where these temporary files (usually in a dedicated per-task temporary folder, so thread-safe and stuff) will be stored. This is useful if your system-wide temporary directory (usually :file:`/tmp` in Linux) is low on free space - then you may use some larger volume for storing temporary files by adding ``_tmpdir=/large/volume/tmp`` to function call.
- Functions that operate on temporary files in temporary folders may be supplied with ``_workdir`` argument - then they will store all temporary files related to their work in this specific folder, and will not remove them afterwards. So e.g. you will be able to directly see e.g what configuration files were created, manually re-run the failed command to experiment with its options (with ``verbose=True`` function call typically prints the complete command line of all calls to external programs, so you may just directly copy and paste it to terminal to repeat its invocation), etc.
- Functions that run external programs (e.g. SExtractor, HOTPANTS, or Astrometry.Net) usually accept ``_exe`` argument to directly specify the path to corresponding executable. If not specified, the code will try to automatically find it for you, so normally you do not need to worry about it.
This diff is collapsed.
......@@ -415,6 +415,10 @@ def reproject_swarp(input=[], wcs=None, shape=None, width=None, height=None, hea
If `use_nans=True`, the regions with zero weights will be filled with NaNs (or 0xFFFF).
Any additional configuration parameter may be passed to SWarp through `extra` argument which
should be the dictionary with parameter names as keys.
# Simple wrapper around print for logging in verbose mode only
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment