Coverage for /wheeldirectory/casa-6.7.0-12-py3.10.el8/lib/py/lib/python3.10/site-packages/casatasks/private/task_plotprofilemap.py: 16%

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