Skip to content

Commit

Permalink
added dark and fringe correction to joinup
Browse files Browse the repository at this point in the history
  • Loading branch information
@trmrsh
trmrsh committed May 26, 2021
1 parent b6f7bc0 commit 688440c
Show file tree
Hide file tree
Showing 3 changed files with 161 additions and 58 deletions.
2 changes: 1 addition & 1 deletion hipercam/scripts/averun.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -306,7 +306,7 @@ def averun(args=None):
resource = hcam.scripts.grab(args)

# at this point 'resource' is a list of files, no matter the input
# method.
# method.

with CleanUp(
resource, server_or_local or bias is not None or dark is not None
Expand Down
179 changes: 141 additions & 38 deletions hipercam/scripts/joinup.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,7 @@
from astropy.io import fits

import hipercam as hcam
from hipercam import cline, utils, spooler, defect
from hipercam import cline, utils, spooler, defect, fringe
from hipercam.cline import Cline

__all__ = [
Expand All @@ -22,7 +22,8 @@

def joinup(args=None):
"""``joinup [source] (run first [twait tmax] | flist) trim ([ncol
nrow]) (ccd) bias flat msub dtype dmax nmax overwrite compress``
nrow]) (ccd) bias dark flat fmap (fpair nhalf rmin rmax) msub
dtype dmax nmax overwrite compress``
Converts a run or a list of hcm images into as near as possible
"standard" FITS files with one image in the primary HDU per file,
Expand Down Expand Up @@ -56,23 +57,26 @@ def joinup(args=None):
run : string [if source ends 's' or 'l']
run number to access, e.g. 'run034'
flist : string [if source ends 'f']
name of file list
first : int [if source ends 's' or 'l']
exposure number to start from. 1 = first frame. For data
from the |hiper| server, a negative number tries to get a frame not
quite at the end. i.e. -10 will try to get 10 from the last
frame. This is mainly to sidestep a difficult bug with the
acquisition system.
last : int [if source ends 's' or 'l']
Last frame to access, 0 for the lot
twait : float [if source ends 's' or 'l'; hidden]
time to wait between attempts to find a new exposure, seconds.
tmax : float [if source ends 's' or 'l'; hidden]
maximum time to wait between attempts to find a new exposure,
seconds.
flist : string [if source ends 'f']
name of file list
trim : bool [if source starts with 'u']
True to trim columns and/or rows off the edges of windows nearest
the readout which can sometimes contain bad data.
Expand All @@ -90,69 +94,97 @@ def joinup(args=None):
bias : str
Name of bias frame to subtract, 'none' to ignore.
dark : str
Name of dark frame to correct for dark counts. 'none' to
ignore.
flat : str
Name of flat field to divide by, 'none' to ignore. Should normally
only be used in conjunction with a bias, although it does allow you
to specify a flat even if you haven't specified a bias.
fmap : str
Fringe map to remove fringes, 'none' to ignore.
fpair : str [if fmap is not == 'none']
File of peak/trough pairs for fringe amplitude measurement.
nhalf : int [if fmap is not == 'none', hidden]
Half-size of regions around each point for measuring intensity
differences from peak/trough pairs.
rmin : float [if fmap is not == 'none', hidden]
Minimum individual ratio for pruning peak/tough ratios prior to
taking their median.
rmin : float [if fmap is not == 'none', hidden]
Maximum individual ratio for pruning peak/trough ratios prior to
taking their median.
msub : bool
subtract the median from each window. If set this happens after any
bias subtraction.
bias subtraction etc.
ndigit : int
number of digits to be used in the frame count attached to the
output file names.
number of digits to be used in the frame counter attached
to the output file names. These are zero-padded so that the
frames order alphabetically. Thus 'run0002_ccd1_0001.fits',
'run0002_ccd1_0002.fits', 'run0002_ccd1_0003.fits' ... for
instance.
dtype : str
output data type. 'unit16', 'float32', 'float64'. The first
of these (2-byte unsigned) is only probably a good idea if
no bias, flat field or median subtraction has been applied because
it involves rounding and it will fail if any data are out of the range
0 to 65535. 32-bit (4 byte) floats should be OK for most purposes,
but require twice the space of uint16.
no bias, flat field or median subtraction has been applied
because it involves rounding and it will fail if any data
are out of the range 0 to 65535. 32-bit (4 byte) floats
should be OK for most purposes, but require twice the space
of uint16.
dmax : float
Maximum amount of data in GB to write out. Just a safety device
to avoid disaster by applying this script to highly windowed data
where you can end up expanding the total amount of data by a large
factor.
Maximum amount of data in GB to write out. A safety device
to avoid disaster in case this script was applied to highly
windowed data where you can end up expanding the total
amount of data by a large factor.
nmax : int
Maximum number of frames. A similar safety device to dmax to
avoid inadvertent application of this script to a million+ frame run.
Maximum number of frames. A similar safety device to dmax
to avoid inadvertent application of this script to a
million+ frame run. File systems tend not to behave well
with vast numbers of files.
overwrite : bool
overwrite any pre-existing files.
compress : str
allows data to be compressed with FITS's internal lossless
compression mechanisms. The file will still end in ".fits"
compression mechanisms. The file will still end as ".fits"
but has a different internal format; 'ds9' copes seamlessly
with all of them. The options are: 'none', 'rice', 'gzip1',
'gzip2'. 'rice' gave about a factor of 2 compression in a
short test I ran, and was as fast as gzip2, but it may
depend upon the nature of the data. 'none' is fastest. See
depend upon the nature of the data. 'none' is fastest. See
astropy.io.fits for further documentation.
.. Note::
Be careful of running this on highly windowed data since it
could end up expanding the total amount of "data" hugely.
It's really aimed at full frame runs above all. The "dmax"
and "nmax" parameters are aimed at heading off disaster.
could end up expanding the total amount of "data" hugely. It's
really aimed at full frame runs above all. The "dmax" and
"nmax" parameters are aimed at heading off disaster.
This routine will fail if windows have been binned but are out
of step (not "in sync") with each other because there is no
way to register such data within a single image.
The routine only creates a window big enough to contain all
the windows. Thus it might end up representing a sub-array of
the CCD as opposed to all of it. The location can be
determined from the 'LLX' and 'LLY' parameters that are
written to the header. These represent the location of the
lowest and leftmost unbinned pixel that is contained within
the data array. The bottom-left pixel of the CCD is considered
to be (1,1), so full frame images have LLX=LLY=1.
of step (not "in sync") with each other because there is no way
to register such data within a single image.
The routine only creates a window big enough to contain all the
windows. Thus it might end up representing a sub-array of the
CCD as opposed to all of it. The location can be determined
from the 'LLX' and 'LLY' parameters that are written to the
header. These represent the location of the lowest and leftmost
unbinned pixel that is contained within the data array. The
bottom-left pixel of the CCD is considered to be (1,1), so full
frame images have LLX=LLY=1.
A HipercamError will be raised if an attempt is made to write out
data outside the range 0 to 65535 into unit16 format and nothing
Expand All @@ -169,6 +201,7 @@ def joinup(args=None):
cl.register("source", Cline.GLOBAL, Cline.HIDE)
cl.register("run", Cline.GLOBAL, Cline.PROMPT)
cl.register("first", Cline.LOCAL, Cline.PROMPT)
cl.register("last", Cline.LOCAL, Cline.PROMPT)
cl.register("trim", Cline.GLOBAL, Cline.PROMPT)
cl.register("ncol", Cline.GLOBAL, Cline.HIDE)
cl.register("nrow", Cline.GLOBAL, Cline.HIDE)
Expand All @@ -177,7 +210,13 @@ def joinup(args=None):
cl.register("flist", Cline.LOCAL, Cline.PROMPT)
cl.register("ccd", Cline.LOCAL, Cline.PROMPT)
cl.register("bias", Cline.GLOBAL, Cline.PROMPT)
cl.register("dark", Cline.GLOBAL, Cline.PROMPT)
cl.register("flat", Cline.GLOBAL, Cline.PROMPT)
cl.register("fmap", Cline.GLOBAL, Cline.PROMPT)
cl.register("fpair", Cline.GLOBAL, Cline.PROMPT)
cl.register("nhalf", Cline.GLOBAL, Cline.HIDE)
cl.register("rmin", Cline.GLOBAL, Cline.HIDE)
cl.register("rmax", Cline.GLOBAL, Cline.HIDE)
cl.register("msub", Cline.GLOBAL, Cline.PROMPT)
cl.register("ndigit", Cline.LOCAL, Cline.PROMPT)
cl.register("dtype", Cline.LOCAL, Cline.PROMPT)
Expand Down Expand Up @@ -205,6 +244,11 @@ def joinup(args=None):
else:
first = cl.get_value("first", "first frame to process", 1, 1)

last = cl.get_value("last", "last frame to grab", 0)
if last and last < first:
sys.stderr.write("last must be >= first or 0")
sys.exit(1)

twait = cl.get_value(
"twait", "time to wait for a new frame [secs]", 1.0, 0.0
)
Expand Down Expand Up @@ -253,6 +297,15 @@ def joinup(args=None):
else:
fprompt = "flat frame ['none' is normal choice with no bias]"

# dark
dark = cl.get_value(
"dark", "dark frame ['none' to ignore]",
cline.Fname("dark", hcam.HCAM), ignore="none"
)
if dark is not None:
# read the dark frame
dark = hcam.MCCD.read(dark)

# flat (if any)
flat = cl.get_value(
"flat", fprompt, cline.Fname("flat", hcam.HCAM), ignore="none"
Expand All @@ -261,6 +314,33 @@ def joinup(args=None):
# read the flat frame
flat = hcam.MCCD.read(flat)

# fringe file (if any)
fmap = cl.get_value(
"fmap",
"fringe map ['none' to ignore]",
cline.Fname("fmap", hcam.HCAM),
ignore="none",
)
if fmap is not None:
# read the fringe map
fmap = hcam.MCCD.read(fmap)
fpair = cl.get_value(
"fpair", "fringe pair file",
cline.Fname("fpair", hcam.FRNG)
)
fpair = fringe.MccdFringePair.read(fpair)

nhalf = cl.get_value(
"nhalf", "half-size of fringe measurement regions",
2, 0
)
rmin = cl.get_value(
"rmin", "minimum fringe pair ratio", -0.2
)
rmax = cl.get_value(
"rmax", "maximum fringe pair ratio", 1.0
)

msub = cl.get_value("msub", "subtract median from each window?", True)
ndigit = cl.get_value(
"ndigit", "number of digits to use for frame numbers in output names",
Expand Down Expand Up @@ -307,7 +387,7 @@ def joinup(args=None):
with spooler.data_source(source, resource, first) as spool:

# 'spool' is an iterable source of MCCDs
nframe = 0
nframe = first
for mccd in spool:

if server_or_local:
Expand All @@ -331,18 +411,25 @@ def joinup(args=None):

# indicate progress
tstamp = Time(mccd.head["TIMSTAMP"], format="isot", precision=3)
nf = mccd.head.get("NFRAME",nframe+1)
nf = mccd.head.get("NFRAME",nframe)
tok = "ok" if mccd.head.get("GOODTIME", True) else "nok"
print(f"{nf}, utc= {tstamp.iso} ({tok})")

if nframe == 0:
if nframe == first:
if bias is not None:
# crop the bias on the first frame only
bias = bias.crop(mccd)

bexpose = bias.head.get("EXPTIME", 0.0)
else:
bexpose = 0.
if dark is not None:
dark = dark.crop(mccd)
if flat is not None:
# crop the flat on the first frame only
flat = flat.crop(mccd)
if fmap is not None:
fmap = fmap.crop(mccd)
fpair = fpair.crop(mccd, nhalf)

# OK here goes
for nc, cnam in enumerate(ccds):
Expand All @@ -356,10 +443,24 @@ def joinup(args=None):
if bias is not None:
ccd -= bias[cnam]

if dark is not None:
# subtract dark, CCD by CCD
dexpose = dark.head["EXPTIME"]
cexpose = ccd.head["EXPTIME"]
scale = (cexpose - bexpose) / dexpose
ccd -= scale * dark[cnam]

# divide out the flat
if flat is not None:
ccd /= flat[cnam]

if fmap is not None:
if cnam in fmap and cnam in fpair:
fscale = fpair[cnam].scale(
ccd, fmap[cnam], nhalf, rmin, rmax
)
ccd -= fscale*fmap[cnam]

if msub:
# subtract median from each window
for wind in ccd.values():
Expand Down Expand Up @@ -454,5 +555,7 @@ def joinup(args=None):

# update the frame number
nframe += 1
if server_or_local and last and nframe > last:
break

print(f'Written {nfile} FITS files to disk.')
38 changes: 19 additions & 19 deletions hipercam/scripts/makefringe.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -33,28 +33,28 @@ def makefringe(args=None):
patterns you see when there is a thin layer of oil on water. Both
patterns are caused by interference. In the case of CCDs this is
caused by sky emission lines and is variable in strength, and
additive in nature.
additive in nature, but seems mostly fairly fixed in position.
De-fringing in |hiper| is implemented according to the method
suggested by Snodgrass & Carry (2013Msngr.152...14S). The idea is
to create a set of pairs of points marking the peak and troughs of
fringes. The difference in intensity between these in the data to
be corrected is compared with the difference in a reference
"fringe map" by taking their ratio. This is done many times so
that the median can be taken to eliminate ratios affected by
celestial targets or cosmic rays. the ratios of which will be used
to estimate the level of fringing compared to a reference
frame. ``makefringe`` is used to make the reference fringe map. It
works as follows: given an input list of files (or optionally a
single run), it reads them all in, optionally debiases,
dark-subtracts and flat-fields them, calculates a median count
level of each one which is subtracted from each CCD
individually. The pixel-by-pixel median of all frames is then
calculated. It is also possible to apply the fringe pair ratio
measurement to scale each frame before combining which allows
frames of similar pattern but varying amplitude to be
combined. Finally, the output can be smoothed because fringes are
usually a medium to large scale pattern.
fringes. The differences in intensity between these point pairs in
the data to be corrected is compared with the differences for the
same pairs in a reference "fringe map" by taking their
ratios. Many pairs are used so that the median can be taken to
eliminate ratios affected by celestial targets or cosmic rays. The
median ratio is used to scale the reference data and thereby
subtract the fringes. ``makefringe`` is used to make the reference
fringe map. It works as follows: given a single run or a list of
files, it reads them all in, optionally debiases, dark-subtracts
and flat-fields them, calculates a median count level of each one
which is subtracted from each CCD individually. The pixel-by-pixel
median of all frames is then calculated. It is also possible (and
a good idea) to apply the fringe pair ratio measurements to scale
each frame before combining them, which allows frames of similar
pattern but varying amplitude to be combined. Finally, the output
can be smoothed because fringes are usually a medium- to
large-scale pattern.
Parameters:
Expand Down Expand Up @@ -123,7 +123,7 @@ def makefringe(args=None):
output : str
output file. Set by default to match the last part of "run"
(but it will have a different extension so they won't
clash)
clash). The default is
.. Note::
Expand Down

0 comments on commit 688440c

Please sign in to comment.