RomanDatamodelImage

class snappl.image.RomanDatamodelImage(*args, **kwargs)

Bases: Image

An image read from a roman datamodel ASDF file.

See Issue #46 for concerns about performance/memory and imlementation of this object.

Instantiate an image. You probably don’t want to do that.

This is an abstract base class that has limited functionality. You probably want to instantiate a subclass if you’re creating a new image.

If you’re trying to pull an image out of the database, then probably what you really want is to use the Image.get_image or Image.find_images class methods.

If you’re working with non-database images and are trying to get a pre-existing image, then probably what you really want to do is call the get_image() method of an ImageCollection object.

Only instantiate an image directly if you’re creating something yourself that you know you want to write in a specific format, or if you’re trying to read a file that’s not covered by an ImageCollection. When you do this, talk to the photometry working group and find out if this image should be covered by an ImageCollection. (Note that there is an ImageCollectionManualFITS collection for reading loose FITS files.)

Parameters:

• filepath (str or Path, default None) – Path of the image relative to the base path for images, unless less no_base_path is True, in which case this is the full absolute path to the image. For datbase images, you do not want to create a path yourself, but leave it at None and let the class create the filepath. See PathedObject.

• full_filepath (str or Path, default None) – The full path to the image. If you’re using an Image subclass to deal with an image that’s not in the database, you probably want to set this to the absolute path of the image, and you probably want to set no_base_path to True, but you might also set base_path yourself and leave no_base_path at False.

• base_path (str or Path, default None) – Always leave this at None for images associated with database, and the default will be used. Otherwise, the absolute path of the image is base_path / filepath (which should be exactly the same as full_filepath). Must be None if no_base_path is True.

• base_dir (str or Path, default None) – Synonym for base_path

• no_base_path (bool, default False) – For images associated with the database, leave this at False, and make filepath relative to the base path (which may be system dependent). For images that aren’t associated with the database, you can make this True and set filepath to be just the path to the image.

• id (UUID or str that can be converted to UUID, default None) – Database ID of the image. This is only relevant if the image is in the l2image table of the Roman SNPIT internal database (but is required in that case).

• provenance_id (UUID or str that can be converted to UUID, default NOne) – The id of the provenance of the image. Only relevant if the image is in the l2image table of the Roman SNPIT internal database (but is required in that case).

• width (int, default None) – The width and height of the image in pixels if known.

• height (int, default None) – The width and height of the image in pixels if known.

• format (int, default -1) – Index into the table Image._format_def at the bottom of this file.

• is_superclass – Used internally, should ONLY ever be set in the super().__init__(…) lines in subclass constructors. All subclasses should set this to True when calling super().__init__(…). If you aren’t writing an Image subclass, ignore this.

Attributes Summary

data	image data in DAV.
dm	This property should usually not be used outside of this class.
flags	An integer 2d numpy array of pixel masks / flags TBD
noise	1σ pixel noise.
sb_data	NOT A STANDARD Image PROPERTY! Surface-brightness units data array.
sb_noise	NOT A STANDARD Image PROPERTY! Surface-brightness units noise array.

Methods Summary

free()	Try to free memory.
get_cutout(x, y, xsize[, ysize, mode, ...])	See Image.get_cutout The mode and fill_value parameters are passed directly to astropy.nddata.Cutout2D for FITSImage.
get_data([which, always_reload, cache])	Read the data from disk and return one or more 2d numpy arrays of data.
get_ra_dec_cutout(ra, dec, xsize[, ysize, ...])	See Image.get_ra_dec_cutout
get_wcs([wcsclass])	Get image WCS.
get_zeropoint([x, y])	Return the Image zeropoint for AB magnitudes.

Attributes Documentation

data

dm

This property should usually not be used outside of this class.

flags

noise

sb_data

NOT A STANDARD Image PROPERTY! Surface-brightness units data array.

This is the native data array straight out of the roman_datamodel L2 asdf files.

sb_noise

NOT A STANDARD Image PROPERTY! Surface-brightness units noise array.

This is the native noise array straight out of the roman_datamodel L2 asdf files.

Methods Documentation

free()

Try to free memory.

get_cutout(x, y, xsize, ysize=None, mode='strict', fill_value=nan, return_FITS=True)

See Image.get_cutout The mode and fill_value parameters are passed directly to astropy.nddata.Cutout2D for FITSImage.

Inputs

return_FITS: bool, default True

If True, the cutout will be returned as a snappl.image.FITSImage. If False, the cutout will be returned as a snappl.image.RomanDatamodelImage.

get_data(which='all', always_reload=False, cache=False)

Read the data from disk and return one or more 2d numpy arrays of data.

See Image.get_data for definition of parameters.

get_ra_dec_cutout(ra, dec, xsize, ysize=None, mode='strict', fill_value=nan)

See Image.get_ra_dec_cutout

The mode and fill_value parameters are passed directly to astropy.nddata.Cutout2D for FITSImage.

get_wcs(wcsclass=None)

Get image WCS. Will be an object of type BaseWCS (from wcs.py) (really likely a subclass).

Parameters:

wcsclass (str or None) – By default, the subclass of BaseWCS you get back will be defined by the Image subclass of the object you call this on. If you want a specific subclass of BaseWCS, you can put the name of that class here. It may not always work; not all types of images are able to return all types of wcses.

Return type:

object of a subclass of snappl.wcs.BaseWCS

get_zeropoint(x=None, y=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).