Coverage for /wheeldirectory/casa-6.7.0-12-py3.10.el8/lib/py/lib/python3.10/site-packages/casatasks/private/task_plotants.py: 85%
201 statements
« prev ^ index » next coverage.py v7.6.4, created at 2024-11-01 07:19 +0000
2import os
3import sys
4import numpy as np
5import pylab as pl
6from textwrap import wrap
8from casatools import table, msmetadata, quanta, ms, measures, ctsys
9from casatasks import casalog
11def plotants( vis=None, figfile=None,
12 antindex=None, logpos=None,
13 exclude=None, checkbaselines=None,
14 title=None, showgui=None ):
15 """Plot the antenna distribution in the local reference frame:
16 The location of the antennas in the MS will be plotted with
17 X-toward local east; Y-toward local north. The name of each
18 antenna is shown next to its respective location.
20 Keyword arguments:
21 vis -- Name of input visibility file.
22 default: none. example: vis='ngc5921.ms'
24 figfile -- Save the plotted figure in this file.
25 default: ''. example: figfile='myFigure.png'
27 antindex -- Label antennas with name and antenna ID
28 default: False. example: antindex=True
30 logpos -- Produce a logarithmic position plot.
31 default: False. example: logpos=True
33 exclude -- antenna IDs or names to exclude from plotting
34 default: []. example: exclude=[2,3,4], exclude='DV15'
36 checkbaselines -- Only plot antennas in the MAIN table.
37 This can be useful after a split. WARNING: Setting
38 checkbaselines to True will add to runtime in
39 proportion to the number of rows in the dataset.
40 default: False. example: checkbaselines=True
42 title -- Title written along top of plot
43 default: ''
45 You can zoom in by pressing the magnifier button (bottom,
46 third from right) and making a rectangular region with
47 the mouse. Press the home button (left most button) to
48 remove zoom.
50 A hard-copy of this plot can be obtained by pressing the
51 button on the right at the bottom of the display. A file
52 dialog will allow you to choose the directory, filename,
53 and format of the export.
54 """
56 # use ctsys values
57 showplot = showgui and not ctsys.getnogui() and not ctsys.getagg() and not ctsys.getpipeline()
59 if not showplot:
60 pl.close()
61 pl.ioff()
62 else:
63 pl.show()
64 pl.ion()
65 pl.clf()
67 # remove trailing / for title basename
68 if vis[-1]=='/':
69 vis = vis[:-1]
70 myms = ms( )
71 try:
72 exclude = myms.msseltoindex(vis, baseline=exclude)['antenna1'].tolist()
73 except RuntimeError as rterr: # MSSelection failed
74 errmsg = str(rterr)
75 errmsg = errmsg.replace('specificion', 'specification')
76 errmsg = errmsg.replace('Antenna Expression: ', '')
77 raise RuntimeError("Exclude selection error: " + errmsg)
79 telescope, names, ids, xpos, ypos, stations = getPlotantsAntennaInfo(vis,
80 logpos, exclude, checkbaselines)
81 if not names:
82 raise ValueError("No antennas selected. Exiting plotants.")
84 if not title:
85 msname = os.path.basename(vis)
86 title = "Antenna Positions for "
87 if len(msname) > 55:
88 title += '\n'
89 title += msname
91 if logpos:
92 plotAntennasLog(telescope, names, ids, xpos, ypos, antindex, stations)
93 else:
94 plotAntennas(telescope, names, ids, xpos, ypos, antindex, stations, showplot)
95 pl.title(title, {'fontsize':12})
97 if figfile:
98 pl.savefig(figfile)
100def getPlotantsAntennaInfo(msname, log, exclude, checkbaselines):
102 tb = table( )
103 me = measures( )
104 qa = quanta( )
106 telescope, arrayPos = getPlotantsObservatoryInfo(msname)
107 arrayWgs84 = me.measure(arrayPos, 'WGS84')
108 arrayLon, arrayLat, arrayAlt = [arrayWgs84[i]['value']
109 for i in ['m0','m1','m2']]
111 # Open the ANTENNA subtable to get the names of the antennas in this MS and
112 # their positions. Note that the entries in the ANTENNA subtable are pretty
113 # much in random order, so antNames translates between their index and name
114 # (e.g., index 11 = STD155). We'll need these indices for later, since the
115 # main data table refers to the antennas by their indices, not names.
117 anttabname = msname + '/ANTENNA'
118 tb.open(anttabname)
119 # Get antenna names from antenna table
120 antNames = np.array(tb.getcol("NAME")).tolist()
121 stationNames = np.array(tb.getcol("STATION")).tolist()
122 if telescope == 'VLBA': # names = ant@station
123 antNames = ['@'.join(antsta) for antsta in zip(antNames,stationNames)]
124 # Get antenna positions from antenna table
125 antPositions = np.array([me.position('ITRF', qa.quantity(x, 'm'),
126 qa.quantity(y, 'm'), qa.quantity(z, 'm'))
127 for (x, y, z) in tb.getcol('POSITION').transpose()])
128 tb.close()
130 allAntIds = list(range(len(antNames)))
131 if checkbaselines:
132 # Get antenna ids from main table; this will add to runtime
133 tb.open(msname)
134 ants1 = tb.getcol('ANTENNA1')
135 ants2 = tb.getcol('ANTENNA2')
136 tb.close()
137 antIdsUsed = list(set(np.append(ants1, ants2)))
138 else:
139 # use them all!
140 antIdsUsed = allAntIds
142 # handle exclude -- remove from antIdsUsed
143 for antId in exclude:
144 try:
145 antNameId = antNames[antId] + " (id " + str(antId) + ")"
146 antIdsUsed.remove(antId)
147 casalog.post("Exclude antenna " + antNameId)
148 except ValueError:
149 casalog.post("Cannot exclude antenna " + antNameId + ": not in main table", "WARN")
151 # apply antIdsUsed mask
152 antNames = [antNames[i] for i in antIdsUsed]
153 antPositions = [antPositions[i] for i in antIdsUsed]
154 stationNames = [stationNames[i] for i in antIdsUsed]
156 nAnts = len(antIdsUsed)
157 casalog.post("Number of points being plotted: " + str(nAnts))
158 if nAnts == 0: # excluded all antennas
159 return telescope, antNames, [], [], []
161 # Get the names, indices, and lat/lon/alt coords of "good" antennas.
162 antWgs84s = np.array([me.measure(pos, 'WGS84') for pos in antPositions])
164 # Convert from lat, lon, alt to X, Y, Z (unless VLBA)
165 # where X is east, Y is north, Z is up,
166 # and 0, 0, 0 is the center
167 # Note: this conversion is NOT exact, since it doesn't take into account
168 # Earth's ellipticity! But it's close enough.
169 if telescope == 'VLBA' and not log:
170 antLons, antLats = [[pos[i] for pos in antWgs84s] for i in ['m0','m1']]
171 antXs = [qa.convert(lon, 'deg')['value'] for lon in antLons]
172 antYs = [qa.convert(lat, 'deg')['value'] for lat in antLats]
173 else:
174 antLons, antLats = [np.array( [pos[i]['value']
175 for pos in antWgs84s]) for i in ['m0','m1']]
176 radE = 6370000.
177 antXs = (antLons - arrayLon) * radE * np.cos(arrayLat)
178 antYs = (antLats - arrayLat) * radE
179 return telescope, antNames, antIdsUsed, antXs, antYs, stationNames
181def getPlotantsObservatoryInfo(msname):
182 metadata = msmetadata()
183 metadata.open(msname)
184 telescope = metadata.observatorynames()[0]
185 arrayPos = metadata.observatoryposition()
186 metadata.close()
187 return telescope, arrayPos
189def getAntennaLabelProps(telescope, station, log=False):
190 # CAS-7120 make plots more readable (rotate labels)
191 vAlign = 'center'
192 hAlign = 'left'
193 rotAngle = 0
194 if station and "VLA" in telescope:
195 # these have non-standard format:
196 # strip off VLA: or VLA:_ prefix if any
197 if 'VLA:' in station:
198 station = station[4:]
199 if station[0] == '_':
200 station = station[1:]
201 if station in ['W01', 'E01', 'W1', 'E1']:
202 vAlign = 'top'
203 hAlign = 'center'
204 else:
205 vAlign = 'bottom'
206 if 'W' in station or 'MAS' in station:
207 hAlign = 'right'
208 rotAngle = -35
209 elif 'E' in station or ('N' in station and not log):
210 rotAngle = 35
211 return vAlign, hAlign, rotAngle
213def plotAntennasLog(telescope, names, ids, xpos, ypos, antindex, stations):
214 fig = pl.figure(1)
215 # Set up subplot.
216 ax = fig.add_subplot(1, 1, 1, polar=True, projection='polar')
217 ax.set_theta_zero_location('N')
218 ax.set_theta_direction(-1)
219 # Do not show azimuth labels.
220 ax.set_xticklabels([])
221 ax.set_yticklabels([])
222 # Do not show grid.
223 ax.grid(False)
225 # code from pipeline summary.py
226 # PlotAntsChart draw_polarlog_ant_map_in_subplot
227 if 'VLA' in telescope:
228 # For (E)VLA, set a fixed local center position that has been
229 # tuned to work well for its array configurations (CAS-7479).
230 xcenter, ycenter = -32, 0
231 rmin_min, rmin_max = 12.5, 350
232 else:
233 # For non-(E)VLA, take the median of antenna offsets as the
234 # center for the plot.
235 xcenter = np.median(xpos)
236 ycenter = np.median(ypos)
237 rmin_min, rmin_max = 3, 350
239 # Derive radial offset w.r.t. center position.
240 r = ((xpos-xcenter)**2 + (ypos-ycenter)**2)**0.5
241 # Set rmin, clamp between a min and max value, ignore station
242 # at r=0 if one is there.
243 rmin = min(rmin_max, max(rmin_min, 0.8*np.min(r[r > 0])))
244 # Update r to move any points below rmin to r=rmin.
245 r[r <= rmin] = rmin
246 rmin = np.log(rmin)
247 # Set rmax.
248 rmax = np.log(1.5*np.max(r))
249 # Derive angle of offset w.r.t. center position.
250 theta = np.arctan2(xpos-xcenter, ypos-ycenter)
252 # Draw circles at specific distances from the center.
253 angles = np.arange(0, 2.01*np.pi, 0.01*np.pi)
254 show_circle = True
255 circles = [30, 100, 300, 1000, 3000, 10000]
256 if telescope == "VLBA":
257 circles = [1e5, 3e5, 1e6, 3e6, 1e7]
258 for cr in circles:
260 # Only draw circles outside rmin.
261 if cr > np.min(r) and show_circle:
263 # Draw the circle.
264 radius = np.ones(len(angles))*np.log(cr)
265 ax.plot(angles, radius, 'k:')
267 # Draw tick marks on the circle at 1 km intervals.
268 inc = 0.1*10000/cr
269 if telescope == "VLBA":
270 inc = 0.1*1e8/cr
271 if cr > 100:
272 for angle in np.arange(inc/2., 2*np.pi+0.05, inc):
273 ax.plot([angle, angle],
274 [np.log(0.95*cr), np.log(1.05*cr)], 'k-')
276 # Add text label to circle to denote distance from center.
277 va = 'top'
278 circle_label_angle = -20.0 * np.pi / 180.
279 if cr >= 1000:
280 if np.log(cr) < rmax:
281 ax.text(circle_label_angle, np.log(cr),
282 '%d km' % (cr/1000), size=8, va=va)
283 ax.text(circle_label_angle + np.pi, np.log(cr),
284 '%d km' % (cr / 1000), size=8, va=va)
285 else:
286 ax.text(circle_label_angle, np.log(cr), '%dm' % (cr),
287 size=8, va=va)
288 ax.text(circle_label_angle + np.pi, np.log(cr),
289 '%dm' % (cr), size=8, va=va)
291 # Find out if most recently drawn circle was outside all antennas,
292 # if so, no more circles will be drawn.
293 if np.log(cr) > rmax:
294 show_circle = False
296 # plot points and antenna names/ids
297 for i, (name, station) in enumerate(zip(names, stations)):
298 if station and 'OUT' not in station:
299 ax.plot(theta[i], np.log(r[i]), 'ko', ms=5, mfc='r')
300 if antindex:
301 name += ' (' + str(ids[i]) + ')'
302 # set alignment and rotation angle (for VLA)
303 valign, halign, angle = getAntennaLabelProps(telescope, station, log=True)
304 # adjust so text is not on the circle:
305 yoffset = 0
306 if halign == 'center':
307 yoffset = 0.1
308 ax.text(theta[i], np.log(r[i])+yoffset, ' '+name, size=8, va=valign, ha=halign,
309 rotation=angle, weight='semibold')
311 # Set minimum and maximum radius.
312 ax.set_rmax(rmax)
313 ax.set_rmin(rmin)
315 # Make room for 2-line title
316 pl.subplots_adjust(top=0.88)
318def plotAntennas(telescope, names, ids, xpos, ypos, antindex, stations, showplot):
319 fig = pl.figure(1)
320 ax = fig.add_subplot(111)
322 if telescope == 'VLBA':
323 labelx = 'Longitude (deg)'
324 labely = 'Latitude (deg)'
325 else:
326 # use m or km units
327 units = ' (m)'
328 if np.median(xpos) > 1e6 or np.median(ypos) > 1e6:
329 xpos /= 1e3
330 ypos /= 1e3
331 units = ' (km)'
332 labelx = 'X' + units
333 labely = 'Y' + units
334 if "VLA" in telescope:
335 spacer = ' '
336 else:
337 spacer = ' '
339 # plot points and antenna names/ids
340 for i, (x, y, name, station) in enumerate(zip(xpos, ypos, names, stations)):
341 if station and 'OUT' not in station:
342 ax.plot(x, y, 'ro')
343 if antindex:
344 name += ' (' + str(ids[i]) + ')'
345 # set alignment and rotation angle (for VLA)
346 valign, halign, angle = getAntennaLabelProps(telescope, station)
347 # adjust so text is not on the circle:
348 if halign == 'center':
349 y -= 10
350 ax.text(x, y, ' '+name, size=8, va=valign, ha=halign, rotation=angle,
351 weight='semibold')
352 if showplot:
353 fig.show()
355 pl.xlabel(labelx)
356 pl.ylabel(labely)
357 pl.margins(0.1, 0.1)