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

1import os 

2 

3import matplotlib.pyplot as plt 

4import numpy as np 

5 

6from casatasks import casalog 

7from casatools import image, quanta 

8 

9from . import sdutil 

10 

11qa = quanta() 

12 

13 

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

21 

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

24 

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

29 

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) 

38 

39 

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 

45 

46dsyb = '$^\circ$' 

47hsyb = ':' 

48msyb = ':' 

49 

50 

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) 

58 

59 

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) 

65 

66 

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) 

79 

80 

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) 

88 

89 

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 

106 

107 

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 

118 

119 

120class ProfileMapAxesManager(object): 

121 label_map = {'Right Ascension': 'RA', 

122 'Declination': 'Dec'} 

123 

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

146 

147 self.normalization_factor = 1 

148 

149 self._axes_spmap = None 

150 

151 # to resize matplotlib window to specified size 

152 plt.figure(self.MATPLOTLIB_FIGURE_ID) 

153 plt.close() 

154 

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

161 

162 @property 

163 def MATPLOTLIB_FIGURE_ID(self): 

164 return "ProfileMap" 

165 

166 @property 

167 def axes_spmap(self): 

168 if self._axes_spmap is None: 

169 self._axes_spmap = list(self.__axes_spmap()) 

170 

171 return self._axes_spmap 

172 

173 @property 

174 def nrow(self): 

175 return self.nv 

176 

177 @property 

178 def ncolumn(self): 

179 return self.nh + 1 

180 

181 @property 

182 def left_margin(self): 

183 return 0.01 + 0.2 / self.ncolumn 

184 

185 @property 

186 def right_margin(self): 

187 return 0.01 

188 

189 @property 

190 def bottom_margin(self): 

191 return 0.01 + 0.5 / self.nrow 

192 

193 @property 

194 def top_margin(self): 

195 return 0.01 

196 

197 @property 

198 def horizontal_space(self): 

199 if self.separatepanel: 

200 return self.horizontal_subplot_size * 0.1 

201 else: 

202 return 0. 

203 

204 @property 

205 def vertical_space(self): 

206 if self.separatepanel: 

207 return self.vertical_subplot_size * 0.1 

208 else: 

209 return 0. 

210 

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. 

217 

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. 

224 

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. 

231 

232 @property 

233 def horizontal_subplot_size(self): 

234 return (1.0 - self.left_margin - self.right_margin - self.ylabel_area) / self.ncolumn 

235 

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 

239 

240 def set_normalization_factor(self, factor): 

241 self.normalization_factor = factor 

242 

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

287 

288 else: 

289 axes.yaxis.set_major_locator(plt.NullLocator()) 

290 axes.xaxis.set_major_locator(plt.NullLocator()) 

291 

292 yield axes 

293 

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

319 

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

332 

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

343 

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) 

356 

357 

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 

396 

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

406 

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 

418 

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) 

425 

426 casalog.post('num_panel=%s, xSTEP=%s, ySTEP=%s, NH=%s, NV=%s' % (num_panel, xSTEP, ySTEP, NH, NV)) 

427 

428 chan0 = 0 

429 chan1 = image_data.nchan 

430 

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) 

460 

461 masked_data = image_data.data * image_data.mask 

462 masked_data[np.logical_not(np.isfinite(masked_data))] = 0.0 

463 

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

469 

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

479 

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) 

493 

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

502 

503 Plot[plot_mask] /= normalization_factor 

504 

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) 

515 

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) 

529 

530 finally: 

531 plotter.done() 

532 

533 

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 

562 

563 @property 

564 def nh(self): 

565 return self.axes.nh 

566 

567 @property 

568 def nv(self): 

569 return self.axes.nv 

570 

571 @property 

572 def TickSize(self): 

573 return self.axes.ticksize 

574 

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) 

597 

598 def setup_reference_level(self, level=0.0): 

599 self.reference_level = level 

600 

601 def set_global_scaling(self): 

602 self.global_scaling = True 

603 

604 def unset_global_scaling(self): 

605 self.global_scaling = False 

606 

607 def set_normalization_factor(self, factor): 

608 self.axes.set_normalization_factor(factor) 

609 

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

623 

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 

630 

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

635 

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 

652 

653 casalog.post('global_ymin=%s, global_ymax=%s' % (global_ymin, global_ymax)) 

654 

655 plt.ioff() 

656 

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) 

703 

704 plt.axis((xmin, xmax, ymin, ymax)) 

705 

706 plt.ion() 

707 plt.draw() 

708 

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) 

713 

714 return True 

715 

716 def done(self): 

717 #plt.close() 

718 del self.axes 

719 

720 

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

783 

784 @property 

785 def nx(self): 

786 return self.image_shape[self.id_direction[0]] 

787 

788 @property 

789 def ny(self): 

790 return self.image_shape[self.id_direction[1]] 

791 

792 @property 

793 def nchan(self): 

794 return self.image_shape[self.id_spectral] 

795 

796 @property 

797 def npol(self): 

798 return self.image_shape[self.id_stokes] 

799 

800 @property 

801 def brightnessunit(self): 

802 return self._brightnessunit 

803 

804 @property 

805 def direction_label(self): 

806 return [self.names[i] for i in self.id_direction] 

807 

808 @property 

809 def spectral_label(self): 

810 return self.names[self.id_spectral] 

811 

812 @property 

813 def spectral_unit(self): 

814 return self.units[self.id_spectral] 

815 

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 

826 

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 

834 

835 def spectral_axis(self, unit='GHz'): 

836 return self.__axis(self.id_spectral, unit=unit) 

837 

838 def direction_axis(self, idx, unit='deg'): 

839 return self.__axis(self.id_direction[idx], unit=unit) 

840 

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)