Coverage for /wheeldirectory/casa-6.7.0-12-py3.10.el8/lib/py/lib/python3.10/site-packages/casatasks/private/task

1########################################################################3

2# task_immath.py

5# Associated Universities, Inc. Washington DC, USA.

7# This script is free software; you can redistribute it and/or modify it

8# under the terms of the GNU Library General Public License as published by

9# the Free Software Foundation; either version 2 of the License, or (at your

10# option) any later version.

11#

12# This library is distributed in the hope that it will be useful, but WITHOUT

13# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or

14# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public

15# License for more details.

16#

17# You should have received a copy of the GNU Library General Public License

18# along with this library; if not, write to the Free Software Foundation,

19# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.

20#

21# Correspondence concerning AIPS++ should be adressed as follows:

22# Internet email: casa-feedback@nrao.edu.

23# Postal address: AIPS++ Project Office

24# National Radio Astronomy Observatory

25# 520 Edgemont Road

26# Charlottesville, VA 22903-2475 USA

27#

28# <summary>

29# CASA task for smoothing an image, by doing Forier-based convolution

30# on a CASA image file.

31# </summary>

32#

33# <reviewed reviwer="" date="" tests="" demos="">

34# </reviewed

35#

36# <author>

37# Shannon Jaeger, University of Calgary (image math)

38# Takeshi Nakazato, National Radio Astronomy Obaservatory (polarization)

39# </author>

40#

41# <prerequisite>

42# </prerequisite>

43#

44# <etymology>

45# immath stands for image mathematics

46# </etymology>

47#

48# <synopsis>

49# This task evaluates mathematical expressions involving existing

50# image files. The results of the calculations are stored in the

51# designated output file. Options are available to specify mathematical

52# expression directly or pre-defined expression for calculation of

53# spectral index image, and polarization intensity and position angle

54# images are available. The image file names imbedded in the expression or

55# specified in the imagename parameter for the pre-defined calculations may

56# be CASA images or FITS images.

57#

58#

59# NOTE: Index values for axes start at 0 for the box and chans

60# parameters, but at 1 when used with the indexin function

61# in expression. Use the imhead task to see the range of

62# values for each axes.

63#

64#

65# Keyword arguments:

66# outfile -- The file where the results of the image calculations

67# are stored. Overwriting an existing outfile is not permitted.

68# Default: none; Example: outfile='results.im'

69# mode -- mode for mathematical operation

70# Default: evalexpr

71# Options: 'evalexpr' : evalulate a mathematical expression defined in 'expr'

72# 'spix' : spectalindex image

73# 'pola' : polarization position angle image

74# 'poli' : polarization intesity image

75# mode expandable parameters

76# expr -- (for mode='evalexpr') A mathematical expression, with image file names.

77# Image file names MUST be enclosed in double quotes (")

78# Default: none

79# Examples:

80# Make an image that is image1.im - image2.im

81# expr=' ("image1.im" - "image2.im" )'

82# Clip an image below a value (0.5 in this case)

83# expr = ' iif("image1.im">=0.5, "image1.im", 0.0) '

84# Note: iif (a, b, c) a is the boolian expression

85# b is the value if true

86# c is the value if false

87# Take the rms value of two images

88# expr = ' sqrt("image1.im" * "image1.im" + "image2.im" * "image2.im") '

89# Note: No exponentiaion available?

90# Build an image pixel by pixel from the minimum of (image2.im, 2*image1.im)

91# expr='min("image2.im",2*max("image1.im"))'

92# imagename -- (for mode='spix','pola','poli') input image names

93# Default: none;

94# Examples: mode='spix'; imagename=['image1.im','image2.im'] will calculate

95# an image of log(S1/S2)/log(f1/f2), where S1 and S2 are fluxes and

96# f1 and f2 are frequencies

97# mode='pola'; imagename=['imageQ.im','imageU.im'] will calculate

98# an image of polarization angle distribution, where imageQ.im and

99# imageU.im are Stokes Q and U images, respectively. Calculate 0.5*arctan(U/Q).

100# mode='poli'; imagename=['imageQ.im','imageU.im','imageV.im'] will calculate

101# total polarization intensity image, where imageQ.im, imageU.im, imageV.im

102# are Stokes Q, U, and V images, respectively.

103# sigma - (for mode='poli') standard deviation of noise of Stokes images with unit such as

104# Jy/beam to correct for bias

105# Default: '0.0Jy/beam' (= no debiasing)

106# mask -- Name of mask applied to each image in the calculation

107# Default '' means no mask; Example: mask='orion.mask'.

108# region -- File path to an ImageRegion file.

109# An ImageRegion file can be created with the CASA

110# viewer's region manager. Typically ImageRegion files

111# will have the suffix '.rgn'. If a region file is given

112# then the box, chans, and stokes selections whill be

113# ignored.

114# Default: none

115# Example: region='myimage.im.rgn'

116# box -- A box region on the directional plane

117# Only pixel values acceptable at this time.

118# Default: none (whole 2-D plane); Example: box='10,10,50,50'

119# chans -- channel numbers, velocity, and/or frequency

120# Only channel numbers acceptable at this time.

121# Default: none (all); Example: chans='3~20'

122# stokes -- Stokes parameters to image, may or may not be separated

123# by commas but best if you use commas.

124# Default: none (all); Example: stokes='IQUV';

125# Options: 'I','Q','U','V','RR','RL','LR','LL','XX','YX','XY','YY', ...

126#

127# Available functions in the <i>expr</i> and <i>mask</i> paramters:

128# pi(), e(), sin(), sinh(), asinh(), cos(), cosh(), tan(), tanh(),

129# atan(), exp(), log(), log10(), pow(), sqrt(), complex(), conj()

130# real(), imag(), abs(), arg(), phase(), aplitude(), min(), max()

131# round(), isgn(), floor(), ceil(), rebin(), spectralindex(), pa(),

132# iif(), indexin(), replace(), ...

133#

134# For a full description of the allowed syntax see the

135# Lattice Expression Language (LEL) documentation on the at:

136# http://aips2.nrao.edu/docs/notes/223/223.html

137#

138# NOTE: where indexing and axis numbering are used in the above

139# functions they are 1-based, ie. numbering starts at 1.

140#

141# </synopsis>

142#

143# <example>

144# <srcblock>

145# </srcblock

146#

147# </example>

148#

149# <motivation>

150# To provide a user-friendly task interface to imagecalc and ???

151# as well as an more user-friendling math syntax then what is

152# provided by the CASA Lattice Exprssion Language.

153# </motivation>

154#

155# <todo>

156# Crystal wanted different masks for different inputs

157# but unlikely that its really needed.

158#

159# Add an "overwrite" output file parameter

160#

161# Add polygon and circle region selection

162# </todo>

163########################################################################3

164

165import os

166import re

167import shutil

168import sys

169

170from casatools import image, imagepol, regionmanager, coordsys, quanta

171from casatasks import casalog

172from .ialib import write_image_history

173

174def immath(

175 imagename, mode, outfile, expr, varnames, sigma,

176 polithresh, mask, region, box, chans, stokes, stretch,

177 imagemd, prec

178):

179 # Tell CASA who will be reporting

180 casalog.origin('immath')

181 tmpFilePrefix='_immath_tmp' + str(os.getpid()) + '_'

182 try:

183 myia = image()

184 myia.dohistory(False)

185 outia = None

186 _immath_initial_cleanup(tmpFilePrefix)

187 outfile = _immath_check_outfile(outfile)

188 # Find the list of filenames in the expression

189 # also do a quick check to see if all of the files

190 # exist

191 tmpfilenames = ''

192 filenames = imagename

193 if mode=='evalexpr':

194 tmpfilenames = _immath_parse(expr)

195 if isinstance(filenames, str):

196 filenames = [filenames]

197 varnames = _immath_varnames(varnames, filenames, tmpfilenames)

198 filenames = _immath_filenames(filenames, tmpfilenames, varnames, mode)

199 expr = expr.replace(' ', '')

200 if mode == 'spix':

201 expr = _immath_dospix(len(filenames), varnames)

202 if mode == 'pola':

203 _immath_new_pola(

204 filenames, outfile, tmpFilePrefix, mask, region,

205 box, chans, stokes, stretch, polithresh, myia

206 )

207 return True

208 elif mode == 'poli' or mode == 'lpoli' or mode == 'tpoli':

209 _immath_new_poli(

210 filenames, outfile, tmpFilePrefix, mask, region,

211 box, chans, stokes, stretch, sigma, myia, mode

212 )

213 return True

214 elif mode == 'evalexpr' or mode == 'spix':

215 if box or chans or stokes or region or mask:

216 (subImages, file_map) = _immath_createsubimages(

217 box, chans, stokes, region, mask,

218 stretch, filenames, myia, tmpFilePrefix

219 )

220 if imagemd:

221 casalog.post(

222 "Specifying region, box, chan, or stokes will "

223 + "create smaller sub-images. The image "

224 + "metadata specified in imagemd will have to "

225 + "conform to the output, not the input image "

226 + "dimensions. Please check your output image "

227 + "for accurate header definition.", 'WARN'

228 )

229 (expr, varnames, subImages) = _immath_updateexpr(

230 expr, varnames, subImages, filenames, file_map

231 )

232 outia = _immath_compute(

233 imagename, expr, outfile, imagemd, myia, prec

234 )

235 else:

236 # If the user didn't give any region or mask information

237 # then just evaluated the expression with the filenames in it.

238 outia = _immath_dofull(

239 imagename, imagemd, outfile, mode, expr,

240 varnames, filenames, myia, prec

241 )

242 else:

243 raise(Exception, "Unsupported mode " + str(mode))

244 try:

245 vars = locals()

246 param_names = immath.__code__.co_varnames[:immath.__code__.co_argcount]

247 param_vals = [vars[p] for p in param_names]

248 write_image_history(

249 outia, sys._getframe().f_code.co_name,

250 param_names, param_vals, casalog

251 )

252 except Exception as instance:

253 casalog.post("*** Error \'%s\' updating HISTORY" % (instance), 'WARN')

254 return True

255 except Exception as error:

256 if mode == 'evalexpr':

257 casalog.post("Unable to process expression " + expr, 'SEVERE')

258 else:

259 casalog.post("Error running immath", 'SEVERE')

260 casalog.post("Exception caught was: " + str(error), 'SEVERE')

261 raise

262 finally:

263 if myia:

264 myia.done()

265 if outia:

266 outia.done()

267 _immath_cleanup(tmpFilePrefix)

268

269def _immath_concat_stokes(filenames, target, myia):

270 myia.open(filenames[0])

271 stokes_axis = myia.coordsys().findaxisbyname("stokes")

272 myia.done()

273 casalog.post("Concatenating images along stokes axis")

274 myia = myia.imageconcat(

275 outfile=target, infiles=filenames, axis=stokes_axis

276 )

277 myia.done()

278

279def _immath_getregion(region, box, chans, stokes, mode, myia, target):

280 myreg = region

281 if (type(region) != type({})):

282 if stokes:

283 casalog.post(

284 "Ignoring stokes parameters selection for mode='"

285 + mode + "'."

286 ,'WARN'

287 )

288 stokes=''

289 myia.open(target)

290 myrg = regionmanager()

291 myreg = myrg.frombcs(

292 csys=myia.coordsys().torecord(), shape=myia.shape(), box=box,

293 chans=chans, stokes=stokes, stokescontrol="a", region=region

294 )

295 myia.done()

296 myrg.done()

297 return myreg

298

299def _immath_new_pola(

300 filenames, outfile, tmpFilePrefix, mask, region,

301 box, chans, stokes, stretch, polithresh, myia

302):

303 target = filenames[0]

304 if len(filenames) > 1:

305 target = tmpFilePrefix + "_concat_for_pola"

306 _immath_concat_stokes(filenames, target, myia)

307 myreg = _immath_getregion(region, box, chans, stokes, "pola", myia, target)

308 mypo = imagepol()

309 myqa = quanta()

310 if (polithresh):

311 if (mask != ""):

312 mask = ""

313 casalog.post("Ignoring mask parameter in favor of polithresh parameter", 'WARN')

314 if (myqa.getunit(polithresh) != ""):

315 initUnit = myqa.getunit(polithresh)

316 myia.dohistory(False)

317 myia.open(filenames[0])

318 bunit = myia.brightnessunit()

319 polithresh = myqa.convert(polithresh, bunit)

320 myia.done()

321 if (myqa.getunit(polithresh) != bunit):

322 raise Exception("Units of polithresh " + initUnit \

323 + " do not conform to input image units of " + bunit \

324 + " so cannot perform thresholding. Please correct units and try again.")

325 polithresh = myqa.getvalue(polithresh)[0]

326 lpol = tmpFilePrefix + "_lpol"

327 mypo.open(target)

328 myia = mypo.linpolint(debias=False, outfile=lpol, region=myreg)

329 myia.done()

330 mypo.done()

331 mypo.open(target)

332 myia = mypo.linpolposang(

333 outfile=outfile, region=myreg, mask=mask, stretch=stretch

334 )

335 mypo.done()

336 if (polithresh):

337 myexpr = "'" + lpol + "' >= " + str(polithresh)

338 myia.dohistory(False)

339 myia.calcmask(name='mask0', mask=myexpr)

340 casalog.post(

341 'Calculated mask based on linear polarization threshold '

342 + str(polithresh),

343 'INFO'

344 )

345 myia.done()

346

347def _immath_new_poli(

348 filenames, outfile, tmpFilePrefix, mask, region,

349 box, chans, stokes, stretch, sigma, myia, mode

350):

351 target = filenames[0]

352 if len(filenames) > 1:

353 target = tmpFilePrefix + "_concat_for_poli"

354 _immath_concat_stokes(filenames, target, myia)

355 if mode == 'tpoli':

356 myia.open(target)

357 csys = myia.coordsys()

358 myia.done()

359 stokes = csys.stokes()

360 csys.done()

361 for p in ["Q", "U", "V"]:

362 if stokes.count(p) == 0:

363 raise Exception(

364 "Stokes " + p + " is required for mode tpoli but was not found"

365 )

366 debias = False

367 newsigma = 0

368 myqa = quanta()

369 if sigma:

370 qsigma = myqa.quantity(sigma)

371 if myqa.getvalue(qsigma)[0] > 0:

372 debias = True

373 sigmaunit = myqa.getunit(qsigma)

374 try:

375 myia.open(filenames[0])

376 iunit = myia.brightnessunit()

377 myia.done()

378 except:

379 raise Exception('Unable to get brightness unit from image file ' + filenames[0])

380 if sigmaunit != iunit:

381 newsigma = myqa.convert(qsigma,iunit)

382 else:

383 newsigma = sigma

384 myreg = _immath_getregion(region, box, chans, stokes, "poli", myia, target)

385 mypo = imagepol()

386 mypo.open(target)

387 # for some annoying reason, qa.getvalue() returns an array in this context

388 numeric_sigma = myqa.getvalue(myqa.quantity(newsigma))[0]

389

390 if mode == 'tpoli' or mode == 'poli':

391 myia = mypo.totpolint(

392 debias=debias, sigma=numeric_sigma, outfile=outfile,

393 region=myreg, mask=mask, stretch=stretch

394 )

395 elif mode == 'lpoli':

396 myia = mypo.linpolint(

397 debias=debias, sigma=numeric_sigma, outfile=outfile,

398 region=myreg, mask=mask, stretch=stretch

399 )

400 else:

401 raise Exception("Logic Error: Unhandled mode " + mode)

402 myia.done()

403 mypo.done()

404

405def _immath_compute(

406 imagename, expr, outfile, imagemd, myia, prec

407):

408 # Do the calculation

409 res = myia.imagecalc(

410 pixels=expr, outfile=outfile,

411 imagemd=_immath_translate_imagemd(imagename, imagemd), prec=prec

412 )

413 res.dohistory(False)

414 # modify stokes type for polarization intensity image

415 return res

416

417def _immath_updateexpr(expr, varnames, subImages, filenames, file_map):

418 # Make sure no problems happened

419 if len(filenames) != len(subImages):

420 raise Exception(

421 'Unable to create subimages for all image names given'

422 )

423 # because real file names also have to be mapped to a corresponding subimage, CAS-1830

424 for k in file_map.keys():

425 # we require actual image names to be in quotes when used in the expression

426 varnames.extend(["'" + k + "'", '"' + k + '"'])

427 subImages.extend(2 * [file_map[k]])

428 # Put the subimage names into the expression

429 try:

430 expr = _immath_expr_from_varnames(expr, varnames, subImages)

431 except Exception as e:

432 casalog.post(

433 "Unable to construct pixel expression aborting immath: " + str(e),

434 'SEVERE'

435 )

436 raise

437 return (expr, varnames, subImages)

438

439def _immath_createsubimages(

440 box, chans, stokes, region, mask,

441 stretch, filenames, myia, tmpFilePrefix

442):

443 subImages = []

444 file_map = {}

445 i = 0

446 for image in filenames:

447 try:

448 myia.open(image)

449 myrg = regionmanager()

450 reg = myrg.frombcs(csys=myia.coordsys().torecord(),

451 shape=myia.shape(), box=box, chans=chans, stokes=stokes,

452 stokescontrol="a", region=region

453 )

454 tmpFile = tmpFilePrefix + str(i)

455 subim = myia.subimage(

456 region=reg, mask=mask, outfile=tmpFile, stretch=stretch

457 )

458 subim.done()

459 file_map[image] = tmpFile

460 subImages.append(tmpFile)

461 myia.done()

462 i = i + 1

463 except Exception as e:

464 raise Exception(

465 'Unable to apply region to file: ' + image

466 )

467 finally:

468 myia.done()

469 return (subImages, file_map)

470

471def _immath_dofull(

472 imagename, imagemd, outfile, mode, expr,

473 varnames, filenames, myia, prec

474):

475 expr = _immath_expr_from_varnames(expr, varnames, filenames)

476 return _immath_compute(

477 imagename, expr, outfile, imagemd, myia, prec

478 )

479

480def _immath_dospix(nfiles, varnames):

481 if nfiles != 2:

482 raise Exception('Requires two images at different frequencies')

483 return 'spectralindex(' + varnames[0] + ', ' + varnames[1] + ')'

484

485def _immath_filenames(filenames, tmpfilenames, varnames, mode):

486 ignoreimagename = False

487 if mode=='evalexpr':

488 varnamesSet = set(varnames)

489 count = 0

490 for imname in tmpfilenames:

491 # check if it is one of varnames, if not check the files in expr exist

492 if(not varnamesSet.issuperset(imname)):

493 if( not os.path.exists(imname)):

494 raise Exception('Image data set not found - please verify ' + imname)

495 else:

496 count = count + 1

497 if len(tmpfilenames) == count:

498 ignoreimagename = True

499 filenames = tmpfilenames

500 if not ignoreimagename:

501 for i in range(len(filenames)):

502 if not os.path.exists(filenames[i]):

503 casalog.post("Image data set not found - please verify " +filenames[i], "SEVERE")

504 raise Exception('Image data set not found - please verify '+filenames[i])

505 return filenames

506

507def _immath_varnames(varnames, filenames, tmpfilenames):

508 # Construct the variable name list. We append to the list the

509 # default variable names if the user hasn't supplied a full suite.

510 if not isinstance(varnames, list):

511 name0 = varnames

512 varnames = []

513 if name0:

514 varnames.append(name0)

515 nfile = max(len(filenames),len(tmpfilenames))

516 for i in range(len(varnames), nfile):

517 varnames.append('IM'+str(i))

518 casalog.post( 'Variable name list is: '+str(varnames), 'DEBUG1' )

519 return varnames

520

521def _immath_initial_cleanup(tmpFilePrefix):

522 try:

523 _immath_cleanup(tmpFilePrefix)

524 except Exception as e:

525 casalog.post( 'Unable to cleanup working directory '+os.getcwd()+'\n'+str(e), 'SEVERE' )

526 raise

527

528def _immath_check_outfile(outfile):

529 if not outfile:

530 outfile = 'immath_results.im'

531 casalog.post(

532 "The outfile parameter is empty, consequently the "

533 + "resultant image will be saved on disk in an image named "

534 + outfile, 'WARN'

535 )

536 if (os.path.exists(outfile)):

537 raise Exception(

538 'Output file '+ outfile

539 + ' exists. immath can not proceed, please '

540 + 'remove it or change the output file name.'

541 )

542 return outfile

543

544def _immath_cleanup(filePrefix):

545 # Remove any old tmp files that may be left from

546 # a previous run of immath

547 fileList=os.listdir( os.getcwd() )

548 for file in fileList:

549 if ( file.startswith( filePrefix ) ):

550 shutil.rmtree( file )

551

552def _immath_parse( expr='' ):

553 retValue=[]

554

555 # Find out if the names are surrounded by single or double quotes

556 quote=''

557 if ( expr.find('"') > -1 ):

558 quote='"'

559 if ( expr.find("'") > -1 ):

560 quote="'"

561

562 current=expr;

563 while( current.count(quote) > 1 ):

564 start = current.index(quote)+1

565 end = current[start:].index(quote)+start

566 if ( retValue.count( current[start:end] ) > 0 ) :

567 # We already have this file name so we won't add it

568 # to the list again. This saves us work and disk space.

569 current=current[end+1:]

570 continue;

571

572 retValue.append( current[start:end] )

573 current=current[end+1:]

574

575 return retValue

576

577# it is important to sort the varnames in reverse order before doing

578# the substitutions to assure the substitution set is performed correctly

579# CAS-1678

580def _immath_expr_from_varnames(expr, varnames, filenames):

581 tmpfiles = {}

582 for i in range(len(filenames)):

583 tmpfiles[varnames[i]] = filenames[i]

584 # python 3 requires explicit list conversion

585 tmpvars = list(tmpfiles.keys())

586

587 tmpvars.sort()

588 tmpvars.reverse()

589 for varname in tmpvars:

590 expr = expr.replace(varname, '"' + tmpfiles[varname] + '"')

591 return(expr)

592

593def _immath_translate_imagemd(imagename, imagemd):

594 # may IM0 etc is a real file name

595 if os.path.exists(imagemd):

596 return imagemd

597 # match IM0, IM1, ... etc

598 m = re.match("^IM(\d+)$", imagemd)

599 if not m:

600 return imagemd

601 idx = int(m.groups()[0])

602 if idx == 0 and type(imagename) == str:

603 return imagename

604 # if here, then imagename is an array of strings

605 if idx >= len(imagename):

606 # out of range

607 return imagemd

608 return imagename[idx]

609

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

288 statements