Pretty Plots - 2D Histogram with 1D Histograms on Axes

Note: This post has been cross-posted to AstroBetter with a slightly more readable version of the code.  You might want to check it out there.

In a long list of plots that Martin wants me to add to the paper draft I sent out last week, the following:

Is it possible to show the histograms projected along each axis in addition to the 2D density? I know people do this in IDL frequently, I'm not sure how to do this in matplotlib. If it's possible we could play a similar game with the L-z figure. The 1D histograms contain lots of useful information, including how significant our clustering detection is in each bin.

Ask and you shall receive Martin. Below is the "money plot" (the temperature 2D histogram which shows amplitudes the bright versus dim correlation function) with histograms on the side of either axis.

I based this plot on code from here.

There are some cool features that I'll describe in the comments below. I think this plot rocks.

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.ticker import NullFormatter
def makeTempHistogramPlot(xdata,ydata,rexp,filename=None,xlims=-99, ylims =-99 , \
nxbins = 50,nybins=50, bw=0, nbins=100,contours=1,sigma=1,line=1):
#bw = 0 for color, = 1 for black and white
#line = 0 for no line, =1 for line
#sigma = 1 for display % below line, =0 for not
#contours = 1 for display 1,2,3 sigma contours, = 0 for not.

# Define the x and y data
x = xdata
y = ydata

# Set up default x and y limits
if (xlims == -99): xlims = [0,max(x)]
if (ylims == -99): ylims = [0,max(y)]

# Set up your x and y labels
xlabel = '$\mathrm{Your\\ X\\ Label}$'
ylabel = '$\mathrm{Your\\ X\\ Label}$'
mtitle = ''

# Define the locations for the axes
left, width = 0.12, 0.55
bottom, height = 0.12, 0.55
bottom_h = left_h = left+width+0.02

# Set up the geometry of the three plots
rect_temperature = [left, bottom, width, height] # dimensions of temp plot
rect_histx = [left, bottom_h, width, 0.25] # dimensions of x-histogram
rect_histy = [left_h, bottom, 0.25, height] # dimensions of y-histogram

# Set up the size of the figure
fig = plt.figure(1, figsize=(9.5,9))

# Make the three plots
axTemperature = plt.axes(rect_temperature) # temperature plot
axHistx = plt.axes(rect_histx) # x histogram
axHisty = plt.axes(rect_histy) # y histogram

# Remove the inner axes numbers of the histograms
nullfmt = NullFormatter()
axHistx.xaxis.set_major_formatter(nullfmt)
axHisty.yaxis.set_major_formatter(nullfmt)

# Find the min/max of the data
xmin = min(xlims)
xmax = max(xlims)
ymin = min(ylims)
ymax = max(y)

# Make the 'main' temperature plot
xbins = linspace(start = 0, stop = xmax, num = nxbins)
ybins = linspace(start = 0, stop = ymax, num = nybins)
xcenter = (xbins[0:-1]+xbins[1:])/2.0
ycenter = (ybins[0:-1]+ybins[1:])/2.0
aspectratio = 1.0*(xmax - 0)/(1.0*ymax - 0)
H, xedges,yedges = N.histogram2d(y,x,bins=(ybins,xbins))
X = xcenter
Y = ycenter
Z = H

# Plot the temperature data
if(bw): cax = axTemperature.imshow(H, extent=[xmin,xmax,ymin,ymax], \
interpolation='nearest', origin='lower',aspect=aspectratio, cmap=cm.gist_yarg)
else : cax = axTemperature.imshow(H, extent=[xmin,xmax,ymin,ymax], \
interpolation='nearest', origin='lower',aspect=aspectratio)

# Plot the temperature plot contours
if(bw): contourcolor = 'black'
else: contourcolor = 'white'

if (contours==0):
print ''
elif (contours==1):
xcenter = N.mean(x)
ycenter = N.mean(y)
ra = N.std(x)
rb = N.std(y)
ang = 0
X,Y=ellipse(ra,rb,ang,xcenter,ycenter)
axTemperature.plot(X,Y,"k:",ms=1,linewidth=2.0)
axTemperature.annotate('$1\\sigma$', xy=(X[15], Y[15]), xycoords='data',xytext=(10, 10), textcoords='offset points',horizontalalignment='right', verticalalignment='bottom',fontsize=25)
X,Y=ellipse(2*ra,2*rb,ang,xcenter,ycenter)
axTemperature.plot(X,Y,"k:",color = contourcolor,ms=1,linewidth=2.0)
axTemperature.annotate('$2\\sigma$', xy=(X[15], Y[15]), xycoords='data',xytext=(10, 10), textcoords='offset points',horizontalalignment='right', verticalalignment='bottom',fontsize=25, color = contourcolor)
X,Y=ellipse(3*ra,3*rb,ang,xcenter,ycenter)
axTemperature.plot(X,Y,"k:",color = contourcolor, ms=1,linewidth=2.0)
axTemperature.annotate('$3\\sigma$', xy=(X[15], Y[15]), xycoords='data',xytext=(10, 10), textcoords='offset points',horizontalalignment='right', verticalalignment='bottom',fontsize=25, color = contourcolor)
else:
xcenter = N.mean(x)
ycenter = N.mean(y)
ra = N.std(x)
rb = N.std(y)
ang = contours*N.pi/180.0
X,Y=ellipse(ra,rb,ang,xcenter,ycenter)
axTemperature.plot(X,Y,"k:",ms=1,linewidth=2.0)
axTemperature.annotate('$1\\sigma$', xy=(X[15], Y[15]), xycoords='data', xytext=(10, 10), textcoords='offset points',horizontalalignment='right', verticalalignment='bottom',fontsize=25)
X,Y=ellipse(2*ra,2*rb,ang,xcenter,ycenter)
axTemperature.plot(X,Y,"k:",ms=1,linewidth=2.0, color = contourcolor)
axTemperature.annotate('$2\\sigma$', xy=(X[15], Y[15]), xycoords='data', xytext=(10, 10), textcoords='offset points',horizontalalignment='right', verticalalignment='bottom',fontsize=25, color = contourcolor)
X,Y=ellipse(3*ra,3*rb,ang,xcenter,ycenter)
axTemperature.plot(X,Y,"k:",ms=1,linewidth=2.0, color = contourcolor)
axTemperature.annotate('$3\\sigma$', xy=(X[15], Y[15]), xycoords='data', xytext=(10, 10), textcoords='offset points',horizontalalignment='right', verticalalignment='bottom',fontsize=25, color = contourcolor)

#Plot the % below line
belowline = 1.0*size(where((x - y) > 0.0))/size(x)*1.0*100
if(sigma): axTemperature.annotate('$%.2f\%%\mathrm{\\ Below\\ Line}$'%(belowline), xy=(xmax-100, ymin+3),fontsize=20, color = contourcolor)

#Plot the axes labels
axTemperature.set_xlabel(xlabel,fontsize=25)
axTemperature.set_ylabel(ylabel,fontsize=25)

#Make the tickmarks pretty
ticklabels = axTemperature.get_xticklabels()
for label in ticklabels:
label.set_fontsize(18)
label.set_family('serif')

ticklabels = axTemperature.get_yticklabels()
for label in ticklabels:
label.set_fontsize(18)
label.set_family('serif')

#Plot the line on the temperature plot
if(line): axTemperature.plot([-1000,1000], [-1000,1000], 'k-', linewidth=2.0, color = contourcolor)

#Set up the plot limits
axTemperature.set_xlim(xlims)
axTemperature.set_ylim(ylims)

#Set up the histogram bins
xbins = N.arange(xmin, xmax, (xmax-xmin)/nbins)
ybins = N.arange(ymin, ymax, (ymax-ymin)/nbins)

#Plot the histograms
if (bw):
axHistx.hist(x, bins=xbins, color = 'silver')
axHisty.hist(y, bins=ybins, orientation='horizontal', color = 'dimgray')
else:
axHistx.hist(x, bins=xbins, color = 'blue')
axHisty.hist(y, bins=ybins, orientation='horizontal', color = 'red')

#Set up the histogram limits
axHistx.set_xlim( 0, max(x) )
axHisty.set_ylim( 0, max(y))

#Make the tickmarks pretty
ticklabels = axHistx.get_yticklabels()
for label in ticklabels:
label.set_fontsize(12)
label.set_family('serif')

#Make the tickmarks pretty
ticklabels = axHisty.get_xticklabels()
for label in ticklabels:
label.set_fontsize(12)
label.set_family('serif')

#Cool trick that changes the number of tickmarks for the histogram axes
axHisty.xaxis.set_major_locator(MaxNLocator(4))
axHistx.yaxis.set_major_locator(MaxNLocator(4))

if(filename):
savefig(filename + '.eps',format = 'eps', transparent=True)
savefig(filename + '.pdf',format = 'pdf', transparent=True)
savefig(filename + '.png',format = 'png', transparent=True)

return 0

Temperature Plots in Python

I figured out how to do temperature plots in python. The only problem with below is that the numbers on the axis are wrong. I if I set 'extent' to be the actual extent of the axis (in this case x = [0.5,1], y = [-28,-18]) then it plots a very skinny, tall plot. So I had to correct the axis by hand later.

import matplotlib
import numpy as np
import matplotlib.cm as cm
import matplotlib.mlab as mlab
import matplotlib.pyplot as plt
import matplotlib.pyplot as plt

redbins = linspace(start = 0.5, stop = 1.0, num = 51)
rcenter = (redbins[0:-1]+redbins[1:])/2.0
mbins = linspace(start = -28, stop = -18, num = 51)
mcenter = (mbins[0:-1]+mbins[1:])/2.0


H, xedges,yedges = N.histogram2d(qsoM[dr7cut],qsored[dr7cut],bins=(mbins,redbins))
X = rcenter
Y = mcenter
Z = H

#plt.imshow(H, extent=[-3,3,-3,3], interpolation='bilinear', origin='lower')
#plt.colorbar()


plt.imshow(H, extent=[-3,3,-3,3], interpolation='nearest', origin='lower')
plt.colorbar()

xlabel('redshift (z)')
ylabel('absolute magnitude (i-band)')
title('SDSS DR7 Quasar Density')

The above code is also in the following log file: ../logs/110614log.py

Changing QSO Catalog Inputs

I'm trying to figure out how to modify the inputs (currently SDSS DR5 quasars with good photometry) into Joe's Monte Carlo. Here is what I have figured out so far (all the below files are on riemann):

hiz_kde_numerator.pro ;main program to generate the QSO Catalog
(calls) →
qso_fakephoto.pro ;generates simulated redshift and i-magnitude based on luminosity function and DR5 inputs
(calls) →
qso_photosamp.pro ;This generates the training set photometry
(calls) →
sdss_read_data.pro ;This reads in SDSS data it currently has the following 'DR5' call
qsos = sdss_read_data('DR5', Z_MIN = Z_MIN, Z_MAX = Z_MAX)

Here is a color-color redshift temperature plot of these DR5 QSOs:


I'm going to try to add in the BOSS QSOs with co-added photometry, especially in the redshift range z > 2.0 where we seem to be missing QSOs:


I found that there are 1,973 BOSS quasars that have a redshift above 2.0 for which we also have coadded photometry. Here they are plotted:


While it is a little more filled in in the u-g = 0.8 space, comparing it to all the BOSS quasars (regardless of how good their photometry is) below:

Comparing Quasars

Here are the color-color diagrams with redshift temperature plots for the current QSO Catalog (in the BOSS redshift range) and the BOSS 3PC Quasars (same temperature-z map as before). I think this shows pretty well that something is wrong with the QSO Catalog's color distribution. These both have the same number of quasars, in the same redshift range:

Their redshift distributions are not that different so this is really a matter of the input QSOs not properly representing the spread of the possible colors, I believe:

White is BOSS QSOs, Green is the QSO Catalog.
The histograms are scaled as a percentage in each bin.

I want to add in BOSS QSOs to Joe's Monte Carlo. We can do this by just adding in the "chunk 1" objects where we have co-added photometry from Stripe-82, or we can apply the same cuts in terms of brightness as we did to the DR5 catalog.

Plots for Joe (and Brandon)

Aside - A couple weeks ago Brandon Basso complained that my buzzing wasn't frequent enough. He requested more plots and perhaps an update on Adam's research. Well ask and you shall receive! Brandon this plot-heavy post is dedicated to you.

I am trying to understand why we seem to be missing objects around u-g = 0.8. Joe suggested I make the color plots in different redshift bins. So here we go. Warning, there are a lot of plots.

First, let me show you a temperature plots of the different luminosity functions. Below are color-color diagrams where the color of the points changes as a function of redshift in the following way:

dark red = 1.0 < z < 1.2
red = 1.2 < z < 1.4
orange red = 1.4 < z < 1.6
orange = 1.6 < z < 1.8
gold = 1.8 < z < 2.0
lawn green = 2.0 < z < 2.2
lime green =2.2 < z < 2.4
dark green = 2.4 < z < 2.6
teal = 2.6 < z < 2.8
dodger blue = 2.8 < z < 3.0
royal blue = 3.0 < z < 3.2
blue = 3.2 < z < 3.4
navy= 3.4 < z < 3.6
dark slate blue = 3.6 < z < 3.8
dark orchid = 3.8 < z < 4.0

As a reminder here are the luminosity functions I am using:

The above plot shows several luminosity functions.

Green is Richards 06, cyan is Jiang, purple is Jiang Combined (what we used in previous QSO Catalog, white is a Richard/Jiang average (suggested by Myers).

Here are the color-color redshift temperature plots for the QSO Catalogs generated by the above luminosity functions (click plots on below to make larger):

Also I have made color-color plots binned by redshift separately. I've done them for all of the luminosity functions, but they all look pretty much the same, so I'll just post them here for the Richards function (green luminosity function above), the redshift range is in the main title of the plots:

Here are some histograms of various quantities for the different luminosity functions the color scheme is same as above, Green is Richards 06, cyan is Jiang, purple is Jiang Combined (what we used in previous QSO Catalog, white is a Richard/Jiang average:

redshift

Color-band psffluxes


Color-band magnitudes

Color-color plot u-g/g-r

Hopefully that is enough plots for Joe to see what is going on with these QSO catalogs and Brandon to be satisfied.

Oh and here is an update on Adam's research (check out his shirt).

Temperature Guide

dark red = 1.0 < z < 1.2
red = 1.2 < z < 1.4
orange red = 1.4 < z < 1.6
orange = 1.6 < z < 1.8
gold = 1.8 < z < 2.0
lawn green = 2.0 < z < 2.2
lime green =2.2 < z < 2.4
dark green = 2.4 < z < 2.6
teal = 2.6 < z < 2.8
dodger blue = 2.8 < z < 3.0
royal blue = 3.0 < z < 3.2
blue = 3.2 < z < 3.4
navy= 3.4 < z < 3.6
dark slate blue = 3.6 < z < 3.8
dark orchid = 3.8 < z < 4.0