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
« prev ^ index » next coverage.py v7.6.4, created at 2024-11-01 07:19 +0000
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']
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]
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()
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
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]
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