Skip to content
This repository was archived by the owner on May 26, 2022. It is now read-only.
Open
Changes from 5 commits
Commits
Show all changes
28 commits
Select commit Hold shift + click to select a range
4669b09
adding ccdproc
crawfordsm Mar 5, 2013
8ecb546
update ccdproc api
crawfordsm Oct 7, 2013
eff0904
Whitespace/formatting changes for PEP8 compliance
mwcraig Jan 31, 2014
b743529
Fix imports/function namespace so that API is closer to ready to execute
mwcraig Jan 31, 2014
bccfe95
Merge pull request #1 from mwcraig/pep8
crawfordsm Feb 5, 2014
62b672c
Add API for a Keyword class used to access image metadata and a coupl…
mwcraig Feb 6, 2014
446a434
Merge pull request #2 from mwcraig/add_keywords
crawfordsm Feb 14, 2014
23a8721
Add logging proposals
mwcraig Feb 11, 2014
7e77d23
Incorporate suggestions by @crawfordsm
mwcraig Feb 15, 2014
fe818eb
Update user-specified keyword option to allow string or ccdproc.Keyword
mwcraig Feb 18, 2014
ac603c5
Fix assert statements for first two options and logging so they will …
mwcraig Feb 18, 2014
839022f
Initial attempt at image combination API
mwcraig Feb 12, 2014
3e2b0dd
Clarify handling of offsets/transforms
mwcraig Feb 15, 2014
f288e32
Perform image combination with a Combiner class
mwcraig Feb 18, 2014
9913a19
Merge pull request #4 from mwcraig/image-combine
crawfordsm Feb 18, 2014
d6fc09b
Update the reader information to use the astropy.io.registry
crawfordsm Mar 11, 2014
16b1878
Simplify behavior of Keyword class
mwcraig Apr 4, 2014
3b5b776
Clarify handling of units
mwcraig Apr 4, 2014
2a5fb64
Explicitly state that units are required for initialization
mwcraig Apr 9, 2014
d2292c3
Add detail to i/o portion of api
mwcraig Apr 9, 2014
0bdcb50
Pass read keywords to FITS reader
mwcraig Apr 9, 2014
8ff8175
User may not disable scaling when reading data from a FITS file.
mwcraig Apr 9, 2014
f31e5db
Describe method for converting CCDData object to an hdu.
mwcraig Apr 10, 2014
2ec4280
Simplify handling of keywords that are strings
mwcraig Apr 10, 2014
5171b3a
Fix typo
mwcraig Apr 11, 2014
d25b877
Merge pull request #5 from mwcraig/ccdproc
crawfordsm Apr 15, 2014
26e5493
simplified and clarified combiner
crawfordsm Apr 22, 2014
9cf829c
updates for current status of ccdproc
crawfordsm May 6, 2014
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
222 changes: 222 additions & 0 deletions ccdproc_api.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,222 @@
#!/usr/bin/python
# -*- coding: utf-8 -*-

from __future__ import print_function

import numpy as np

from astropy.nddata import NDData
from astropy.io import fits
from astropy import units as u

import ccdproc
'''
The ccdproc package provides tools for the reduction and
analysis of optical data captured with a CCD. The package
is built around the CCDData class, which has built into
it all of the functions to process data. The CCDData object
contains all the information to describe the 2-D readout
from a single amplifier/detector.

The CCDData class inherits from the NDData class as its base object
and the object on which all actions will be performed. By
inheriting from the CCD data class, basic array manipulation
and error handling are already built into the object.

The CCDData task should be able to properly deal with the
propogation of errors and propogate bad pixel frames
through each of the tasks. It should also update the meta
data, units, and WCS information as the data are processed
through each step.

The following functions are required for performing basic CCD correction:
-creation of variance frame
-overscan subtraction
-bias subtraction
-trimming the data
-gain correction
-xtalk correction
-dark frames correction
-flat field correction
-illumination correction
-fringe correction
-scattered light correction
-cosmic ray cleaning
-distortion correction

In addition to the CCDData and CCDList class, the ccdproc does
require some additional features in order to properly
reduce CCD data. The following features are required
for basic processing of CCD data:
-fitting data
-combining data
-re-sampling data
-transforming data

All actions of ccdproc should be logged and recorded.

Multi-Extension FITS files can be handled by treating
each extension as a CCDData object and

'''

# ============
# Base Objects
# ============
'''
CCDData is an object that inherits from NDData class and specifically
describes an object created from the single readout of a CCD.

Users should be able to create a CCDData object from scratch, from
an existing NDData object, or a single extension from a FITS file.

In the case of the CCDData, the parameter 'uncertainty' will
be mapped to variance as that will be more explicitly describing
the information that will be kept for the processing of the

'''
data = 100 + 10 * np.random.random((110, 100))
ccddata = ccdproc.CCDData(data=data)
ccddata = ccdproc.CCDData(NDData.NDData(data))
ccddata = ccdproc.CCDData(fits.ImageHDU(data))

#Setting basic properties of the object
# ----------------------
ccddata.variance = data**0.5
ccddata.mask = np.ones(110, 100)
ccddata.flags = np.zeros(110, 100)
ccddata.wcs = None
ccddata.meta = {}
ccddata.units = u.adu # is this valid?

#The ccddata class should have a functional form to create a CCDData
#object directory from a fits file
ccddata = ccdproc.CCDData.fromFITS('img.fits')

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Could this be called from_fits instead?

@astrofrog Might want to chime in here it but might also be better if it registered a FITS reader via the I/O registry; though it might not hurt to have an explicit from_fits method as well, which the I/O handler would simply wrap.

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@embray - I agree, if CCDData is an NDData sub-class, it should define it's own reader/writer and then use CCDData.read

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@astrofrog can you point me to an example of this so I know what needs to be added?

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@crawfordsm - yes, this is described here: http://docs.astropy.org/en/stable/io/registry.html

We're planning on refactoring how this is done a little, but it won't change too much.


# Functional Requirements
# ----------------------
# A number of these different fucntions are convenient functions that
# just outline the process that is needed. The important thing is that
# the process is being logged and that a clear process is being handled
# by each step to make building a pipeline easy. Then again, it might
# not be easy to handle all possible steps which are needed, and the more
# important steps will be things which aren't already handled by NDData.

#All functions should propogate throught to the variance frame and
#bad pixel mask

#convenience function based on a given value for
#the readnoise and gain. Units should be specified
#for these values though.
#Question: Do we have an object that is just a scalar

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This may be a dated question, but the answer is: astropy.units.Quantity - it's either a scalar + unit or an array + unit, depending on what you feed it at the time of creation. You definitely want to use those wherever it makes sense, as that's one of the ways a lot of astropy is "glued together"

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In a first bit of code that I was writing on this, that is what I did although I did have to create some of the units like adu.

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ah, ok, cool. There is an adu unit in astropy.units.astrophys, though (http://docs.astropy.org/en/stable/units/index.html#module-astropy.units.astrophys) - or did you add that for this purpose?

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Well, that's one less line of code that is needed now. Thanks for pointing that out!

#and a unit? Or a scalar, unit and an error? ie, This
#could actually be handled by the gain and readnoise being
#specified as an NDData object
ccddata = ccdproc.createvariance(ccddata, gain=1.0, readnoise=5.0)

#Overscan subtract the data
#Should be able to provide the meta data for
#the keyworkd or provide a section to define and
#possible an axis to specify the oritation of the
#Question: Best way to specify the section? Should it be given

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think what you have here might work, but another option might be a dictionary with (only) keywords xmin, xmax, ymin, ymax.

An (easy to implement) alternative would be to allow slice objects or, more generally, anything that you can do somearray[section] on. That would be trivial to implement, but would allow for situations like weird overscans that aren't rectangular or something (not sure if that exists in the wild, but you never know...).

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@eteq @crawfordsm -- I ended up implementing this as written, i.e. section is a string like "[:, 100:110]". This can, of course be changed, but there were a few reasons:

  • It should be clear to both new users of python and users of IRAF what "[:, 100:110]" means...something like a dict requires explanation and takes more typing.
  • To the extent that keywords like BIASSEC are used in FITS files they seem to have a format like this, and this approach something like section=ccddata.meta['BIASSEC'] to be used.
  • It turned out to be easy to implement a function to translate these strings into slice objects.

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm not attached to the exact approach I described above, @mwcraig, but I'm concerned about the "[:, 100:110]" approach for two reasons:

  1. I'm pretty sure the IRAF version of this is one-indexed, not zero-indexed, as python is. So people will certainly try to use it like IRAF when they should be assuming python, or vice vera, and this is certain to lead to massive confusion/headaches. (I dealt with a similar indexing problem just a few weeks ago, so it's fresh in my mind!)
  2. We really don't want to encourage the use of older idioms when one already exists in python. There are certainly situations where a user will want to pass in something that they used somewhere else on a numpy array, and we want that to work.

So here's an intermediate solution: accept both. That is, have the docstrings encourage the use of slice objects, and specifically address the zero-vs-one index problem. But have them accept the slice strings and interpret them in the IRAF fashion to be compatible with FITS keywords that use that convention.

How does that sound? Also, @mwcraig, when you say "I ended up", do you just mean you implemented this directly in ccdproc? Do you want some code review of that at some point?

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@eteq - in no particular order:

do you just mean you implemented this directly in ccdproc? Do you want some code review of that at some point?

Yes and yes; @crawfordsm and are making progress on ccdproc and will be talking on skype tomorrow. I'd guess we are 1-2 weeks from being ready for a review, with the understanding that some of what has been implemented will need to change.

I tend towards getting a functional draft out, knowing it will need revision, over trying to get all potential revisions worked before having something that works (hopefully that made sense).

specifically address the zero-vs-one index problem

shudder...it has been a few years since I've worked in a language that wasn't zero-indexed.

Points well taken, though -- it would be nice (though probably impossible) to have some way of catching these index errors, and documentation needs to emphasize this issue.

We really don't want to encourage the use of older idioms when one already exists in python....So here's an intermediate solution: accept both. That is, have the docstrings encourage the use of slice objects

Oddly, I see "[:, 100:110]" as more pythonic than (slice(None), slice(100, 110)) but there may be a better way of doing rectangular slices than that...I think that was one of the other things that led me down the string route. slice can only do one dimension, so the check would need to be for a tuple (or list?) of slice objects, each of which is a bit of a wordy to generate.

Another alternative is to go with the implementation @crawfordsm originally had for this:

def subtract_overscan(ccd, oscan)

where both ccd and oscan are CCDData objects, so the call becomes
ccddata = subtract_overscan(ccddata, ccddata[:, 100:100]). This has the virtue of being clear and compact!

It is maybe worth re-emphasizing that something is getting implemented so that there will be a concrete implementation to react to, with the expectation that the implementation will be re-worked, potentially substantially.

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@mwcraig - ok, now that I understand you're doing sort of a "prototype implementation", that makes sense. So do whatever you think best given the 1-index IRAF issue. Then we can revisit it once there's something concrete to look at.

Looking forward to it!

#Error Checks: That the section is within the image
ccddata = ccdproc.subtract_overscan(ccddata, section='[:,100:110]',
function='polynomial', order=3)

#trim the images--the section gives the part of the image to keep
#That the trim section is within the image
ccddata = ccdproc.trim_image(ccddata, section='[0:100,0:100]')

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I guess here you can't have the arbitrary indexing like numpy bool arrays or whatever (that I mentioned above) because you can't necessarily trim them. I suppose it could just mask those if its not a rectangular slice sort of thing, though.


#subtract the master bias. Although this is a convenience function as
#subtracting the two arrays will do the same thing. This should be able
#to handle logging of subtracting it off (logging should be added to NDData

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

by "logging", you just mean noting that it happened, right? That's what NDData.meta is for - you could just add a key 'bias_subtracted' and set it to True or the name of the file or something.

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Actually, that is probably the best way to handle the logging. These tasks should be pretty low level and different pipelines will have different desired types of logging. And so including all the information in the meta-data makes it easy to pull it out either for a log or for a later time.

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm coming into this quite late, but I wonder if there is any interest in adapting our pipeline framework we are developing for JWST to use for things like ccdproc. It would need some changes, but it's intended to handle things like logging, configuration, omitting or injecting new steps and stuff like that so that each calibration step doesn't have to deal much with that. Here's a link to the documentation:

http://ssb.stsci.edu/doc/jwst_dev/jwstlib.stpipe.doc/html/index.html

We've also been using data model classes for these steps that isolate the details of the files from the calibration code itself (the model class handles mapping the relevant information in the file into the data object). For general use, it would be up to the user to supply the appropriate keyword mapping for things needed by ccdproc (far fewer things needed for that than for JWST pipelines). Presumably the nddata object would serve this role with appropriate adapters for the data a user has to read the data in from FITS files.

If there is any interest, I can see about making the necessary changes to use for astropy.

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

On that last note, I would add at some point astropy.io.fits is going to start returning some kind of NDData object rather than just plain Numpy arrays. It will probably be a subclass specially designed to meet the special needs of FITS. But part of this interface might include ways to provide mappings from FITS headers to data models accessible through the NDData object's meta attribute. I hesitate to bring this up because I can't make any promises right now as to when that would happen, but it's just some kind of big picture stuff to keep in mind.

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hi Perry--Thanks for reminding me about this. I had taken a look at this previously and my feeling was stpipe was at a higher level than what we want to be developing for ccdproc. I believe we are going to follow something closer to @eteq suggestion of updating the meta-data with logging information and leave recording of the processing in a different format up to the user as everyone will have slightly different needs and preferences. It would be my hope that ccdproc would be able to be integrated into something like stpipe if ccdproc does prove useful and my hope is that how we are setting things up wouldn't prevent that from happening.

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ok, thanks. I just wanted to raise the issue. I'm happy to leave it the way it is.

#and then this is really just a convenience function
#Error checks: the masterbias and image are the same shape
masterbias = NDData.NDData(np.zeros(100, 100))
ccddata = ccdproc.subtract_bias(ccddata, masterbias)

#correct for dark frames
#Error checks: the masterbias and image are the same shape
masterdark = NDData.NDData(np.zeros(100, 100))
ccddata = ccdproc.subtract_dark(ccddata, masterdark)

#correct for gain--once again gain should have a unit and even an error
#associated with it.
ccddata = ccdproc.gain_correct(ccddata, gain=1.0)
#Also the gain may be non-linear
ccddata = ccdproc.gain_correct(ccddata, gain=np.array([1.0, 0.5e-3]))
#although then this step should be apply before any other corrections
#if it is non-linear, but that is more up to the person processing their
#own data.

#crosstalk corrections--also potential a convenience function, but basically
#multiples the xtalkimage by the coeffient and then subtracts it. It is kept
#general because this will be dependent on the CCD and the set up of the CCD.
#Not applicable for a single CCD situation
#Error checks: the xtalkimage and image are the same shape
xtalkimage = NDData.NDData(np.zeros(100, 100))
ccddata = ccdproc.xtalk_correct(ccddata, xtalkimage, coef=1e-3)

#flat field correction--this can either be a dome flat, sky flat, or an
#illumination corrected image. This step should normalize by the value of the
#flatfield after dividing by it.
#Error checks: the flatimage and image are the same shape
#Error checks: check for divive by zero
#Features: If the flat is less than minvalue, minvalue is used
flatimage = NDData.NDData(np.ones(100, 100))
ccddata = ccdproc.flat_correct(ccddata, flatimage, minvalue=1)

#fringe correction or any correction that requires subtracting
#off a potentially scaled image
#Error checks: the flatimage and image are the same shape
fringeimage = NDData.NDData(np.ones(100, 100))
ccddata = ccdproc.fringe_correct(ccddata, fringeimage, scale=1,
operation='multiple')

#cosmic ray cleaning step--this should have options for different
#ways to do it with their associated steps. We also might want to
#implement this as a slightly different step. The cosmic ray cleaning
#step should update the mask and flags. So the user could have options
#to replace the cosmic rays, only flag the cosmic rays, or flag and
#mask the cosmic rays, or all of the above.
ccddata = ccdproc.cosmicray_laplace(ccddata, method='laplace')
ccddata = ccdproc.cosmicray_median(ccddata, method='laplace')

#Apply distortion corrections
#Either update the WCS or transform the frame
ccddata = ccdproc.distortion_correct(ccddata, distortion)


# ================
# Helper Functions
# ================

#fit a 1-D function with iterative rejections and the ability to
#select different functions to fit.
#other options are reject parameters, number of iteractions
#and/or convergernce limit
coef = ccdproc.iterfit(x, y, function='polynomial', order=3)

#fit a 2-D function with iterative rejections and the ability to
#select different functions to fit.
#other options are reject parameters, number of iteractions
#and/or convergernce limit
coef = ccdproc.iterfit(data, function='polynomial', order=3)

#in addition to these operations, basic addition, subtraction
# multiplication, and division should work for CCDDATA objects
ccddata = ccddata + ccddata
ccddata2 = ccddata * 2

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yeah, definitely - you may know this already, but NDData has methods like add subtract multiply, etc. that are intended to make this easy for subclasses to implement.

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I was looking at that the other day -- I gather these examples should be written as:

ccddata = ccddata.add(ccddata, ccddata) 
ccddata2 = ccddata.multiply(2)

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Actually, no, I'd say they should be as written. To implement this, you'd put this in the CCDImage class:

def __add__(self, other):
    if isinstance(other, CCDImage):
        return self.add(other)
    else:
        raise TypeError("can't add a CCDImage to something that's not a CCDImage")

The logic behind this is that for the generic NDData, you do not want to make the assumption that standard operations make sense. For example, we certainly don't want ccdimage + datacube to work, because that would have to make assumptions about how to align the axes that may not be correct. But ccdimage + ccdimage is a sensible operation, so it should work.

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@eteq @crawfordsm -- FYI, I am holding off on implementing __mul__ and friends for the moment.

The problem isn't with how to do it -- it is easy enough to write the method to allow for multiplication by a scalar or another CCDData and propagate the error correctly in both cases...and in principle multiplication by a Quantity can be handled too.

However, until NDData and Quantity better coexist it is impossible to make sure that something like (1.5 * u.ph/u.adu) * ccd_data_object returns a CCDData object. In this case the fact that Quantity defines __mul__ means that its multiplication method gets called and the result is a Quantity. Multiplying in the opposite order returns a CCDData object, rendering multiplication non-commutative.

What I've done for now is overridden the multiply, etc., methods from NDData. They check for the case of multiplication by a scalar or Quantity and handle that case, with multiplication by another CCDData passed off to super.

Just to be clear, I am not at all opposed to implementing __mul__ and friends -- in fact, I would prefer it and went down that road first! But for now I can't control how __mul__ behaves.

Actually, I could, I suppose, by removing the __array__ method from CCDData. Then Quantity's __mul__ ought to fail.

Speak up if that would be the preferred behavior -- doing so would not, I think, interfere with slicing, which is implemented in NDData.__getitem__.



#combine a set of NDData objects
alldata = ccdproc.combine([ccddata, ccddata2], method='average',
reject=None)

#re-sample the data to different binnings (either larger or smaller)
ccddata = ccdproc.rebin(ccddata, binning=(2, 2))

#tranform the data--ie shift, rotate, etc
#question--add convenience functions for image shifting and rotation?
#should udpate WCS although that would actually be the prefered method
ccddata = ccdproc.transform(ccddata, transform, conserve_flux=True)