This notebook is a case study in working with python and several 3rd-party modules. There are many ways to attack a problem such as this; this is simply one way. The point is to illustrate how you can get existing modules to do the heavy-lifting for you and that visualization is a powerful diagnostic tool. Try not to get caught up in the details of the model; it's quite complex and the point is not to understand all the equations, but the procedure for fitting data to the model (read the citation if you're really interested).
This notebook requires the following modules:
numpy: dealing with arrays of numbers and mathematicsscipy: collection of scientific algorithmsmatplotlib: de-facto plotting modulepandas: module for organizing arrays of number into tablesbokeh: another module for plotting, with emphasis on interactive visualizationThe problem I needed to solve: predict the background sky brightness caused by the moon at a given location in the sky on a given date. This is to help plan observations at the telescope. As with all problems of this type, we need to do several things:
In this case, the data to model is roughly 10 years of photometry from the Carnegie Supernova Project (CSP). Each and every measurement of the flux from a standard star has an associated estimate of the sky background (which must be subtracted from the raw counts of the star). These data were taken over many different times of the month and a many different altitudes, so are ideal for this problem.
Let's start by getting the data. For convenience, this has been included in the data folder and so we can load it up immediately into a pandas dataframe.
In [ ]:
import pandas as pd
data = pd.read_csv('data/skyfit.dat')
We can take a quick look at what's in this DataFrame by printing out the first few rows.
In [ ]:
print(data[0:10])
Let's have a look at the distribution of sky brightnesses to make sure they "make sense". The units should be magnitudes per square-arc-second and be on order of 22 or so, but should be smaller for bright time (full moon). Since we're just doing a quick-look, we can use pandas' built-in histogram plotter.
In [ ]:
%matplotlib inline
data.hist('magsky', bins=50)
As you can see, there is peak near 22 mag/square-arc-sec, as expected, but a broader peak at brighter backgrounds. We expect this is due to moonlight. Something to think about: why would this be bi-modal?
Whatever model we use is going to require knowledge of the moon's position and phase. There are mathematical formulae for this, but we'll use the handy astropy.coordinates module to do all the work for us. First, let's compute the lunar phase for each date in our table. To do this, we need the position of the moon and the sun at these times.
In [ ]:
from astropy.coordinates import get_moon, get_sun
from astropy.time import Time
times = Time(data['jd'], format='jd') # makes an array of time objects
moon = get_moon(times) # makes an array of moon positions
sun = get_sun(times) # makes an array of sun positions
Currently, astropy.coordinates does not have a lunar phase function, so we'll just use the angular separation as a proxy. If the angular separation is 0 degrees, that's new moon, whereas an angular separation of 180 degrees is full moon. Other phases lie in between.
In [ ]:
seps = moon.separation(sun) # angular separation from moon to sun
data['phase'] = pd.Series(seps, index=data.index) # Add this new parameter to the data frame
Now that we have the phase information, let's see if our earlier hypothesis about the moon being a source of background light is valid. We'll plot one versus the other, again using the pandas built-in plotting functionality.
In [ ]:
data.plot.scatter('phase','magsky')
Great! There's a definite trend there, but also some interesting patterns. Remember these are magnitudes per square arc-second, so brighter sky is down, not up. We can also split up the data based on the phase and plot the resulting histograms together. You can run this next snippet of code with different phasecut values to see how they separate out. We use matplotlib's gca function to "get the current axis", allowing us to over-plot two histograms.
In [ ]:
from matplotlib.pyplot import gca,legend
phasecut = 90.
res = data[data.phase>phasecut].hist('magsky', bins=50, label='> %.2f degrees' % phasecut, alpha=0.7)
ax = gca()
res = data[data.phase<phasecut].hist('magsky', ax=ax, bins=50, label='< %.2f degrees' % phasecut, alpha=0.7)
legend(loc='upper left')
Success! It definitely looksl like scattered moonlight is responsible for the bulk of the added sky brightness. But there's also a portion of data where the moon was bright but the sky was still dark. There's more to it than just phase. Now we turn to the task of fitting a model to this.
Turns out that the definitive reference for this was authored by a colleague of mine: Kevin Krisciunas at Texas A&M. His paper can be found at the ADS abstract service: http://adsabs.harvard.edu/abs/1991PASP..103.1033K
You can read the details (lots of empirical formulas, scattering theory, and unit conversions), but the short of it is that we get a predictive model of the sky-brightness as a function of the following variables:
The following diagram shows some of these variables:
Actually, $\alpha$, $\rho$, $Z$, and $Z_m$ are all functions of the date of observations and sky coordinates, which we have already. That leaves $k_x$ and $m_{dark}$ to be determined. Given these variables, the flux from the moon is given by an empirically-determined function that takes into account the fact that the moon is not a perfect sphere:
$$I^* = 10^{-0.4(3.84 + 0.026|\alpha | + 4\times 10^{-9}\alpha^4)}$$This flux is then scattered by angle $\rho$ into our line of sight, contributing to the sky background. The fraction of light scattered into angle $\rho$ is given empirically by:
$$f(\rho) = 10^{5.36}\left[1.06 + \cos^2\rho\right] + 10^{6.15 - \rho/40} $$This just tells us how quickly the sky brightness falls off as we look further away from the moon. We can visualize this by making a 2D array of angles from the center of an image ($\rho$) and comptuing $f(\rho)$.
In [ ]:
from numpy import sqrt, power, cos, sin, indices, pi
from matplotlib.pyplot import imshow, xlabel, ylabel
jj,ii = indices((1024,1024))/1024 # 2D index arrays scaled 0->1
rho = sqrt((ii-0.5)**2 + (jj-0.5)**2)*45.0 # 2D array of angles from center in degrees
f = 10**5.36*(1.06 + (cos(rho*pi/180)**2)) + power(10, 6.15-rho/40)
imshow(f, origin='lower', extent=(-45,45,-45,45), )
xlabel('X angular separation')
ylabel('Y angular separation')
So there's less and less scattered light farther from the moon (at the center). This scattered light is then attenuated by the atmosphere. This attenuation is parametrized by the airmass $X$, the relative amount of atmosphere the light has to penetrate (with $X=1$ for the zenith). Krisciunas & Schaefer (1991) present this formula for the airmass: $X(Z) = \left(1 - 0.96 \sin^2 Z\right)^{-1/2}$. We'll come back to this later. Suffice it to say for the moment that this is an approximation very close to the "infinite slab" model of the atmosphere. Putting it all together, the surface brigthness (in the interesting units of nanoLamberts) from the moon will be:
$$ B_{moon} = f(\rho)I^*10^{-0.4 k_X X(Z_m)}\left[1 - 10^{-0.4k_X X(Z)}\right] $$Let's visualize that first factor, which attenuates the light from the moon. I'll just set $I^*=1$ and $k_X=5$ to make the effect obvious. We'll define the airmass function for later use as well. Let's assume the top of the graph is the zenith ($Z=0$) and the bottom is the horizon ($Z=90$).
In [ ]:
def X(Z):
'''Airmass as a function zenith angle Z in radians'''
return 1./sqrt(1 - 0.96*power(sin(Z),2))
Z = (45 - jj*45)*pi/180. # need in radians
imshow(f*power(10, -0.4*5*X(Z)), origin='lower', extent=(-45,45,-45,45))
So as we get closer to the horizon, there's less moonlight, as it's been attenuated by the larger amount of atmosphere. Lastly, to convert these nanoLamberts into magnitudes per square arc-second, we need the dark (no moon) sky brightness at the zenith, $m_{dark}$, and convert that to nanoLamberts using this formula:
$$ B_{dark} = 34.08\exp (20.7233 - 0.92104 m_{dark})10^{-0.4 k_X (X(Z)-1)}X(Z) $$where we have also corrected for attenuation by the atmosphere and air-glow (which increases with airmass). The final model for observed sky brightness $m_{sky}$ is:
$$ m_{sky} = m_{dark} - 2.5 \log_{10}\left(\frac{B_{moon} + B_{dark}}{B_{dark}}\right) $$Whew! That's a lot of math. But that's all it is, and we can make a python function that will do it all for us.
In [ ]:
from numpy import absolute,exp,log10,arccos,log,where
def modelsky(alpha, rho, kx, Z, Zm, mdark):
Istar = power(10, -0.4*(3.84+0.026*absolute(alpha)+4e-9*power(alpha,4)))
frho = power(10, 5.36)*(1.06 + power(cos(rho),2))+power(10, 6.15-rho*180./pi/40)
Bmoon = frho*Istar*power(10,-0.4*kx*X(Zm))*(1-power(10,-0.4*kx*X(Z)))
Bdark = 34.08*exp(20.723 - 0.92104*mdark)*power(10,-0.4*kx*(X(Z)-1))*X(Z)
return mdark - 2.5*log10((Bmoon+Bdark)/Bdark)
Note that all angles should be entered in radians to work with numpy trig functions.
Now, we just need the final ingredients: $\alpha$, $\rho$, $Z$, and $Z_m$, all of which are computed using astropy.coordinates. The lunar phase angle $\alpha$ is defined as the angular separation between the Earth and Sun as observed on the moon. Alas, astropy can't compute this directly (guess they never thought lunar astronauts would use the software). But since the Earth-moon distance is much less than the Earth-sun distance, this is close enough to 180 degrees minus the angular separation between the moon and sun as observed on Earth (call it $\beta$, which we already computed). See diaram below.
In [ ]:
alpha = (180. - data['phase']) # Note: these need to be in degrees
data['alpha'] = pd.Series(alpha, index=data.index)
Next, in order to compute zenith angles and azimuths, we need to tell the astropy functions where on Earth we are located, since these quantities depend on our local horizon. Luckily, Las Campanas Observatory (LCO) is in astropy's database of locations. We'll also need to create locations on the sky for all our background observations.
In [ ]:
from astropy.coordinates import EarthLocation, SkyCoord, AltAz
from astropy import units as u
lco = EarthLocation.of_site('lco')
fields = SkyCoord(data['RA']*u.degree, data['Decl']*u.degree)
f_altaz = fields.transform_to(AltAz(obstime=times, location=lco))
m_altaz = moon.transform_to(AltAz(obstime=times, location=lco))
rho = moon.separation(fields)*pi/180.0 # angular distance between moon and all fields
Z = (90. - f_altaz.alt.value)*pi/180.0
Zm = (90. - m_altaz.alt.value)*pi/180.0
skyaz = f_altaz.az.value
data['rho'] = pd.Series(rho, index=data.index)
data['Z'] = pd.Series(Z, index=data.index) # radians
data['Zm'] = pd.Series(Zm, index=data.index)
data['skyaz'] = pd.Series(skyaz, index=data.index)
I've added the variables to the Pandas dataFrame as it will help with plotting later. We can try plotting some of these variables against others to see how things look. Let's try a scatter plot of moon/sky separation vs. sky brightness and color the points according to lunar phase. I tried this with the Pandas scatter() and it didn't look that great, so we'll do it with the matplotlib functions directly.
In [ ]:
from matplotlib.pyplot import scatter,colorbar,xlabel, ylabel
scatter(data['rho'], data['magsky'], marker='.', c=data['alpha'], cmap='viridis_r')
xlabel(r'$\rho$', fontsize=16)
ylabel('Sky brightness (mag/sq-arc-sec)', fontsize=12)
colorbar()
There certainly seems to be a trend that the closer to full ($\alpha = 0$, yellow), the brighter the background and the closer the moon is to the field (lower $\rho$), the higher the background. Looks good.
Let's try and fit this data with our model and solve for $m_{dark}$, and $k_x$, the only unknowns in the problem. For this we need to create a dummy function that we can use with scipy's leastsq function. It needs to take a list of parameters (p) as its first argument, followed by any other arguments and return the weighted difference between the model and data. We don't have any weights (uncertainties), so it will just return the differences.
In [ ]:
from scipy.optimize import leastsq
def func(p, alpha, rho, Z, Zm, magsky):
mdark,kx = p
return magsky - modelsky(alpha, rho, kx, Z, Zm, mdark)
We now run the least-squares function, which will find the parameters p which minimize the squared sum of the residuals (i.e. $\chi^2$). leastsq takes as arguments the function we wrote above, func, an initial guess of the parameters, and a tuple of extra arguments needed by our function. It returns the best-fit parameters and a status code. We can print these out, but also use them in our modelsky function to get the prediction that we can compare to the observed data.
In [ ]:
pars,stat = leastsq(func, [22, 0.2], args=(data['alpha'],data['rho'],data['Z'],data['Zm'],data['magsky']))
print(pars)
# save the best-fit model and residuals
data['modelsky']=pd.Series(modelsky(data['alpha'],data['rho'],pars[1],data['Z'],data['Zm'],pars[0]), index=data.index)
Now we want to compare the model to the data. The typical way to do this when you have many variables is to plot residuals versus the variables and see how things look. But a cool package called bokeh gives a very powerful diagnostic tool: linking graphs so that selecting objects in one will select the corresponding objects in all other graphs that share the same dataset. This is why we've been adding our variables to the pandas dataFrame, data. Try selecting different regions of the upper-left panel (which compares the model with the observations).
In [ ]:
from bokeh.plotting import figure
from bokeh.layouts import gridplot
from bokeh.io import show,output_notebook
from bokeh.models import ColumnDataSource
output_notebook()
source = ColumnDataSource(data)
TOOLS = ['box_select','lasso_select','reset','box_zoom','help']
vars = [('modelsky','magsky'),('alpha','rho'),('alpha','Zm'),
('jd','alpha'),('Z','Zm'),('RA','Decl')]
plots = []
for var in vars:
s = figure(tools=TOOLS, plot_width=300, plot_height=300)
s.circle(*var, source=source, selection_color='red')
s.xaxis.axis_label = var[0]
s.yaxis.axis_label = var[1]
plots.append(s)
plots[0].line([17.8,22.3],[17.8,22.3], line_color='orangered')
p = gridplot([plots[0:3],plots[3:]])
show(p)
With a little data exploring, it's pretty obvious that the majority of the outlying points comes from observations when the moon is relatively full and very low (or even below) the horizon. The reason is that the airmass formula that we implemented above has a problem with $Zm > \pi/2$. To see this, we can simply plot X(Z) as a function of 'Z':
In [ ]:
from numpy import linspace
from matplotlib.pyplot import plot, xlabel, ylabel,ylim
Z = linspace(0, 3*pi/4, 100) # make a range of Zenith angles
plot(Z*180/pi, X(Z), '-')
xlabel('Zenith angle (degrees)')
ylabel('Airmass')
So the airmass (amount of air the light travels through) increases as you get to the horizon ($Z=90^\circ$), but then decreases. That's not the right behaviour. Can you think of a way to fix this problem? Try it out. Just go back to the cell above where X(Z) is defined and change it. Then select Cell -> Run All Below from the menu so see how the results change.
At this point you might be feeling overwhelmed. How did I know which modules to use? How did I know how to use them? The answer: Google, ADS, and 20+ years (eek!) of experience coding in Python. I also neglected to show all the dead-ends and mistakes I made on the way to getting the final solution, all the emails I sent to Kevin asking about the details of his paper, and advice I got from Shannon about using Bokeh.
Before you start tackling a particular problem it's well worth your time to research whether there is already a solution "out there" that you can use or modify for your use. It has never been so easy to do this, thanks to search engines (Google, et al.), data/software catalogs (PyPI, et al.), discussion groups (Stackoverflow, et al.) and even social media (python users in astronomy facebook group, etc). And your friendly neighborhood python experts are there to make helpful suggestions.
Don't re-invent the wheel, but improve it by all means.