Coverage for /wheeldirectory/casa-6.7.0-12-py3.10.el8/lib/py/lib/python3.10/site-packages/casatasks/private/task_plotprofilemap.py: 81%
601 statements
« prev ^ index » next coverage.py v7.6.4, created at 2024-11-01 07:19 +0000
« prev ^ index » next coverage.py v7.6.4, created at 2024-11-01 07:19 +0000
1import os
3import matplotlib.pyplot as plt
4import numpy as np
6from casatasks import casalog
7from casatools import image, quanta
9from . import sdutil
11qa = quanta()
14@sdutil.sdtask_decorator
15def plotprofilemap(imagename=None, figfile=None, overwrite=None, transparent=None,
16 pol=None, spectralaxis=None, restfreq=None, plotrange=None, title=None,
17 linecolor=None, linestyle=None, linewidth=None,
18 separatepanel=None, plotmasked=None, maskedcolor=None,
19 showaxislabel=None, showtick=None, showticklabel=None,
20 figsize=None, numpanels=None):
22 if len(figfile) > 0 and os.path.exists(figfile) and not overwrite:
23 raise RuntimeError('overwrite is False and output file exists: \'%s\''%(figfile))
25 image_data = SpectralImage(imagename)
26 npol = len(image_data.stokes)
27 if pol < 0 or pol > npol - 1:
28 raise RuntimeError('pol {pol} is out of range (Stokes axis {stokes})'.format(pol=pol,stokes=image_data.stokes))
30 parsed_size = parse_figsize(figsize)
31 nx, ny = parse_numpanels(numpanels)
32 if (not isinstance(restfreq, str)) or len(restfreq) == 0:
33 restfreq = None
34 plot_profile_map(image_data, figfile, pol, spectralaxis, restfreq, title, linecolor, linestyle, linewidth,
35 separatepanel, plotmasked, maskedcolor,
36 showaxislabel, showtick, showticklabel, parsed_size,
37 nx, ny, transparent, plotrange)
40NoData = -32767.0
41NoDataThreshold = NoData + 10000.0
42LightSpeedQuantity = qa.constants('c')
43LightSpeed = qa.convert(LightSpeedQuantity, 'km/s')['value'] # speed of light in km/s
44DPIDetail = 130
46dsyb = '$^\circ$'
47hsyb = ':'
48msyb = ':'
51def Deg2HMS(x, arrowance):
52 # Transform degree to HHMMSS.sss format
53 xx = x % 360 + arrowance
54 h = int(xx / 15)
55 m = int((xx % 15) * 4)
56 s = ((xx % 15) * 4 - m) * 60.0
57 return (h, m, s)
60def HHMMSSss(x, pos):
61 # HHMMSS.ss format
62 (h, m, s) = Deg2HMS(x, 1 / 240000.0)
63 #return '%02dh%02dm%05.2fs' % (h, m, s)
64 return '%02d%s%02d%s%05.2f' % (h, hsyb, m, msyb, s)
67def Deg2DMS(x, arrowance):
68 # Transform degree to +ddmmss.ss format
69 xxx = (x + 90) % 180 - 90
70 xx = abs(xxx) + arrowance
71 if xxx < 0:
72 sign = -1
73 else:
74 sign = 1
75 d = int(xx * sign)
76 m = int((xx % 1) * 60)
77 s = ((xx % 1) * 60 - m) * 60.0
78 return (d, m, s)
81def DDMMSSs(x, pos):
82 # +DDMMSS.s format
83 (d, m, s) = Deg2DMS(x, 1 / 360000.0)
84 #return '%+02dd%02dm%04.1fs' % (d, m, s)
85 sint = int(s)
86 sstr = ('%3.1f' % (s - int(s))).lstrip('0')
87 return '%+02d%s%02d\'%02d\"%s' % (d, dsyb, m, sint, sstr)
90def parse_figsize(figsize):
91 """
92 return figsize in inches
93 """
94 casalog.post('parse_figsize input: {input}'.format(input=figsize), priority='DEBUG')
95 parsed = None
96 if figsize is not None and isinstance(figsize, str) and len(figsize) > 0:
97 size_list = figsize.split(',')
98 size_inch_list = [x / 25.4 for s in size_list for x in qa.getvalue(qa.convert(qa.quantity(s),outunit='mm'))]
99 if len(size_inch_list) == 1:
100 s = size_inch_list[0]
101 parsed = (s, s)
102 else:
103 parsed = tuple(size_inch_list[:2])
104 casalog.post('parse_figsize output: {output}'.format(output=parsed), priority='DEBUG')
105 return parsed
108def parse_numpanels(numpanels):
109 parsed = (-1, -1)
110 if isinstance(numpanels, str) and len(numpanels) > 0:
111 n = list(map(int, numpanels.split(',')))
112 if len(n) == 1:
113 parsed = (n[0], n[0])
114 else:
115 parsed = (n[0], n[1])
116 casalog.post('parse_numpanels output: {output}'.format(output=parsed), priority='DEBUG')
117 return parsed
120class ProfileMapAxesManager(object):
121 label_map = {'Right Ascension': 'RA',
122 'Declination': 'Dec'}
124 def __init__(self, nh, nv, brightnessunit, direction_label, direction_reference,
125 spectral_label, spectral_unit, ticksize, title='',
126 separatepanel=True,
127 showaxislabel=False, showtick=False, showticklabel=False,
128 figsize=None,
129 clearpanel=True):
130 self.nh = nh
131 self.nv = nv
132 self.ticksize = ticksize
133 self.brightnessunit = brightnessunit
134 self.numeric_formatter = plt.FormatStrFormatter('%.2f')
135 self.direction_label = direction_label
136 self.direction_reference = direction_reference
137 self.separatepanel = separatepanel
138 self.spectral_label = spectral_label
139 self.spectral_unit = spectral_unit
140 self.showaxislabel = showaxislabel
141 self.showtick = showtick
142 self.showticklabel = showticklabel
143 self.title = title
144 self.figsize = figsize
145 casalog.post('figsize={figsize}'.format(figsize=self.figsize), priority='DEBUG')
147 self.normalization_factor = 1
149 self._axes_spmap = None
151 # to resize matplotlib window to specified size
152 plt.figure(self.MATPLOTLIB_FIGURE_ID)
153 plt.close()
155 if self.figsize is None:
156 plt.figure(self.MATPLOTLIB_FIGURE_ID)
157 else:
158 plt.figure(self.MATPLOTLIB_FIGURE_ID, figsize=self.figsize)
159 if clearpanel:
160 plt.clf()
162 @property
163 def MATPLOTLIB_FIGURE_ID(self):
164 return "ProfileMap"
166 @property
167 def axes_spmap(self):
168 if self._axes_spmap is None:
169 self._axes_spmap = list(self.__axes_spmap())
171 return self._axes_spmap
173 @property
174 def nrow(self):
175 return self.nv
177 @property
178 def ncolumn(self):
179 return self.nh + 1
181 @property
182 def left_margin(self):
183 return 0.01 + 0.2 / self.ncolumn
185 @property
186 def right_margin(self):
187 return 0.01
189 @property
190 def bottom_margin(self):
191 return 0.01 + 0.5 / self.nrow
193 @property
194 def top_margin(self):
195 return 0.01
197 @property
198 def horizontal_space(self):
199 if self.separatepanel:
200 return self.horizontal_subplot_size * 0.1
201 else:
202 return 0.
204 @property
205 def vertical_space(self):
206 if self.separatepanel:
207 return self.vertical_subplot_size * 0.1
208 else:
209 return 0.
211 @property
212 def xlabel_area(self):
213 if self.showaxislabel or (self.showtick and self.showticklabel):
214 return 0.02
215 else:
216 return 0.
218 @property
219 def ylabel_area(self):
220 if self.showaxislabel or (self.showtick and self.showticklabel):
221 return 0.03
222 else:
223 return 0.
225 @property
226 def title_area(self):
227 if isinstance(self.title, str) and len(self.title) > 0:
228 return 0.04 * (self.title.count('\n') + 1)
229 else:
230 return 0.
232 @property
233 def horizontal_subplot_size(self):
234 return (1.0 - self.left_margin - self.right_margin - self.ylabel_area) / self.ncolumn
236 @property
237 def vertical_subplot_size(self):
238 return (1.0 - self.bottom_margin - self.top_margin - self.xlabel_area - self.title_area) / self.nrow
240 def set_normalization_factor(self, factor):
241 self.normalization_factor = factor
243 def __axes_spmap(self):
244 for x in range(self.nh):
245 for y in range(self.nv):
246 width = self.horizontal_subplot_size
247 height = self.vertical_subplot_size
248 left = 1.0 - self.right_margin - width * (x + 1) + 0.5 * self.horizontal_space
249 bottom = self.bottom_margin + self.ylabel_area + height * y + 0.5 * self.vertical_space
250 axes = plt.axes([left, bottom, width - self.horizontal_space, height - self.vertical_space])
251 axes.cla()
252 if self.showaxislabel and y == 0 and x == self.nh - 1:
253 casalog.post('label "{label}" unit "{unit}"'.format(label=self.spectral_label, unit=self.spectral_unit), priority='DEBUG')
254 if self.spectral_unit is not None and len(self.spectral_unit) > 0:
255 spectral_label = '%s [%s]' % (self.spectral_label, self.spectral_unit)
256 else:
257 spectral_label = self.spectral_label
258 axes.xaxis.set_label_text(spectral_label,
259 size=self.ticksize)
260 if self.normalization_factor < 100 and self.normalization_factor > 0.01:
261 label_text = 'Intensity [%s]' % self.brightnessunit
262 else:
263 label_text = 'Intensity [1e%d x %s]' % (int(np.log10(self.normalization_factor)),
264 self.brightnessunit)
265 axes.yaxis.set_label_text(label_text,
266 size=self.ticksize, rotation='vertical')
267 if self.showtick:
268 axes.xaxis.tick_bottom()
269 axes.yaxis.tick_left()
270 if self.showticklabel:
271 xlocator = plt.AutoLocator()
272 xlocator.set_params(nbins=6, prune='upper')
273 axes.xaxis.set_major_locator(xlocator)
274 ylocator = plt.AutoLocator()
275 ylocator.set_params(nbins=4)
276 axes.yaxis.set_major_locator(ylocator)
277 xformatter = plt.ScalarFormatter(useOffset=False)
278 axes.xaxis.set_major_formatter(xformatter)
279 axes.xaxis.set_tick_params(labelsize=max(self.ticksize - 1, 1))
280 axes.yaxis.set_tick_params(labelsize=max(self.ticksize - 1, 1))
281 if y != 0 or x != self.nh - 1:
282 axes.xaxis.set_major_formatter(plt.NullFormatter())
283 axes.yaxis.set_major_formatter(plt.NullFormatter())
284 else:
285 axes.xaxis.set_major_formatter(plt.NullFormatter())
286 axes.yaxis.set_major_formatter(plt.NullFormatter())
288 else:
289 axes.yaxis.set_major_locator(plt.NullLocator())
290 axes.xaxis.set_major_locator(plt.NullLocator())
292 yield axes
294 def setup_labels(self, label_ra, label_dec):
295 for x in range(self.nh):
296 width = self.horizontal_subplot_size
297 left = 1.0 - self.right_margin - width * (x + 1)
298 height = self.bottom_margin * 0.5
299 bottom = self.bottom_margin - height
300 a1 = plt.axes([left, bottom, width, height])
301 a1.set_axis_off()
302 if len(a1.texts) == 0:
303 plt.text(0.5, 0.2, HHMMSSss((label_ra[x][0] + label_ra[x][1]) / 2.0, 0),
304 horizontalalignment='center', verticalalignment='center', size=self.ticksize)
305 else:
306 a1.texts[0].set_text(HHMMSSss((label_ra[x][0] + label_ra[x][1]) / 2.0, 0))
307 for y in range(self.nv):
308 left = self.left_margin
309 width = self.horizontal_subplot_size
310 height = self.vertical_subplot_size
311 bottom = self.bottom_margin + y * height
312 a1 = plt.axes([left, bottom, width, height])
313 a1.set_axis_off()
314 if len(a1.texts) == 0:
315 plt.text(0.5, 0.5, DDMMSSs((label_dec[y][0] + label_dec[y][1]) / 2.0, 0),
316 horizontalalignment='center', verticalalignment='center', size=self.ticksize)
317 else:
318 a1.texts[0].set_text(DDMMSSs((label_dec[y][0] + label_dec[y][1]) / 2.0, 0))
320 # longitude label
321 left = self.left_margin + self.xlabel_area
322 height = self.bottom_margin * 0.5
323 bottom = 0.
324 width = 1.0 - left - self.right_margin
325 a1 = plt.axes([left, bottom, width, height])
326 a1.set_axis_off()
327 xpos = (1.0 + 0.5 * self.nh) / self.ncolumn
328 casalog.post('xpos=%s' % (xpos), priority='DEBUG')
329 plt.text(xpos, 0.5, '%s (%s)' % (self.direction_label[0], self.direction_reference),
330 horizontalalignment='center', verticalalignment='center',
331 size=(self.ticksize + 2))
333 # latitude label
334 left = 0.0
335 width = self.left_margin
336 height = self.vertical_subplot_size
337 bottom = self.bottom_margin + 0.5 * (height * self.nrow - self.vertical_subplot_size)
338 a1 = plt.axes([left, bottom, width, height])
339 a1.set_axis_off()
340 plt.text(1.0, 0.5, '%s (%s)' % (self.direction_label[1], self.direction_reference),
341 horizontalalignment='right', verticalalignment='center',
342 rotation='vertical', size=(self.ticksize + 2))
344 # title
345 if self.title_area > 0.:
346 left = self.left_margin + self.xlabel_area
347 bottom = 1.0 - self.title_area - self.top_margin
348 width = 1.0 - left - self.right_margin
349 height = self.title_area
350 a1 = plt.axes([left, bottom, width, height])
351 a1.set_axis_off()
352 xpos = (1.0 + 0.5 * self.nh) / self.ncolumn
353 plt.text(xpos, 0.1, self.title,
354 horizontalalignment='center', verticalalignment='bottom',
355 size=self.ticksize + 4)
358def plot_profile_map(image_data, figfile, pol, spectralaxis='', restfreq=None, title=None,
359 linecolor='b', linestyle='-', linewidth=0.2,
360 separatepanel=True, plotmasked=None, maskedcolor=None,
361 showaxislabel=False, showtick=False, showticklabel=False,
362 figsize=None, nx=-1, ny=-1, transparent=False, user_xlim=None):
363 """
364 image_data
365 figfile
366 linecolor
367 linestyle
368 linewidth
369 separatepanel
370 plotmasked -- 'none': show neither lines nor panels
371 'empty': show empty panel
372 'text': show 'NO DATA' symbol
373 'zero': plot zero level
374 'plot': plot masked data with different color
375 """
376 if linecolor is None:
377 linecolor = 'b'
378 if linestyle is None:
379 linestyle = '-'
380 if linewidth is None:
381 linewidth = 0.2
382 if separatepanel is None:
383 separatepanel = True
384 if plotmasked is None:
385 plotmasked = 'none'
386 if maskedcolor is None:
387 maskedcolor = 'gray' if linecolor != 'gray' else 'black'
388 if showaxislabel is None:
389 showaxislabel = False
390 if showtick is None:
391 showtick = False
392 if showticklabel is None:
393 showticklabel = False
394 if transparent is None:
395 transparent = False
397 # user-specified range for horizontal (spectral) axis
398 user_xmin = None
399 user_xmax = None
400 if isinstance(user_xlim, str) and user_xlim.find('~') != -1:
401 user_range = user_xlim.split('~')
402 if len(user_range[0]) > 0:
403 user_xmin = float(user_range[0])
404 if len(user_range[1]) > 0:
405 user_xmax = float(user_range[1])
407 x_max = image_data.nx - 1
408 x_min = 0
409 y_max = image_data.ny - 1
410 y_min = 0
411 MaxPanel = 8
412 num_panel = min(max(x_max - x_min + 1, y_max - y_min + 1), MaxPanel)
413 STEP = int((max(x_max - x_min + 1, y_max - y_min + 1) - 1) / num_panel) + 1
414 NH = (x_max - x_min) // STEP + 1
415 NV = (y_max - y_min) // STEP + 1
416 xSTEP = STEP
417 ySTEP = STEP
419 if nx > 0:
420 NH = nx
421 xSTEP = image_data.nx // NH + (1 if image_data.nx % NH > 0 else 0)
422 if ny > 0:
423 NV = ny
424 ySTEP = image_data.ny // NV + (1 if image_data.ny % NV > 0 else 0)
426 casalog.post('num_panel=%s, xSTEP=%s, ySTEP=%s, NH=%s, NV=%s' % (num_panel, xSTEP, ySTEP, NH, NV))
428 chan0 = 0
429 chan1 = image_data.nchan
431 direction_label = image_data.direction_label
432 direction_reference = image_data.direction_reference
433 default_spectral_label = image_data.spectral_label
434 default_spectral_unit = image_data.spectral_unit
435 if default_spectral_label == 'Frequency':
436 default_spectral_unit = 'GHz'
437 default_spectral_data = image_data.spectral_data[chan0:chan1]
438 if spectralaxis is None or spectralaxis == '' or spectralaxis == default_spectral_label.lower():
439 spectral_label = default_spectral_label
440 spectral_unit = default_spectral_unit
441 spectral_data = default_spectral_data
442 elif spectralaxis == 'channel':
443 spectral_label = 'Channel'
444 spectral_unit = ''
445 spectral_data = np.arange(chan0, chan1, dtype=np.int32)
446 else:
447 spectral_label = spectralaxis.capitalize()
448 if spectralaxis == 'frequency':
449 spectral_unit = 'GHz'
450 spectral_data = np.fromiter((image_data.to_frequency(v, freq_unit='GHz')
451 for v in default_spectral_data), dtype=np.float64)
452 elif spectralaxis == 'velocity':
453 if restfreq is not None:
454 casalog.post('User-specified rest frequency %s' % (restfreq))
455 else:
456 casalog.post('Default rest frequency from input image %s' % (image_data.coordsys.restfrequency()))
457 spectral_unit = 'km/s'
458 spectral_data = np.fromiter((image_data.to_velocity(f, freq_unit='GHz', restfreq=restfreq)
459 for f in default_spectral_data), dtype=np.float64)
461 masked_data = image_data.data * image_data.mask
462 masked_data[np.logical_not(np.isfinite(masked_data))] = 0.0
464 refpix = [0, 0]
465 refval = [0, 0]
466 increment = [0, 0]
467 refpix[0], refval[0], increment[0] = image_data.direction_axis(0, unit='deg')
468 refpix[1], refval[1], increment[1] = image_data.direction_axis(1, unit='deg')
470 stokes = image_data.stokes[pol]
471 casalog.post('Generate profile map for pol {stokes}'.format(stokes=stokes))
472 casalog.post('masked_data.shape=%s id_stokes=%s' % (list(masked_data.shape), image_data.id_stokes), priority='DEBUG')
473 masked_data_p = masked_data.take([pol], axis=image_data.id_stokes).squeeze(axis=image_data.id_stokes)
474 Plot = np.zeros((NH, NV, (chan1 - chan0)), np.float32) + NoData
475 mask_p = image_data.mask.take([pol], axis=image_data.id_stokes).squeeze(axis=image_data.id_stokes)
476 isvalid = np.any(mask_p, axis=2)
477 Nsp = sum(isvalid.flatten())
478 casalog.post('Nsp=%s' % (Nsp))
480 for x in range(NH):
481 x0 = x * xSTEP
482 x1 = (x + 1) * xSTEP
483 for y in range(NV):
484 y0 = y * ySTEP
485 y1 = (y + 1) * ySTEP
486 valid_index = isvalid[x0:x1, y0:y1].nonzero()
487 chunk = masked_data_p[x0:x1, y0:y1]
488 valid_sp = chunk[valid_index[0], valid_index[1], :]
489 if len(valid_sp) == 0:
490 Plot[x][y] = NoData
491 else:
492 Plot[x][y] = valid_sp.mean(axis=0)
494 # normalize plot data
495 plot_mask = np.logical_and(np.isfinite(Plot), Plot != NoData)
496 max_data = np.abs(Plot[plot_mask]).max()
497 casalog.post('max_data = %s' % (max_data), priority='DEBUG')
498 normalization_factor = np.power(10.0, int(np.log10(max_data)))
499 if normalization_factor < 1.0:
500 normalization_factor /= 10.
501 casalog.post('normalization_factor = %s' % (normalization_factor), priority='DEBUG')
503 Plot[plot_mask] /= normalization_factor
505 plotter = SDProfileMapPlotter(NH, NV, xSTEP, ySTEP, image_data.brightnessunit,
506 direction_label, direction_reference,
507 spectral_label, spectral_unit,
508 title=title,
509 separatepanel=separatepanel,
510 showaxislabel=showaxislabel,
511 showtick=showtick,
512 showticklabel=showticklabel,
513 figsize=figsize,
514 clearpanel=True)
516 try:
517 plotter.setup_labels(refpix, refval, increment)
518 plotter.set_normalization_factor(normalization_factor)
519 status = plotter.plot(figfile, Plot,
520 spectral_data,
521 linecolor=linecolor,
522 linestyle=linestyle,
523 linewidth=linewidth,
524 plotmasked=plotmasked,
525 maskedcolor=maskedcolor,
526 transparent=transparent,
527 user_xmin=user_xmin,
528 user_xmax=user_xmax)
530 finally:
531 plotter.done()
534class SDProfileMapPlotter(object):
535 def __init__(self, nh, nv, xstep, ystep, brightnessunit, direction_label, direction_reference,
536 spectral_label, spectral_unit, title=None, separatepanel=True,
537 showaxislabel=False, showtick=False, showticklabel=False,
538 figsize=None,
539 clearpanel=True):
540 self.xstep = xstep
541 self.ystep = ystep
542 if self.xstep > 1 or self.ystep > 1:
543 step = max(self.xstep, self.ystep)
544 ticksize = 10 - int(max(nh, nv) * step / (step - 1)) / 2
545 else:
546 ticksize = 10 - int(max(nh, nv)) / 2
547 self.axes = ProfileMapAxesManager(nh, nv, brightnessunit,
548 direction_label, direction_reference,
549 spectral_label, spectral_unit,
550 ticksize, title=title,
551 separatepanel=separatepanel,
552 showaxislabel=showaxislabel,
553 showtick=showtick,
554 showticklabel=showticklabel,
555 figsize=figsize,
556 clearpanel=clearpanel)
557 self.lines_averaged = None
558 self.lines_map = None
559 self.reference_level = None
560 self.global_scaling = True
561 self.deviation_mask = None
563 @property
564 def nh(self):
565 return self.axes.nh
567 @property
568 def nv(self):
569 return self.axes.nv
571 @property
572 def TickSize(self):
573 return self.axes.ticksize
575 def setup_labels(self, refpix_list, refval_list, increment_list):
576 LabelRA = np.zeros((self.nh, 2), np.float32) + NoData
577 LabelDEC = np.zeros((self.nv, 2), np.float32) + NoData
578 refpix = refpix_list[0]
579 refval = refval_list[0]
580 increment = increment_list[0]
581 #casalog.post('axis 0: refpix,refval,increment=%s,%s,%s'%(refpix,refval,increment))
582 for x in range(self.nh):
583 x0 = (self.nh - x - 1) * self.xstep
584 x1 = (self.nh - x - 2) * self.xstep + 1
585 LabelRA[x][0] = refval + (x0 - refpix) * increment
586 LabelRA[x][1] = refval + (x1 - refpix) * increment
587 refpix = refpix_list[1]
588 refval = refval_list[1]
589 increment = increment_list[1]
590 #casalog.post('axis 1: refpix,refval,increment=%s,%s,%s'%(refpix,refval,increment))
591 for y in range(self.nv):
592 y0 = y * self.ystep
593 y1 = (y + 1) * self.ystep - 1
594 LabelDEC[y][0] = refval + (y0 - refpix) * increment
595 LabelDEC[y][1] = refval + (y1 - refpix) * increment
596 self.axes.setup_labels(LabelRA, LabelDEC)
598 def setup_reference_level(self, level=0.0):
599 self.reference_level = level
601 def set_global_scaling(self):
602 self.global_scaling = True
604 def unset_global_scaling(self):
605 self.global_scaling = False
607 def set_normalization_factor(self, factor):
608 self.axes.set_normalization_factor(factor)
610 def plot(self, figfile, map_data, frequency,
611 linecolor='b', linestyle='-', linewidth=0.2,
612 plotmasked='none', maskedcolor='gray', transparent=False,
613 user_xmin=None, user_xmax=None):
614 if len(frequency) == 1:
615 casalog.post('Number of spectral channels is 1. Resulting profile map will be useless.',
616 priority='WARN')
617 global_xmin = frequency[0] - 1.
618 global_xmax = frequency[0] + 1.
619 else:
620 global_xmin = min(frequency[0], frequency[-1])
621 global_xmax = max(frequency[0], frequency[-1])
622 casalog.post('full spectral range: global_xmin=%s, global_xmax=%s' % (global_xmin, global_xmax))
624 if user_xmin is not None:
625 casalog.post('user-specified global_xmin %s' % (user_xmin))
626 global_xmin = user_xmin
627 if user_xmax is not None:
628 casalog.post('user-specified global_xmax %s' % (user_xmax))
629 global_xmax = user_xmax
631 if global_xmin == global_xmax:
632 raise RuntimeError('Wrong spectral axis range: no range to plot.')
633 elif global_xmin > global_xmax:
634 raise RuntimeError('Wrong spectral axis range: minimum > maximum')
636 # Auto scaling
637 # to eliminate max/min value due to bad pixel or bad fitting,
638 # 1/10-th value from max and min are used instead
639 valid_index = np.where(map_data.min(axis=2) > NoDataThreshold)
640 valid_data = map_data[valid_index[0], valid_index[1], :]
641 ListMax = valid_data.max(axis=1)
642 ListMin = valid_data.min(axis=1)
643 casalog.post('ListMax=%s' % (list(ListMax)))
644 casalog.post('ListMin=%s' % (list(ListMin)))
645 if len(ListMax) == 0:
646 return False
647 global_ymax = np.sort(ListMax)[len(ListMax) - len(ListMax) // 10 - 1]
648 global_ymin = np.sort(ListMin)[len(ListMin) // 10]
649 global_ymax = global_ymax + (global_ymax - global_ymin) * 0.2
650 global_ymin = global_ymin - (global_ymax - global_ymin) * 0.1
651 del ListMax, ListMin
653 casalog.post('global_ymin=%s, global_ymax=%s' % (global_ymin, global_ymax))
655 plt.ioff()
657 no_data = np.zeros(len(frequency), dtype=np.float32)
658 for x in range(self.nh):
659 for y in range(self.nv):
660 if self.global_scaling is True:
661 xmin = global_xmin
662 xmax = global_xmax
663 ymin = global_ymin
664 ymax = global_ymax
665 else:
666 xmin = global_xmin
667 xmax = global_xmax
668 if map_data[x][y].min() > NoDataThreshold:
669 median = np.median(map_data[x][y])
670 mad = np.median(map_data[x][y] - median)
671 sigma = map_data[x][y].std()
672 ymin = median - 2.0 * sigma
673 ymax = median + 5.0 * sigma
674 else:
675 ymin = global_ymin
676 ymax = global_ymax
677 casalog.post('Per panel scaling turned on: ymin=%s, ymax=%s (global ymin=%s, ymax=%s)' % (ymin, ymax, global_ymin, global_ymax))
678 plt.gcf().sca(self.axes.axes_spmap[y + (self.nh - x - 1) * self.nv])
679 if map_data[x][y].min() > NoDataThreshold:
680 plt.plot(frequency, map_data[x][y], color=linecolor, linestyle=linestyle,
681 linewidth=linewidth)
682 else:
683 if plotmasked == 'empty':
684 pass
685 elif plotmasked == 'text':
686 plt.text((xmin + xmax) / 2.0, (ymin + ymax) / 2.0, 'NO DATA', ha='center', va='center',
687 size=(self.TickSize + 1))
688 elif plotmasked == 'zero':
689 plt.plot(frequency, no_data,
690 color=maskedcolor, linestyle=linestyle, linewidth=linewidth)
691 elif plotmasked == 'none':
692 a = plt.gcf().gca()
693 if self.axes.xlabel_area > 0. and y == 0 and x == 0:
694 pass
695 else:
696 a.set_axis_off()
697 elif plotmasked == 'plot':
698 m = map_data[x][y] == NoDataThreshold
699 if not all(m):
700 ma = np.ma.masked_array(map_data[x][y], m)
701 plt.plot(frequency, ma,
702 color=maskedcolor, linestyle=linestyle, linewidth=linewidth)
704 plt.axis((xmin, xmax, ymin, ymax))
706 plt.ion()
707 plt.draw()
709 casalog.post('figfile=\'%s\'' % (figfile), priority='DEBUG')
710 if figfile is not None and len(figfile) > 0:
711 casalog.post('Output profile map to %s' % (figfile))
712 plt.savefig(figfile, format='png', dpi=DPIDetail, transparent=transparent)
714 return True
716 def done(self):
717 #plt.close()
718 del self.axes
721class SpectralImage(object):
722 def __init__(self, imagename):
723 self.imagename = imagename
724 # read data to storage
725 ia = image()
726 ia.open(self.imagename)
727 try:
728 self.image_shape = ia.shape()
729 self.coordsys = ia.coordsys()
730 coord_types = self.coordsys.axiscoordinatetypes()
731 self.units = self.coordsys.units()
732 self.names = self.coordsys.names()
733 self.direction_reference = self.coordsys.referencecode('dir')[0]
734 self.spectral_reference = self.coordsys.referencecode('spectral')[0]
735 self.id_direction = coord_types.index('Direction')
736 self.id_direction = [self.id_direction, self.id_direction + 1]
737 self.id_spectral = coord_types.index('Spectral')
738 try:
739 self.id_stokes = coord_types.index('Stokes')
740 stokes_axis_exists = True
741 except Exception:
742 # if no Stokes axis exists, set dummy axis at the end
743 self.id_stokes = len(coord_types)
744 stokes_axis_exists = False
745 casalog.post('id_direction=%s' % (self.id_direction), priority='DEBUG')
746 casalog.post('id_spectral=%s' % (self.id_spectral), priority='DEBUG')
747 casalog.post('id_stokes=%s (%s stokes axis)' % (self.id_stokes, ('real' if stokes_axis_exists else 'dummy')),
748 priority='DEBUG')
749 self.data = ia.getchunk()
750 self.mask = ia.getchunk(getmask=True)
751 if not stokes_axis_exists:
752 # put degenerate Stokes axis at the end
753 self.data = np.expand_dims(self.data, axis=-1)
754 self.mask = np.expand_dims(self.mask, axis=-1)
755 casalog.post('add dummy Stokes axis to data and mask since input image doesn\'t have it.')
756 casalog.post('data.shape={}'.format(self.data.shape), priority='DEBUG')
757 casalog.post('mask.shape={}'.format(self.mask.shape), priority='DEBUG')
758 bottom = ia.toworld(np.zeros(len(self.image_shape), dtype=int), 'q')['quantity']
759 top = ia.toworld(self.image_shape - 1, 'q')['quantity']
760 key = lambda x: '*%s' % (x + 1)
761 ra_min = bottom[key(self.id_direction[0])]
762 ra_max = top[key(self.id_direction[0])]
763 if qa.gt(ra_min, ra_max):
764 ra_min, ra_max = ra_max, ra_min
765 self.ra_min = ra_min
766 self.ra_max = ra_max
767 self.dec_min = bottom[key(self.id_direction[1])]
768 self.dec_max = top[key(self.id_direction[1])]
769 self._brightnessunit = ia.brightnessunit()
770 if self.spectral_label == 'Frequency':
771 refpix, refval, increment = self.spectral_axis(unit='GHz')
772 self.spectral_data = np.array([refval + increment * (i - refpix) for i in range(self.nchan)])
773 elif self.spectral_label == 'Velocity':
774 refpix, refval, increment = self.spectral_axis(unit='km/s')
775 self.spectral_data = np.array([refval + increment * (i - refpix) for i in range(self.nchan)])
776 if stokes_axis_exists:
777 self.stokes = self.coordsys.stokes()
778 else:
779 # if no Stokes axis exists, set ['I']
780 self.stokes = ['I']
781 finally:
782 ia.close()
784 @property
785 def nx(self):
786 return self.image_shape[self.id_direction[0]]
788 @property
789 def ny(self):
790 return self.image_shape[self.id_direction[1]]
792 @property
793 def nchan(self):
794 return self.image_shape[self.id_spectral]
796 @property
797 def npol(self):
798 return self.image_shape[self.id_stokes]
800 @property
801 def brightnessunit(self):
802 return self._brightnessunit
804 @property
805 def direction_label(self):
806 return [self.names[i] for i in self.id_direction]
808 @property
809 def spectral_label(self):
810 return self.names[self.id_spectral]
812 @property
813 def spectral_unit(self):
814 return self.units[self.id_spectral]
816 def to_velocity(self, frequency, freq_unit='GHz', restfreq=None):
817 rest_frequency = self.coordsys.restfrequency()
818 # user-defined rest frequency takes priority
819 if restfreq is not None:
820 vrf = qa.convert(qa.quantity(restfreq), freq_unit)['value']
821 elif rest_frequency['unit'] != freq_unit:
822 vrf = qa.convert(rest_frequency, freq_unit)['value']
823 else:
824 vrf = rest_frequency['value']
825 return (1.0 - (frequency / vrf)) * LightSpeed
827 def to_frequency(self, velocity, freq_unit='km/s'):
828 rest_frequency = self.coordsys.restfrequency()
829 if rest_frequency['unit'] != freq_unit:
830 vrf = qa.convert(rest_frequency, freq_unit)['value']
831 else:
832 vrf = rest_frequency['value']
833 return (1.0 - velocity / LightSpeed) * vrf
835 def spectral_axis(self, unit='GHz'):
836 return self.__axis(self.id_spectral, unit=unit)
838 def direction_axis(self, idx, unit='deg'):
839 return self.__axis(self.id_direction[idx], unit=unit)
841 def __axis(self, idx, unit):
842 refpix = self.coordsys.referencepixel()['numeric'][idx]
843 refval = self.coordsys.referencevalue()['numeric'][idx]
844 increment = self.coordsys.increment()['numeric'][idx]
845 _unit = self.units[idx]
846 if _unit != unit:
847 refval = qa.convert(qa.quantity(refval, _unit), unit)['value']
848 increment = qa.convert(qa.quantity(increment, _unit), unit)['value']
849 #return np.array([refval+increment*(i-refpix) for i in range(self.nchan)])
850 return (refpix, refval, increment)