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