python-script
python-script
— 94.9 KB
File contents
"""
Version: 1.1
Author
John Carpenter (jcarpent @ alma.cl)
Change history:
March 18, 2016 : First releaset of script (JMC)
Purpose
plotobs_cycle4.py prints and plots the pointing positions of previous
ALMA observations (Cycle 1, 2, or 3, but not Cycle 0) and active
Grade A projects that are scheduled to be observed.
The script provides three basic function:
1) Plot observations around a specified source name or coordinate,
and overlay a proposed observation (single pointing or
rectangular mosaic) on the plot.
2) List the sources observed for a project
3) List the details (sensitivity, angular resolution, etc...)
of an observations
Support
plotobs_cycle4.py is user-contributed software that is distributed
"as-is". plotobs_cycle4.py is not an offical ALMA product and is not
supported by the ALMA Regional Centers (ARCs). Please report any
bugs to John Carpenter.
Dependencies
1) The plotobs_cycle4.py script requires a data file containing
the list of ALMA observations. The root name of the file is
contained in the variable LIST_OF_OBSERVATIONS, and the data file
can be found on the ALMA Science Portal.
The science portal has two versions of the data: an excel spreadsheet
(.xls extension) and a comma-separated-variable (.csv extension)
format. Both formats work with this script. To read the excel
file, the python package "xlrd" must be installed. My experience
is that this package is often not installed, and using the .csv
file is more likely to be successful.
2) Various python packages must be installed, which are listed below
under "Required software".
To run the program
# Start python with graphics enabled by the --pylab option
ipython --pylab
# Import program
import plotobs_cycle4 as p4
# Main routines
p4.plot(...) # Plots ALMA observations and a proposed observation
p4.project(...) # Show observations for ALMA projects
p4.row(...) # List observing parameters for rows in spreadsheet
Examples
# Plot observations within 480" x 480" of M83 along with a proposed
# pointing at the position of the galaxy
p4.plot('M83', plotsize=480)
# Plot observations within 480" x 480" of M83, along with a
# proposed mosaic at 345 GHz
p4.plot('M83', plotsize=480, length=240, width=120,
pa=60, freq=345., mosaic=True)
# Plot observations within 480" of (ra, dec) = (17h20m42s, -35d48m04s)
p4.plot('17h20m42s', '-35d48m04s', plotsize=280, frame='icrs', unit='deg')
# Plot 1 deg region around ra, dec = (150, 2.2) J2000
p4.plot(150, 2.2, plotsize=3600)
# Plot 1 deg around galactic center, but mosaics only
p4.plot(0, 0, frame='galactic', plotsize=3600., mosonly=True)
# Plot observations within a 60" x 60" of HL Tau and DG Tau
p4.plot('HL Tau, DG Tau', plotsize=60)
# Plot observations near coordinates in row 272 of the spreadsheet
p4.plot(272)
# Plot observations only for row 272 in the spreadsheet,
# but plot only row 272
p4.plot(272, include=272)
# Plot observations for sources listed in a text file
# Contents of input file "input.txt":
HL Tau
DG Tau
plotsize = 480
COSMOS-16199
frame = galactic
plotsize = 60
178.8590 -19.9724
p4.plot(inputFile='input.txt')
# Print observations of project 2015.1.01602.S
p4.project('2015.1.01602.S')
# Print the observational parameters for row 5719 in Excel spreadsheet
p4.row(5719)
Notes
1) After running "import plotobs as p4", additional help can
be obtained by typing in python:
help(p4)
help(p4.plot)
help(p4.row)
help(p4.project)
2) When plotobs_cycle4 is run the first time in a python session, the
spreadsheet containing the observations is read into memory and
stored in the global variable OBSERVATIONS. The data stored in
memory is then used on subsequent calls.
3) Target of Opportunity or solar system observations will likely (but
not always) have a RA/Dec coordinate of (0 deg, 0 deg). Such
observations can be identified by running the script as:
p4.plot(0, 0)
Cautions
1) The spreadsheet used by plotobs_cycle4.py contains the sensitivity
and angular resolution requested by the principal investigator.
The achieved sensitivity and angular resolution of the actual
observations will differ.
2) The sensitivity per spectral window is computed using the reference
frequency and reference bandwidth, but does not take into account
the variation in the system temperature amongst the various
spectral windows.
Software
Platforms
1) plotobs_cycle4.py program has been successfully tested on
a MacBook Pro and an iMac running OS X El Capitan, and
on linux clusters running Redhat 6.6.
The program has successfully worked with the following:
1) python 2.7
2) numpy 1.10.2
3) pandas version 0.17.1; see http://pandas.pydata.org .
4) astropy version 1.0.5; see http://www.astropy.org
5) If you want to read the excel (.xls) instead of the csv version
of the spreadsheet, the xlrd python package needs to be installed.
Internet connection
plotobs_cycle4.py does not require a internet connection to run.
plotobs_cycle4.py does use the Sesame database to resolve the
coordinates for source names. If an internet connection is
not available, the script cannot use Sesame, and will to determine
the coordinates by looking for the source name in the spreadsheet.
"""
# Load packages
import matplotlib.pylab as py
import numpy as np
import pandas as pd
import re
import os
import astropy
from astropy.coordinates import SkyCoord
import astropy.constants as G
from matplotlib.patches import Circle, Polygon
# Seaborn package is not required, but it makes the plot looks nicer
try:
import seaborn as sns
except:
pass
# Name of columns in observations
ANG_RES = 'Req. Ang. Res.'
BAND = 'Band'
DEC_DEG = 'Dec'
DEC_DMS = 'Dec_DMS'
FWHM_PB = 'FWHM PB' # Created in readObservations
IS_MOSAIC = 'Mosaic'
IS_SPECTRAL = 'isSpectralScan' # Created in readObservations
IS_SPW_SKY_FREQ = 'Is Sky Freq?'
LAS = 'Req. LAS'
LAT_OFFSET = 'Lat Offset'
LON_OFFSET = 'Long Offset'
MOS_AREA = 'mosaic area' # Created in readObservations
MOS_LENGTH = 'Mos. Length'
MOS_PA = 'Mos. PA'
MOS_WIDTH = 'Mos. Width'
MAX_SIZE = 'Max size' # Created in readObservations
MOS_SPACING = 'Mos. Spacing'
MOS_COORD = 'Mos. Coord.'
POLARIZATION = 'Polarization'
PROJECT = 'Project Code'
RA_DEG = 'RA'
RA_HMS = 'RA_HMS'
REF_FREQ = 'Ref.Frequency'
REF_BW = 'Ref.Freq.Width'
REQ_SENS = 'Req.Sensitivity'
SPS_BW = 'SPS Bandwidth'
SPS_END = 'SPS End Freq.'
SPS_SPW_RES = 'SPS Spec. Res.'
SPS_START = 'SPS Start Freq.'
SPW_FREQ = {}
SPW_BW = {}
SPW_RES = {}
TARGET = 'Target Name'
USE_ACA = 'Use 7-m?'
USE_TPA = 'Use TP?'
VELOCITY = 'Velocity'
VELOCITY_UNIT = 'Velocity unit' # added by readObservations()
VELOCITY_REF = 'Vel. Frame'
VELOCITY_DOP = 'Vel. Convention'
# Set array keywords
NWINDOWS = 16
for win in range(1, NWINDOWS+1):
SPW_FREQ[win] = 'Freq SPW %d' % win
SPW_BW[win] = 'Bandwidth SPW %d' % win
SPW_RES[win] = 'Spec.Res. SPW %d' % win
# List of Cycle 1/2/3 observations
LIST_OF_OBSERVATIONS = 'duplication_cycle4_march18'
# Rows to skip in LIST_OF_OBSERVAIONS file since they contain header information
SKIPROWS = range(0,41)
SKIPROWS.append(max(SKIPROWS)+2)
# Store observations in the spreadsheet in global variable
OBSERVATIONS = None
# List of entries in the observations, their label, and unit
OBS_ITEMS = dict()
OBS_ITEMS[ANG_RES] = ['Angular resolution', 'arcsec']
OBS_ITEMS[BAND] = ['Band', None]
OBS_ITEMS[DEC_DEG] = ['Declination', 'dms']
OBS_ITEMS[FWHM_PB] = ['12m primary beam size', 'arcsec']
OBS_ITEMS[IS_MOSAIC] = ['Is mosaic?', None]
OBS_ITEMS[IS_SPECTRAL] = ['Is spectral scan?', None]
OBS_ITEMS[LAS] = ['Largest angular scale', 'arcsec']
OBS_ITEMS[LAT_OFFSET] = ['Latitude offset', 'arcsec']
OBS_ITEMS[LON_OFFSET] = ['Latitude offset', 'arcsec']
OBS_ITEMS[MAX_SIZE] = ['Maximum field of view', 'arcsec']
OBS_ITEMS[MOS_AREA] = ['Mosaic area', 'sq. arcmin']
OBS_ITEMS[MOS_COORD] = ['Coordinate system', None]
OBS_ITEMS[MOS_LENGTH] = ['Length', 'arcsec']
OBS_ITEMS[MOS_PA] = ['Position angle', 'deg']
OBS_ITEMS[MOS_SPACING] = ['Pointing spacings', 'arcsec']
OBS_ITEMS[MOS_WIDTH] = ['Width', 'arcsec']
OBS_ITEMS[POLARIZATION] = ['Polarization', None]
OBS_ITEMS[PROJECT] = ['Project code', None]
OBS_ITEMS[RA_DEG] = ['Right ascension', 'hms']
OBS_ITEMS[REF_FREQ] = ['Reference frequency', 'GHz']
OBS_ITEMS[REF_BW] = ['Reference bandwidth', 'MHz']
OBS_ITEMS[REQ_SENS] = ['Reference sensitivity', 'mJy']
OBS_ITEMS[SPS_START] = ['Spectral scan start', 'GHz']
OBS_ITEMS[SPS_END] = ['Spectral scan end', 'GHz']
OBS_ITEMS[SPS_BW] = ['Spectral scan bandwidth', 'MHz']
OBS_ITEMS[SPS_SPW_RES] = ['Spectral scan resolution', 'MHz']
OBS_ITEMS[TARGET] = ['Target name', None]
OBS_ITEMS[VELOCITY] = ['Center velocity', None]
OBS_ITEMS[VELOCITY_UNIT] = ['Unit for velocity', None]
OBS_ITEMS[VELOCITY_REF] = ['Velocity frame', None]
OBS_ITEMS[VELOCITY_DOP] = ['Velocity convention', None]
OBS_ITEMS[USE_ACA] = ['Use ACA?', None]
OBS_ITEMS[USE_TPA] = ['Use Total Power Array?', None]
# Keywords for dictionary in readObservations only
DATA = 'data'
# Keywords for dictionary in getSourceCoordinates()
RA = 'ra'
DEC = 'dec'
ORIGINAL = 'original'
PLOTSIZE = 'plotsize'
ISCOORD = 'iscoord'
NOCOORD = 'nocoord'
# Miscellaneous keywords
COORDS = 'c'
MOSAIC_TYPE_RECTANGLE = 'Rectangle'
MOSAIC_TYPE_CUSTOM = 'Custom'
MOSAIC_TYPES = [MOSAIC_TYPE_RECTANGLE, MOSAIC_TYPE_CUSTOM]
# Miscellaneous parameters
EPSILON = 1e-4
# Set band frequencies
# http://www.almaobservatory.org/en/about-alma/how-does-alma-work/technology/front-end
BAND_UNKNOWN = 'unknown'
BAND_FREQ = dict()
BAND_FREQ['ALMA_RB_01'] = [31, 50]
BAND_FREQ['ALMA_RB_02'] = [67, 90]
BAND_FREQ['ALMA_RB_03'] = [84, 116]
BAND_FREQ['ALMA_RB_04'] = [125, 163]
BAND_FREQ['ALMA_RB_05'] = [162, 211]
BAND_FREQ['ALMA_RB_06'] = [211, 275]
BAND_FREQ['ALMA_RB_07'] = [275, 373]
BAND_FREQ['ALMA_RB_08'] = [385, 500]
BAND_FREQ['ALMA_RB_09'] = [602, 720]
BAND_FREQ['ALMA_RB_10'] = [787, 950]
BAND_FREQ[BAND_UNKNOWN] = [0, 0]
def getUsableBandwidth(bw, indx=None, throwException=True):
"""
Returns the usable bandwidth in MHz.
Inputs:
bw : nominal bandwidth in MHz.
indx: Row number in the data structure. This is used for informational
purposes only if there is an error in the input bw.
throwException : If True, then throw an exception if bandwidth is invalid.
Otherwise, return None as the usable bandwidth.
Output:
usable bandwidth in MHz
Notes:
The usable bandwidths are given in Table 5.3 of the
Cycle 3 ALMA Technical Handbook.
"""
msg = None
if abs(bw - 62.5) < EPSILON:
use_bw = 58.6
elif abs(bw - 125.) < EPSILON:
use_bw = 117.2
elif abs(bw - 250.) < EPSILON:
use_bw = 234.4
elif abs(bw - 500.) < EPSILON:
use_bw = 468.8
elif abs(bw - 1000.) < EPSILON:
use_bw = 937.5
elif abs(bw-1875) < EPSILON or abs(bw-2000) < EPSILON:
use_bw = 1875.
else:
if throwException:
msg = 'Do not recognize bandwidth value: %.1f' % bw
if indx is not None:
msg = 'Do not recognize bandwidth value (%.1f) on row %d' % (bw, getExcelFromIndex(indx))
else:
use_bw = None
if msg is not None:
raise Exception, msg
else:
return use_bw
def checkEphemerisNames(data):
"""
Check the names for sources that have a name "ephemeris". They
should have been replaced with the actual source names from
the proposals.
"""
# Set codes/sources that are moving objects and names need to be modified.
tuples = [ ]
# Separate tuples
found_tuples = [False] * len(tuples)
codes = []
sgnames = []
for i in range(len(tuples)):
codes.append(tuples[i][0])
sgnames.append(tuples[i][1])
codes = np.array(codes)
sgnames = np.array(sgnames)
# Loop over rows in data structure
for indx in range(data.shape[0]):
# Get information
code = data[PROJECT][indx]
name = str(data[TARGET][indx]) # str is needed since some source names are numeric
# Loop over cases
if name.find('Ephemeris') == 0 and \
abs(data[RA_DEG][indx]) < EPSILON and \
abs(data[DEC_DEG][indx]) < EPSILON:
# Find entry
j = np.where( (codes == code) & (sgnames == name) )[0]
if len(j) == 0:
raise Exception,'Could not find project code %s and source %s: see row=%d' % (code, name, getExcelFromIndex(indx))
elif len(j) > 1:
raise Exception,'Found multiple entries for project code %s and source %s' % (code, name)
# Modify name
data.loc[indx, (TARGET)] = tuples[j][2]
found_tuples[j] = True
# All rows should have been identified
if found_tuples.count(False) > 0:
print 'Did not find the following ephemeris sources:'
for j, found in enumerate(found_tuples):
if not found:
print '%s %s' % (tuples[j][0], tuples[j][1])
raise Exception,'Not all ephemeris sources were identified.'
def getIndexFromExcel(excel, check=True):
"""
Return the index number of the data structure from the excel row number.
If check = True, the resultant values is checked for errors.
"""
# Set index to account for header row and the index=1 in excel
indx = excel - 2
# Account for skipped rows. Assumes all skipped rows are before data
if SKIPROWS is not None:
indx -= len(SKIPROWS)
# Check for errors
if check and OBSERVATIONS is not None and \
(indx < 0 or indx >= OBSERVATIONS[DATA].shape[0]):
minRow = getExcelFromIndex(0, check=False)
maxRow = getExcelFromIndex(OBSERVATIONS[DATA].shape[0]-1, check=False)
raise Exception,'Excel row (%d) is out of range. Allowed range is between %d and %d' % \
(excel, minRow, maxRow)
return indx
def getExcelFromIndex(indx, check=True):
"""
Return the row number in the excel spreadsheet from the index number
in the data structure.
If check = True, the resultant valueis checked for errors.
"""
# Set index to account for header row and the index=1 in excel
excel = indx + 2
# Account for skipped rows. Assumes all skipped rows are before data
if SKIPROWS is not None:
excel += len(SKIPROWS)
# Check for errors
if check and OBSERVATIONS is not None:
minIndex = 0
maxIndex = OBSERVATIONS[DATA].shape[0]-1
minExcel = getExcelFromIndex(minIndex, check=False)
maxExcel = getExcelFromIndex(maxIndex, check=False)
if excel < minExcel or excel > maxExcel:
raise Exception,'Index (%d) is out of range. Allowed range is between %d and %d' % \
(indx, minIndex, maxIndex)
return excel
def getBandNumber(band):
"""
Returns band number from the band string in the spreadsheet.
This assumes the band format is in ALMA_RB_NN.
If there is an error, then 0 is return.
"""
try :
bn = int(band.split('_')[-1])
except:
bn = 0
return bn
def makeList(sources, parse=True, delimiter=','):
""" Convert the variable "sources" into a list.
Input : sources - a variable or list
parse - if True, parse elements with delimiter separator
delimiter - the delimiter when parsing strings
Output: Return value is a list.
If input is a single value, then output is [sources]
If input is a list, the result is copied
Examples:
(1) l = makeList('3c273')
(2) l = makeList('3c273, 3c279')
(3) l = makeList(['3c273','3c279'])
"""
# Return if sources is empty
if sources == None: return None
# Convert to arrays
l = list()
t = [type(sources)]
if str in t or unicode in t:
l = [sources]
else:
try:
i = len(sources)
l = np.array(sources)
except:
l = np.array([sources])
if len(l) == 0: return None
# Parse string
if parse:
t = list()
for z in l:
if str in [type(z)]:
for zz in z.split(delimiter):
t.append(re.sub(' +',' ', zz.strip()))
else:
t.append(z)
l = np.array(t[:])
return l
def convertHmsString(value, ndec=0, showSeconds=True, delimiter=':'):
"""
Converts floating point to HH:MM:[SS.S]
Inputs
value : floating point number
ndec : number of decimal points to print seconds
showSeconds : If True, show seconds, otherwise just
use HH:MM. Default:True
delimiter : the characters to be used between the numbers.
Output
a string of format hh:mm:[ss.s]
Examples:
convertHmsString(15.0)
convertHmsString(15.0, ndec=2)
convertHmsString(15.0, ndec=2, delimiter='hms')
convertHmsString(15.0, ndec=2, delimiter='dms')
"""
# Set delimiter
spacing = delimiter
if len(spacing) == 1: spacing = [spacing] * 3
# Construct string
t = value
st = str(type(value))
if st.find('int') >= 0 or st.find('float') >= 0:
x = abs(value)
h = int(x)
m = int(60 * (x - h))
sec = 3600.0 * (x - h - m/60.0)
t = str("%.2d" % h) + spacing[0] + str('%.2d' % m)
if showSeconds :
# Add seconds
t += spacing[1]
format = '%0' + str(ndec+3) + '.' + str(ndec) + 'f'
if ndec <= 0: format = '%.2d'
t += str(format % sec)
if spacing[2] not in [' ', ':']: t += spacing[2]
if value < 0.0: t = '-' + t
return t
def computeZ(data, indx):
"""
Computes redshift based on velocity for row indx in data.
"""
# Get velocity
velocity = data[VELOCITY][indx]
velocity_unit = data[VELOCITY_UNIT][indx]
doppler = data[VELOCITY_DOP][indx]
# Convert velocity to km/s
if velocity_unit == 'km/s':
v_kms = velocity
elif velocity_unit == 'm/s':
v_kms = velocity / 1e3
else:
raise Exception,'Velocity unit not recognized: %s' % velocity_unit
# Compute redshift
c_kms = G.c.value / 1e3
if doppler == 'OPTICAL':
z = v_kms / c_kms
elif doppler == 'RADIO':
z = v_kms/c_kms / (1.0 - v_kms/c_kms)
elif doppler == 'RELATIVISTIC':
z = np.sqrt( (1. + v_kms/c_kms) / (1.0 - v_kms/c_kms) ) - 1.0
else:
raise Exception,'Doppler frame not recognized: %s' % doppler
# Done
return z
def fwhmPB(freqGHz, diameter):
"""
Returns primary FWHM in arcsec
freqGHz : frequency in GHz
diameter: telescope diameter in meters
See the knowledgebase article for more information:
https://help.almascience.org/index.php?/Knowledgebase/Article/View/90/0/how-do-i-model-the-alma-primary-beam-and-how-can-i-use-that-model-to-obtain-the-sensitivity-profile-for-an-image-mosaic
"""
return 1.13 * G.c.value / (freqGHz*1e9) / diameter / np.pi * 180.0 * 3600.0
def plotMosaic(ax, corners, fc='black', ec='None',
linewidth=2, alpha=0.5, hatch=None):
"""
Plot a rectangular region for a mosaic.
alpha : transperancy of the plot rectangle
ax : pylab plot handle
corners : dictionary containing the corners of the mosaic, specific
in RA,Dec offsets in arcseconds relative to the (0,0). It
is best to use getMosaicCorners() to compute corner
positions.
ec : edge color of the plot rectangle
fc : solid face color of the rectangle
hatch : hatch pattern for the rectangle
linewidth : linewidth of the plot edge
"""
UL = corners['UL']
UR = corners['UR']
BL = corners['BL']
BR = corners['BR']
xypolygon = np.zeros( (4, 2) )
xypolygon[:,0] = np.array([UL[0], UR[0], BR[0], BL[0]])
xypolygon[:,1] = np.array([UL[1], UR[1], BR[1], BL[1]])
p = Polygon(xypolygon, closed=True, fc=fc, ec=ec, alpha=alpha, linewidth=linewidth, hatch=hatch)
result = ax.add_artist(p)
# ax.plot([UL[0], UR[0], BR[0], BL[0], UL[0]],
# [UL[1], UR[1], BR[1], BL[1], UL[1]],
# color=ec, linewidth=linewidth)
return result
def plotMosaicOld(ax, xy, width, height, angle, color='black', ec='None',
linewidth=2, alpha=0.5, hatch=None):
"""
Plot a rectangular region for a mosaic.
alpha : transperancy of the plot rectangle
angle : position angle of the rectangle in degrees
ax : pylab plot handle
color : solid face color of the rectangle
ec : edge color of the plot rectangle
hatch : hatch pattern for the rectangle
height : height of the rectangle. Note this is not the same
definition as the OT
linewidth : linewidth of the plot edge
width : width of the rectangle. Note this is not the same
definition as the OT
xy : center of the mosaic
"""
x = xy[0]
y = xy[1]
# A = (-90+angle) / 180.0 * np.pi
A = -angle / 180.0 * np.pi
UL = (x + ( width / 2. ) * np.cos(A) - ( height / 2. ) * np.sin(A), y + ( height / 2. ) * np.cos(A) + ( width / 2. ) * np.sin(A))
UR = (x - ( width / 2. ) * np.cos(A) - ( height / 2. ) * np.sin(A), y + ( height / 2. ) * np.cos(A) - ( width / 2. ) * np.sin(A))
BL = (x + ( width / 2. ) * np.cos(A) + ( height / 2. ) * np.sin(A), y - ( height / 2. ) * np.cos(A) + ( width / 2. ) * np.sin(A))
BR = (x - ( width / 2. ) * np.cos(A)+ ( height / 2. ) * np.sin(A), y - ( height / 2. ) * np.cos(A) - ( width / 2. ) * np.sin(A))
xypolygon = np.zeros( (4, 2) )
xypolygon[:,0] = np.array([UL[0], UR[0], BR[0], BL[0]])
xypolygon[:,1] = np.array([UL[1], UR[1], BR[1], BL[1]])
p = Polygon(xypolygon, closed=True, color=color, ec=ec, alpha=alpha, linewidth=linewidth, hatch=hatch)
result = ax.add_artist(p)
# ax.plot([UL[0], UR[0], BR[0], BL[0], UL[0]],
# [UL[1], UR[1], BR[1], BL[1], UL[1]],
# color=ec, linewidth=linewidth)
return result
def plotPrimaryBeam(ax, xy, freq, diameter, fc='black', ec='black',
alpha=1, linewidth=2):
"""
Plot the primary beam FWHM for an observation at frequency freq in GHz
for a telescope diameter in meters.
ax : pylab plot handle
xy : center of the primary beam
fc : face color of the plot circle
ec : edge color of the plot circle
alpha : transperancy of the plot circle
linewidth : linewidth of the plot edge
"""
radius = 0.5 * fwhmPB(freq, diameter)
c = Circle( xy, radius=radius, fc=fc, ec=ec, alpha=alpha, linewidth=linewidth)
result = ax.add_artist(c)
# color = 'yellow'
# e = Ellipse( (0, 0), width=2*radius, height=2*radius, angle=0,
# color=color, fc=color, ec='black')
# result = ax.add_artist(e)
return result
def getBandColor(freq):
"""
Set a plot color based on the input frequency (GHz),
which a different color is chosen per ALMA band.
"""
# Set band colors for plotting
band_colors = dict()
band_colors['ALMA_RB_01'] = 'gray'
band_colors['ALMA_RB_02'] = 'silver'
band_colors['ALMA_RB_03'] = 'blue'
band_colors['ALMA_RB_04'] = 'peru'
band_colors['ALMA_RB_05'] = 'yellow'
band_colors['ALMA_RB_06'] = 'black'
band_colors['ALMA_RB_07'] = 'green'
band_colors['ALMA_RB_08'] = 'orange'
band_colors['ALMA_RB_09'] = 'magenta'
band_colors['ALMA_RB_10'] = 'red'
band_colors[BAND_UNKNOWN] = 'tan'
# Find color
band_input = None
for key, freq_range in BAND_FREQ.iteritems():
if freq >= freq_range[0] and freq <= freq_range[1]:
band_input = key
break
if band_input is None:
print 'Warning: Frequency outside of ALMA bands.'
band_input = BAND_UNKNOWN
# Done
return band_colors[band_input]
def checkObservations(obsdata, isdata=False):
"""
Check that the user-supplied observation seems reasonable.
Only a light check is done.
obsdata is either "observations" from readObservations (isdata=False)
or observations[DATA] (if isdata=True)
"""
# Checks if isdata=False
if isdata:
# Must be panda structure
if pd.core.frame.DataFrame not in [type(obsdata)]:
raise Exception,'observations must be a dictionary'
# Get keys
keys = obsdata.keys().tolist()
msg = 'data'
else:
# Must have two keywords, if observations
if dict not in [type(obsdata)]:
raise Exception,'observations must be a dictionary'
# Check keys
for key in [DATA, COORDS]:
if not obsdata.has_key(key):
raise Exception,'%s keyword not found in observations' % key
# Get keys
keys = obsdata[DATA].keys().tolist()
msg = 'observations[%s]' % DATA
# Check entries in OBS_ITEMS
for key in OBS_ITEMS.keys():
if keys.count(key) != 1:
raise Exception,'Keyword "%s" not found in %s' % (key, msg)
# Check spectral keywords
for w in SPW_FREQ.keys():
if SPW_FREQ[w] not in keys:
raise Exception,'Keyword "%s" not found in %s' % (SPW_FREQ[w], msg)
if SPW_BW[w] not in keys:
raise Exception,'Keyword "%s" not found in %s' % (SPW_BW[w], msg)
if SPW_RES[w] not in keys:
raise Exception,'Keyword "%s" not found in %s' % (SPW_RES[w], msg)
# Done
return True
def checkData(data, verbose=True, spreadsheet=None):
"""
Run checks on the data read from the Excel spreadsheet to
catch any obvious anomalies.
"""
# Message
if verbose:
print 'Running checks on data'
# Initialize
error = False
# Check keywords
if verbose:
print ' ... checking keywords in data structure'
checkObservations(data, isdata=True)
# Check RA
if verbose:
print ' ... checking right ascensions'
j = np.where( (data[RA_DEG] < 0) | (data[RA_DEG] > 360.0) )[0]
if len(j) > 0:
error = True
if verbose: print ' --- invalid right ascension in %d rows' % len(j)
# Check Declination
if verbose:
if verbose: print ' ... checking declinations'
j = np.where( (data[DEC_DEG] < -90.) | (data[DEC_DEG] > 90.0) )[0]
if len(j) > 0:
error = True
if verbose:
print ' --- invalid declination in %d rows' % len(j)
# Check mosaic offsets
if verbose:
print ' ... checking mosaic parameters are present for rectangular mosaics'
j = np.where(data[IS_MOSAIC] == MOSAIC_TYPE_RECTANGLE)[0]
# if verbose:
# print ' --- found %d rectangular mosaics' % len(j)
if len(j) > 0:
# Check mosaic system
k = np.where(data[MOS_COORD][j] != data[MOS_COORD][j])[0]
if len(k) > 0:
error = True
if verbose:
print ' --- found %d mosaics where mosaic system was not set' % len(k)
# Check PA
k = np.where(np.isnan(data[MOS_PA][j]) == True)[0]
if len(k) > 0:
error = True
if verbose:
print ' --- found %d mosaics where PA was not set' % len(k)
# Check length
k = np.where(np.isnan(data[MOS_LENGTH][j]) == True)[0]
if len(k) > 0:
error = True
if verbose:
print ' --- found %d mosaics where mosaic length was not set' % len(k)
k = np.where(np.abs(data[MOS_LENGTH][j]) < EPSILON)[0]
if len(k) > 0:
error = True
if verbose:
print ' --- found %d mosaics where mosaic length was zero' % len(k)
# Check width
k = np.where(np.isnan(data[MOS_WIDTH][j]) == True)[0]
if len(k) > 0:
error = True
if verbose:
print ' --- found %d mosaics where mosaic width was not set' % len(k)
k = np.where(np.abs(data[MOS_WIDTH][j]) < EPSILON)[0]
if len(k) > 0:
error = True
if verbose:
print ' --- found %d mosaics where mosaic width was zero' % len(k)
# Check if reference sensitivity is valid
if verbose:
print ' ... checking values of requested sensitivity are non-zero'
j = np.where( (data[REQ_SENS] <= 0.) | (np.isnan(data[REQ_SENS]) == True))[0]
if len(j) > 0:
error = True
if verbose:
print ' --- invalid sensitivity values in %d rows' % len(j)
# Check if reference bandwidth is valid
if verbose:
print ' ... checking values of reference bandwidth are non-zero'
j = np.where( (data[REF_BW] <= 0.) | (np.isnan(data[REF_BW]) == True)) [0]
if len(j) > 0:
error = True
if verbose:
print ' --- invalid reference bandwidth values in %d rows' % len(j)
# Check spectral scans
if verbose:
print ' ... checking spectral scan frequencies and bandwidths'
k = np.where(data[IS_SPECTRAL] == True)[0]
for keys in [ [SPS_START, 'starting frequency'],
[SPS_END, 'ending frequency'],
[SPS_BW, 'bandwidth'],
[SPS_SPW_RES, 'spectral resolution'],
]:
j = np.where( (np.isnan(data[keys[0]][k]) == True) | (data[keys[0]][k] <= 0) )[0]
if len(j) > 0:
error = True
if verbose:
print ' --- spectral scan %s is not set in %d rows' % (keys[1], len(j))
# Check that the frequencies match the bands
if verbose:
print ' ... checking values of frequencies and bandwidths'
for indx in range(data.shape[0]):
# Get range of frequencies defined this ALMA band
freq_range = BAND_FREQ[data[BAND][indx]]
# Check reference frequency
if data[REF_FREQ][indx] < freq_range[0] or \
data[REF_FREQ][indx] > freq_range[1]:
error = True
if verbose:
print ' --- reference frequency out of range in row %d' % getExcelFromIndex(indx, check=False)
# Check spectral scan frequencies
if data[IS_SPECTRAL][indx]:
# Starting frequency
if data[SPS_START][indx] < freq_range[0] or \
data[SPS_START][indx] > freq_range[1]:
error = True
if verbose:
print ' --- spectral scan starting frequency out of range in row %d' % getExcelFromIndex(indx, check=False)
# Ending frequency
if data[SPS_END][indx] < freq_range[0] or \
data[SPS_END][indx] > freq_range[1]:
error = True
if verbose:
print ' --- spectral scan ending frequency out of range in row %d' % getExcelFromIndex(indx)
else:
# Check each window
for w in SPW_FREQ.keys():
# Check frequencies
nu = data[SPW_FREQ[w]][indx]
if np.isfinite(nu) and (nu < freq_range[0] or nu > freq_range[1]):
error = True
if verbose:
print ' --- spectral window frequency out of range in window %d on row %d' % \
(w, getExcelFromIndex(indx, check=False))
# Both frequencies and bandwidths need to be set
bw = data[SPW_BW[w]][indx]
if (np.isfinite(nu) and np.isnan(bw)) or \
(np.isfinite(bw) and np.isnan(nu)):
error = True
if verbose:
print ' --- frequencies and bandwidths are not both set in window %d on row %d' % \
(w, getExcelFromIndex(indx, check=False))
# Check usable bandwidth
if np.isfinite(nu) and getUsableBandwidth(bw, indx=indx) is None:
error = True
if verbose:
print ' --- bandwidth is not recognized in window %d on row %d' % \
(w, getExcelFromIndex(indx, check=False))
# Exit if there were any errors
if verbose:
print ' ... done'
if error:
if spreadsheet is None:
msg = 'The spreadsheet contains unexpected errors.'
else:
msg = 'The spreadsheet %s contains unexpected errors.' % spreadsheet
if not verbose:
msg += '\nTry with verbose=True to see further information'
raise Exception, msg
# Done
return
def readObservations(input=LIST_OF_OBSERVATIONS, verbose=True):
"""
Read in observations using pandas.
A global variable called SKIPROWS indicate which rows
to skip in reading the data file.
Inputs:
input : root name for the file containing the existing and
scheduled observations. The .csv or .xls extension can
be omitted.
Output:
A python dictionary with the following entries:
DATA : a panda data set containing the spreadsheet
COORDS : astropy sky coordinates for each source in DATA
"""
# Read in observations. Input list may either be an excel spreadsheet or
# a csv file. Both are tried, with the .csv first and then the .xls.
data = None
input_without_extension = input.replace('.xls', '').replace('.csv', '')
if verbose:
print 'Reading observations'
for ext in ['csv', 'xls']:
# Set file name
inputFile = '%s.%s' % (input_without_extension, ext)
# Try reading file
if os.path.exists(inputFile):
if ext == 'csv':
data = pd.read_csv(inputFile, skiprows=SKIPROWS, header=0, low_memory=False)
elif ext == 'xls':
data = pd.read_excel(inputFile, skiprows=SKIPROWS, header=0)
else:
raise Exception,'Do not recognize extension for file %s' % input
# Did we succeed?
if data is not None:
if verbose:
print ' ... read %d rows in %s' % (data.shape[0], inputFile)
break
# Check if data were read
if data is None:
raise Exception,'Could not read data. Check that the input file %s exists with either a .csv for .xls extension.' % input
# Delete column containing RA/DEC in string form since RA/DEC may
# be modified below.
del data[RA_HMS]
del data[DEC_DMS]
# Correct the RA/DEC for sources with an RA/DEC offset. Note the offsets
# may be given in galactic coordinates for mosaics.
if verbose:
print ' ... correcting centroid RA/DEC in mosaics for offsets'
offset_type = data[MOS_COORD]
k = np.where(offset_type != offset_type)[0]
offset_type = offset_type.values
offset_type[k] = 'J2000'
j = np.where( (np.isfinite(data[LAT_OFFSET]) == True) & \
(np.isfinite(data[LON_OFFSET]) == True) & \
( (np.abs(data[LON_OFFSET]) > EPSILON) |
(np.abs(data[LAT_OFFSET]) > EPSILON)) ) [0]
if len(j) > 0:
# Get system
coord_system = np.unique(offset_type[j])
# Correction depends on the coordinate system
for cs in coord_system:
# Find observations with this system
k = np.where(offset_type[j] == cs)[0]
# Get data
x = data.loc[j[k], (RA_DEG)].values
y = data.loc[j[k], (DEC_DEG)].values
dx = data.loc[j[k], (LON_OFFSET)].values
dy = data.loc[j[k], (LAT_OFFSET)].values
# Correct coordinates
if cs == 'J2000':
# Message
# if verbose:
# print ' --- correcting %4d positions for equatorial offsets' % len(k)
# Correct RA/Dec
new_ra = x + dx / 3600.0 / np.cos(y / 180.0 * np.pi)
new_dec = y + dy / 3600.0
pass
elif cs == 'galactic':
# Message
# if verbose:
# print ' --- correcting %4d mosaic positions for galactic offsets' % len(k)
# Convert mosaic center to galactic coordinates
c = SkyCoord(ra=x, dec=y, frame='icrs', unit='degree')
glon = c.galactic.l.deg
glat = c.galactic.b.deg
# Add offset
glon += dx / 3600.0 / np.cos(glat / 180.0 * np.pi)
glat += dy / 3600.0
# Convert back to RA/DEC
c = SkyCoord(glon, glat, unit='degree', frame='galactic')
new_ra = c.icrs.ra.deg
new_dec = c.icrs.dec.deg
else:
raise Exception,'Mosaic coordinate system not implemented: %s' % cs
# Check right ascension
l = np.where(new_ra < 0)[0]
if len(l) > 0:
new_ra[l] += 360.0
l = np.where(new_ra > 360)[0]
if len(l) > 0:
new_ra[l] -= 360.0
# Check declination
if np.any(np.abs(new_dec) > 90.0):
j = np.where(np.abs(new_dec) > 90)
raise Exception,'Error setting mosaic declination'
# Save
data.loc[j[k], (RA_DEG)] = new_ra
data.loc[j[k], (DEC_DEG)] = new_dec
# Convert RA/DEC to degree in sky coordinates
if verbose:
print ' ... converting RA/DEC to astropy sky coordinates'
coords = SkyCoord(data[RA_DEG]*astropy.units.degree, data[DEC_DEG]*astropy.units.degree, frame='icrs')
# Set the velocity unit.
# Until now, all SG have used km/s as the velocity unit. but I want to keep
# the flexibility in case m/s is used in the future.
data[VELOCITY_UNIT] = np.array(['km/s'] * data.shape[0])
# Compute redshift
if verbose:
print ' ... computing redshifts'
redshifts = np.zeros(data.shape[0])
for n in range(redshifts.size):
redshifts[n] = computeZ(data, n)
# Set variable indicating if observation is a spectral scan
data[IS_SPECTRAL] = np.isfinite(data[SPS_START])
# Correct frequencies for spectral windows that have rest frequencies
if verbose:
print ' ... correcting rest frequencies to sky frequencies'
jndx = dict()
znu = dict()
for w in SPW_FREQ.keys():
jndx[w] = []
znu[w] = []
for indx in range(data.shape[0]):
# Determine if windows are doppler corrected or not
if data[IS_SPECTRAL][indx]:
if data[IS_SPW_SKY_FREQ][indx] == True:
pass
else:
raise Exception,'Not expecting spectral scan to be rest frequencies: excel line = %d' % getExcelFromIndex(indx, check=False)
else:
# Determine if windows are sky or rest frequencies
if data[IS_SPW_SKY_FREQ][indx] not in [0, 1, False, True]:
raise Exception,'Error reading IS_SPW_SKY_FREQ on row %d' % \
getExcelFromIndex(indx)
elif not data[IS_SPW_SKY_FREQ][indx]:
# Loop over windows and correct frequencies
for w in SPW_FREQ.keys():
nu = data[SPW_FREQ[w]][indx]
if np.isfinite(nu):
jndx[w].append(indx)
znu[w].append(nu / (1.0 + redshifts[indx]))
for w in SPW_FREQ.keys():
if len(jndx[w]) > 0:
data.loc[jndx[w],(SPW_FREQ[w])] = znu[w]
# Since all windows are now doppler corrected, delete IS_SKY_FREQ columns
del data[IS_SPW_SKY_FREQ]
# Get FWHM primary beam size
freq = data[REF_FREQ].values # GHz
j = np.where( (freq < 0) | (np.isnan(freq) == True) )[0]
if len(j) < 0:
raise Exception,'Invalid frequencies in spreadsheet'
diameter = 12.0
data[FWHM_PB] = fwhmPB(freq, diameter)
# Set largest size in arcseconds, including mosaics and FWHM of primary beam
data[MAX_SIZE] = data[ [MOS_LENGTH, MOS_WIDTH, FWHM_PB] ].max(axis=1)
# Compute mosaic area
if verbose:
print ' ... computing area of rectangular mosaics'
mosarea = np.zeros_like(data[MOS_LENGTH])
j = np.where(data[IS_MOSAIC] == MOSAIC_TYPE_RECTANGLE)[0]
if len(j) > 0:
mosarea[j] = data[MOS_LENGTH][j] * data[MOS_WIDTH][j] / 3600.0
data[MOS_AREA] = mosarea
# Modify ephemeris names
if verbose:
print ' ... checking names for ephemeris sources'
checkEphemerisNames(data)
# Message
if verbose:
print ' ... done'
# Run checks on the data structure
checkData(data, verbose=verbose, spreadsheet=input)
# Done
return {DATA: data, COORDS:coords}
def computeContinuumSensitivity(data, indx):
"""
Compute the continuum sensitivity for a single pointing
For spectral scans, the sensitivity is computed for one tuning,
not the full spectral scan.
Inputs:
data : data structure from readObservations()
indx : Index number in data
Output:
A dictionary with the following keys:
'bw' : aggregrate bandwidth in MHz
'rmsmJy': continuum rms in mJy
'rmsmK' : continuum rms in mK for the specified angular
resolution and reference refrequencin data
"""
# Gather sensitivity information
req_sen = data[REQ_SENS][indx]
ref_freq = data[REF_FREQ][indx]
ref_bw = data[REF_BW][indx]
angres = data[ANG_RES][indx]
# Compute bandwidth if single tuning or spectral scans
if data[IS_SPECTRAL][indx]:
# Get spectral scan settings
sps_bw = getUsableBandwidth(data[SPS_BW][indx])
nwin = 4
aggregateBandwidth = nwin * sps_bw
else:
# Get frequencies of each window
windows = SPW_FREQ.keys()
nu1 = np.zeros(len(windows))
nu2 = np.zeros(len(windows))
use = np.zeros(len(windows), dtype=bool)
for i, w in enumerate(windows):
# Is this a valid window?
if np.isnan(data[SPW_FREQ[w]][indx]):
use[i] = False
continue
freq = data[SPW_FREQ[w]][indx] * 1e3
bw = getUsableBandwidth(data[SPW_BW[w]][indx])
use[i] = True
f1 = freq - 0.5 * bw
f2 = freq + 0.5 * bw
nu1[i] = min([f1, f2])
nu2[i] = max([f1, f2])
# Sort by decreasing bandwidth
j = np.where(nu1 > 0)
nu1 = nu1[j]
nu2 = nu2[j]
bw = nu2 - nu1
iarg = np.argsort(bw)[::-1]
nu1 = nu1[iarg]
nu2 = nu2[iarg]
# Compute aggregate bandwidth
j = np.where(nu1 > 0)
for iw in range(nu1.size):
# Skip if not using the window
if not use[iw]:
continue
# Check other windows
for jw in range(iw+1, nu1.size):
# Skip name windows
if not use[jw]:
continue
# Check frequency range
if nu1[jw] >= nu1[iw] and nu2[jw] <= nu2[iw]:
nu1[jw] = 0
nu2[jw] = 0
use[jw] = False
elif nu1[jw] < nu1[iw] and nu2[jw] > nu2[iw]:
raise Exception,'Should not be here since I sorted by decreasing bandwidth. Indx=%d' % indx
elif nu1[jw] >= nu2[iw] or nu2[jw] <= nu1[iw]:
pass
elif nu1[jw] >= nu1[iw] and nu1[jw] <= nu2[iw] and nu2[jw] >= nu2[iw]:
nu1[jw] = nu2[iw]
if nu1[jw] == nu2[jw]:
use[jw] = False
elif nu2[jw] < nu1[jw]:
raise Exception,'Error setting bandwidth'
elif nu1[jw] <= nu1[iw] and nu2[jw] >= nu1[iw] and nu2[jw] <= nu2[iw]:
nu2[jw] = nu1[iw]
if nu1[jw] == nu2[jw]:
use[jw] = False
elif nu2[jw] < nu1[jw]:
raise Exception,'Error setting bandwidth'
else:
raise Exception,'Error computing aggregate bandwidth'
# Compute aggregate bandwidth in MHz
aggregateBandwidth = np.sum(nu2 - nu1)
# Compute continuum sensitivity
if ref_bw <= 0 or aggregateBandwidth <= 0:
rmsContinuum = np.nan
raise Exception,'Error computing continuum sensitivity'
else:
rmsContinuum_mJy = req_sen * np.sqrt(ref_bw / aggregateBandwidth)
rmsContinuum_mK = rmsContinuum_mJy / mjypermk(ref_freq, angres, ref_freq)
# Done
return {'bw': aggregateBandwidth, 'rmsmJy': rmsContinuum_mJy,
'rmsmK': rmsContinuum_mK}
def printRowProject(data, indx, spaces=''):
"""
Print project information for entry indx in "data".
"""
# Functions to print data
def printEntries(obsdata, indx, items):
for key in items:
a = OBS_ITEMS[key]
printEntry(obsdata, indx, key, a[0], a[1])
def printEntry(obsdata, indx, key, label, unit):
# Set unit
sunit = unit
if sunit is None: sunit = ''
# Set value - there must be a better way!
if key in [FWHM_PB, REF_FREQ, REF_BW]:
svalue = '%.1f' % obsdata[key][indx]
elif key in [MOS_LENGTH, MOS_WIDTH]:
svalue = '%.1f' % obsdata[key][indx]
elif key in [MOS_AREA]:
svalue = '%.3f' % obsdata[key][indx]
elif key in [REQ_SENS, ANG_RES]:
svalue = '%.3f' % obsdata[key][indx]
elif key in [RA_DEG]:
svalue = convertHmsString(obsdata[key][indx] / 15.0, delimiter=':', ndec=2)
elif key in [DEC_DEG]:
svalue = convertHmsString(obsdata[key][indx], delimiter=':', ndec=2)
elif key == BAND:
svalue = '%d' % getBandNumber(obsdata[key][indx])
elif key in [POLARIZATION]:
svalue = obsdata[key][indx].lower()
elif key in [MOS_SPACING]:
svalue = '%.1f' % obsdata[key][indx]
elif key == IS_MOSAIC:
if isObsMosaic(obsdata, indx, mtype=MOSAIC_TYPE_RECTANGLE):
svalue = 'Rectangular mosaic'
elif isObsMosaic(obsdata, indx, mtype=MOSAIC_TYPE_CUSTOM):
svalue = 'Custom mosaic'
elif obsdata[IS_MOSAIC][indx] == 'N/A' or np.isnan(obsdata[IS_MOSAIC][indx]):
svalue = 'False'
else:
raise Exception,'Unknown type of mosaic: %s' % data[IS_MOSAIC][indx]
else:
svalue = '%s' % obsdata[key][indx]
# Print item
print '%s%-22s %-13s %s' % (spaces, label, svalue, sunit)
# Add some space
print ''
print ''
# Print header
lineExcel = getExcelFromIndex(indx)
print '%sSource information for spreadsheet line %d (project = %s)' % \
(spaces, lineExcel, data[PROJECT][indx])
# Print source information
items = [TARGET,
RA_DEG,
DEC_DEG,
]
printEntries(data, indx, items)
# Print observing parameters
print ''
print '%sObserving parameters' % spaces
items = [BAND,
FWHM_PB,
ANG_RES,
LAS,
]
printEntries(data, indx, items)
# Print observing modes
print ''
print '%sObserving modes' % spaces
items = [POLARIZATION, USE_ACA, USE_TPA, IS_SPECTRAL]
printEntries(data, indx, items)
# Print mosaic information
print ''
print '%sMosaic information' % spaces
items = [IS_MOSAIC]
if isObsMosaic(data, indx, mtype=MOSAIC_TYPE_RECTANGLE):
items.extend([MOS_AREA, MOS_LENGTH, MOS_WIDTH, MOS_PA, MOS_SPACING, MOS_COORD])
printEntries(data, indx, items)
# Print continuum sensitivity if not a spectral scan
if not data[IS_SPECTRAL][indx]:
print ''
print '%sEstimated continuum sensitivity' % spaces
result = computeContinuumSensitivity(data, indx)
print '%s%-22s %-13.0f %s' % \
(spaces, 'Aggregate bandwidth', result['bw'], 'MHz')
if np.isnan(result['rmsmJy']):
print '%s%-22s %-13s %s' % \
(spaces, 'Continuum RMS', 'N/A', '')
else:
print '%s%-22s %-13.3f %s' % \
(spaces, 'Continuum RMS', result['rmsmJy'], 'mJy')
print '%s%-22s %-13.1f %s' % \
(spaces, 'Continuum RMS', result['rmsmK'], 'mK')
def mjypermk(freq_ghz, angres_ref, freq_ref_ghz):
"""
Returns the conversion factor to convert mK to mJy.
"""
omega_ref = np.pi / 4.0 / np.log(2.0) * (angres_ref / 3600.0 / 180.0 * np.pi)**2
omega = omega_ref * (freq_ref_ghz / freq_ghz)**2
wavelength = G.c.value / (freq_ghz * 1e9)
constant = 2.0 * G.k_B.value * 1e-3 / wavelength**2 * omega * 1e29
return constant
def printRowCorrelator(data, indx, spaces=''):
"""
Print correlator configuration for a row entry observations.
The output will be different for spectral scans and single tunings.
Inputs:
data : Contains data from readObservations() in the DATA key
indx : Entry number in observations
spaces : Spaces to format output
"""
if data[IS_SPECTRAL][indx]:
printRowCorrelatorSpectralScan(data, indx, spaces=spaces)
else:
printRowCorrelatorSingle(data, indx, spaces=spaces)
def printRowCorrelatorSingle(data, indx, spaces=''):
"""
Print project information for entry indx in "observations" if it
is a single tuning.
Inputs:
data : Contains data from readObservations() in the DATA key
indx : Entry number in observations
spaces : Spaces to format output
"""
# Add some space
print ''
# Print header
lineExcel = getExcelFromIndex(indx)
print '%sCorrelator setup' % spaces
# Print header
print '%s%3s %8s %17s %19s %17s %17s' % \
(spaces, 'Win', 'Sky Freq', 'Usable Bandwidth', 'Spectral resolution',
'RMS/bandwidth', 'RMS/resolution')
print '%s%3s %8s %17s %19s %17s %17s' % \
(spaces, '---', '--------', '-----------------', '-------------------',
'----------------', '-----------------')
print '%s%3s %8s %8s %8s %9s %9s %8s %8s %8s %8s' % \
(spaces, '', '(GHz)', '(MHz)', '(km/s)', '(MHz)', '(km/s)',
'mJy', 'mK', 'mJy', 'K'
)
# Gather sensitivity information
req_sen = data[REQ_SENS][indx]
ref_freq = data[REF_FREQ][indx]
ref_bw = data[REF_BW][indx]
# Get frequencies
windows = SPW_FREQ.keys()
winfreq = []
winnumber = []
for w in windows:
if np.isfinite(data[SPW_FREQ[w]][indx]):
winnumber.append(w)
winfreq.append(data[SPW_FREQ[w]][indx])
# Sort by frequency
iarg = np.argsort(winfreq)
winfreq = np.array(winfreq)[iarg]
winnumber = np.array(winnumber)[iarg]
# Loop over windows
for w in winnumber:
# Is this a valid window?
if np.isnan(data[SPW_FREQ[w]][indx]):
continue
# Gather spectral window information
freq = data[SPW_FREQ[w]][indx]
bw = getUsableBandwidth(data[SPW_BW[w]][indx])
res = data[SPW_RES[w]][indx]
# Print window
print '%s%3d' % (spaces, w),
# Frequency
print '%8.3f' % freq,
# Bandwidth in MHz
print ' %8.1f' % bw,
# Bandwidth in km/s
bw_kms = bw * 1e-3 / freq * G.c.value / 1e3
print '%8.1f' % bw_kms,
# Channel resolution in MHz
print ' %9.3f' % res,
# Channel resolution in km/s
res_kms = res * 1e-3 / freq * G.c.value / 1e3
print '%9.3f' % res_kms,
# Sensitivity per bandwidth and channel
if ref_bw == 0:
srms_bw_mjy = '%8s' % 'N/A'
srms_bw_mk = '%8s' % 'N/A'
srms_res_mjy = '%8s' % 'N/A'
srms_res_k = '%8s' % 'N/A'
else:
# In mJy
rms_bandwidth_mJy = req_sen * np.sqrt(ref_bw / bw)
rms_resolution_mJy = req_sen * np.sqrt(ref_bw / res)
# In mK
angres = data[ANG_RES][indx]
rms_bandwidth_mK = rms_bandwidth_mJy / mjypermk(freq, angres, ref_freq)
rms_resolution_K = rms_resolution_mJy / mjypermk(freq, angres, ref_freq) / 1e3
# Save as strings
srms_bw_mjy = '%8.3f' % rms_bandwidth_mJy
srms_bw_mk = '%8.3f' % rms_bandwidth_mK
srms_res_mjy = '%8.3f' % rms_resolution_mJy
srms_res_k = '%8.3f' % rms_resolution_K
print ' %s' % srms_bw_mjy,
print '%s' % srms_bw_mk,
print ' %s' % srms_res_mjy,
print '%s' % srms_res_k,
# Done width entry
print ''
def printRowCorrelatorSpectralScan(data, indx, spaces=''):
"""
Print project information for entry indx in "observations"
that is a spectral scan.
Inputs:
data : Contains data from readObservations() in the DATA key
indx : Entry number in data
spaces : Spaces to format output
"""
# Add some space
print ''
# Print header
lineExcel = getExcelFromIndex(indx)
print '%sCorrelator information for spectral scan' % spaces
# Gather sensitivity information
req_sen = data[REQ_SENS][indx]
ref_freq = data[REF_FREQ][indx]
ref_bw = data[REF_BW][indx]
# Get spectral window information
sps_start = data[SPS_START][indx]
sps_end = data[SPS_END][indx]
sps_bw = getUsableBandwidth(data[SPS_BW][indx])
sps_res = data[SPS_SPW_RES][indx]
# Compute quantities
sps_res_kms = sps_res/1e3 / sps_start * G.c.value / 1e3
# Sensitivity per bandwidth and resolution element
if ref_bw == 0:
srms_bw_mjy = '%8s' % 'N/A'
srms_bw_mk = '%8s' % 'N/A'
srms_res_mjy = '%8s' % 'N/A'
srms_res_k = '%8s' % 'N/A'
else:
# In mJy
rms_bandwidth_mJy = req_sen * np.sqrt(ref_bw / sps_bw)
rms_resolution_mJy = req_sen * np.sqrt(ref_bw / sps_res)
# In mK
angres = data[ANG_RES][indx]
rms_bandwidth_mK = rms_bandwidth_mJy / mjypermk(sps_start, angres, ref_freq)
rms_resolution_K = rms_resolution_mJy / mjypermk(sps_start, angres, ref_freq) / 1e3
# Save as strings
srms_bw_mjy = '%8.2f' % rms_bandwidth_mJy
srms_bw_mk = '%8.2f' % rms_bandwidth_mK
srms_res_mjy = '%8.3f' % rms_resolution_mJy
srms_res_k = '%8.3f' % rms_resolution_K
# Print spectral scale information
print '%sStarting frequency : %8.3f GHz' % (spaces, sps_start)
print '%sEnding frequency : %8.3f GHz' % (spaces, sps_end)
print '%sUsable bandwidth per window : %8.3f MHz' % (spaces, sps_bw)
print '%sSpectral resolution : %8.3f MHz' % (spaces, sps_res)
print '%sSpectral resolution : %8.3f km/s @ %.1f GHz' % (spaces, sps_res_kms, sps_start)
print '%sRMS per spectral resolution : %s mJy' % (spaces, srms_res_mjy)
print '%sRMS per spectral resolution : %s K ' % (spaces, srms_res_k)
def isObsMosaic(data, indx, mtype=MOSAIC_TYPES):
"""
Returns True or False if an observation is a mosaic of a
type listed in MOSAIC_TYPES.
data : Contains data from readObservations() in the DATA key
indx : Row number within data
mtype : A list of mosaic types
The type of mosaics are currently MOSAIC_TYPE_CUSTOM and
MOSAIC_TYPE_RECTANGLE.
"""
ismosaic = False
if data[IS_MOSAIC][indx] in MOSAIC_TYPES:
ismosaic = data[IS_MOSAIC][indx] in makeList(mtype)
elif data[IS_MOSAIC][indx] == 'N/A' or np.isnan(data[IS_MOSAIC][indx]):
ismosaic = False
else:
raise Exception,'Unexpected value of IS_MOSAIC on row %d' % getExcelFromIndex(indx)
return ismosaic
def printSummaryHeader(label, spaces=''):
"""
Print header for summary table.
Inputs:
label : label to print in the header row
spaces : spaces used for formatting
"""
print ''
print 'Summary information for %s' % label
print '%s%6s %6s %-14s %-23s %12s %12s %8s %8s %8s %7s %7s %4s %4s %5s' % \
(spaces, 'N', 'Excel', 'Project code', 'Target name', 'RA', 'Dec', 'Sky Freq', 'Ang.Res.', 'L.A.S.',
'Polar-', 'MosArea', 'ACA?', 'TPA?', 'Spec.')
print '%s%6s %6s %-14s %-23s %12s %12s %8s %8s %8s %7s %7s %4s %4s %5s' % \
(spaces, '', 'row', '', '', 'J2000', 'J2000', '(GHz)', '(arcsec)', '(arcsec)', 'ization', 'amin^2', '', '', 'Scan?')
def printSummarySource(number, observations, indx, spaces=''):
"""
Print summary information for sources with indx of the observations.
Inputs:
number : Index to keep track of number of sources
observations : The observations data from readObservations()
indx : The index number if observations
Output:
A printed summary of the sources
"""
project = observations[DATA][PROJECT][indx]
target = observations[DATA][TARGET][indx]
scoords = observations[COORDS][indx].to_string(style='hmsdms', precision=1).split()
sra = scoords[0]
sdec = scoords[1]
sfreq = '%.1f' % observations[DATA][REF_FREQ][indx]
sangular = '%.2f' % observations[DATA][ANG_RES][indx]
slas = '%.1f' % observations[DATA][LAS][indx]
smosaic = '-'
saca = '-'
saca = '-'
stpa = '-'
sspectral = '-'
spol = observations[DATA][POLARIZATION][indx].lower()
if observations[DATA][IS_SPECTRAL][indx]:
sspectral = 'Yes'
if isObsMosaic(observations[DATA], indx):
if isObsMosaic(observations[DATA], indx, mtype=MOSAIC_TYPE_RECTANGLE):
smosaic = '%.1f' % observations[DATA][MOS_AREA][indx]
else:
smosaic = 'custom'
if observations[DATA][USE_ACA][indx] == True:
saca = 'Yes'
if observations[DATA][USE_TPA][indx] == True:
sata = 'Yes'
print '%s%6d %6d %-14s %-23s %12s %12s %8s %8s %8s %7s %7s %4s %4s %5s' % \
(spaces, number, getExcelFromIndex(indx), project, target, sra, sdec,
sfreq, sangular, slas, spol, smosaic, saca, stpa, sspectral)
def getMosaicCorners(ra_mosaic, dec_mosaic, width, height, angle, frame,
center=None, indx=None):
"""
Returns the corners of the mosaic in RA/Dec offsets from the mosaic
center. The function takes into account if the mosaic is specified
in galactic coordinates.
Inputs
angle : position angle of the rectangle in degrees.
center : The center of the plot in RA/Dec coordinates in degrees.
Enter as a tuple. e.g., center=[180.0, -10.0]
dec_mosaic : Declination center of the mosaic
frame : coordinate frame of the mosaic (J2000, icrs, or galactic)
heighto : height of the rectangle in arcseconds. Note this
corresponds to "width" in the spreadsheet.
ra_mosaic : RA center of the mosaic
widtho : width of the rectangle in arcseconds. Note this
corresponds to "height" in the spreadsheet.
Output:
A dictionary containing the RA/Dec offsets of the corners of the
mosaic relative to coords. If center, is none, then the mosaic
center is used. The dictionary keywords are UL, UR, BL, and BR.
"""
# Get the center of the mosaic in the native frame
if frame in ['J2000', 'icrs', 'ICRS']:
xcen_mosaic = ra_mosaic
ycen_mosaic = dec_mosaic
elif frame.lower() == 'galactic':
# Convert to ra/dec
c = SkyCoord(ra=ra_mosaic, dec=dec_mosaic, frame='icrs', unit='degree')
xcen_mosaic = c.galactic.l.deg
ycen_mosaic = c.galactic.b.deg
else:
msg = 'Unknown mosaic frame (%s)' % (frame)
if indx is not None:
msg += ' on row %d' % (getExcelFromIndex(indx))
raise Exception, msg
# Compute vertices relative to mosaic center
x = 0
y = 0
A = -angle / 180.0 * np.pi
results = dict()
results['UL'] = (x + ( width / 2. ) * np.cos(A) - ( height / 2. ) * np.sin(A), y + ( height / 2. ) * np.cos(A) + ( width / 2. ) * np.sin(A))
results['UR'] = (x - ( width / 2. ) * np.cos(A) - ( height / 2. ) * np.sin(A), y + ( height / 2. ) * np.cos(A) - ( width / 2. ) * np.sin(A))
results['BL'] = (x + ( width / 2. ) * np.cos(A) + ( height / 2. ) * np.sin(A), y - ( height / 2. ) * np.cos(A) + ( width / 2. ) * np.sin(A))
results['BR'] = (x - ( width / 2. ) * np.cos(A)+ ( height / 2. ) * np.sin(A), y - ( height / 2. ) * np.cos(A) - ( width / 2. ) * np.sin(A))
# Set plot center.
# The plot center is either specified in center, or if center=None,
# the plot center is the mosaic itself.
ra_center = ra_mosaic
dec_center = dec_mosaic
if center is not None:
ra_center = center[0]
dec_center = center[1]
# Set mosaic offsets in arcseconds relative to center
corners = dict()
for key, c in results.iteritems():
# Add offsets to central coordinates in native frame of the mosaic
xcorner = xcen_mosaic + c[0] / 3600.0 / np.cos(ycen_mosaic / 180.0 * np.pi)
ycorner = ycen_mosaic + c[1] / 3600.0
# Check limits on x
if xcorner < 0:
xcorner += 360.
if xcorner > 360:
xcorner -= 360.
# Check limits on y
if abs(ycorner) > 90.0:
msg = 'Error setting y-axis mosaic'
if indx is not None:
msg += ' on row %d' % getExcelFromIndex(indx)
raise Exception, msg
# Convert corner coordinates to ra/dec, if needed
if frame.lower() == 'galactic':
c = SkyCoord(xcorner, ycorner, frame=frame, unit='degree')
xcorner = c.icrs.ra.deg
ycorner = c.icrs.dec.deg
# Compute offsets relative to mosaic center
dra = (xcorner - ra_center)
ddec = (ycorner - dec_center) * 3600.0
if abs(dra) > 180.0:
# we are on either side of RA=0
if dra > 0:
dra -= 360.0
else:
dra += 360.0
dalp = dra * 3600.0 * np.cos(dec_center / 180. * np.pi)
# Save corner
corners[key] = (dalp, ddec)
# Done
return corners
def plotSources(coords, observations,
diameter=12,
include=None, exclude=None,
mosaic=False, width=60, length=60, pa=0.,
mframe='icrs', freq=345., mosonly=False,
plotroot='source', plotsizeo=120., plottype='pdf'):
"""
Plot observations that are nearby the specified coordinates.
Inputs:
coords : proposed coordinates from the user set by the
function getSourceCoordinates()
diameter : diameter of the telescope in meters
exclude : if set, it contains a list of spreadsheet row numbers
that should not be plotted.
freq : proposed observing frequency in GHz
include : if set, it contains a list of spreadsheet row
numbers that can be plotted if all other criteria
are set.
length : length of the mosaic in arcseconds
mframe : coordinate system of the mosaic (icrs or galactic)
mosaic : if True, the proposed observations are a mosaic
mosonly : if True, plot/print mosaic observations only
observations : existing observations from readObservations()
pa : position angle of the mosaic in degrees
plotroot : root name for the file containing the plot. The root
name will be appended with a plot number, which is
useful when plotting.
multiple sources) and the plottype.
plotsizeo : plot size in arcseconds
plottype : Type of plot to generate. The plot type must be
supported by your version of python. "pdf" and "png"
are usually safe. If plottype=None, then no plot is
saved.
Default is "pdf".
width : width of the mosaic in arcseconds
"""
# Set spacing for formatting purposes
spaces = ''
# Make a plot for each source
for i in range(coords[ORIGINAL].size):
# Initialize plot
py.figure(i+1, figsize=(9,9))
py.clf()
ax = py.subplot(1, 1, 1, aspect='equal')
# Set plot width
plotsize = coords[PLOTSIZE][i]
if plotsize is None:
plotsize = plotsizeo
# Find sources that overlap
if coords[NOCOORD][i]:
j = np.where(observations[DATA][TARGET].str.lower() == coords[ORIGINAL][i].lower())[0]
else:
# Compute separation in degrees
sep = coords[COORDS][i].separation(observations[COORDS])
sepArcsec = sep.arcsec
# Find entries within plot size
separationArcsec = sepArcsec - observations[DATA][MAX_SIZE]
j = np.where( separationArcsec <= (0.5*plotsize) )[0]
# Print summary of observation
if len(j) == 0:
sra = convertHmsString(float(coords[COORDS][i].ra.deg/15.0), ndec=1, delimiter='hms')
sdec = convertHmsString(float(coords[COORDS][i].dec.deg), ndec=1, delimiter='dms')
print ''
print 'Summary information for %s' % coords[ORIGINAL][i]
print ' No observations found within %g x %g arcsec region centered around (ra, dec) = (%s, %s) J2000' % \
(plotsize, plotsize, sra, sdec)
else:
printSummaryHeader('%g x %g arcsec region around %s' % \
(plotsize, plotsize, coords[ORIGINAL][i]), spaces=spaces)
# Set limits based on plotsize
ax.set_xlim( 0.5*plotsize, -0.5*plotsize)
ax.set_ylim(-0.5*plotsize, 0.5*plotsize)
ax.set_xlabel('$\\Delta\\alpha\ \ \\mathrm{(arcsec)}$')
ax.set_ylabel('$\\Delta\\delta\ \ \\mathrm{(arcsec)}$')
# Get row lists to include/exclude
rows_include = makeList(include)
rows_exclude = makeList(exclude)
# Set plot center.
# Sources will be plotted in offsets relative to this coordinate.
ra_center = coords[COORDS][i].ra.deg
dec_center = coords[COORDS][i].dec.deg
# Loop over observations which overlap the field
legend_symbols = []
legend_labels = []
number = 0
for indx in j:
# Get excel line
excelRow = getExcelFromIndex(indx)
if rows_include is not None and excelRow not in rows_include:
continue
if rows_exclude is not None and excelRow in rows_exclude:
continue
# If not mosaic and mosonly=True, then skip
if not isObsMosaic(observations[DATA], indx) and mosonly:
continue
# Get ra and dec offset from nominal position in arcseconds
if coords[NOCOORD][i]:
dalp = 0
ddec = 0
else:
# Compute offset
ddec = (observations[DATA][DEC_DEG][indx] - dec_center) * 3600.0
dra = (observations[DATA][RA_DEG][indx] - ra_center)
if abs(dra) > 180.0:
# we are on either side of RA=0
if dra > 0:
dra -= 360.0
else:
dra += 360.0
dalp = dra * 3600.0 * np.cos(dec_center / 180.0 * np.pi)
# Set center as offset coordinates
xy = (dalp, ddec)
# Print summary of observation
number += 1
printSummarySource(number, observations, indx, spaces=spaces)
# Set plot color based on band
color = getBandColor(observations[DATA][REF_FREQ][indx])
# Plot mosaic or single pointing
label = 'N = %d' % number
if isObsMosaic(observations[DATA], indx, mtype=MOSAIC_TYPE_RECTANGLE):
# Get the coordinates of the rectangle in RA/DEC offset units
mosaicCorners = getMosaicCorners(\
observations[DATA][RA_DEG][indx],
observations[DATA][DEC_DEG][indx],
observations[DATA][MOS_LENGTH][indx],
observations[DATA][MOS_WIDTH][indx],
observations[DATA][MOS_PA][indx],
observations[DATA][MOS_COORD][indx],
center=[ra_center, dec_center])
result = plotMosaic(ax, mosaicCorners, fc=color, ec=color , hatch=None, alpha=0.1)
# result = plotMosaicOld(ax, xy, observations[DATA][MOS_LENGTH][indx],
# observations[DATA][MOS_WIDTH][indx],
# observations[DATA][MOS_PA][indx],
# color='None', ec=color , hatch='/')
else:
result = plotPrimaryBeam(ax, xy, observations[DATA][REF_FREQ][indx], diameter,
fc='None', ec=color)
legend_symbols.append(result)
legend_labels.append(label)
# Loop over observations which overlap the field
# color = getBandColor(freq)
color = 'tan'
if mosaic:
corners = getMosaicCorners(ra_center, dec_center, length, width, pa, mframe)
result = plotMosaic(ax, corners, fc=color, ec='None', alpha=0.5, linewidth=3)
# result = plotMosaicOld(ax, (0,0), length, width, pa, color=color, ec='black', alpha=0.5, linewidth=3)
else:
result = plotPrimaryBeam(ax, (0,0), freq, diameter, fc=color, ec='None', alpha=0.5)
legend_symbols.append(result)
legend_labels.append('Proposed')
# Plot title with original entry and translation
sra = convertHmsString(float(coords[COORDS][i].ra.deg/15.0), ndec=1, delimiter='hms')
sdec = convertHmsString(float(coords[COORDS][i].dec.deg), ndec=1, delimiter='dms')
if coords[NOCOORD][i]:
labelTranslated = ''
else:
labelTranslated = '%s, %s' % (sra, sdec)
labelOriginal = '%s' % coords[ORIGINAL][i]
fs = 14
# ax.set_title(labelOriginal, loc='left', fontsize=fs)
# ax.set_title(labelTranslated, loc='right', fontsize=fs)
ax.annotate(s='Entered:', xy=(0,1.05), xycoords='axes fraction', size=fs)
ax.annotate(s=labelOriginal, xy=(0.2,1.05), xycoords='axes fraction', size=fs)
ax.annotate(s='Translated:', xy=(0,1.01), xycoords='axes fraction', size=fs)
ax.annotate(s=labelTranslated, xy=(0.2,1.01), xycoords='axes fraction', size=fs)
# py.legend(legend_symbols, legend_labels)
# Save plot to a file
py.show()
if plottype is not None:
# Set name
root = plotroot
if root is None:
root = 'source'
# Try creating plot
try:
plotfile = '%s%d.%s' % (root, i+1, plottype)
py.savefig(plotfile)
print ' Plot saved to %s' % plotfile
except:
print ' Warning: Could not create plot %s.' % plotfile
print ' Is that plot type supported by your python extension?'
def isMovingObject(data, indx):
"""
Returns True/False that object is a moving object.
"""
return (data[RA_DEG][indx]==0 and data[DEC_DEG][index]==0)
def getCoordinatesFromName(name, data=None, lookup=True):
"""
Try to get coordinates for the source name.
First, it tries using Sesame. If that fails, then we try using
finding the source in the observations spreadsheet.
name : Name of the source
data : Contains data from readObservations() in the DATA key
"""
# Try resolving source name
found = False
if lookup:
try:
result = astropy.coordinates.name_resolve.get_icrs_coordinates(name)
found = True
except:
pass
if not found:
# Initialize
result = None
# Is source in observation list?
if data is not None:
j = np.where(data[TARGET].str.lower() == name.lower())[0]
if len(j) > 0 and not isMovingObject(data, j[0]):
result = SkyCoord(data[RA_DEG][j[0]],
data[DEC_DEG][j[0]], unit='deg', frame='icrs')
# Done
return result
def getSourceCoordinates(longitude, latitude, sources, inputFile, rows,
unit='deg', frame='icrs', data=None,
spreadsheet=None, lookup=True):
"""
Get list of coordinates to search.
longitude : longitude (equatorial or galactic coordinates. Examples:
longitude=180.0
longitude='180.0d'
longitude='15h00m00s'
longitude=['10h12m13s', '180.0d', 180.0]
latitude : longitude (equatorial or galactic coordinates. Examples:
latitude=-4.0
latitude='-4.0d'
latitude=['-21d11m12s', '-4d', -4.0]
sources : List of source names. Examples:
sources='DG Tau'
sources='DG Tau, HL Tau, DM Tau'
sources=['DG Tau', 'HL Tau', 'DM Tau']
inputFile : Name of text file containing sources to search. In addition
to source names, the file may contain commands to change
plot sizes. The pound symbol (#) is used as a comment marker.
Example input for the data file.
# Specify sources to sea
HL Tau
DG Tau
# Change plot size from default for remaining sources
plotsize = 480
# Additional source name
COSMOS-16199
# Change coordinate frames and plot size
frame = galactic
plotsize = 60
# Search by coordinate
178.8590 -19.9724
rows : List of row numbers in excel spreadsheet
rows=3502
rows=[3502, 3510]
unit : Unit for longitude/latitude if they are entered as floating
points. Most likely values are
a) unit='deg' , which means both longitude and latitude
are in degrees.
b) unit='hour,deg', which means longitude is in
hours and latitude is in degrees.
frame : Coordinate frame for longitude/latitude. Accepted values
are 'icrs' for equatorial or 'galactic'.
data : Contains data from readObservations() in the DATA key.
This is used to search for individual sources if
a source name is not name resolved.
lookup : If True, resolve the source name using Sesame.
lookup=False can be useful if the name listed in the
spreadsheet is desired.
spreadsheet: Name of the spreadsheet containing the observations.
This parameter is optional, and only used to
customize the output.
There are four option to search:
1) Enter list of numpy array of RA/DEC.
The format is flexible in that RA/Dec be in string, hex, or
decimal format. If the units are not specified, the default
units are given by the 'unit' variable.
For example:
longitude=['01h05m11s', 1.32]
latitude=['-10d05m11s', -20.0]
2) Enter of a list of source names that will be name resolved
3) Enter a file with coordinates/names, one entry per line
The comment symbol per line is the pound symbol.
The entry will name resolved first, and if that is not successful,
it will be coordinate resolved. If that is not successful, then
the source is skipped.
4) Enter list of row numbers
"""
# Initialize ra/dec list
ra = []
dec = []
original = []
iscoord = []
nocoord = []
plotsize = []
# Add coordinates
def addCoordinates(c, psize=None):
# Save in degrees
if c is None:
ra.append(0)
dec.append(0)
nocoord.append(True)
plotsize.append(psize)
else:
if type(c.icrs.ra.deg) in [list, np.ndarray]:
zra = np.array(c.icrs.ra.deg)
else:
zra = np.array([c.icrs.ra.deg])
if type(c.icrs.dec.deg) in [list, np.ndarray]:
zdec = np.array(c.icrs.dec.deg)
else:
zdec = np.array([c.icrs.dec.deg])
ra.extend(zra)
dec.extend(zdec)
nocoord.extend([False] * zra.size)
plotsize.extend([psize] * zra.size)
# Process RA/Dec
if longitude is not None or latitude is not None:
# Both must be set
if longitude is None or latitude is None:
raise Exception, 'Both longitude/latitude must be set if you are entering coordinates'
# Convert entries into list so we can loop over them
llon = makeList(longitude)
llat = makeList(latitude)
# Must have same size
if llon.size != llat.size:
raise Exception, 'longitude/latitude must have the same size'
if llon.ndim != 1:
raise Exception, 'Only 1D arrays can be processed for longitude'
if llat.ndim != 1:
raise Exception, 'Only 1D arrays can be processed for latitude'
# Get coordinates
c = SkyCoord(llon, llat, unit=unit, frame=frame)
# Add coordinates
addCoordinates(c)
# Save original entries
iscoord.extend([True] * c.icrs.ra.size)
s = []
for t in zip(llon, llat):
s.append('%s, %s %s' % (t[0], t[1], frame))
original.extend(s)
# Process row numbers
if rows is not None:
# Convert to list
lrows = getIndexFromExcel(makeList(rows))
# Get coordinates
c = SkyCoord(data[RA_DEG][lrows], data[DEC_DEG][lrows], unit='deg', frame='icrs')
# Add coordinates
addCoordinates(c)
# Save original entries
iscoord.extend([True] * c.icrs.ra.size)
s = []
for t in lrows:
ss = 'Row %d (%s)' % (getExcelFromIndex(t), str(data[TARGET][t]))
if spreadsheet is not None:
ss += ' in %s' % spreadsheet
s.append(ss)
original.extend(s)
# Process source names
if sources is not None:
# Convert into list
names = makeList(sources)
# Loop over names
for name in names:
# Get coordinates
c = getCoordinatesFromName(name, data=data, lookup=lookup)
# Add to coordinate list
if c is None:
print 'Warning: Could not find coordinates for %s. Searching by object name.' % name
addCoordinates(c)
iscoord.append(False)
original.append(name)
# Process input file
psize_new = None
unit_new = unit
frame_new = frame
if inputFile is not None:
# Read file one line at a time
finp = open(inputFile, 'r')
# Process one line at a time
nlines = 0
for line in finp:
# Strip leading and trailing spaces, as well as multiple spaces
line = re.sub(' +',' ', line.strip())
# Remove comment portion of line
j = line.find('#')
if j >= 0: line = line[0:j]
# Skip blank lines
if line.strip() == '':
continue
# Is this the plot size?
if line.lower().find('plotsize') == 0:
value = line.split('=')[1]
try:
psize_new = float(value)
except:
pass
continue
# Is this the unit?
if line.lower().find('unit') == 0:
value = line.split('=')[1]
try:
unit_new = value.strip()
except:
pass
continue
# Is this the frame?
if line.lower().find('frame') == 0:
value = line.split('=')[1]
try:
frame_new = value.strip()
except:
pass
continue
# First, try if parse line as coordinates
try:
nlines += 1
c = SkyCoord(line, unit=unit_new, frame=frame_new)
addCoordinates(c, psize=psize_new)
iscoord.extend([True] * c.icrs.ra.size)
original.extend(['%s %s' % (line, frame_new)])
except:
# Coordinates was not successful. Try resolving it as a name
c = getCoordinatesFromName(line, data=data)
if c is not None:
nlines += 1
addCoordinates(c, psize=psize_new)
iscoord.append(False)
original.append(line)
else:
print 'Warning: Could not find coordinates for %s. Searching by object name.' % line
addCoordinates(None, psize=psize_new)
iscoord.append(False)
original.append(line)
# Message
print 'Read in %d sources from %s' % (nlines, inputFile)
# Prepare results
results = dict()
results[COORDS] = SkyCoord(ra, dec, unit='deg')
results[RA] = makeList(ra)
results[DEC] = makeList(dec)
results[ORIGINAL] = makeList(original)
results[ISCOORD] = makeList(iscoord)
results[NOCOORD] = makeList(nocoord)
results[PLOTSIZE] = makeList(plotsize)
# Check that all elements have the same size
keys = results.keys()
for i, k in enumerate(keys):
if len(results[k]) != len(results[keys[0]]):
raise Exception,'Arrays do not have the same size: %s and %s' % (keys[0], k)
# Done
return results
def row(rows, showProject=True, showCorrelator=True,
inputObs=LIST_OF_OBSERVATIONS, verbose=True, refresh=False):
"""
Print detailed information about the project and correlator for a
a line in the spreadsheet.
Inputs:
rows : The row numbers in the excel spreadsheet that
should be printed. An integer or a list can be
entered.
showProject : If True, then print the correlator information for
each line
showCorrelator: If True, then print the project information for
each line
inputObs : The name of the excel spreadsheet containing the
observations.
verbose : If True, print out messages while reading spreadsheet
refresh : If True, re-read the spreadsheet
Output:
A text listing of the project and correlator information.
Example usage:
import plotobs as p4
p4.showRows(3329)
p4.showRows([3329, 3330])
"""
# Read list of previous or scheduled observations
global OBSERVATIONS
if OBSERVATIONS is None or refresh:
OBSERVATIONS = readObservations(inputObs, verbose=verbose)
observations = OBSERVATIONS
# Loop over input excel sheet spread lines
for excel in makeList(rows):
# Convert to index number if observation list
indx = getIndexFromExcel(excel)
# Print
if showProject: printRowProject(observations[DATA], indx)
if showCorrelator: printRowCorrelator(observations[DATA], indx)
def project(projects, inputObs=LIST_OF_OBSERVATIONS, verbose=True, refresh=False):
"""
Print information about a project in the spreadsheet.
Inputs:
projects : A string or list containing the projects to list.
inputObs : The name of the excel spreadsheet or csv file
containing the observations.
verbose : If True, print out messages while reading spreadsheet
refresh : If True, re-read the spreadsheet
Output:
A text listing of the project and correlator information.
Example usage:
import plotobs as p4
p4.project('2015.1.01289.S')
p4.project(['2015.1.01289.S', '2013.1.00385.S'])
"""
# Read list of previous or scheduled observations
global OBSERVATIONS
if OBSERVATIONS is None or refresh:
OBSERVATIONS = readObservations(inputObs, verbose=verbose)
observations = OBSERVATIONS
# Get projects
obsProjects = np.char.lower(observations[DATA][PROJECT].values.tolist())
# Loop over projects
for p in makeList(projects):
# Print header
printSummaryHeader(p)
# Find indices
indices = np.where(obsProjects == p.lower())[0]
# Print indices
if len(indices) == 0:
print ' No observations found for project %s' % p
else:
for n, indx in enumerate(indices):
printSummarySource(n+1, observations, indx)
def test(projects=True, rows=True, plots=True, refresh=True,
inputObs=LIST_OF_OBSERVATIONS, verbose=True):
"""
Run through all sources and through the projects to test the scripts
to make all sources can be processed without run-time errors.
projects: If True, process all projects
rows : If True, process all rows
plots : Run the plot() command on each row (very slow)
refresh : Reload spreadsheet instead of using data in memory
"""
# Load data
global OBSERVATIONS
if OBSERVATIONS is None or refresh:
OBSERVATIONS = readObservations(inputObs, verbose=verbose)
observations = OBSERVATIONS
# Check projects
if projects:
# Get unique projects
projectList = np.unique(observations[DATA][PROJECT])
# Run each project
project(projectList)
# Check rows
if rows:
for indx in np.arange(observations[DATA].shape[0]):
row(getExcelFromIndex(indx))
# Plots
if plots:
for indx in np.arange(observations[DATA].shape[0]):
plot(rows=getExcelFromIndex(indx))
def plot(longitude=None, latitude=None, frame='icrs', unit='deg',
sources=None, lookup=True, inputFile=None,
rows=None, include=None, exclude=None,
freq=345.0, width=30., length=60, pa=0.,
mframe='icrs', mosaic=False,
refresh=False, verbose=True, inputObs=LIST_OF_OBSERVATIONS,
mosonly=False,
plotroot='source', plotsize=120., plottype="pdf"):
"""
Purpose:
plot() lists and displays existing observations from Cycle 1, 2, and
3 and any scheduled Cycle 3 Grade A observations.
Each observation is plotted as a single pointing with a circle the
size of the primary beam, or a rectangle for rectangular mosaics. For
custom mosaics, each individual pointing is plotted.
One can specify an observing frequency and mosaic parameters
(optional), which will also be plotted with a filled circle or
rectangle.
Usage:
1) The user enters coordinates or source names, the observed
frequency, the proposed mosaic parameters, and the size of the
plot region (PLOTSIZE)
2) plot() will list all observation within PLOTSIZE size, as well
as plot the observations.
How to input source name or coordinates:
Parameters:
exclude : If set, it contains a list of spreadsheet row numbers
that should not be plotted.
frame : reference frame for longitude, either "icrs"
(i.e., equatorial) or "galactic"
include : If set, it contains a list of spreadsheet row numbers
that can be plotted if all other criteria are set
inputFile : Input list sources or positions
inputObs : Input excel file containing previous observations
and scheduled observations
latitude : Latitude coordinate with reference frame specified
by "frame". Examples:
latitude=-20.0, unit='deg'
latitude='-20 deg'
latitude=['-20 deg', '20d00m00s']
longitude : Longitude coordinate with reference frame specified
by "frame". Examples:
longitude=100.0, unit='deg'
longitude='100 deg'
longitude=['100 deg', '17h15m:10s']
lookup : If True, resolve the source name using Sesame.
lookup=False can be useful if the name listed in the
spreadsheet is desired.
mosonly : If True, print/display only mosaic observations
refresh : If True, re-read excel spreadsheet into memory
rows : List of row numbers in excel spreadsheet
sources : List of source names. The coordinates of the source
names are obtain from the Sesame database, and if the
name is not resolved, the name is searched
(case-insensitive) in the input spreadsheet. If
lookup=False, the search in the Sesame database is
skipped.
unit : units of longitude and latitude if not specified by
in longitude/latitude. Examples:
unit='deg' implies longitude/latitude are in degrees.
unit='hour,deg' implies longitude is in hours
and latitude is in degrees.
verbose : If True, print out messages when reading the excel
spreadsheet
Plot parameters
plotroot : Root name for the file containing the plot. The root
name will be appended with a plot number (useful when
plotting multiple sources) and the plottype.
plotsize : Size of the square region to plot in arcseconds
plottype : Type of plot to generate. The plot type must be
supported by your version of python. "pdf" and "png"
are usually safe. Default is "pdf". Set plottype=None
to not save the plots to a file.
Proposed observed parameters
freq : Proposed observing frequency in GHz
length : Proposed length of the mosaic in arcseconds
mframe : Coordinate system of the mosaic (icrs or galactic)
mosaic : If True, the proposed observations are a mosaic
pa : Proposed position angle of the mosaic in degrees
width : Proposed width of the mosaic in arcseconds
Usage:
There are three options to search for source names:
1) Enter list of coordinates.
The format is flexible in that coordinates may be a string,
hex, or decimal format. If the units are not specified, the
default units are given by the 'unit' variable.
For example:
longitude=['01h05m11s', 1.32]
latitude=['-10d05m11s', -20.0]
2) Enter of a list of source names that will be name resolved
by querying the Sesame database. An internet connection
must be available for the coordinate lookup
3) Enter a file with coordinates/names, one entry per line
The comment symbol per line is the pound symbol.
The entry will be name resolved first, and if that is not
successful, it will be coordinate resolved. If that is not
successful, then the source is skipped.
4) The source name will be matched to the OBSERVATION list.
The search will be case insensitive.
Examples:
# Plot observations within 480" x 480" of M83 along with a proposed
# pointing at the position of the galaxy
p4.plot('M83', plotsize=480)
# Plot observations within 480" x 480" of M83, along with a
# proposed mosaic at 345 GHz
p4.plot('M83', plotsize=480, length=240, width=120,
pa=60, freq=345., mosaic=True)
# Plot observations within 480" of (ra,dec) = (17h20m42s,-35d48m04s)
p4.plot('17h20m42s', '-35d48m04s', plotsize=280,
frame='icrs', unit='deg')
# Plot 1 deg region around ra, dec = (150, 2.2) J2000
p4.plot(150, 2.2, plotsize=3600)
# Plot 1 deg around galactic center
p4.plot(0, 0, frame='galactic', plotsize=3600., mosonly=True)
# Plot observations within a 60" x 60" of HL Tau and DG Tau
p4.plot('HL Tau, DG Tau', plotsize=60)
# Plot observations near coordinates in row 272 of the spreadsheet
p4.plot(272)
# Plot observations only for row 272 in the spreadsheet,
# but plot only row 272
p4.plot(272, include=272)
# Plot observations for sources listed in a text file
p4.plot(inputFile='input.txt')
# Contents of input file "input.txt":
HL Tau
DG Tau
plotsize = 480
COSMOS-16199
frame = galactic
plotsize = 60
178.8590 -19.9724
"""
# If longitude is set but no other source parameters are set, then
# assume the following:
# 1) If longitude is a string, then assume it is a source name
# 2) If it is an integer, then it is a row number
if longitude is not None and latitude is None and \
rows is None and sources is None and inputFile is None:
if str in [type(longitude)]:
sources = longitude
longitude = None
elif int in [type(longitude)]:
rows = longitude
longitude = None
# Check inputs
error = False
# Longitude/latitude must both be set or not at all
if (longitude is not None and latitude is None) or \
(longitude is None and latitude is not None):
error = True
print 'Error: Longitude and latitude must both be set to enter coordinates'
# Some coordinates must be set
if (longitude is None or latitude is None) and \
rows is None and \
inputFile is None and \
sources is None:
error = True
print 'Error: longitude/latitude, sources, rows, or inputFile must be set'
# Check frame of coordinates
if frame.lower() not in ['icrs', 'galactic']:
error = True
print 'Error: frame must be set to "icrs" or "galactic"'
# Check plotsize
if plotsize <= 0:
error = True
print 'Error: plotsize must be > 0'
# Check frequency
if freq <= 0:
error = True
print 'Error: frequency must be > 0'
# Check mosaic parameters, if mosaic is set
if mosaic:
if length <= 0:
error = True
print 'Error: length must be > 0'
if width <= 0:
error = True
print 'Error: width must be > 0'
if mframe.lower() not in ['icrs', 'galactic']:
error = True
print 'Error: mframe must be set to "icrs" or "galactic"'
if error:
print ''
print 'Error starting plotobs_cycle4'
return
# Read list of previous or scheduled observations, if not entered
# on command line
global OBSERVATIONS
if OBSERVATIONS is None or refresh:
OBSERVATIONS = readObservations(input=inputObs, verbose=verbose)
observations = OBSERVATIONS
# Get coordinates of source(s)
coords = getSourceCoordinates(longitude, latitude, sources, inputFile, rows,
unit=unit, frame=frame, data=observations[DATA],
spreadsheet=inputObs, lookup=lookup)
# Plot results
plotSources(coords, observations,
include=include, exclude=exclude,
freq=freq, mosonly=mosonly,
mosaic=mosaic, mframe=mframe, width=width, length=length, pa=pa,
plotroot=plotroot, plottype=plottype, plotsizeo=plotsize)