Skip to content

Instantly share code, notes, and snippets.

@thespacedoctor thespacedoctor/README.md Secret
Last active Dec 8, 2017

Embed
What would you like to do?
[IceCube FOV Plot] Script to plot location of IceCube 160427A #icecube #plot #160427A

Here is the plot produced by this script:

Icecube FOV 160427A

#!/usr/local/bin/python
# encoding: utf-8
"""
icecube-160427A-plot.py
=======================
:Summary:
Script to plot location of IceCube 160427A
:Author:
David Young
:Date Created:
August 23, 2016
Usage:
icecube-160427A-plot
-h, --help show this help message
-v, --version show version
-s, --settings the settings file
..todo::
@review: when complete pull all general functions and classes into dryxPython
"""
## USER DEFINED VARIABLES - ONLY PLACE THIS SCRIPT NEEDS MODIFIED ##
################# GLOBAL IMPORTS ####################
import sys
import os
os.environ['TERM'] = 'vt100'
import readline
import glob
import docopt
# SUPPRESS MATPLOTLIB WARNINGS
import warnings
warnings.filterwarnings("ignore")
import matplotlib as mpl
from fundamentals import tools, times
from astropy import wcs as awcs
import numpy as np
from matplotlib.pyplot import savefig
from crowdedText import adjust_text
def main(arguments=None):
"""
The main function used when ``icecube-160427A-plot.py`` is run as a single script from the cl, or when installed as a cl command
"""
# setup the command-line util settings
su = tools(
arguments=arguments,
docString=__doc__,
logLevel="DEBUG",
options_first=False,
projectName="lv",
tunnel=False
)
arguments, settings, log, dbConn = su.setup()
plot(
log=log
)
return
def plot(
log,
fileFormats=["pdf"],
folderName="",
projection="tan",
outputDirectory=False):
"""
"""
log.info('starting the ``plot`` method')
import healpy as hp
from matplotlib.font_manager import FontProperties
font = FontProperties()
import matplotlib.pyplot as plt
# HEALPY REQUIRES RA, DEC IN RADIANS AND AS TWO SEPERATE ARRAYS
import math
pi = (4 * math.atan(1.0))
DEG_TO_RAD_FACTOR = pi / 180.0
RAD_TO_DEG_FACTOR = 180.0 / pi
# VARIABLES
font.set_family("Arial")
pixelSizeDeg = 0.066667
unit = "likelihood"
cmap = "YlOrRd"
colorBar = False
# PARAMETERS
raCentre = 240.57
decCentre = 9.34
centralCoordinate = [raCentre, decCentre]
plotParameters = {}
plotParameters["raRange"] = 1.5
plotParameters["decRange"] = 1.5
plotTitle = "IceCube-160427A-QSOs-PS16cgx-Location"
# INITIALISE FIGURE
fig = plt.figure()
# pixelSizeDeg = 1.0 / 60.
# CREATE A NEW WCS OBJECT
w = awcs.WCS(naxis=2)
# SET THE REQUIRED PIXEL SIZE
w.wcs.cdelt = np.array([pixelSizeDeg, pixelSizeDeg])
# WORLD COORDINATES AT REFERENCE PIXEL
w.wcs.crval = centralCoordinate
projectionDict = {
"mollweide": "MOL",
"aitoff": "AIT",
"hammer": "AIT",
"lambert": "ZEA",
"polar": "TAN",
"rectilinear": "MER"
}
if projection in ["mollweide", "aitoff", "hammer", "lambert", "polar", "rectilinear"]:
# MAP VISULISATION RATIO IS ALWAYS 1/2
xRange = 2000
yRange = xRange / 2.
# SET THE REFERENCE PIXEL TO THE CENTRE PIXEL
w.wcs.crpix = [xRange / 2., yRange / 2.]
# FULL-SKY MAP SO PLOT FULL RA AND DEC RANGES
# DEC FROM 180 to 0
theta = np.linspace(np.pi, 0, yRange)
latitude = np.radians(np.linspace(-90, 90, yRange))
# RA FROM -180 to +180
phi = np.linspace(-np.pi, np.pi, xRange)
longitude = np.radians(np.linspace(-180, 180, xRange))
X, Y = np.meshgrid(longitude, latitude)
# PROJECT THE MAP TO A RECTANGULAR MATRIX xRange X yRange
PHI, THETA = np.meshgrid(phi, theta)
healpixIds = hp.ang2pix(nside, THETA, PHI)
probs = aMap[healpixIds]
# healpixIds = np.reshape(healpixIds, (1, -1))[0]
# CTYPE FOR THE FITS HEADER
thisctype = projectionDict[projection]
w.wcs.ctype = ["RA---%(thisctype)s" %
locals(), "DEC--%(thisctype)s" % locals()]
stampProb = np.sum(aMap)
print "Probability for the plot stamp is %(stampProb)s" % locals()
# MATPLOTLIB IS DOING THE PROJECTION
# ax = fig.add_subplot(111, projection=projection)
# USE WCS AS THE PROJECTION
if projection == "mollweide":
ax = plt.axes(projection=ccrs.Mollweide())
ax.set_extent([80, 170, -45, 30])
else:
ax = fig.add_axes([0.15, 0.1, 0.8, 0.8], projection=projection)
print ax.get_extent()
# RASTERIZED MAKES THE MAP BITMAP WHILE THE LABELS REMAIN VECTORIAL
# FLIP LONGITUDE TO THE ASTRO CONVENTION
ax.contourf(longitude[
::-1], latitude, probs,
transform=ccrs.PlateCarree(),
cmap=cmap)
# ax.coastlines()
ax.set_global()
plt.show()
sys.exit(0)
image = ax.pcolormesh(longitude[
::-1], latitude, probs, rasterized=True, cmap=cmap)
plt.show()
# GRATICULE
if projection not in ["polar", "rectilinear", "mollweide"]:
ax.set_longitude_grid(60)
ax.set_latitude_grid(45)
ax.xaxis.set_major_formatter(ThetaFormatterShiftPi(60))
ax.set_longitude_grid_ends(90)
# CONTOURS - NEED TO ADD THE CUMMULATIVE PROBABILITY
i = np.flipud(np.argsort(aMap))
cumsum = np.cumsum(aMap[i])
cls = np.empty_like(aMap)
cls[i] = cumsum * 99.99999999 * stampProb
# EXTRACT CONTOUR VALUES AT HEALPIX INDICES
contours = []
contours[:] = [cls[i] for i in healpixIds]
# contours = np.reshape(np.array(contours), (yRange, xRange))
CS = ax.contour(longitude[::-1], latitude,
contours, linewidths=.5, alpha=0.7, zorder=2, origin='image')
if projection not in ["lambert", "aitoff"]:
CS.set_alpha(0.5)
CS.clabel(fontsize=10, inline=True,
fmt='%2.1f', fontproperties=font, alpha=0.0)
# COLORBAR
if colorBar:
cb = fig.colorbar(image, orientation='horizontal',
shrink=.6, pad=0.05, ticks=[0, 1])
cb.ax.xaxis.set_label_text("likelihood")
cb.ax.xaxis.labelpad = -8
# WORKAROUND FOR ISSUE WITH VIEWERS, SEE COLORBAR DOCSTRING
cb.solids.set_edgecolor("face")
ax.tick_params(axis='x', labelsize=16)
ax.tick_params(axis='y', labelsize=16)
# lon.set_ticks_position('bt')
# lon.set_ticklabel_position('b')
# lon.set_ticklabel(size=20)
# lat.set_ticklabel(size=20)
# lon.set_axislabel_position('b')
# lat.set_ticks_position('lr')
# lat.set_ticklabel_position('l')
# lat.set_axislabel_position('l')
# # REMOVE TICK LABELS
# ax.xaxis.set_ticklabels([])
# ax.yaxis.set_ticklabels([])
# # REMOVE GRID
# ax.xaxis.set_ticks([])
# ax.yaxis.set_ticks([])
# REMOVE WHITE SPACE AROUND FIGURE
spacing = 0.01
plt.subplots_adjust(bottom=spacing, top=1 - spacing,
left=spacing, right=1 - spacing)
plt.grid(True)
elif projection in ["tan"]:
# UNPACK THE PLOT PARAMETERS
raRange = plotParameters["raRange"]
decRange = plotParameters["decRange"]
raMax = centralCoordinate[0] + raRange / 2.
raMin = centralCoordinate[0] - raRange / 2.
decMax = centralCoordinate[1] + decRange / 2.
decMin = centralCoordinate[1] - decRange / 2.
# DETERMINE THE PIXEL GRID X,Y RANGES
xRange = int(raRange / pixelSizeDeg)
yRange = int(decRange / pixelSizeDeg)
# SET THE REFERENCE PIXEL TO THE CENTRE PIXEL
w.wcs.crpix = [xRange / 2., yRange / 2.]
# USE THE "GNOMONIC" PROJECTION ("COORDINATESYS---PROJECTION")
w.wcs.ctype = ["RA---TAN", "DEC--TAN"]
# CREATE A PIXEL GRID - 2 ARRAYS OF X, Y
columns = []
px = np.tile(np.arange(0, xRange), yRange)
py = np.repeat(np.arange(0, yRange), xRange)
# CONVERT THE PIXELS TO WORLD COORDINATES
wr, wd = w.wcs_pix2world(px, py, 1)
# MAKE SURE RA IS +VE
nr = []
nr[:] = [r if r > 0 else r + 360. for r in wr]
wr = np.array(nr)
# CREATE THE FITS HEADER WITH WCS
header = w.to_header()
# CREATE THE FITS FILE
# GRAB THE WCS FROM HEADER GENERATED EARLIER
from wcsaxes import datasets, WCS
from astropy.wcs import WCS
from wcsaxes import WCSAxes
wcs = WCS(header)
# USE WCS AS THE PROJECTION
ax = WCSAxes(fig, [0.15, 0.1, 0.8, 0.8], wcs=wcs)
# note that the axes have to be explicitly added to the figure
ax = fig.add_axes(ax)
# RESET THE AXES TO THE FRAME OF THE FITS FILE
ax.set_xlim(-0.5, xRange - 0.5)
ax.set_ylim(-0.5, yRange - 0.5)
# THE COORDINATES USED IN THE PLOT CAN BE ACCESSED USING THE COORDS
# ATTRIBUTE (NOT X AND Y)
lon = ax.coords[0]
lat = ax.coords[1]
lon.set_axislabel('RA (deg)', minpad=0.87,
size=20)
lat.set_axislabel('DEC (deg)', minpad=0.87,
size=20)
lon.set_major_formatter('d.d')
lat.set_major_formatter('d.d')
# THE SEPARATORS FOR ANGULAR COORDINATE TICK LABELS CAN ALSO BE SET BY
# SPECIFYING A STRING
lat.set_separator(':-s')
# SET THE APPROXIMATE NUMBER OF TICKS, WITH COLOR & PREVENT OVERLAPPING
# TICK LABELS FROM BEING DISPLAYED.
lon.set_ticks(number=4, color='#657b83',
exclude_overlapping=True, size=10)
lat.set_ticks(number=10, color='#657b83',
exclude_overlapping=True, size=10)
# MINOR TICKS NOT SHOWN BY DEFAULT
lon.display_minor_ticks(True)
lat.display_minor_ticks(True)
lat.set_minor_frequency(2)
# CUSTOMISE TICK POSITIONS (l, b, r, t == left, bottom, right, or
# top)
lon.set_ticks_position('bt')
lon.set_ticklabel_position('b')
lon.set_ticklabel(size=20)
lat.set_ticklabel(size=20)
lon.set_axislabel_position('b')
lat.set_ticks_position('lr')
lat.set_ticklabel_position('l')
lat.set_axislabel_position('l')
# HIDE AXES
# lon.set_ticklabel_position('')
# lat.set_ticklabel_position('')
# lon.set_axislabel('', minpad=0.5, fontsize=12)
# lat.set_axislabel('', minpad=0.5, fontsize=12)
# ADD A GRID
ax.coords.grid(color='#657b83', alpha=0.5, linestyle='dashed')
plt.gca().invert_xaxis()
else:
log.error(
'please give a valid projection. The projection given was `%(projection)s`.' % locals())
header = w.to_header()
# CREATE THE FITS FILE
subTitle = ""
raDeg = raCentre
decDeg = decCentre
# PLOT ICECUBE FOOTPRINT
from matplotlib.patches import Circle
from matplotlib.patches import Ellipse
height = 36 * 2 / 60.
width = height / math.cos(decDeg * DEG_TO_RAD_FACTOR)
# MULTIPLE CIRCLES
if projection in ["tan"]:
circ = Ellipse(
(raDeg, decDeg), width=width, height=height, alpha=0.2, color='#268bd2', fill=True, transform=ax.get_transform('fk5'), zorder=3)
else:
if raDeg > 180.:
raDeg = raDeg - 360.
circ = Ellipse(
(-raDeg * DEG_TO_RAD_FACTOR, decDeg * DEG_TO_RAD_FACTOR), width=width * DEG_TO_RAD_FACTOR, height=height * DEG_TO_RAD_FACTOR, alpha=0.2, color='#268bd2', fill=True, zorder=3)
ax.add_patch(circ)
# # LEGEND FOR PS1
# if len(ps1Pointings):
# raDeg = 90.
# decDeg = 10.
# height = 3.5
# width = height / math.cos(decDeg * DEG_TO_RAD_FACTOR)
# if projection in ["tan"]:
# circ = Ellipse(
# (raDeg, decDeg), width=width, height=height, alpha=0.2, color='#859900', fill=True, transform=ax.get_transform('fk5'), zorder=3)
# ax.text(
# raDeg - 3,
# decDeg - 1.,
# "PS1",
# fontsize=14,
# zorder=4,
# family='monospace',
# transform=ax.get_transform('fk5')
# )
# else:
# height = 6.
# width = height / math.cos(decDeg * DEG_TO_RAD_FACTOR)
# raDeg = 285.
# decDeg = 38.
# if raDeg > 180.:
# raDeg = raDeg - 360.
# circ = Ellipse(
# (-raDeg * DEG_TO_RAD_FACTOR, decDeg * DEG_TO_RAD_FACTOR), width=width * DEG_TO_RAD_FACTOR, height=height * DEG_TO_RAD_FACTOR, alpha=0.2, color='#859900', fill=True, zorder=3)
# ax.text(
# (-raDeg + 8.) * DEG_TO_RAD_FACTOR,
# (decDeg - 3.) * DEG_TO_RAD_FACTOR,
# "PS1",
# fontsize=14,
# zorder=4,
# family='monospace'
# )
# ax.add_patch(circ)
# ADD DATA POINTS FOR TRANSIENTS
names = []
ra = []
dec = []
raRad = []
decRad = []
texts = []
ps1QSOs = []
ps1SN = [
{"ps1_designation": "PS16cgv",
"ra_psf": "15 59 53.84", "dec_psf": "+09 11 08.4"},
{"ps1_designation": "PS16cgw",
"ra_psf": "16 00 39.69", "dec_psf": "+09 39 21.2"},
{"ps1_designation": "PS16cfu",
"ra_psf": "16 01 15.66", "dec_psf": "+09 25 04.7"},
{"ps1_designation": "PS16cfz",
"ra_psf": "16 02 11.96", "dec_psf": "+09 54 07.9"},
{"ps1_designation": "PS16cgb",
"ra_psf": "16 02 19.12", "dec_psf": "+09 34 50.1"},
{"ps1_designation": "PS16cgi", "ra_psf": "16 03 27.72", "dec_psf": "+08 54 05.2"}
]
ps1BigSN = [{"ps1_designation": "PS16cgx",
"ra_psf": "16 01 18.60", "dec_psf": "+09 51 53.1"}]
# ps1QSOs = [
# {"ps1_designation": ""}
# ]
from astrocalc.coords import unit_conversion
converter = unit_conversion(
log=log
)
for trans in ps1SN:
# if trans["ps1_designation"] in ["PS15dpg", "PS15dpp", "PS15dpq", "PS15don", "PS15dpa", "PS15dom"]:
# continue
name = trans["ps1_designation"]
names.append(name)
decDeg = converter.dec_sexegesimal_to_decimal(
dec=trans["dec_psf"]
)
raDeg = converter.ra_sexegesimal_to_decimal(
ra=trans["ra_psf"]
)
ra.append(raDeg)
dec.append(decDeg)
raRad.append(-raDeg * DEG_TO_RAD_FACTOR)
decRad.append(decDeg * DEG_TO_RAD_FACTOR)
if len(ra) > 0:
# MULTIPLE CIRCLES
if projection in ["tan"]:
ax.scatter(
x=np.array(ra),
y=np.array(dec),
transform=ax.get_transform('fk5'),
s=15,
c='black',
edgecolor='black',
alpha=1,
zorder=4
)
xx, yy = w.wcs_world2pix(np.array(ra), np.array(dec), 0)
print xx, yy
# ADD TRANSIENT LABELS
for r, d, n in zip(xx, yy, names):
# if n == "PS16cgv":
# r = r + 0.5
texts.append(ax.text(
r,
d,
n,
fontsize=16,
zorder=4,
family='monospace'
))
else:
ax.scatter(
x=np.array(raRad),
y=np.array(decRad),
s=15,
c='#dc322f',
edgecolor='#dc322f',
alpha=1,
zorder=4
)
ps1QSOs = [
{"ps1_designation": "PS16cfw",
"ra_psf": "16:01:39.10", "dec_psf": "+09:21:59.5"},
{"ps1_designation": "PS16cfx",
"ra_psf": "16:01:49.65", "dec_psf": "+08:50:29.7"},
{"ps1_designation": "PS16cfy",
"ra_psf": "16:02:05.20", "dec_psf": "+09:46:30.4"},
{"ps1_designation": "PS16cgg",
"ra_psf": "16:03:09.19", "dec_psf": "+09:20:19.5"},
{"ps1_designation": "PS16cgh",
"ra_psf": "16:03:17.92", "dec_psf": "+09:00:38.0"},
{"ps1_designation": "PS16cgm",
"ra_psf": "15:59:47.90", "dec_psf": "+09:10:22.4"},
{"ps1_designation": "PS16cgn",
"ra_psf": "15:59:54.76", "dec_psf": "+09:14:26.3"},
{"ps1_designation": "PS16cgo",
"ra_psf": "16:01:06.27", "dec_psf": "+08:46:05.8"}
]
# ADD DATA POINTS FOR TRANSIENTS
names = []
ra = []
dec = []
raRad = []
decRad = []
for trans in ps1QSOs:
# if trans["ps1_designation"] in ["PS15dpg", "PS15dpp", "PS15dpq", "PS15don", "PS15dpa", "PS15dom"]:
# continue
name = trans["ps1_designation"]
names.append(name)
decDeg = converter.dec_sexegesimal_to_decimal(
dec=trans["dec_psf"]
)
raDeg = converter.ra_sexegesimal_to_decimal(
ra=trans["ra_psf"]
)
ra.append(raDeg)
dec.append(decDeg)
raRad.append(-raDeg * DEG_TO_RAD_FACTOR)
decRad.append(decDeg * DEG_TO_RAD_FACTOR)
if len(ra) > 0:
# MULTIPLE CIRCLES
if projection in ["tan"]:
ax.scatter(
x=np.array(ra),
y=np.array(dec),
transform=ax.get_transform('fk5'),
s=15,
c='#268bd2',
edgecolor='#268bd2',
alpha=1,
zorder=4,
marker="D"
)
xx, yy = w.wcs_world2pix(np.array(ra), np.array(dec), 0)
print xx, yy
# ADD TRANSIENT LABELS
for r, d, n in zip(xx, yy, names):
# if n in ["PS16cgn", "PS16cgm"]:
# r = r + 0.8
texts.append(ax.text(
r,
d,
n,
color='#268bd2',
fontsize=16,
zorder=4,
family='monospace'
))
else:
ax.scatter(
x=np.array(raRad),
y=np.array(decRad),
s=15,
c='#268bd2',
edgecolor='#268bd2',
alpha=1,
zorder=4
)
# ADD DATA POINTS FOR TRANSIENTS
names = []
ra = []
dec = []
raRad = []
decRad = []
for trans in ps1BigSN:
# if trans["ps1_designation"] in ["PS15dpg", "PS15dpp", "PS15dpq", "PS15don", "PS15dpa", "PS15dom"]:
# continue
name = trans["ps1_designation"]
names.append(name)
decDeg = converter.dec_sexegesimal_to_decimal(
dec=trans["dec_psf"]
)
raDeg = converter.ra_sexegesimal_to_decimal(
ra=trans["ra_psf"]
)
ra.append(raDeg)
dec.append(decDeg)
raRad.append(-raDeg * DEG_TO_RAD_FACTOR)
decRad.append(decDeg * DEG_TO_RAD_FACTOR)
if len(ra) > 0:
# MULTIPLE CIRCLES
if projection in ["tan"]:
ax.scatter(
x=np.array(ra),
y=np.array(dec),
transform=ax.get_transform('fk5'),
s=50,
facecolor='none',
edgecolor='#dc322f',
alpha=1,
zorder=4
)
ax.scatter(
x=np.array(ra),
y=np.array(dec),
transform=ax.get_transform('fk5'),
s=15,
c='#dc322f',
edgecolor='#dc322f',
alpha=1,
zorder=4
)
xx, yy = w.wcs_world2pix(np.array(ra), np.array(dec), 0)
# ADD TRANSIENT LABELS
for r, d, n in zip(xx, yy, names):
# if n == "PS16cgv":
# r = r + 0.5
texts.append(ax.text(
r,
d,
n,
color='#dc322f',
fontsize=16,
zorder=4,
family='monospace'
))
else:
ax.scatter(
x=np.array(raRad),
y=np.array(decRad),
s=15,
c='#dc322f',
edgecolor='#dc322f',
alpha=1,
zorder=4
)
from random import shuffle
shuffle(texts)
adjust_text(
xx,
yy,
texts,
expand_text=(1.0, 1.0),
expand_points=(1.0, 1.0),
va='center',
ha='right',
autoalign=False,
force_text=1.,
force_points=1.,
lim=100,
precision=0,
only_move={},
text_from_text=True,
text_from_points=True,
save_steps=False,
save_prefix='',
save_format='png',
add_step_numbers=True,
min_arrow_sep=3.,
draggable=True,
arrowprops=dict(arrowstyle="->", facecolor='black', lw=1.2,
patchB=None, shrinkB=0, connectionstyle="arc3,rad=0.1", zorder=3, alpha=0.5),
fontsize=16,
family='monospace',
repel_from_axes=True
)
# TIME-RANGE LABEL
fig = plt.gcf()
fig.set_size_inches(8.0, 8.0)
fWidth, fHeight = fig.get_size_inches()
print fWidth, fHeight
if projection == "tan":
plt.text(
xRange * 0.25,
# xRange * 0.95,
yRange * 0.93,
"",
fontsize=16,
zorder=4,
color="#dc322f",
fontproperties=font
)
plt.text(
xRange * 0.1,
# xRange * 0.95,
yRange * 0.93,
"",
fontsize=20,
zorder=4,
color="black",
fontproperties=font
)
else:
plt.text(
fWidth * 0.7,
# xRange * 0.95,
fHeight * 0.95,
timeRangeLabel,
fontsize=16,
zorder=4,
color="#dc322f",
fontproperties=font
)
plt.text(
fWidth * 0.87,
# xRange * 0.95,
fHeight * 0.315,
"",
fontsize=20,
zorder=4,
color="black",
fontproperties=font
)
# Recursively create missing directories
plotDir = "./"
plotTitle = plotTitle.replace(" ", "_").replace(
"<", "lt").replace(">", "gt").replace(",", "").replace("\n", "_").replace("&", "and")
figureName = """%(plotTitle)s""" % locals(
)
for f in fileFormats:
figurePath = "%(plotDir)s/%(figureName)s_%(projection)s.%(f)s" % locals()
savefig(figurePath, bbox_inches='tight', dpi=300)
log.info('completed the ``plot`` method')
return None
if __name__ == '__main__':
main()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.