Final GSoC Post

Wow, what a journey! It’s hard to believe that GSoC is coming to an end. This project has taken quite a few twists and turns, which I will attempt to lead you through here, but ultimately I think it has all come together into a product that will be useful for anyone in the astronomical community who is planning an observation.

Looking back on the proposal to Google Summer of Code, our original plan was to build “a simple API for retrieving signal-to-noise estimates for an arbitrary astronomical observation” – in other words, we wanted to create something which could simulate an observation given a celestial object and some telescope specifications, then return a signal to noise ratio one could expect from such an exposure. I’m happy to say that we have succeeded, though not exactly in the way we had originally imagined.

In the beginning we had envisioned this functionality taking the form of a whole new package called telescopy, but it quickly became clear that a good chunk of the work behind this already existed as an astropy affiliated package called synphot. Specifically, synphot provided the means to create the synthetic photometry we needed to evaluate the expected count rate of an observation, and therefore the exposure time/signal to noise.

Much of the first half of GSoC therefore revolved around unveiling what exactly we could use from synphot and what more we needed to code. To do this, we decided that I would write a suite of Jupyter tutorials showing how to use synphot to obtain a count rate, signal to noise ratio (SNR), and exposure time for a few example observations, which can be seen at the following links:

The conversations between me and my mentors on the development of these tutorials can be seen by browsing these pull requests (which are just made to my own github simply for the sake of following the open source philosophy).

Parts of these tutorials now exist as a pull request to astropy-tutorials, which I expect to be merged soon! All that needs to be done in that regard is to get the tutorials reviewed by one more astropy-tutorials collaborator (and for me to respond to those suggestions). When that is done, our tutorials will exist on the astropy-tutorials website to help others learn how to use synphot.

The original tutorials I wrote include SNR and exposure time calculations, which I had integrated into the synphot code on my own branch. My additions to synphot (which will not be merged – more on that later) can be seen at the PR here, where the relevant discussion on the changes can be seen by navigating to the “Conversation” tab on the same page.

Once I had written these tutorials (and after much head scratching), it became evident that my additions would live much more comfortably inside of astroplan instead of synphot, since the exposure time and SNR calculators I had written are much more relevant to planning observations (indeed, the express purpose of astroplan), than to creating synthetic photometry. So at this point I shifted gears and started thinking about how what I had written could fit into the astroplan API.

The code (and relevant tests) that I’ve written and reference to in the remainder of the write-up can be viewed by clicking this link. To be clear, everything highlighted in green/red has been added/removed by me.

In brief, the additions I’ve made to astroplan now allow it to not only give when it’s best to observe, but also for how long to observe, given a desired SNR. It sounds simple at face value, but to do this, astroplan needed to learn how to:

  1. Model an observation, which means convolving a model spectrum with:
    • a bandpass
    • CCD response function
    • atmospheric transmission
    • etc…
  2. Measure the count rate of the observation
  3. Calculate the required exposure time for the given SNR

The bandpass and CCD response function are all inherent to whatever telescope/instrument is in use, so it seemed appropriate to write a Telescope class to the new astroplan.telescope module. The Telescope object serves as a container for the specs of the observer’s instrument, where the diameter, bandpass, and gain are given by the user and set as attributes (optional attributes include the read noise and CCD response function).

The atmospheric transmission is instead inherent to the observer’s location. astroplan already has a module for the observer’s location (astroplan.observer), but it didn’t yet include an option to model the atmospheric transmission. Enter astroplan.skymodel, a new module I’ve written which queries an atmospheric transmission model from the SkyCalc Sky Model Calculator. To obtain a sky model, the user just has to set the skymodel attribute of their Observer object, where all of the variable parameters seen on the SkyCalc page can be set as so (if none are set, skymodel just queries the default parameters seen on the SkyCalc page):

>>> from astroplan import Observer, skymodel

>>> observer = Observer.at_site('apo')
>>> observer.skymodel = skymodel(airmass=2)

Finally, I wrote the exptime module and added it to astroplan for calculating an exposure time using the Observer and Telescope objects. This module (and astroplan.telescope) require that synphot be installed, since exptime calls on synphot to convolve the spectral elements given and to calculate the subsequent count rate of the simulated observation. The expressions for actually calculating the exposure times given an idealized SNR are taken from pgs 57-58 of the trusty “Handbook of CCD Astronomy” by Howell (2000).

To see how all of this is used on the user-side of things, I have also written a tutorial to the astroplan documentation (unfortunately this can only really be seen once it’s been merged, but the code can be viewed in the /docs/tutorials/exptime.rst file).

What’s left to do:

  • Respond to suggestions on my astropy-tutorials PR
  • Submit a PR to astroquery for the SkyCalc query I wrote (permission has been granted by the maintainers of the SkyCalc Calculator)
  • Finish writing tests for astroplan.exptime
  • Get suggestions from collaborators on any changes that should be made to my astroplan additions
  • Get merged!!!

I have learned SO MUCH as a GSoC student, and am really grateful to have been a part of it. I can confidently say that I’ve become a better programmer thanks to all of the excellent advice from my GSoC mentors. You all are the best! I look forward to contributing to astropy in the future and hope to run into you soon!

Week 13: Hello from Reykjavik!

What I completed this week:

  • Moved basically all of the stuff I’ve been exploring with synphot over to astroplan, which in hindsight makes more sense! In brief, this is what I did:
    • Added a new module exptime.py to astroplan, which uses synphot to calculate the exposure time needed to obtain a given signal to noise ratio. Even though it uses synphot, it is not a required package of astroplan unless the user wants to use exptime_from_ccd_snr
    • Added a new class called Telescope which carries information about the instrument the observer is using, such as the bandpass, gain, read noise, and CCD response function.
    • Added a new module skycalc.py which has a function to query the SKYCALC Sky Calculator for optional modeling of atmospheric transmission
  • Wrote a tutorial for using astroplan to calculate exposure time, which you can see here
I’m in Iceland this week for the 4th Extreme Solar Systems conference! Which is appropriate, because Iceland has what is closest to an alien landscape that I can imagine.

This week’s goals:

  • Write tests for the new functionality in astroplan
  • Add examples to the astroplan documentation
    • I think I’ll just copy and paste my Jupyter notebook tutorial into something Sphinx will display nicely
    • Add a reference to the exptime example mentioned above to the “getting started” part of a documentation. I originally wanted to add a simple 1-2 line example to “getting started”, but since you have to model the spectrum of the target and all that I think it would be more appropriate to just keep the example on its own page
  • Write my final blog post, which will summarize everything I’ve done and serve as my “code submission” in the final evaluation!

Week 12: astroplanning

What I completed this week:

  • Submitted my synphot tutorials as a pull request to astropy/astropy-tutorials. Just waiting for someone (outside of my GSoC mentors) to review!
  • Moved my synphot work (including the wrapper and signal-to-noise functions) into an experimental folder. I’m still uncertain whether this folder will be merged as an experimental folder or if its contents will be merged back into the synphot API, but I guess we will see..! Also, I’m still open to suggestions on any other functionality I should include.
  • Emailed ESO inquiring about any plans to merge skycalc_cli into astroquery, and if they would be open to me submitting a PR (no response yet).
  • Thought more about how my wrapper, signal-to-noise functions, and “synspec” can be implemented into already-existing astronomical observing tools. Simply put, I think my GSoC work (and synphot in general) would serve as valuable additions to astroplan. As it stands, astroplan is a great tool for scheduling observations, but its definition leaves its functionality open to aspects of planning which are outside of the scheduling realm:

The goal of astroplan is to make a flexible toolbox for observation planning and scheduling. 

The astroplan Documentation

Specifically, I think a count rate function and exposure time calculator would be of greatest use. Then not only could astroplan help determine when to observe, it could also suggest exposure durations given model photometric/spectroscopic observations generated by synphot/synphot tools.

I also took a weekend trip with my family to go rafting on the Yakima river (pictured here). And I fed a cow. Its tongue was disturbingly pokey.

This week’s goals:

  • If everyone is open to the idea, I would like to submit a PR to astroplan for adding the exposure time calculator (supported by synphot and its tools)

Week 11: A wrapper says its first word

What I completed this week:

  • Continued working on the synphot wrapper class, which we’ve renamed “Exposure” (since “Observe” is a little too close to the other synphot class “Observation”). So far this class works pretty much like how I imagined it last week, where the most simple (read: default) synthetic observation can be initiated given a source spectrum, telescope diameter, and exposure time:
>>> from synphot.spectrum import SourceSpectrum
>>> from synphot.exposure import Exposure
>>> import astropy.units as u

>>> source_spec = SourceSpectrum(...)
>>> tele_diameter = 1 * u.m
>>> exptime = 1 * u.s

>>> some_observation = Exposure(source_spec, tele_diameter, exptime)

Then attributes can be called and methods applied to view some characteristics of this observation, such as the count rate, signal to noise ratio (SNR), and the exposure time needed to obtain a specified SNR:

>>> some_observation.aperture_area
<Quantity 0.78539816 m2>

>>> some_observation.countrate
<Quantity 126398.52488845 ct / s>

>>> some_observation.snr()
<Quantity 355.52570215>

>>> desired_snr = 100
>>> some_observation.required_exptime(desired_snr)
<Quantity 0.07911485 s>
I also took a little break to see family and go “tubing”, as seen above (wee)

This week’s goals:

  • Finish the wrapper! Now taking suggestions for anything else I should include
  • Take a look at telescopy and astroplan and think how they, along with synspec and any higher order functions we come up with, can be merged (okay, I really need to do this one this week)

Longer term goals:

Week 10: A wrapper is born

What I completed this week:

  • Continued working on the signal to noise functions pull request to synphot
    • Changed function names from howell_snr to ccd_snr and exptime_from_howell_snr to exptime_from_ccd_snr, which are more descriptive/informative
    • The units returned by ccd_snr are now dimensionless
  • Began writing a synphot wrapper. This work in progress PR is a bit of an embarrassing mess right now, but I’m envisioning something where you can initiate a class with the source spectrum of an object, give a telescope diameter, and boom, you can get attributes like SNR, countrate, exptime, etc… :
>>> from synphot.spectrum import SourceSpectrum
>>> from synphot.wrapper import Observe
>>> import astropy.units as u

>>> tele_diameter = 1 * u.m
>>> source_spec = SourceSpectrum(...)
>>> some_observation = Observe(source_spec, tele_diameter)

>>> some_observation.countrate
1e6 ct/s
>>> some_observation.snr
100
...

The observation can be customized by setting a bandpass, giving a desired SNR (to get a required exposure time), providing the quantum efficiency of the CCD or an atmospheric transmission model, etc.

But I’m not really sure if attributes are the best way to get these, er, “attributes” of the observation. Maybe argument-less methods would be better, e.g. some_observation.snr()? That might make it easier to alter some parameters like exposure time without having to create a new object… Anyway, this week has been full of thoughts like these.

How I kind of feel writing this. Something good will come together soon, I’m sure

This week’s goals:

  • Take a look at telescopy and astroplan and think how they, along with synspec and any higher order functions we come up with, can be merged

Longer term goals:

Week 9: synspec and other thoughts

What I completed this week:

  • Traveled to Seattle, where I’ll be for the rest of the summer
  • Added documentation on the signal to noise ratio (SNR) functions to the synphot docs
  • Integrated the SNR functions into the synphot tutorials I’ve been working on
  • Checked that the binned photometry (aka spectroscopy) is functioning correctly by comparing the count rate for the HAT-P-11 source spectrum and to the countrate for the binned spectrum (i.e. the “synspec” spectrum).
    The count rate for the binned spectrum is only %0.06 greater than the source spectrum’s, so I think it’s safe to say that synphot‘s binned photometry is behaving how we want it to.
  • Wrote a rough draft of what could be the synspec aspect of a higher function (see “Longer term goals” below)
I also took my dog to the groomers and I haven’t stopped squealing since

This week’s goals:

  • (minor) Fix units in SNR functions
  • Spend time with the notebook tutorials I’ve written and think about any higher order functions that could make their functionalities more streamlined
  • Take a look at telescopy and astroplan and think how they, along with synspec and any higher order functions we come up with, can be merged

Longer term goals:

  • Some thoughts on a possible deliverable, copied from my Slack message to the team:

So I’ve been thinking a bit more about how to incorporate synthetic spectroscopy into a package for synthetic photometry, and I’ve come full circle back to “telescopy”. At our last meeting we talked about writing a kind of wrapper for the synphot concepts covered in the tutorials, and I’m wondering if telescopy can be that wrapper.

Where synphot has a wide range of functionality, perhaps telescopy could be a package which utilizes synphot for a more “streamlined” purpose, for example “telescopy: a package for estimating the exposure time needed to obtain x-SNR with y-instrument”, where y-instrument could be either photometric or spectroscopic in nature.

The problem of this which we’ve discussed before is that we would be creating “yet another package”. But it would be a solution which provides both a wrapper and a synspec resource while avoiding the whole “why is there a synthetic spectroscopy generator inside of a synthetic photometry package??” confusion.

Week 8: More synspec and SNR

What I completed this week:

  • Made a pull request to synphot for implementing the signal-to-noise calculators that I wrote last week. I forgot to mention in my blog post last week that writing this module was a really useful exercise for me because I was introduced to so many new things! Like:
    • writing tests for pytest (and using the assert_quantity_allclose function from astropy.tests.helper)
    • python decorators (and how to ensure parameters are given in compatible units with astropy.units.quantity_input)
    • “duck typing” and the hasatrr() function
    • private functions
    • citing references in the docstring in the “astropythonic” way
    • a little bit of Sphinx documentation syntax
Me every day as a GSoC student
  • Wrote the first-draft of a “synspec” (synthetic spectroscopy) example. To simulate spectroscopy, I essentially create a custom waveset derived from some spectral resolution R and use synphot to bin the flux accordingly. There are a couple of things I’m unsure about though:
    • Does the binning done by synphot do (for all intents and purposes) the same thing as a real spectrometer? I haven’t done much work with spectroscopy so the mechanics are a little mysterious to me…
    • The spectral resolution R is defined to be R = λ / Δ λ , but what is the “correct” way of constructing an array of wavelengths according to this resolution? Do we keep R constant over the range of wavelengths (in which case Δ λ changes for each wavelength), or do we keep Δ λ constant and have a range of R over the wavelength range? For example, if we have an instrument which ranges from 0.3 – 1.1 microns with a spectral resolution of R = 1000, then either:
      • We let Δ λ change from 0.0003 to 0.0011 microns over the wavelength range
      • We define R = 1000 to be the resolution at the center of the waveset, (0.7 microns), which gives Δ λ = 0.0007 microns. Then for a constant Δ λ, the range of R is [430, 1600] (how it’s done in the current version of synspec)
  • Started working on the astroquery pull request for querying the SVO filter service, borrowing heavily from tynt (no need to reinvent the wheel).

This week’s goals:

  • Finish the synphot SNR pull request by adding to the docs
  • Adding SNR functionality into the other examples
  • Continue working on synspec
    • To check that the binned photometry (aka spectroscopy) is functioning correctly, we can compute the counts for a source spectrum and compare to the counts for the spectroscopized(?) spectrum. Then if they are the same within some error we can rest easier.
    • As it is now, we set a constant Δ λ and allow R to vary accordingly. Add functionality such that the user can either set a constant Δ λ and not worry about what R is, or they can give a specific R and not worry about what Δ λ is.
  • Submit the astroquery PR once it has been reviewed by mentors

Longer term goals:

  • Make 4-5 notebooks which explore different use cases in order to get an idea of how we want to implement any changes or enhancements to synphot:
    1. Model spectrum -> Synthetic photometry
      • Ground-based example: existing APO notebook
      • Space-based example: existing APO notebook + Kepler
    2. Empirical spectrum -> Synthetic photometry
      • Example: Galaxy with high H-alpha emission observed by SDSS
    3. Model spectrum -> synthetic spectroscopy
      • Example: Observations of a G dwarf with CHEOPS at R~1k vs 100k – what count rates do you get?
    4. Empirical spectrum -> synthetic spectroscopy
      • Possible example: Some HII region spectrum -> how many hours to a S/N of X
  • Make a pull request to synphot for the signal-to-noise predictor
  • Create a pull request to Astroquery to enable queries to the filter VO service for filter transmittance curves
  • Think about whether we want to add a higher-level function which takes the telescope params and target and spits out a SNR

Weeks 6 and 7: Synspec and SNR

What I completed these last two weeks:

  • Trying to solve the offset problem for example 2:
    • Do the same analysis for multiple objects. Is the systematic error still there?
      The systematic error persists strongly for both objects, as you can see in the Figure 1 below. So it’s not object-specific, and I think this also means it’s not weather related (assuming the two objects were observed on different nights).
    • Plot bandpasses on top of spectrum
      Figure 2 shows spectra of both objects ontop of the bandpasses, where we’ve removed the u ad z bands because they don’t fully overlap with the spectrum, which resulted in underestimated magnitudes for the u and z bands since synphot has to guess what the spectrum looks like to fill in the gaps. Other than that, I can’t see a reason why the filter models would be responsible for the systematic offset.
    • How systematic is this offset? Plot the fractional error in each band to see.
      Figure 3 shows the fractional error between the synphot fluxes and the observed fluxes. The offset isn’t quite as systematic as I thought, since the three objects are affected to different degrees. Could this mean that the error *is* weather related…?
    • Add a “boring” elliptical galaxy for comparison.
      The offset persists for the elliptical galaxy, and is in the same ballpark as the others in terms of fractional error. I’m still unsure why synphot would be overestimating the photometric flux…
Figure 1. Photometry of two different H-alpha emission galaxies and one elliptical galaxy observed by SDSS, with a persistant systematic offset between synphot fluxes and observed fluxes.
Figure 2. Spectra of three objects in relation to the bandpass models.
Figure 3. Fractional error between the synphot fluxes and the observed fluxes in each band.
  • Wrote a signal-to-noise predictor based on “the CCD equation” in Steve Howell’s “Handbook of CCD Astronomy”, 2000. I would like to make a pull request to implement this into synphot, possibly as a method on Observation, e.g.:
>>> from synphot.observation import Observation
>>> import astropy.units as u
...
>>> observation = Observation(spectrum, bandpass)
>>> countrate = observation.countrate()
>>> counts = countrate * 10 * u.s
>>> snr = observation.snr(counts)
>>> snr
42
  • Also wrote a function which returns the exposure time needed to reach a desired signal to noise ratio.
    • I’m not really sure what to name it, exposure_time_from_snr() seems too long…
    • In the case where the error due to the background estimation and gain are significant, solving for t can not be done analytically, so I use scipy.optimize.fsolve to solve it numerically. However, fsolve doesn’t preserve quantity types, so I have to add some annoying conditional statements to convert t back to a quantity… is there a better way to do this?

This week’s goals:

  • Make a pull request to synphot for the signal-to-noise predictor
  • Begin working on “synspec” (synthetic spectroscopy)
    • Is there an easy way to use synphot for synthetic spectroscopy?
    • See how long each run takes – ideally it should be fast, for people who want to run this on thousands of spectra

Longer term goals:

  • Make 4-5 notebooks which explore different use cases in order to get an idea of how we want to implement any changes or enhancements to synphot:
    1. Model spectrum -> Synthetic photometry
      • Ground-based example: existing APO notebook
      • Space-based example: existing APO notebook + Kepler
    2. Empirical spectrum -> Synthetic photometry
      • Example: Galaxy with high H-alpha emission observed by SDSS
    3. Model spectrum -> synthetic spectroscopy
      • Example: Observations of a G dwarf with CHEOPS at R~1k vs 100k – what count rates do you get?
    4. Empirical spectrum -> synthetic spectroscopy
      • Possible example: Some HII region spectrum -> how many hours to a S/N of X
  • Make a pull request to synphot for the signal-to-noise predictor
  • Create a pull request to Astroquery to enable queries to the filter VO service for filter transmittance curves Is this still on the table given tynt‘s existence?
Oh and it was also the 4th of July sometime in between the coding

Week 5: This time with galaxies

What I completed this last week:

  • Nearly finished the empirical spectrum -> synthetic photometry tutorial for a high H-alpha emission galaxy. This example can be seen here.
    • As you can see from the 1-1 plot at the end of the tutorial, the predicted photometry does not perform as well compared to the 1-1 plot for the first use-case. As Erik (one of my mentors, for anyone reading this who is not already one of my mentors) mentioned, even having a percentage error within 30% is not terrible if you just want a ballpark estimate of the required exposure time.
    • This systematic offset still bugs me though, because it seems like it could be easily solved:
I suppose the mystery continues, for now…
  • Took part in an astropy hack day at the Center for Computational Astrophysics in NYC! Erik and I got to work together in person, I met some astronomers I knew of but never actually met, and I got to witness a bit more of the inner workings of open source development. Also, free food. Lots of free food.

This week’s goals:

  • Try to solve the offset problem for example 2
    • Do the same analysis for multiple objects. Is the systematic error still there?
    • Plot bandpasses on top of spectrum
  • Begin working on synthetic *spectroscopy* (synspec)?
    • Is there an easy way to use synphot for synthetic spectroscopy?
    • See how long each run takes – ideally it should be fast, for people who want to run this on thousands of spectra
  • Begin working on a signal-to-noise predictor (maybe as an eventual method on Observation?)
Me when there’s free food

Longer term goals:

  • Make 4-5 notebooks which explore different use cases in order to get an idea of how we want to implement any changes or enhancements to synphot:
    1. Model spectrum -> Synthetic photometry
      • Ground-based example: existing APO notebook
      • Space-based example: existing APO notebook + Kepler
    2. Empirical spectrum -> Synthetic photometry
      • Example: Galaxy with high H-alpha emission observed by SDSS
    3. Model spectrum -> synthetic spectroscopy
      • Example: Observations of a G dwarf with CHEOPS at R~1k vs 100k – what count rates do you get?
    4. Empirical spectrum -> synthetic spectroscopy
      • Possible example: Some HII region spectrum -> how many hours to a S/N of X
  • Implement a signal-to-noise predictor
  • Create a pull request to Astroquery to enable queries to the filter VO service for filter transmittance curves

Week 4: SkyCalc and SVO and tynt oh my

What I completed this last week:

  • Successfully queried from the SVO filter database for our tutorial’s throughput models. This inspired Brett to write a package called tynt, a “super lightweight package containing approximate transmittance curves for more than five hundred astronomical filters”, which I’ve implemented into the synphot examples.
  • Added a Kepler (i.e. space-based) example to our tutorial. The errors between the synphot counts and the empirical counts are less than 15%.
  • Wrote atmospheric_transmittance.py to model the atmospheric transmittance with the SkyCalc Sky Model Calculator. Parameters can be set when calling get() for a more precise model of the sky, otherwise get() will use the default parameters provided by SkyCalc.

This week’s goals:

  • Begin working on Example 2: Empirical spectrum (like from SDSS/Hubble website) -> Synthetic photometry
  • Begin working on the Astroquery pull request? (mentioned below)

Longer term goals:

  • Make 4-5 notebooks which explore different use cases in order to get an idea of how we want to implement any changes or enhancements to synphot:
    1. Model spectrum -> Synthetic photometry
      • Ground-based example: existing APO notebook
      • Space-based example: existing APO notebook + Kepler
    2. Empirical spectrum (like from SDSS/Hubble website) -> Synthetic photometry
      • Example: Erik’s palomar spectrum + MDM Halpha observations
    3. Model spectrum -> synthetic spectroscopy
      • Example: Observations of a G dwarf with [the space-based mission Brett mentioned called CHEOPS, but don’t worry about that] at R~1k vs 100k – what count rates do you get?
    4. Empirical spectrum -> synthetic spectroscopy
      • Possible example: Some HII region spectrum -> how many hours to a S/N of X
  • Implement a signal-to-noise predictor
  • Create a pull request to Astroquery to enable queries to the filter VO service for filter transmittance curves