Calculating WFC3 Zeropoints with STSynphot

Calculating WFC3 Zeropoints with STSynphot#

Learning Goals#

By the end of this tutorial, you will:

  • Calculate zeropoints and other photometric properties using stsynphot.

  • Create, plot, and save ‘total system throughput’ tables.

Table of Contents#

Introduction
1. Imports
2. Download throughput tables and define variables
3. Set up the ‘obsmode’ string
4. Basic usage for a single ‘obsmode’
5. Compute zeropoints and other photometric properties
6. Iterate over multiple ‘obsmodes’
7. Create and plot ‘total system throughput’ tables
8. Conclusions
Additional Resources
About the Notebook
Citations

Introduction#

This notebook shows how to calculate photometric zeropoints using the Python package stsynphot for any WFC3 detector, filter, date, or aperture. This tutorial is especially useful for calculating Vegamag zeropoints, which require an input spectrum. The notebook is also useful for computing time-dependent WFC3/UVIS zeropoints for any observation date, as the values listed in WFC3 ISR 2021-04 are defined for the reference epoch. As of mid-2021, the WFC3/IR zeropoints are not time-dependent.

More documentation on stsynphot is available here. Using stsynphot requires downloading the throughput curves for the HST instruments and optical path. One method of doing this is shown in Section 2. More information on the throughput tables can be found here.

1. Imports#

This notebook assumes you have created the virtual environment in WFC3 notebooks’ installation instructions.

We import:

  • os for setting environment variables

  • tarfile for extracting a .tar archive

  • numpy for handling array functions

  • matplotlib.pyplot for plotting data

  • astropy for astronomy related functions

  • synphot and stsynphot for evaluating synthetic photometry

We will need to set the PYSYN_CDBS environment variable before importing stsynphot. We will also create a custom Vega spectrum, as the stsynphot will supercede the usual synphot functionality regarding the Vega spectrum and would otherwise require a downloaded copy of the spectrum to be provided.

import os
import tarfile

import numpy as np
import matplotlib.pyplot as plt

from astropy.table import Table

from synphot import Observation

2. Download throughput tables and define variables#

This section obtains the WFC3 throughput component tables for use with stsynphot. This step only needs to be done once. If these reference files have already been downloaded, this section can be skipped.

cmd_input = 'curl -O https://archive.stsci.edu/hlsps/reference-atlases/hlsp_reference-atlases_hst_multi_everything_multi_v11_sed.tar'
os.system(cmd_input)

Once the downloaded is complete, extract the file and set the environment variable PYSYN_CDBS to the path of the trds subdirectory. The next cell will do this for you, as long as the .tar file downloaded above has not been moved.

tar_archive = 'hlsp_reference-atlases_hst_multi_everything_multi_v11_sed.tar'
extract_to = 'hlsp_reference-atlases_hst_multi_everything_multi_v11_sed'
with tarfile.open(tar_archive, 'r') as tar:
    tar.extractall(path=extract_to)

os.environ['PYSYN_CDBS'] = 'hlsp_reference-atlases_hst_multi_everything_multi_v11_sed/grp/redcat/trds/'

Now, after having set up PYSYN_CDBS, we import stsynphot. A warning regarding the Vega spectrum is expected here.

import stsynphot as stsyn

Rather than downloading the entire calspec database (synphot6.tar.gz), we can point directly to the latest Vega spectrum which is required for computing VEGAMAG.

vega_url = 'https://ssb.stsci.edu/trds/calspec/alpha_lyr_stis_010.fits'
stsyn.Vega = stsyn.spectrum.SourceSpectrum.from_file(vega_url)

3. Set up the ‘obsmode’ string#

Parameters to set in the obsmode string include:

  1. detector,

  2. filter,

  3. observation date (WFC3/UVIS only), and

  4. aperture size (in arcsec).

Note that a 6.0” aperture is considered to be “infinite”, thus containing all of the flux. The zeropoints posted on the WFC3 website are calculated for an infinite aperture, so when computing photometry for smaller radii, aperture corrections must be applied.

The inputs below can be changed to any desired obsmode, with examples of alternate parameters shown as commented lines.

First, here are some detector examples with WFC3/UVIS1 as the default, and other options including both WFC3/UVIS chips or the WFC3/IR detector.

Note: if the IR detector is chosen, the filtnames below must be updated.

detectors = ['uvis1']
# detectors = ['uvis1', 'uvis2']
# detectors = ['ir']

Next, here are some filter examples with all WFC3/UVIS filters as the default, and other options including just F606W and the WFC3/IR filters.

Note: if WFC3/IR filters is chosen, the detectors above must be set to [‘ir’].

filtnames = ['f200lp', 'f218w', 'f225w', 'f275w', 'f280n', 'f300x', 'f336w', 'f343n', 'f350lp', 
             'f373n', 'f390m', 'f390w', 'f395n', 'f410m', 'f438w', 'f467m', 'f469n', 'f475w', 
             'f475x', 'f487n', 'f502n', 'f547m', 'f555w', 'f600lp', 'f606w', 'f621m', 'f625w', 
             'f631n', 'f645n', 'f656n', 'f657n', 'f658n', 'f665n', 'f673n', 'f680n', 'f689m', 
             'f763m', 'f775w', 'f814w', 'f845m', 'f850lp', 'f953n']
# filtnames = ['f606w']   
# filtnames = ['f098m', 'f105w', 'f110w', 'f125w', 'f126n', 'f127m', 'f128n', 'f130n', 
#              'f132n', 'f139m', 'f140w', 'f153m', 'f160w', 'f164n', 'f167n']

Now, here are some date examples with the WFC3/UVIS reference epoch (55008 in MJD; 2009-06-26) as the default, and the other option being the time right now.

mjd = '55008'
# mjd = str(Time.now().mjd)

Finally, here are some aperture radius examples with 6.0” (151 pixels; “infinity”) as the default, and the other options including 0.396” (10 pixels for WFC3/UVIS) and 0.385” (3 pixels for WFC3/IR).

aper = '6.0'
# aper = '0.396'
# aper = '0.385'

4. Basic usage for a single ‘obsmode’#

The calculation of the zeropoints starts with creating a specific bandpass object. Bandpasses generally consist of at least an instrument name, detector name, and filter name, though other parameters (such as the MJD and aperture radius shown above) are optional.

The cell below defines obsmode and creates a bandpass object.

obsmode = 'wfc3,uvis1,f200lp'
bp = stsyn.band(obsmode)

Optional parameters are supplied on the end of the basic bandpass:

obsmode = 'wfc3,uvis1,f200lp,mjd#55008,aper#6.0'
bp = stsyn.band(obsmode)

In addition, we can use the parameters defined in Section 3.

obsmode = f'wfc3, {detectors[0]}, {filtnames[0]}, mjd#{mjd}, aper#{aper}'
bp = stsyn.band(obsmode)

5. Compute zeropoints and other photometric properties#

With the bandpass objects, we can now calculate zeropoints, pivot wavelengths, and photometric bandwidths. To calculate Vegamag zeropoints, we use the Vega spectrum to calculate the flux in a given bandpass.

def calculate_values(detector, filt, mjd, aper):
    # parameters can be removed from obsmode as needed
    obsmode = f'wfc3, {detector}, {filt}, mjd#{mjd}, aper#{aper}'
    bp = stsyn.band(obsmode)  
    
    # STMag
    photflam = bp.unit_response(stsyn.conf.area)  # inverse sensitivity in flam
    stmag = -21.1 - 2.5 * np.log10(photflam.value)
    
    # Pivot Wavelength and bandwidth
    photplam = bp.pivot() # pivot wavelength in angstroms
    bandwidth = bp.photbw() # bandwidth in angstroms
    
    # ABMag
    abmag = stmag - 5 * np.log10(photplam.value) + 18.6921
    
    # Vegamag
    obs = Observation(stsyn.Vega, bp, binset=bp.binset)  # synthetic observation of vega in bandpass using vega spectrum
    vegamag = -1 * obs.effstim(flux_unit='obmag', area=stsyn.conf.area)
    
    return obsmode, photplam.value, bandwidth.value, photflam.value, stmag, abmag, vegamag.value

obsmode, photplam, bandwidth, photflam, stmag, abmag, vegamag = calculate_values(detectors[0], filtnames[0], mjd, aper)

# print values
print('Obsmode                              PivotWave Photflam   STMAG   ABMAG   VEGAMAG')
print(f'{obsmode}, {photplam:.1f}, {photflam:.4e}, {stmag:.3f}, {abmag:.3f}, {vegamag:.3f}')

6. Iterate over multiple ‘obsmodes’#

To calculate zeropoints for multiple detectors and/or filters, we can use the function defined above and loop through detectors and filters defined in Section 3.

oms, pivots, bws, pfs, st, ab, vm = [], [], [], [], [], [], []

print('Obsmode                              PivotWave Photflam   STMAG   ABMAG   VEGAMAG')
for detector in detectors:
    for filt in filtnames:
        res = calculate_values(detector, filt, mjd, aper)
        obsmode, photplam, bandwidth, photflam, stmag, abmag, vegamag = res # solely for readability
        
        # print values
        print(f'{obsmode}, {photplam:.1f}, {photflam:.4e}, {stmag:.3f}, {abmag:.3f}, {vegamag:.3f}')
        
        oms.append(obsmode)
        pivots.append(photplam)
        bws.append(bandwidth)
        pfs.append(photflam)
        st.append(stmag)
        ab.append(abmag)
        vm.append(vegamag)

Values can also be written into an astropy table.

tbl = Table([oms, pivots, bws, pfs, st, ab, vm], 
            names=['Obsmode', 'Pivot Wave', 'Bandwidth', 'Photflam', 'STMag', 'ABMag', 'VegaMag'])

We’ll also round columns to a smaller number of decimals.

for col in tbl.itercols():
    if col.name == 'Photflam':
        col.info.format = '.4e'
    elif col.info.dtype.kind == 'f':        
        col.info.format = '.3f'

Let’s view our astropy table:

tbl

We can finally save the table as a .txt file.

if not os.path.exists('./uvis_zp_tbl.txt'):
    tbl.write('uvis_zp_tbl.txt', format='ascii.commented_header')

7. Create and plot ‘total system throughput’ tables#

The function below returns a tuple containing two objects, the first being an array of wavelengths, and the second being the throughput at each of those wavelengths.

def calculate_bands(bp, save=False, overwrite=True):
    # Pass in bandpass object as bp
    waves = bp.waveset
    throughput = bp(waves)
    
    if save:
        tmp = Table([waves, throughput], names=['WAVELENGTH', 'THROUGHPUT'])
        tmp.write(', '.join(bp.obsmode.modes)+'.txt', format='ascii.commented_header', overwrite=overwrite)
        
    return (waves, throughput)

We’ll calculate the throughput table for WFC3/UVIS1 in F200LP.

obsmode = 'wfc3,uvis1,f200lp'
bp = stsyn.band(obsmode)
wl, tp = calculate_bands(bp)

Now, let’s plot our results.

fig = plt.figure(figsize=(10, 5))
plt.plot(wl, tp)
plt.xlim(1500, 11000) 
plt.xlabel('Wavelength [Angstroms]')
plt.ylabel('Throughput')
plt.title('WFC3,UVIS1,F200LP')

To save the curve in an ascii table, simply pass the argument save=True:

calculate_bands(bp, save=True)

To save curves for all obsmodes defined in Section 3 in the input list, we can loop through detectors and filters.

for det in detectors:
    for filt in filtnames:
        obsmode = f'wfc3, {det}, {filt}'
        bp = stsyn.band(obsmode)
        calculate_bands(bp, save=True)

In addition, we’ll create a directory called obsmodes_curves and move all the saved files to that directory.

! mkdir obsmodes_curves
! mv wfc3*txt obsmodes_curves
! ls obsmodes_curves

8. Conclusions#

Thank you for walking through this notebook. Now using WFC3 data, you should be more familiar with:

  • Calculating zeropoints and other photometric properties using stsynphot.

  • Creating, plotting, and saving ‘total system throughput’ tables.

Congratulations, you have completed the notebook!#

Additional Resources#

Below are some additional resources that may be helpful. Please send any questions through the HST Helpdesk.

About this Notebook#

Authors: Varun Bajaj, Jennifer Mack; WFC3 Instrument Team

Updated on: 2024-03-18

Citations#

If you use numpy, astropy, synphot, or stsynphot for published research, please cite the authors. Follow these links for more information about citing the libraries below:

Top of Page Space Telescope Logo

Contents