""" 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)