FITSImageStdHeaders

class snappl.image.FITSImageStdHeaders(*args, zeropoint=None, is_superclass=False, header_kws={'band': 'BAND', 'dec': 'DEC', 'exptime': 'EXPTIME', 'mjd': 'MJD', 'observation_id': 'POINTING', 'position_angle': 'POSANG', 'ra': 'RA', 'sca': 'SCA', 'sky_level': 'SKYLEVEL', 'zeropoint': 'ZPT'}, **kwargs)

Bases: FITSImage

A FITS Image that has standardized header keywords corresponding to the properties defined in Image.

Setting a property also updates the internally stored header. When you construct an object, there is an optional argument header_kws (that has some sane defaults). The keys of this dictionary are (approximately) names of internal properties of the Image object. Allowed values include many of the arguments to Image.__init__. If you look at that docstring, the allwowed

WARNING: THIS CLASS (and subclasses) IS NOT THREAD SAFE. If you are using multithreading or multiprocessing, make sure that the same object of this class is not accessed by more than one thread or process at a time.

Construct a FITSImageStdHeaders model.

As with all these constructors, only use this if you really know what you’re doing. If you’re reading an image from the database, you will never use a Image constructor. Even if you’re not reading an image from the database, most of the time you’re going to use an ImageCollection, and not use an Image directly.

Parameters:

• header_kws (dict) –

  A dictionary of internal_property: header_keyword.

  internal_property is a standard internal property of Image. Those properties are mostly defined in the Image docstring. They are also many of the arguments to Image::__init__, starting with observation_id and going through sky_level. If a internal_property is not present in the header_kws dict, it implies that it can’t be found in the header, which may or may not break things. (Special case: there may also be a “zeropoint” entry in this dictionary, which isn’t a standard Image internal property, but is handled specially in this class.)

  header_keyword is, of course, the FITS header keyword to find the property in the FITS header. This means that it’s no longer than 8 characters, must be ASCII (no é or ξ or 💩) and is almost certainly ALL CAPS.

  HANDLING OF THIS IS SUBTLE; READ AND BE CAREFUL. If you pass the property directly in the constructor (so, for instance, if you make an image with:

  im = FITSImageStdHeaders( full_filepath="/path/to/file.fits", ra=42. )
  

  then later when you access im.ra, you will get 42., NOT what was in the header. However, if you access im.dec, you will get whatever was found in the header, or an exception if nothing was found in the header. The standard headers, thus, really are a fallback, but when you’re using this class, the fallback is probably what you’re really after.

  Be careful: this class blindly sets object properties based on the keys of this dictionary. If you pass the wrong things, you could break its functionality.

• zeropoint (float, default None) – Unlike most image classes (which explicitly make getting the zeropoint a function, because exactly how it’s done will be different for different kinds of images, and because it will depend on sed and maybe position), this class lets you pass one at construction time. The reason is because you might want to be setting the zeropoint in the header. (However, you still can’t set it after the object is constructed… if we need that functionality, we should add it, BUT we may have more complicated zeropoint handling anyway in the future.)

• **kwargs (Everything else is passed to parent class) – constructors (FITSImage and its parent(s)).

Methods Summary

get_fits_header()	This particular subclass will make a new header if it can't read it for the image.
get_zeropoint([x, y, sed])	Return the Image zeropoint for AB magnitudes.

Methods Documentation

get_fits_header()

This particular subclass will make a new header if it can’t read it for the image.

It will populate the header based on the header_kws value passed to the constructor, if the corresponding property in the object is not _UnsetProperty. It will then update the header so that the header has all the things in object_properties that show up in the header_kws dictionary where the object property os not _UnsetProperty.

get_zeropoint(x=None, y=None, sed=None)

Return the Image zeropoint for AB magnitudes.

By definition, the zeropoint returned by this image is a “infinite aperture” zeropoint, or one that may be used with a snappl.psf.PSF whose normalization is done right, i.e., if the clip size were infinte, the get_clip() method of the PSF object would return an infinitely-sized numpy 2d array whose sum was 1. (This is also the definition that STPSF uses when returning PSF/PRFs.) See below for much more discussion.

x, y: integers (or, I guess, floats); optional.

Pixel position on the image. Ideally, given our definition of zeropoint, these aren’t used, because the units of the .data array for a properly flatfielded and illumination-corrected image makes the the image zeropoint constant across the image. The parameters are here to hedge our bets in case a future subclass needs it. To be safe, always pass in the x, y of the position on the image where you need the zeropoint. If you don’t pass anything, and if it matters, a properly-implemented Image subclass will assume the center of the image.

sedDEFINITION STILL INCOMING; optional.

DON’T USE THIS RIGHT NOW. The interface may well change. It’s here as a placeholder to remind us we need it, and also for the docstring below.

The SED of the object for which you want a zeropoint. Exactly how we specify SEDs is not yet known, but hopefully it will be a subclass of something we define in snappl/sed.py. If not given, different subclasses will make different (maybe implicit!) assumptions. It’s possible that the subclass will not be able to take an arbitrary SED. We hope to use snappl.sed for this, but we’re still thinking it through. For now, this parameter is ignored, and you’ll get something that’s for some SED that may be not only subclass dependent, but dependent on execution details (like, for instance, some kind of weighted average of the real SEDs of the stars used to determine the image zeropoint)

zp: float

Can be used in:

m_AB = -2.5 log10(DAV) + zp

for an object with psf-fit or aperture-corrected DAV, if that object has an SED consistent with the sed parameter you passed or that is assumed by the subclass.

So that we are very clear what we mean by zeropoint as returned by the the get_zeropoint() method of a snappl.image.Image or a snappl.image.Image subclass, this is the definition.

First, imagine that you have an Image (i.e., an object of the class defined in snappl/image.py). That image’s data property is a two-dimensional array of floats. Define “DAV” (for “data array value”) as the units of that two dimensional array. To highlight this:

THE DAV IS THE UNIT OF THE NUMBERS WE GET IN THE DATA ARRAY

(This is also what we define in the docstring of the Image class itself.)

Whatever that actually is. Importantly, this definition is agnostic as to whether the data array represents something like “counts” or “counts per second”. However, it still does have opinions about the meaning of the numbers; read on.

Second, imagine that we have an astronomical source (a star, to make it concrete), and we have an image of that star taken by the telescope. (Let’s assume that our thought-experiment stars are not at all variable, so it doesn’t matter if we’re talking about the number of photons that entered the aperture during the time of the exposure, or per second.) For our zeropoint definition, we are going to assume that the number of DAVs in the Image.data array is proportional to the number of photons that entered the telescope’s aperture. [ASIDE 1: this implicitly assumes that something like bias subtraction has already been done, so there isn’t a systematic offset from pure electronic effects.] [ASIDE 2: this defintion means that DAV is NOT a surface-brightness unit! A properly implemented Image subclass is promising to do a conversion when you access the .data and .noise arrays to make sure you aren’t getting something in surface brightness units; see RomanDatamodel Image for example.] In reality, diffraction, quantum efficiency, and electronic effects will mean that some of the light energy that entered the telescope aperture will miss the detector or otherwise not be reflected in the read-out data array, but for now, let’s assume that that is negligible. Also, for definitional purposes, assume that there are absolutely no astronomical sources contributing to the light of hitting the detector than the star we’re currently pointing at.

Third, the star’s SED matters. F_ν(ν), or “flux density”, comes in dimensionality of enery/time/area/frequency. It is defined so that:

dE = A F_ν(dν) dt dν

is the amount of light energy coming from the star at frequency ν within dν that entered a telescope aperture of area A in time dt. Right now, we’re going to assume that the star has a flat spectrum, i.e., F_ν(ν) is constant for all ν. (We will relax this later; see COLOR TERMS below.)

Fourth, when we divide the image into pixels, we want the response of every pixel in the data array to be exacly the same; by “response”, we mean the conversion from number of photons entering the telescope in the sky area subtended by the pixel to DAV of the pixel. (See CORRECTING FOR PIXEL RESPONSE below.) (Also see PIXEL AREA ISSUE below.)

Fifth, let’s assume that all backgrounds (i.e., light from anything other than the one star we’re looking at) has been subtracted from the image.

Under all these assumptions, we can define the flat-spectrum zeropoint zp (which may not be exactly what get_zeropoint returns!) to be:

m = -2.5 * log10( DAVs ) + zp

where DAVs is the sum of the whole data array, and m is an AB magnitude. An AB magnitude is defined by:

m_AB = -2.5 log10( f_ν / (erg s⁻¹ Hz⁻¹ cm⁻¹) ) - 48.60

(at least if Wikipedia can be trusted). This means that a source with a flux density 3631×10⁻²³ erg s⁻¹ Hz⁻¹ cm⁻¹=3631 Jy has m_AB=0. (Closer to 3.63078054770099×10³ Jy assuming 48.60 is a definition (not a measurement with uncertainty), but 4 sig figs is plenty for a docstring.)

CORRECTING FOR PIXEL RESPONSE

For our definition to work, it means that we’re assuming some preprocessing has been done to the image by the time we receive it. Neglecting all issues of pixel area, that means pixel-to-pixel gain variables have been corrected by flatfielding, so the same zeropoint applies to every pixel on the image. It also assumes that if there is any vignetting (e.g., if the “effective telescope aperture” is different for different pixels), an illumination correction has taken all of that out.

PIXEL AREA ISSUE

When we say “pixel area” in this context, we are NOT talking about the physical area of the pixel on the array, but rather than angular area subtended on the sky by a pixel. (Yes, if we’re going to be precise, the existence of diffraction (at the very least) means that there isn’t a hard-edged area on the sky that corresponds exactly to what a given pixel absorbs, but that’s one of the big reasons we talk about PSFs for space-based imaging (on the ground, the atmosphere is usually way more significant). It is still meaningful, by putting in the right kind of delta function or whatever in place of the actual diffraction (and/or atmospheric blurring) function, to map the physical area of a pixel on the array through optics to an angular area on the sky.) This pixel area can come in units like steradian or arcsec².

In the Roman Space Telescope, the pixel area subtended on the sky can vary by ~±2% over a single SCA. The L2 maps provided by the Roman SOC have array values in units of surface brightness, i.e., something like DN/sec/steradian. However, we have defined DAV to be more like DN/sec (though, again, we are explicitly agnostic as to whether DAV is a rate or not).

What this means is that at least for L2 Roman images, the Image.data array will do a pixel-area correction before returning the DAV values; see the docstring on the Image class and on the Image.data property.

As a result of all , pixel area is not an issue for the definition of the zeropoint.

However, that also means that this zeropoint is what you’d use for point-source photometry. It is not the zeropoint you’d use to identify isophots in a galaxy. (Also, Image.data isn’t formally the right thing to use to identify isophots in a galaxy, unless the pixel area really is constant across the array!)

Note that when Image.data corrects the data to give DAV as something proportional to photon counts, not surface brightness, it fixes just purgely optical/geometric effects. For electronic effects, espeically ones that depend on how full the well is, futher corrections that cannot be encapsulated (at least currently) by the Image.get_zeropoint() method will be needed. (Thushara, save us!)

ACTUAL PHOTOMETRY

Importantly, the zeropoint we’ve defined here DOES NOT take into account any aperture size, nor does it take into account any particular realization of a PSF. It is a property of the image, not of the method used to extract photometry. That means to use this zeropoint:

• Aperture photometry values must be properly “aperture corrected” before the DAVs are fed into the zeropoint formula. Ideally, when things aren’t too complicated, this correction is just a single factor that multiplies the number of DAVs in the aperture to give an effective “infinite aperture” number of DAVs. This factor will, of course, be different for apertures of different sizes (and shapes), and will also in principle be different at different positions on a detector array. (For small apertures, it’s also very difficult to do right.) For real images, it’s very difficult to determine this by looking at stars on images; you find yourself stuck between needing a very big aperture to capture, within your precision, “all” the flux, and not wanting your aperture to be too big so that you can find enough isolated stars. If you have a very good estimate of the PSF/PRF, you can determine an aperture correction by integrating that.

• PSF (or PRF) photometry must use PSFs (or PRFs) that are properly normalized to fit the defintion here. “Properly normalized” here means that if you had an infinitely-sized image-scale array of the PSF (really PRF), its sum would be 1; in practice, because you can’t get infinitely-sized data arrays, the sum of the array you get will be less than 1, though for a big enough stamp size it might be very close. The PSFs (really PRFs) returned by snappl.psf.PSF.get_stamp() (and other methods) are supposed to be normalized this way. (Also, as I understand it, the PSFS you get from STPSF are also normalized this way.)

  IT IS POSSIBLE that some further calibration post-processing of photometry after the zeropoint is applied may be entirely convolved with the definition of the PSF. At the moment, snappl’s class structure does not support this, but we will adapt if necessary. However, we should ONLY adapt if it really is necessary! If it’s just a matter of normalizing your PSFs differently, then just normalize them differently to fit our definitions!

FILTERS AND COLOR TERMS

In reality, we never measure something proportional to F_ν(ν) directly. (Spectroscopy gets a lot closer to this than photometry does.) Rather, we’re always measuring some integral of F_ν(ν). There are two things we have to consider.

First, detectors and filters (and the whole telescope system, for that matter) have a different response at different frequencies. Filters, in particular, only transmit light within a finite range of ν, though real detectors are also not sensitive to all frequencies. We will call the system response D(ν), which we will define “the number of DAVs that we get in our data array per frequency bin for light of frequency ν for a source with f(ν)=3631 Jy”, i.e., if we’re looking at that hypotetical f(ν)=3631 Jy star:

DAVs = ∫ D(ν) dν

This means that D(ν) has units of s (or, more clearly, Hz⁻¹) (or, maybe, if you don’t think of DAVs as dimensionless, units of DAV/Hz).

Second, astronomical sources do not have a flat F_ν(ν), as we assumed in our discussion above and in the definiton of the thing we called zp. The actual light source is going to have some SED S(ν) (in units of Energy/Time/Flux Binwidth/Area).

The total number of DAVs detected, therefore, is:

DAVs = ∫ S(ν) D(ν) / (3631Jy) dν

(Presumably D(ν) goes to zero outside some finite range of ν so we don’t have to think about infinite numbers.)

Given this, the flat-spectrum zeropoint (which is what we defiend as zp above) is defined as:

zp = 2.5 log10( ∫ D(ν) dν )

(To see this: consider S(ν) = 3631 Jy for all ν, which is the definition of a m_AB=0 source. In this case:

0  = -2.5 log10(DAVs) + zp
   = -2.5 log10( ∫ (3631Jy) D(ν) / (3631Jy) dν ) + zp
   = -2.5 log10( ∫ D(ν) dν ) + zp
zp = 2.5 log10( ∫ D(ν) dν )

)

A flat-spectrum source with flux density S₀ at all ν has AB magnitude:

m_AB = -2.5 log10( S₀ / 3631Jy ) = -2.5 log10( S₀/Jy ) + 8.900

(Which is where “8.900 is the AB zeropoint” comes from. You will sometimes see people using a zeropoint of 31.4; this is just the zeropoint where the flux density is in nJy rather than Jy, as 2.5log10(10⁹)=22.5.)

The number of DAVs from such a source would be:

DAVs = ∫ S₀ D(ν) / (3631Jy) dν = S₀ / 3631Jy * ∫ D(ν) dν

or:

DAVs / ( ∫ D(ν) dν ) = S₀ / 3631Jy

Taking logs of both sides:

-2.5 log10( DAVs ) + 2.5 log10( ∫ D(ν) dν ) = -2.5 log10( S₀/Jy ) + 2.5 log10( 3631 )
-2.5 log10( DAVs ) + zp = -2.5 log10( S₀ ) + 8.900 = m_AB

Where it gets painful is when S(ν) is not constant with ν. In this case, the magnitude you will calculate from the flat-spectrum zeropoint would be:

m_calc = -2.5 log10( DAVs ) + zp
       = -2.5 log10( ∫ S(ν) D(ν) / (3631Jy) dν ) + 2.5 log10( ∫ D(ν) dν )

However, for a source that doesn’t have a flat S(ν), the true AB magnitude is ill-defined, because it’s different for every ν! So, for a given filter, we have to define a fiducial frequency ν₀ (which corresponds to a fiducial wavelength λ₀ by the usual ν₀=hc/λ₀). We could then define the “true” apparent magnitude of the object with SED S(ν) as:

m = -2.5 log10( S(ν₀) / 3631 Jy )

I think all the Roman filters have a defined fiducial wavelength, so we should just use that (really, hc/that) for ν₀, but we may need to document this somewhere.

We then have a SED correction:

cor_sed ≡ m - m_calc
        = -2.5 log10( S(ν₀) / 3631 Jy ) + 2.5 log10( ∫ S(ν) D(ν) / (3631Jy) dν ) - 2.5 log10( ∫ D(ν) dν )
        = 2.5 log10( ∫ S(ν) D(ν) / S(ν₀) dν ) - 2.5 log10( ∫ D(ν) dν )

cor_sed = 2.5 log10( ∫ S(ν) D(ν) / S(ν₀) dν ) - zp

The magnitude of an object is then:

m = -2.5 log10( DAVs ) + zp + cor_sed

(Do not become confused by the fact that zp is in cor_sed; we’re not subtracting out the zeropoint from the final magnitude formula, because it’s added back, sorta, inside the integral in cor_sed, we just can’t separate it out to another obvious +zp because it’s inside the integral, and while I’ve known named-chair professors of physics (but not astronomy) to claim that we were all doing cosmology wrong and making it too complicated because he freely factored variable things out of integrals, you aren’t really supposed to do that.)

Notice that you don’t need to know the absolute S(ν) to calculate cor_sed, only S(ν)/S(ν₀). This is why we say the “shape” of the SED. The thing passed to the sed parameter of get_zeropoint() is really a SED shape (though a cautiously implemented subclass will not assume that the user is passing in a properly normalized sed).

What get_zeropoint() returns is:

zp + cor_sed

for the sed specified by the sed parameter (with an implicitly assumed ν₀), or for some default sed if you don’t specify one. IMPORTANT, don’t assume this is a flat spectrum, because in practice that may be difficult or impossible to determine. Each subclass may assume a different cor_sed (at least for now).