Clean simple FITS data

This won't be a lengthy discussion about methods to clean amateur astronomical FITS data; instead, since I've written about it before (albeit when I didn't know very much about Python), I'd rather post some new code I wrote to do the same process in a much clearer way.

#! /usr/bin/env python3.6
# -*- coding: utf-8 -*-

""" Reduce a quick batch of FITS data. Currently run with input file list
/home/brandon/Downloads/Astronomy_Data/V1283_Her/fixed_gain_V1283_Her_*.fit
/home/brandon/Downloads/Astronomy_Data/V1283_Her/Dark-60s_*.fit
/home/brandon/Downloads/Astronomy_Data/V1283_Her/Flat-V_*.fit
-o _Finished_Feb_08_2016
"""

from typing import Generator, Tuple
from astropy.io import fits as FITS
from astropy.io.fits import Header as HDU
from argparse import ArgumentParser
from glob import iglob
from datetime import datetime
import numpy as np

_EXT = 'fits'


def get_fits_file(name: str, axis: int =0) -> Tuple[np.ndarray, HDU]:
    """ Get both np.ndarray of data and the HDU (these are simple data) """
    with FITS.open(name) as image:
        file = image[axis]
        return file.data.astype(np.float32), file.header


def _get_median_dark_frame(darks: Generator[str, None, None]) -> np.ndarray:
    collective = list(zip(*map(get_fits_file, darks)))
    return np.median(collective[0], axis=0)


def _get_master_flat_frame(flats: Generator[str, None, None]) -> np.ndarray:
    collective = list(zip(*map(get_fits_file, flats)))

    # Get absolute mean of frames and scale them
    for frame in collective[0]:
        frame /= np.mean(frame)

    return np.median(collective[0], axis=0)


def clean(args) -> None:
    """ Process data, just a straightforward pipeline """
    lights, flats, darks = args.lights, args.flats, args.darks
    master_dark = _get_median_dark_frame(darks)
    master_flat = _get_master_flat_frame(flats)

    for light in lights:
        light_data, light_header = get_fits_file(light)
        light_data -= master_dark    # undo additive noise
        light_data /= master_flat    # undo convolution with PSF

        # Now write these data with the new name
        new_name = ''.join(light.split('.')[:-1]) + args.output + '.' + _EXT
        light_header['HISTORY'] = f'Dark subtracted and flat corrected {datetime.now()}'
        FITS.writeto(new_name, light_data, light_header, overwrite=True)


if __name__ == '__main__':
    parser = ArgumentParser(description=__doc__)
    parser.add_argument('-o', '--output', type=str,
        help='Base string of output file name')

    # Positional argument requiring a list of files
    parser.add_argument('lights', type=iglob,
        help='List of input files')

    # Another requiring darks
    parser.add_argument('darks', type=iglob,
        help='List of input dark frames')

    # Another requiring flats
    parser.add_argument('flats', type=iglob,
        help='List of input flat frames')

    args = parser.parse_args()
    del parser
    clean(args)
It creates master dark and flat frames with get_median_dark_frame and get_master_flat_frame, then computes the following for each light \(F_{i}\): $$ I_{i}=\dfrac{F_{i}-\mathscr{D}}{\mathscr{F}}. $$ I've written it as a command line tool with argparse so it can be called in a Bash script.