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

200 statements  

« prev     ^ index     » next       coverage.py v7.6.4, created at 2024-11-01 07:19 +0000

1import os 

2import sys 

3import shutil 

4 

5from casatools import image, coordsys, quanta 

6from casatasks import casalog 

7from .ialib import write_image_history 

8 

9 

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] 

50 

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