{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Welcome to Advanced Plotting!" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This is not meant to be a comprehensive guide to plotting, but more of a primer to give you a taste of the fancy things you can to with matplotlib. The examples in here were created by IceCube and RNO grad students as part of their research.\n", "\n", "More examples like this can be found in the advanced plotting folder or at https://icecube.wisc.edu/~icecube-bootcamp/bootcamp2019/advanced_plotting/\n", "\n", "This barely scratches the surface of what you can do. The matplotlib documentation has many good examples https://matplotlib.org/gallery/index.html, and in general is a great reference.\n", "\n", "If you're interested in making plots pretty, readable, and colorblind friendly then you might want to check out: https://seaborn.pydata.org/, https://gka.github.io/palettes/, http://colorbrewer2.org/" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## matplotlib\n", "\n", "### `matplotlib.pyplot` is a state machine interface to `matplotlib` similar to MATLAB\n", "In this interface there is a global state that contains information about our plot and subsequent calls to `pyplot` modify that state.\n", "\n", "### `matplotlib` is python the plotting library most often used by IceCubers. Its interface is object oriented.\n", "\n", "In this figure we can see the hierarchy of some of the high level objects that make up a plot. Each of these is represented as a python object that can be manipulated.\n", "\n", "\n", "\n", "There are however many more parts to a plot. These are all represented by python objects\n", "\n", "" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib\n", "import matplotlib.figure\n", "%matplotlib inline\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Plotting with the state machine interface\n", "Many simple plots can be made with the state-machine interface" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Create a little data for us to plot\n", "coin_16p5 = np.array([1.000000, 0.592732, 0.214104, 0.067033, 0.014277])\n", "coin_18p0 = np.array([1.000000, 0.723583, 0.386608, 0.246627, 0.125687])\n", "\n", "###################### Plot surface coincidence fractions #####################\n", "xlabels = [\"Surface\", \"15 m\", \"60 m\", \"100 m\", \"200 m\"]\n", "width = 0.4 \n", "sep = 0.2 \n", "n_bars = len(coin_16p5)\n", "\n", "# Set the size of the plot\n", "plt.figure(figsize=(8, 5)) \n", "\n", "# Set up horizontal dotted grid lines\n", "plt.grid(axis='y', ls=':')\n", "\n", "# Plot two sets of bar graphs, one for each data set\n", "plt.bar(np.arange(n_bars)-sep, coin_16p5,\n", " width=width, alpha=0.8, label='1e16.5 eV')\n", "plt.bar(np.arange(n_bars)+sep, coin_18p0,\n", " width=width, alpha=0.8, label='1e18.0 eV')\n", "# Note that the first argument of plt.bar specifies the x-position of the center of the bar\n", "# and that sep == width/2.0. This is why these two sets of bars are plotted side-by-side\n", "\n", "# Add a legend to the plot with the given labels\n", "plt.legend(fontsize=10)\n", "\n", "# Set the axes ticks labels and fontsize\n", "plt.xticks(ticks=np.arange(n_bars), labels=xlabels, fontsize=12)\n", "plt.yticks(fontsize=12)\n", "\n", "# Give the plot a title\n", "plt.title(\"Surface Coincidence Fraction\", fontsize=14)\n", "# Show the plot\n", "plt.show()\n", "# Or we could save the plot as a PDF so it scales nicely in presentations\n", "# plt.savefig('station_surface_coincidence_fractions.pdf')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Plotting with the object oriented interface\n", "These same plots can of course be created" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Create a little data for us to plot\n", "coin_16p5 = np.array([1.000000, 0.592732, 0.214104, 0.067033, 0.014277])\n", "coin_18p0 = np.array([1.000000, 0.723583, 0.386608, 0.246627, 0.125687])\n", "\n", "###################### Plot surface coincidence fractions #####################\n", "xlabels = [\"Surface\", \"15 m\", \"60 m\", \"100 m\", \"200 m\"]\n", "width = 0.4\n", "sep = 0.2 \n", "n_bars = len(coin_16p5)\n", "# Set the size of the plot\n", "fig = plt.figure(figsize=(8, 5))\n", "# this method of instatiating the figure ensures that it is managed\n", "# by the pyplot interface and things like plt.show() will work\n", "# we also could have used:\n", "# fig = matplotlib.figure.Figure(figsize=(8, 5))\n", "# if we wanted to do things in a purely object oriented way\n", "\n", "# Add a subplot\n", "ax = fig.add_subplot(1, 1, 1)\n", "\n", "# Set up horizontal dotted grid lines\n", "ax.grid(axis='y', ls=':')\n", "\n", "# Plot two sets of bar graphs, one for each data set\n", "ax.bar(np.arange(n_bars)-sep, coin_16p5,\n", " width=width, alpha=0.8, label='1e16.5 eV')\n", "ax.bar(np.arange(n_bars)+sep, coin_18p0,\n", " width=width, alpha=0.8, label='1e18.0 eV')\n", "# Note that the first argument of ax.bar specifies the x-position of the center of the bar\n", "# and that sep == width/2.0. This is why these two sets of bars are plotted side-by-side\n", "\n", "# Add a legend to the plot with the given labels\n", "ax.legend(fontsize=10)\n", "\n", "# Set the x tick positions to 0 - # of bars\n", "ax.set_xticks(np.arange(n_bars)) \n", "\n", "# Set the labels to be the strings proveded above in xlabels\n", "ax.set_xticklabels(xlabels)\n", "\n", "# Set the fontsize of all axis labels\n", "for item in (ax.get_xticklabels() + ax.get_yticklabels()):\n", " item.set_fontsize(12)\n", "\n", "# Give the plot a title\n", "ax.set_title(\"Surface Coincidence Fraction\", fontsize=14)\n", "\n", "# Show the plot\n", "plt.show()\n", "# Or we could save the plot as a PDF so it scales nicely in presentations\n", "# plt.savefig('station_surface_coincidence_fractions.pdf')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Which interface should I use?\n", "Both work perfectly fine for most tasks, but for more complicated plots the object oriented interface is the better choice.\n", "\n", "For the example below it makes the most sense for us to use the object oriented interface." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Example number 1\n", "This example plots the inverse of the bayes factor as a function of the prompt normalization.\n", "\n", "We also plot regions of exclusion at the substantial, strong, very strong, and decisive levels according to Jefferys' scale, and some arrows to indicate which region is excluded vs. allowed.\n", "\n", "The end result should look like this:\n", "\n", "This plot uses a second axis with invisible elements placed on top of the first axis. This allows us to plot using two different scales simultaneously. In particular the arrows we draw on our second axis with linear scaling would be distorted if we plotted them on the first axis which has log scaling on the y-axis." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Set up our imports\n", "We need `numpy` for some data manipulation, `matplotlib` for plotting, and a bunch of `matplotlib` includes for style manipulation" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib\n", "import matplotlib.style\n", "import matplotlib.cm\n", "import matplotlib.pyplot as plt\n", "from matplotlib import rc, rcParams" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Plot settings\n", "Here we change some of the settings to make the plot look nice. These same settings would usually be stored in a style sheet file but instead we modify the settings inline.\n", "\n", "The style sheet file would look like:\n", "\n", "style.mpl\n", "```\n", "figure.figsize : 5, 5 # figure size in inches\n", "savefig.dpi : 600 # figure dots per inch\n", "\n", "font.size: 18\n", "font.family: serif\n", "text.usetex: True\n", "axes.linewidth: 1.0\n", "font.weight: bold\n", "\n", "patch.linewidth: 2\n", "lines.linewidth: 2\n", "xtick.labelsize: medium\n", "ytick.labelsize: medium\n", "xtick.minor.visible: True # visibility of minor ticks on x-axis\n", "ytick.minor.visible: True # visibility of minor ticks on y-axis\n", "xtick.major.size: 8 # major tick size in points\n", "xtick.minor.size: 4 # minor tick size in points\n", "ytick.major.size: 8 # major tick size in points\n", "ytick.minor.size: 4 # minor tick size in points\n", "xtick.major.width: 1.2\n", "xtick.minor.width: 1.2 \n", "ytick.major.width: 1.2 \n", "ytick.minor.width: 1.2\n", "```\n", "\n", "After we finish making the plot we will reset these settings to their defaults" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Choose a whole bunch of settings\n", "rc('figure',figsize=(5,5))\n", "rc('savefig',dpi=600)\n", "rc('font',size=18)\n", "rc('font',weight='bold')\n", "rc('font',family='serif')\n", "rc('text',usetex=True)\n", "rc('patch', linewidth=2)\n", "rc('lines', linewidth=2)\n", "rc('xtick', labelsize='medium')\n", "rc('ytick', labelsize='medium')\n", "rc('xtick.minor', visible=True)\n", "rc('ytick.minor', visible=True)\n", "rc('xtick.major', size=8)\n", "rc('xtick.minor', size=4)\n", "rc('ytick.major', size=8)\n", "rc('ytick.minor', size=4)\n", "rc('xtick.major', width=1.2)\n", "rc('xtick.minor', width=1.2 )\n", "rc('ytick.major', width=1.2 )\n", "rc('ytick.minor', width=1.2 )" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Fake data\n", "We need some data to plot, so let's make some" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#### Prepare the data ####\n", "\n", "# Define the function our fake data will follow\n", "a = 0.0102757\n", "b = -0.120551\n", "c = 0.0\n", "fake_data_fit = lambda x: 10.0**(a*x**2 + b*x + c)\n", "\n", "# Define the points we are going to sample at\n", "prompt_minimum = 0\n", "prompt_maximum = 25\n", "samples = 100+1\n", "\n", "# Get our fake data\n", "fake_data_x = np.linspace(prompt_minimum, prompt_maximum, samples)\n", "fake_data_y = fake_data_fit(fake_data_x)\n", "\n", "prompts = fake_data_x\n", "bayes_factor = fake_data_y\n", "#print('Prompt normalizations:', prompts)\n", "#print('Bayes factors:', bayes_factor)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Check out what our fake data looks like\n", "Before we write the whole plotting script, let's take a look at the data.\n", "\n", "The code below will make a quick plot for us.\n", "```\n", "plt.plot(prompts, bayes_factor)\n", "plt.yscale('log')\n", "plt.show()\n", "```" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Plotting code goes here\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#### Set up some parameters for plotting ####\n", "\n", "# Choose the plasma color map\n", "cm = matplotlib.cm.get_cmap('plasma')\n", "\n", "# Define the bounds of the data and axis\n", "min_x = max(np.amin(prompts), 0)\n", "max_x = min(np.amax(prompts), 25)\n", "min_y = min(np.amin(bayes_factor), 1.0)\n", "if min_y < 1.0:\n", " min_y = min(np.amin(bayes_factor), 0.9)\n", "max_y = min(np.amax(bayes_factor), 10**4)\n", "\n", "# Choose an alpha value for the fill regions\n", "alpha = 0.3\n", "\n", "# Define labels and values for the Bayes factors\n", "factor_names = ['barely worth mentioning', 'substantial', 'strong', 'very strong', 'decisive']\n", "factors = [1, 10.0**0.5, 10.0, 10.0**(3.0/2.0), 10.0**2]\n", "\n", "# Compute the prompt limit for each Bayes factor using interpolation\n", "# in log space for the y-axis and linear space for the x-axis\n", "limits = []\n", "for factor, name in zip(factors[1:], factor_names[1:]):\n", " mask = bayes_factor > factor\n", " i = np.argmax(mask)\n", " y0 = np.log(bayes_factor[i-1])\n", " x0 = prompts[i-1]\n", " y1 = np.log(bayes_factor[i])\n", " x1 = prompts[i]\n", " y2 = np.log(factor)\n", " x2 = (x1*y0-x0*y1+x0*y2-x1*y2) / (y0-y1)\n", " limit = x2\n", " limits.append(limit)\n", "\n", "# Extend the limit list beyond the plotting range for easy plotting of the filled regions\n", "limits = [min_x-abs(max_x - min_x)] + limits + [max_x+abs(max_x - min_x)]\n", "factor_names = [None] + factor_names[1:] + [None]\n", "\n", "# Define functions that transform from the normal data space to the space in newax\n", "x_rescale = lambda x: (x - min_x) / (max_x - min_x)\n", "y_rescale = lambda y: (np.log10(y) - np.log10(min_y)) / (np.log10(may_y) - np.log10(min_y))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Check out what our limits look like on the plot\n", "```\n", "for i in range(len(limits))[1:-1]:\n", " plt.axvline(limits[i], 0, 1, linestyle='solid', color='red')\n", "plt.plot(prompts, bayes_factor)\n", "plt.yscale('log')\n", "plt.show()\n", "```" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Plotting code goes here\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Try to draw an arrow on the same plot\n", "```\n", "plt.arrow(0,10,10,0, head_length=2, head_width=10)\n", "for i in range(len(limits))[1:-1]:\n", " plt.axvline(limits[i], 0, 1, linestyle='solid', color='red')\n", "plt.plot(prompts, bayes_factor)\n", "plt.yscale('log')\n", "plt.show()\n", "print('Notice how the log scale distorts the arrow')```" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Plotting code goes here\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Make the plot ####\n", "```\n", "# Create the figure and axis (single axis)\n", "fig, ax = plt.subplots()\n", "\n", "# Create a second axis with the same position a the previous one\n", "newax = fig.add_axes(ax.get_position())\n", "\n", "# Make its background invisible\n", "newax.patch.set_visible(False)\n", "\n", "# Make each axis itself invisible\n", "newax.yaxis.set_visible(False)\n", "newax.xaxis.set_visible(False)\n", "\n", "# Make the spines invisible\n", "for spinename, spine in newax.spines.items():\n", " spine.set_visible(False)\n", "\n", "# Set the data limits of newax to be [(0,1),(0,1)] for simplicity\n", "newax.set_xlim((0,1))\n", "newax.set_ylim((0,1))\n", "\n", "# Plot the filled regions\n", "fills = []\n", "names = []\n", "for i in range(1, len(limits)-1):\n", " # Get the color to use for the solid lines\n", " color = cm(1.0 - float(i+1)/float(len(limits)))\n", " \n", " # Get the semi-transparent verison of that color to use for the filled region\n", " transparent_color = tuple(list(color)[:-1] + [alpha])\n", " \n", " # Define the fill limits in the newax space\n", " fill_min = x_rescale(limits[i])\n", " fill_max = x_rescale(limits[i+1])\n", " \n", " # Define the arrow length to be 75% of the the filled region width\n", " arrow_length = (fill_max - fill_min) * 0.75\n", " \n", " # Plot the filled region\n", " fill = newax.fill_between([fill_min, fill_max], 2, -1, edgecolor=color, facecolor=transparent_color)\n", " \n", " # Draw the arrow\n", " newax.arrow(fill_min, 0.5,\n", " min(arrow_length, (x_rescale(max_x)-fill_min)*0.75), 0,\n", " linestyle='solid', color=color,\n", " length_includes_head=True, head_length=0.01, head_width=0.01)\n", " \n", " # Draw a vertical line to cover the boundary between the filled regions\n", " newax.axvline(fill_min, 0, 1, linestyle='solid', color=color)\n", " \n", " # Keep these for later when we plot the legend\n", " fills.append(fill)\n", " names.append(factor_names[i])\n", "\n", "# Plot the Bayes factor\n", "ax.plot(prompts, bayes_factor, color=cm(0.0), zorder=4)\n", "\n", "# Draw a dashed line\n", "ax.axhline(1.0, 0, 1, color='#8e8e8e', linestyle='dashed', linewidth=1, zorder=5)\n", "\n", "# Set up the axis\n", "ax.set_xlim(min_x, max_x)\n", "ax.set_ylim(min_y, max_y)\n", "ax.set_xlabel(r'$\\Phi_\\texttt{prompt}[\\textmd{BERSS}]$')\n", "ax.set_ylabel(r'$1/\\mathcal{B}$')\n", "ax.set_yscale('log')\n", "\n", "# Call tight_layout to arrange things nicely\n", "fig.tight_layout()\n", "\n", "# Ensure that newax is in the same place as ax\n", "newax.set_position(ax.get_position())\n", "\n", "# Draw the legend\n", "newax.legend(fills, [f[0].upper() + f[1:] for f in names], frameon=True, title='Evidence for\\nexclusion',fontsize=12)\n", "\n", "# Show the figure\n", "plt.show()\n", "#fig.savefig('prompt_test.pdf')\n", "#fig.savefig('prompt_test.png')\n", "```" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Plotting code goes here\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Reset all those settings we changed before we make any more plots\n", "matplotlib.rcParams.update(matplotlib.rcParamsDefault)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Example number 2" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this next example we are going to make a plot for a gravitational wave neutrino followup analysis.\n", "\n", "This plot will have several components: heatmap on the sky for the GW location probability, a contour showing the 90% containment of the GW directional probability, and neutrino directions with error contours.\n", "\n", "In the end it should look something like this:\n", "" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For this example we're going to need two additional packages to help us with some of the heavy lifting. We'll need to install `healpy` and `meander`\n", "\n", "If you haven't already, you'll need to install these with the following commands:\n", "```\n", "pip3 install --user healpy\n", "pip3 install --user meander```\n", "\n", "### Let's take care of some imports now" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import healpy as hp\n", "import matplotlib\n", "import matplotlib.pyplot as plt \n", "import meander #pip install this to draw contours around skymap\n", "\n", "from matplotlib import cm" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### We're also going to need some data" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Fits file containing Gravitational wave skymap\n", "fitsFile = 'https://gracedb.ligo.org/api/superevents/S190521r/files/bayestar.fits.gz'\n", "\n", "### Read map and get probabilities\n", "probs = hp.read_map(fitsFile)\n", "nside = hp.pixelfunc.get_nside(probs)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Let's also choose some style options" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Choose color map and set background to white\n", "cmap = cm.YlOrRd\n", "cmap.set_under(\"w\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### First we'll try to plot the heatmap by itself\n", "```\n", "# Plot GW skymap in Mollweide projection\n", "hp.mollview(probs,cbar=True,unit=r'Probability',min=0,max=3e-5,rot=180,cmap=cmap)\n", "hp.graticule() # Set grid lines\n", "```" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Plotting code goes here\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Then we'll add some neutrinos in the sky and plot them as well\n", "```\n", "# Plot GW skymap in Mollweide projection\n", "hp.mollview(probs,cbar=True,unit=r'Probability',min=0,max=3e-5,rot=180,cmap=cmap)\n", "hp.graticule() # Set grid lines\n", "\n", "# Let's generate 5 random neutrinos on the sky and plot them\n", "theta = np.random.uniform(0,np.pi,size=5)\n", "phi = np.random.uniform(0,2*np.pi,size=5)\n", "hp.projscatter(theta,phi,c='b',marker='x',label='GFU Event')\n", "```" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Plotting code goes here\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Defining some functions\n", "To plot the angular errors of the neutrino events and the 90% contour of the GW event, we're going to define some convenience functions. We only use these functions once in this example but in general if you are making a lot of plots that contain similar elements it is good to encapsulate repeated tasks into functions.\n", "\n", "Ecapsulating small well-defined tasks into functions is also good programming practice in general.\n", "\n", "#### First we will define a function for computing the angular error contour of the neutrino events.\n", "\n", "Read through the code carefully to understand what the function is doing, in particular the healpy function calls" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def compute_ang_err(ra, dec, sigma):\n", " r''' Compute circular uncertainty contour around point\n", " on a healpix map.\n", " Parameters:\n", " -----------\n", " ra: float\n", " Right ascension of center point in radians\n", " dec: float \n", " declination of center point in radians\n", " sigma:\n", " Angular uncertainty (radius of circle to draw \n", " around (ra,dec)) in radians\n", " Returns:\n", " --------\n", " Theta: array\n", " theta values of contour in radians\n", " Phi: array\n", " phi values of contour in radians\n", " '''\n", "\n", " dec = np.pi/2 - dec\n", " \n", " delta = sigma\n", " step = 1. / np.sin(delta) / 20.\n", " bins = int(360./step)\n", " Theta = np.zeros(bins+1, dtype=np.double)\n", " Phi = np.zeros(bins+1, dtype=np.double)\n", " \n", " # define the contour\n", " for j in range(0,bins) :\n", " phi = j*step/180.*np.pi\n", " vx = np.cos(phi)*np.sin(ra)*np.sin(delta) + np.cos(ra)*(np.cos(delta)*np.sin(dec) + np.cos(dec)*np.sin(delta)*np.sin(phi))\n", " vy = np.cos(delta)*np.sin(dec)*np.sin(ra) + np.sin(delta)*(-np.cos(ra)*np.cos(phi) + np.cos(dec)*np.sin(ra)*np.sin(phi))\n", " vz = np.cos(dec)*np.cos(delta) - np.sin(dec)*np.sin(delta)*np.sin(phi)\n", " \n", " idx = hp.vec2pix(nside, vx, vy, vz)\n", " DEC, RA = hp.pix2ang(nside, idx)\n", " \n", " Theta[j] = DEC\n", " Phi[j] = RA\n", " Theta[bins] = Theta[0]\n", " Phi[bins] = Phi[0]\n", "\n", " return Theta, Phi" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Next we define a function that computes the 90% containment of the GW event\n", "Read through the code carefully and note the calls to `meander` and `healpy`" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def compute_contours(proportions, samples):\n", " r''' Compute containment contour around desired level.\n", " E.g 90% containment of a PDF on a healpix map\n", "\n", " Parameters:\n", " -----------\n", " proportions: list\n", " list of containment level to make contours for.\n", " E.g [0.68,0.9]\n", " samples: array\n", " array of values read in from healpix map\n", " E.g samples = hp.read_map(file)\n", " Returns:\n", " --------\n", " theta_list: list\n", " List of arrays containing theta values for desired contours\n", " phi_list: list\n", " List of arrays containing phi values for desired contours\n", " '''\n", "\n", " levels = []\n", " \n", " # Get the angular coordinates of each pixel\n", " nside = hp.pixelfunc.get_nside(samples)\n", " sample_points = np.array(hp.pix2ang(nside,np.arange(len(samples)))).T\n", " \n", " # Sort the pixels so we can compute the level\n", " sorted_samples = np.flip(np.sort(samples))\n", " for proportion in proportions:\n", " # Find the pixel that brings us over the threshold\n", " level_index = (np.cumsum(sorted_samples) > proportion).tolist().index(True)\n", " \n", " # Use the midpoint of the two pixels as the level\n", " level = (sorted_samples[level_index] + (sorted_samples[level_index-1] if level_index-1 > 0 else proportion*2)) / 2.0\n", " levels.append(level)\n", " \n", " # Compute the contours\n", " contours_by_level = meander.spherical_contours(sample_points, samples, levels)\n", " \n", " # Transform the coordinates to something healpy will take\n", " theta_list = []\n", " phi_list=[]\n", " for contours in contours_by_level:\n", " for contour in contours:\n", " theta, phi = contour.T\n", " phi[phi<0] += 2.0*np.pi\n", " theta_list.append(theta)\n", " phi_list.append(phi)\n", "\n", " return theta_list, phi_list" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Adding some errors to the neutrino events\n", "```\n", "# Plot GW skymap in Mollweide projection\n", "hp.mollview(probs,cbar=True,unit=r'Probability',min=0,max=3e-5,rot=180,cmap=cmap)\n", "hp.graticule() # Set grid lines\n", "\n", "# Let's generate 5 random neutrinos on the sky and plot them\n", "theta= np.random.uniform(0,np.pi,size=5)\n", "phi = np.random.uniform(0,2*np.pi,size=5)\n", "hp.projscatter(theta,phi,c='b',marker='x',label='GFU Event')\n", "\n", "# Now let's assign some angular error to each event. Generally these\n", "# neutrinos have ~1 deg angular error (39% containment) but we will make it \n", "# larger here for visual effect\n", "sigma = np.random.uniform(np.radians(2),np.radians(5),size=5)\n", "\n", "#plot events on sky with error contours\n", "for i in range(theta.size):\n", " my_contour = compute_ang_err(phi[i],np.pi/2-theta[i],sigma[i])\n", " hp.projplot(my_contour[0], my_contour[1], linewidth=1.75, color=\"gray\",\n", " linestyle='solid',coord='C')```" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Plotting code goes here\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Adding the 90% containment contour\n", "```\n", "# Plot GW skymap in Mollweide projection\n", "hp.mollview(probs,cbar=True,unit=r'Probability',min=0,max=3e-5,rot=180,cmap=cmap)\n", "hp.graticule() # Set grid lines\n", "\n", "# Let's generate 5 random neutrinos on the sky and plot them\n", "theta= np.random.uniform(0,np.pi,size=5)\n", "phi = np.random.uniform(0,2*np.pi,size=5)\n", "hp.projscatter(theta,phi,c='b',marker='x',label='GFU Event')\n", "\n", "# Now let's assign some angular error to each event. Generally these\n", "# neutrinos have ~1 deg angular error (39% containment) but we will make it \n", "# larger here for visual effect\n", "sigma = np.random.uniform(np.radians(2),np.radians(5),size=5)\n", "\n", "#plot events on sky with error contours\n", "for i in range(theta.size):\n", " my_contour = compute_ang_err(phi[i],np.pi/2-theta[i],sigma[i])\n", " hp.projplot(my_contour[0], my_contour[1], linewidth=1.75, color=\"gray\",\n", " linestyle='solid',coord='C')\n", "\n", "# Draw containment contour around GW skymap\n", "reduced_probs = hp.pixelfunc.ud_grade(probs, 64) #reduce nside to make it faster\n", "reduced_probs = reduced_probs/np.sum(reduced_probs)\n", "pixels = np.arange(reduced_probs.size)\n", "sample_points = np.array(hp.pix2ang(nside,pixels)).T\n", "\n", "### plot 90% containment contour of GW PDF\n", "levels = [0.90]\n", "theta_contour, phi_contour = compute_contours(levels,reduced_probs)\n", "hp.projplot(theta_contour[0],phi_contour[0],linewidth=1,c='k',label='GW (90% C.L)')\n", "for i in range(1,len(theta_contour)):\n", " hp.projplot(theta_contour[i],phi_contour[i],linewidth=1,c='k')```" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Plotting code goes here\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Finally let's add some labels and a legend\n", "```\n", "# Plot GW skymap in Mollweide projection\n", "hp.mollview(probs,cbar=True,unit=r'Probability',min=0,max=3e-5,rot=180,cmap=cmap)\n", "hp.graticule() # Set grid lines\n", "\n", "# Let's generate 5 random neutrinos on the sky and plot them\n", "theta= np.random.uniform(0,np.pi,size=5)\n", "phi = np.random.uniform(0,2*np.pi,size=5)\n", "hp.projscatter(theta,phi,c='b',marker='x',label='GFU Event')\n", "\n", "# Now let's assign some angular error to each event. Generally these\n", "# neutrinos have ~1 deg angular error (39% containment) but we will make it \n", "# larger here for visual effect\n", "sigma = np.random.uniform(np.radians(2),np.radians(5),size=5)\n", "\n", "#plot events on sky with error contours\n", "for i in range(theta.size):\n", " my_contour = compute_ang_err(phi[i],np.pi/2-theta[i],sigma[i])\n", " hp.projplot(my_contour[0], my_contour[1], linewidth=1.75, color=\"gray\",\n", " linestyle='solid',coord='C')\n", "\n", "# Draw containment contour around GW skymap\n", "reduced_probs = hp.pixelfunc.ud_grade(probs, 64) #reduce nside to make it faster\n", "reduced_probs = reduced_probs/np.sum(reduced_probs)\n", "pixels = np.arange(reduced_probs.size)\n", "sample_points = np.array(hp.pix2ang(nside,pixels)).T\n", "\n", "### plot 90% containment contour of GW PDF\n", "levels = [0.90]\n", "theta_contour, phi_contour = compute_contours(levels,reduced_probs)\n", "hp.projplot(theta_contour[0],phi_contour[0],linewidth=1,c='k',label='GW (90% C.L)')\n", "for i in range(1,len(theta_contour)):\n", " hp.projplot(theta_contour[i],phi_contour[i],linewidth=1,c='k')\n", "\n", "#Add labels to plot\n", "plt.text(2.0,0., r\"$0^\\circ$\", ha=\"left\", va=\"center\")\n", "plt.text(1.9,0.45, r\"$30^\\circ$\", ha=\"left\", va=\"center\")\n", "plt.text(1.4,0.8, r\"$60^\\circ$\", ha=\"left\", va=\"center\")\n", "plt.text(1.9,-0.45, r\"$-30^\\circ$\", ha=\"left\", va=\"center\")\n", "plt.text(1.4,-0.8, r\"$-60^\\circ$\", ha=\"left\", va=\"center\")\n", "plt.text(1.333, -0.15, r\"$60^\\circ$\", ha=\"center\", va=\"center\")\n", "plt.text(.666, -0.15, r\"$120^\\circ$\", ha=\"center\", va=\"center\")\n", "plt.text(0.0, -0.15, r\"$180^\\circ$\", ha=\"center\", va=\"center\")\n", "plt.text(-.666, -0.15, r\"$240^\\circ$\", ha=\"center\", va=\"center\")\n", "plt.text(-1.333, -0.15, r\"$300^\\circ$\", ha=\"center\", va=\"center\")\n", "plt.text(-2.0, -0.15, r\"$360^\\circ$\", ha=\"center\", va=\"center\")\n", "\n", "plt.title('GW Skymap',fontsize=15)\n", "plt.legend(loc=1,bbox_to_anchor=(1.08,1.09),prop={'size': 14})```" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Plotting code goes here\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If you've gotten this far hopefully you've had the opportunity to go through the above examples carefully.\n", "\n", "There are some more examples in the advanced plotting folder or here https://icecube.wisc.edu/~icecube-bootcamp/bootcamp2019/advanced_plotting/ which you can go through in this notebook" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If you're comfortable with the examples in that directory, you can take a crack at a final exercise.\n", "\n", "See if you can recreate the axes used in figure 5 of this paper: https://arxiv.org/pdf/1507.05287.pdf\n", "\n", "You will need to create a new scale and register it with matplotlib: https://matplotlib.org/1.4.0/devel/add_new_projection.html" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.6.7" } }, "nbformat": 4, "nbformat_minor": 2 }