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

1import os

2import sys

3import shutil

5from casatools import image, coordsys, quanta

6from casatasks import casalog

7from .ialib import write_image_history

10def imregrid(

11 imagename, template, output, asvelocity, axes, shape,

12 interpolation, decimate, replicate, overwrite

13):

14 _myia = None

15 outia = None

16 csys = None

17 try:

18 casalog.origin('imregrid')

19 if hasattr(template, 'lower') and not template.lower() == "get":

20 # First check to see if the output file exists. If it

21 # does then we abort. CASA doesn't allow files to be

22 # over-written, just a policy.

23 if len(output) < 1:

24 output = imagename + '.regridded'

25 casalog.post("output was not specified - defaulting to\n\t"

26 + output, 'INFO')

27 _myia = image()

28 _myia.dohistory(False)

29 # Figure out what the user wants.

30 if not isinstance(template, dict):

31 if template.lower() == 'get':

32 _myia.open(imagename)

33 csys = _myia.coordsys()

34 shape = _myia.shape()

35 _myia.done()

36 return {'csys': csys.torecord(), 'shap': shape}

37 elif template.upper() in (

38 'J2000', 'B1950', 'B1950_VLA', 'GALACTIC',

39 'HADEC', 'AZEL', 'AZELSW', 'AZELNE', 'ECLIPTIC',

40 'MECLIPTIC', 'TECLIPTIC', 'SUPERGAL'

41 ):

42 outia = _imregrid_to_new_ref_frame(

43 _myia, imagename, template, output, axes,

44 shape, overwrite, interpolation, decimate

45 )

46 try:

47 param_names = imregrid.__code__.co_varnames[:imregrid.__code__.co_argcount]

48 vars = locals( )

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

51 write_image_history(

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

53 param_names, param_vals, casalog

54 )

55 except Exception as instance:

56 casalog.post(

57 "*** Error \'%s\' updating HISTORY" % (instance), 'WARN'

58 )

59 outia.done()

60 return

61 else:

62 if (

63 not os.path.isdir(template)

64 or not os.access(template, os.R_OK)

65 ):

66 raise TypeError('Cannot read template image %s' % template)

67 template_ia = image()

68 template_ia.open(template)

69 template_csys = template_ia.coordsys()

70 image_ia = image()

71 image_ia.open(imagename)

72 image_csys = image_ia.coordsys()

73 tempshape = template_ia.shape()

74 imshape = image_ia.shape()

75 axestoregrid = axes

76 if (axes[0] < 0):

77 # default value of axes, need to determine actual axes to

78 # send to ia.regrid()

79 axestoregrid = []

80 image_ncoords = image_csys.ncoordinates()

81 for i in range(image_ncoords):

82 ctype = image_csys.coordinatetype(i)[0]

83 template_coord = template_csys.findcoordinate(ctype)

84 if ctype != 'Stokes' and template_coord["return"]:

85 # only regrid if not Stokes axis and coordinate

86 # exists in template

87 for template_pix_axis in template_coord['pixel']:

88 if tempshape[template_pix_axis] > 1:

89 # only regrid if template axis is not

90 # degenerate

91 world_axes = image_csys.findcoordinate(

92 ctype

93 )['pixel']

94 for world_pix_axis in world_axes:

95 if imshape[world_pix_axis] > 1:

96 # only regrid if the world axis is

97 # not degenerate

98 axestoregrid.append(world_pix_axis)

99 # eliminate dups

100 axestoregrid = list(set(axestoregrid))

101 if len(axestoregrid) == 0:

102 raise RuntimeError("Found no axes to regrid!")

103 axestoregrid.sort()

104 if (len(shape) == 1 and shape[0] == -1):

105 shape = _imregrid_handle_default_shape(

106 imshape, image_csys, template_csys,

107 axestoregrid, tempshape, axes

108 )

109 template_ia.done()

110 image_ia.done()

111 csys = template_csys

112 else:

113 # csys and shape specified in dictionary generated by previous run

114 # with template="get"

115 csys = coordsys()

116 csys.fromrecord(template['csys'])

117 shape = template['shap']

118

119 # The actual regridding.

120 _myia.open(imagename)

121 outia = _myia.regrid(

122 outfile=output, shape=shape, csys=csys.torecord(),

123 axes=axes, asvelocity=asvelocity,

124 method=interpolation, decimate=decimate,

125 replicate=replicate, overwrite=overwrite

126 )

127 try:

128 param_names = imregrid.__code__.co_varnames[:imregrid.__code__.co_argcount]

129 vars = locals( )

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

131

132 write_image_history(

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

134 param_names, param_vals, casalog

135 )

136 except Exception as instance:

137 casalog.post(

138 "*** Error \'%s\' updating HISTORY" % (instance), 'WARN'

139 )

140 return

141 finally:

142 if _myia:

143 _myia.done()

144 if outia:

145 outia.done()

146 if csys:

147 csys.done()

148

149def _imregrid_to_new_ref_frame(

150 _myia, imagename, template, output, axes,

151 shape, overwrite, interpolation, decimate

152):

153 _myia.open(imagename)

154 csys = _myia.coordsys()

155 if len(shape) > 0 and shape != [-1]:

156 casalog.post(

157 "Specified shape parameter will be ignored when regridding to a "

158 + "new reference frame", "WARN"

159 )

160 if len(axes) > 0 and axes != [-1]:

161 casalog.post(

162 "Specified axes parameter will be ignored when "

163 + "regridding to a new reference frame",

164 "WARN"

165 )

166 dirinfo = csys.findcoordinate("direction")

167 if not dirinfo['return']:

168 raise Exception("Image does not have a direction coordinate.")

169 newrefcode = template.upper()

170 oldrefcode = csys.referencecode("direction")[0]

171 if oldrefcode == newrefcode:

172 casalog.post(

173 imagename + ' is already in ' + oldrefcode,

174 'INFO'

175 )

176 casalog.post("...making a straight copy...", 'INFO')

177 subi = _myia.subimage(output)

178 _myia.done()

179 csys.done()

180 return subi

181 if (csys.projection()['type'] == 'SFL'):

182 raise Exception(

183 "The direction coordinate of this image has a projection "

184 "of SFL. Because of the constraints of this projection, "

185 "this image cannot be easily rotated. You may wish to "

186 "consider temporarily modifying the projection using "

187 "cs.setprojection() to allow rotation of the image."

188 )

189 casalog.post(

190 "Changing coordinate system from " + oldrefcode

191 + " to " + newrefcode, 'INFO'

192 )

193 diraxes = dirinfo['pixel']

194 if len(diraxes) != 2:

195 raise Exception(

196 "Unsupported number of direction axes. There must be exactly 2."

197 )

198 dirrefpix = csys.referencepixel("direction")["numeric"]

199 shape = _myia.shape()

200 centerpix = [int(shape[diraxes[0]]/2), int(shape[diraxes[1]]/2)]

201 if centerpix[0] != dirrefpix[0] or centerpix[1] != dirrefpix[1]:

202 casalog.post(

203 "Center direction pixel and reference pixel are "

204 + "different, making a temporary image and setting "

205 + "the reference pixel equal to the center pixel. "

206 + "The output image will have this modified coordinate system."

207 )

208 # so toworld() works in the correct frame

209 csys.setconversiontype(oldrefcode)

210 newrefpix = csys.referencepixel()["numeric"]

211 newrefpix[diraxes[0]] = centerpix[0]

212 newrefpix[diraxes[1]] = centerpix[1]

213 newrefval = csys.toworld(newrefpix)["numeric"]

214 csys.setreferencepixel(newrefpix)

215 csys.setreferencevalue(newrefval)

216 tsub = _myia.subimage()

217 _myia.done()

218 _myia = tsub

219 _myia.dohistory(False)

220 _myia.setcoordsys(csys.torecord())

221 doref = (

222 csys.referencecode("direction")[0] == csys.conversiontype("direction")

223 )

224 angle = csys.convertdirection(newrefcode)

225 myqa = quanta()

226 mysin = myqa.getvalue(myqa.sin(angle))

227 mycos = myqa.getvalue(myqa.cos(angle))

228 xnew = 0

229 ynew = 0

230 for xx in [-centerpix[0], centerpix[0]]:

231 for yy in [-centerpix[1], centerpix[1]]:

232 xnew = max(xnew, abs(xx*mycos - yy*mysin + 1))

233 ynew = max(ynew, abs(xx*mysin + yy*mycos + 1))

234 pad = int(max(xnew - shape[0]/2, ynew - shape[1]/2))

235 if pad > 0:

236 casalog.post(

237 "Padding image by " + str(pad)

238 + " pixels so no pixels are cut off in the regridding",

239 "NORMAL"

240 )

241 _myia = _myia.pad("", pad, wantreturn=True)

242 _myia.dohistory(False)

243 shape = _myia.shape()

244 newrefpix = csys.referencepixel()['numeric']

245 newrefpix[diraxes[0]] = newrefpix[diraxes[0]] + pad

246 newrefpix[diraxes[1]] = newrefpix[diraxes[1]] + pad

247 csys.setreferencepixel(newrefpix)

248 regridded = _myia.regrid(

249 outfile="",shape=shape, csys=csys.torecord(), axes=diraxes,

250 method=interpolation, decimate=decimate ,doref=doref

251 )

252 regridded.dohistory(False)

253 # beam is rotated counterclockwise

254 angle_in_deg = format("%7.3f" % myqa.convert(angle, "deg")['value'])

255 casalog.post(

256 "Will rotate beams counterclockwise by " + angle_in_deg + " degrees, "

257 + "if necessary, to account for angle between original and new frame "

258 + "at the reference pixel", 'NORMAL'

259 )

260 regridded.rotatebeam(angle=myqa.mul(-1, angle))

261 # now crop

262 casalog.post("Cropping masked image boundaries", "NORMAL")

263 cropped = regridded.crop(outfile=output, axes=diraxes, overwrite=overwrite)

264 regridded.done()

265 _myia.done()

266 return cropped

267

268def _imregrid_handle_default_shape(

269 imshape, image_csys, template_csys,

270 axestoregrid, tempshape, original_axes

271):

272 # CAS-4959, output shape should have template shape for axes being

273 # regridded, input image shape for axes not being regridded, CAS-4960 in

274 # cases where the input image and template both have multiple stokes, the

275 # number of pixels on the output stokes axis is to be the number of stokes

276 # the input and template have in common

277 shape = imshape

278 targetaxesnames = image_csys.names()

279 template_spectral_info = template_csys.findcoordinate("Spectral")

280 template_has_spectral = template_spectral_info['return']

281 if template_has_spectral:

282 template_spectral_axis = template_spectral_info['pixel'][0]

283 atr = axestoregrid[:]

284 for i in range(len(imshape)):

285 atr_count = 0

286 for j in atr:

287 if i == j:

288 # axis numbers may not correspond so have to look for the

289 # template axis location by the axis name, CAS-4960

290 template_axis = template_csys.findaxisbyname(targetaxesnames[i])

291 template_axis_length = tempshape[template_axis]

292 if (

293 template_has_spectral

294 and template_spectral_axis == template_axis

295 and template_axis_length == 1

296 ):

297 casalog.post(

298 "You've specified that you want to regrid the spectral axis without specifying "

299 "the output shape. Normally the length chosen would be that of the corresponding "

300 "template axis, however, the template spectral axis is degenerate and one cannot regrid "

301 "an axis such that its output length is one. So, removing axis " + str(j)

302 + " from the axes list and just copying the input spectral information to the output image",

303 "WARN"

304 )

305 shape[i] = imshape[i]

306 del axestoregrid[atr_count]

307 else:

308 shape[i] = tempshape[template_axis]

309

310 break

311 atr_count += 1;

312 if (

313 template_csys.findcoordinate("stokes")['return']

314 and image_csys.findcoordinate("stokes")['return']

315 and len(template_csys.stokes()) > 1

316 and len(image_csys.stokes()) > 1

317 ):

318 stokes_axis = image_csys.findaxisbyname("stokes")

319 found = (

320 len(original_axes) == 0

321 or (len(original_axes) == 1 and original_axes[0] < 0)

322 )

323 if not found:

324 for axis in axestoregrid:

325 if axis == stokes_axis:

326 found = True

327 break

328 if found:

329 # adjust stokes length to be equal to number of common stokes

330 common_stokes_count = 0

331 for image_stokes in image_csys.stokes():

332 for template_stokes in template_csys.stokes():

333 if image_stokes == template_stokes:

334 common_stokes_count += 1

335 break

336 shape[stokes_axis] = common_stokes_count

337 else:

338 shape[stokes_axis] = imshape[stokes_axis]

339 return shape

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

200 statements