Coverage for /wheeldirectory/casa-6.7.0-12-py3.10.el8/lib/py/lib/python3.10/site-packages/casatasks/private/task_fixplanets.py: 4%
333 statements
« prev ^ index » next coverage.py v7.6.4, created at 2024-10-31 17:39 +0000
1import shutil
2import os
3import string
5from .parallel.parallel_task_helper import ParallelTaskHelper
6from casatools import ms as mstool
7from casatools import table as tbtool
8from casatools import imager as imtool
9from casatools import measures, quanta
10from casatasks import casalog
11from .mstools import write_history
12_qa = quanta( )
13_me = measures( )
15def _checkinternalephemtab(vis, field):
16 """
17 This function checks if there is an ephemeris table attached in the MS for
18 the field selected. It will returned table names found under FIELD table.
20 """
21 import glob
22 from casatools import table, ms
23 _tb = table()
24 _ms = ms()
25 fids = _ms.msseltoindex(vis=vis,field=field)['field']
26 _tb.open(vis+'/FIELD')
27 ephemnames = []
28 if 'EPHEMERIS_ID' in _tb.colnames():
29 for i in fids:
30 ephemid = _tb.getcell('EPHEMERIS_ID',i)
31 if ephemid != -1:
32 ephemnames.append(glob.glob(f'{vis}/FIELD/EPHEM{ephemid}*/')[0])
33 _tb.close()
34 return list(set(ephemnames))
36def fixplanets(vis, field, fixuvw=False, direction='', refant=0, reftime='first'):
37 """
38 Fix FIELD, SOURCE, and UVW for given fields based on given direction or pointing
39 table information
41 This task's main purpose is to correct observations which were performed
42 with correct pointing and correlation but for which incorrect direction
43 information was entered in the FIELD and SOURCE table of the MS.
44 If you actually want to change the phase center of the visibilties in an MS,
45 you should use task fixvis.
47 Input Parameters
48 vis -- Name of the input visibility set
50 field -- field selection string
52 fixuvw -- recalc uvw coordinates? (default: False)
54 direction -- if set, don't use pointing table but set direction to this value.
55 The direction can either be given explicitly or as the path
56 to a JPL Horizons ephemeris (for an example of the format,
57 ephemerides/JPL-Horizons/ in the local data directory or "External Data" section
58 of the current CASADocs).
59 Alternatively, the ephemeris table can also be obtained using
60 getephemtable task.
62 example: 'J2000 19h30m00 -40d00m00', default= '' (use pointing table)
64 refant -- if using pointing table information, use it from this antenna
65 default: 0 (antenna id 0)
66 examples: 'DV06' (antenna with name DV06)
67 3 (antenna id 3)
69 reftime -- if using pointing table information, use it from this timestamp
70 default: 'first'
71 examples: 'median' will use the median timestamp for the given field
72 using only the unflagged maintable rows
73 '2012/07/11/08:41:32' will use the given timestamp (must be
74 within the observaton time)
75 """
77 casalog.origin('fixplanets')
79 mst = mstool()
80 tbt = tbtool()
81 imt = None
83 try:
84 fields = mst.msseltoindex(vis=vis,field=field)['field']
85 numfields = 0
86 if(len(fields) == 0):
87 casalog.post( "Field selection returned zero results.", 'WARN')
88 return
90 tbt.open(vis+"/FIELD")
91 oldrefcol = []
92 if('PhaseDir_Ref' in tbt.colnames()):
93 oldrefcol = tbt.getcol('PhaseDir_Ref')
94 tbt.close()
96 fixplanetstemp = 'fixplanetstemp-'+os.path.basename(vis) # temp file name
97 fixplanetstemp2 = 'fixplanetstemp2-'+os.path.basename(vis) # temp file name
99 for fld in fields:
100 thenewra_rad = 0.
101 thenewdec_rad = 0.
102 thenewref = -1
103 thenewrefstr = ''
104 theephemeris = ''
106 if(direction==''): # use information from the pointing table
107 # find median timestamp for this field in the main table
108 shutil.rmtree(fixplanetstemp, ignore_errors=True)
109 thetime = 0
110 if(reftime.lower()=='median'):
111 tbt.open(vis)
112 tttb = tbt.query('FIELD_ID=='+str(fld)+' AND FLAG_ROW==False',
113 name=fixplanetstemp, columns='TIME')
114 tttb.close()
115 tttb = None
116 tbt.close()
117 tbt.open(fixplanetstemp)
118 if(tbt.nrows()>0):
119 thetime = tbt.getcell('TIME',tbt.nrows()//2)
120 casalog.post( "MEDIAN TIME "+str(thetime), 'NORMAL')
121 tbt.close()
122 else:
123 casalog.post( "No pointing table rows for field "+field, 'NORMAL')
124 tbt.close()
125 shutil.rmtree(fixplanetstemp, ignore_errors=True)
126 return
127 elif(reftime.lower()=='first'):
128 tbt.open(vis)
129 tttb = tbt.query('FIELD_ID=='+str(fld), name=fixplanetstemp, columns='TIME')
130 tttb.close()
131 tttb = None
132 tbt.close()
133 tbt.open(fixplanetstemp)
134 if(tbt.nrows()>0):
135 thetime = tbt.getcell('TIME',0)
136 casalog.post( "FIRST TIME "+str(thetime), 'NORMAL')
137 tbt.close()
138 else:
139 casalog.post( "No pointing table rows for field "+field, 'NORMAL')
140 tbt.close()
141 shutil.rmtree(fixplanetstemp, ignore_errors=True)
142 return
143 else:
144 try:
145 myqa = _qa.quantity(reftime)
146 thetime = _qa.convert(myqa,'s')['value']
147 except Exception as instance:
148 raise TypeError("reftime parameter is not a valid date (e.g. YYYY/MM/DD/hh:mm:ss): %s" % reftime)
149 tbt.open(vis)
150 tttb = tbt.query('FIELD_ID=='+str(fld), name=fixplanetstemp, columns='TIME')
151 tttb.close()
152 tttb = None
153 tbt.close()
154 tbt.open(fixplanetstemp)
155 if(tbt.nrows()>0):
156 thefirsttime = tbt.getcell('TIME',0)
157 thelasttime = tbt.getcell('TIME',tbt.nrows()-1)
158 tbt.close()
159 else:
160 casalog.post( "No pointing table rows for field "+field, 'NORMAL')
161 tbt.close()
162 shutil.rmtree(fixplanetstemp, ignore_errors=True)
163 return
164 shutil.rmtree(fixplanetstemp, ignore_errors=True)
165 if (thefirsttime<=thetime and thetime<=thelasttime):
166 casalog.post( "GIVEN TIME "+reftime+" == "+str(thetime), 'NORMAL')
167 else:
168 casalog.post( "GIVEN TIME "+reftime+" == "+str(thetime)+" is not within the observation time ("
169 +str(thefirsttime)+'s'+" - "
170 +str(thelasttime)+'s'+")", 'SEVERE')
171 raise TypeError
172 shutil.rmtree(fixplanetstemp, ignore_errors=True)
174 # determine reference antenna
175 antids = mst.msseltoindex(vis=vis,baseline=refant)['antenna1']
176 antid = -1
177 if(len(antids) == 0):
178 msg = "Antenna selection returned zero results."
179 raise RuntimeError(msg)
180 antid = int(antids[0])
181 casalog.post('Using antenna id '+str(antid)+' as reference antenna.', 'NORMAL')
182 tbt.open(vis+'/ANTENNA')
183 flgcol = tbt.getcol('FLAG_ROW')
184 tbt.close()
185 if(flgcol[antid]==True):
186 msg = 'Antenna id '+str(antid)+' is flagged. Please select a different one.'
187 raise RuntimeError(msg)
189 # get direction for the timestamp from pointing table
190 shutil.rmtree(fixplanetstemp2, ignore_errors=True)
191 tbt.open(vis+'/POINTING')
192 ttb = tbt.query('TRACKING==True AND NEARABS(TIME,'+str(thetime)+',INTERVAL/2.) AND ANTENNA_ID=='
193 +str(antid),
194 name=fixplanetstemp2)
195 nr = ttb.nrows()
196 ttb.close()
197 ttb = None
198 if(nr==0):
199 shutil.rmtree(fixplanetstemp2, ignore_errors=True)
200 ttb2 = tbt.query('TRACKING==True AND NEARABS(TIME,'+str(thetime)+',3.) AND ANTENNA_ID=='
201 +str(antid), # search within 3 seconds
202 name=fixplanetstemp2)
203 nr = ttb2.nrows()
204 ttb2.close()
205 ttb2 = None
206 if(nr==0):
207 casalog.post( "Cannot find any POINTING table rows for antenna "+str(antid)
208 +" with TRACKING==True within 3 seconds of TIME "+str(thetime), 'NORMAL')
209 casalog.post( "Will try without requiring TRACKING==True ...", 'NORMAL')
210 shutil.rmtree(fixplanetstemp2, ignore_errors=True)
211 ttb3 = tbt.query('NEARABS(TIME,'+str(thetime)+',INTERVAL/2.) AND ANTENNA_ID=='
212 +str(antid),
213 name=fixplanetstemp2)
214 nr = ttb3.nrows()
215 ttb3.close()
216 ttb3 = None
217 if(nr==0):
218 shutil.rmtree(fixplanetstemp2, ignore_errors=True)
219 ttb4 = tbt.query('NEARABS(TIME,'+str(thetime)+',3.) AND ANTENNA_ID=='
220 +str(antid), # search within 3 seconds
221 name=fixplanetstemp2)
222 nr = ttb4.nrows()
223 ttb4.close()
224 ttb4 = None
225 if(nr==0):
226 tbt.close()
227 msg = ("Cannot find any POINTING table rows for antenna "+str(antid) +
228 " within 3 seconds of TIME "+str(thetime))
229 raise RuntimeError(msg) # give up
230 tbt.close()
231 tbt.open(fixplanetstemp2)
232 thedir = tbt.getcell('DIRECTION',0)
233 tbt.close()
234 casalog.post( ' field id '+str(fld)+ ' AZ EL '+str(thedir[0])+" "+str(thedir[1]), 'NORMAL')
235 thedirme = _me.direction(rf='AZELGEO',v0=_qa.quantity(thedir[0][0], 'rad'),
236 v1=_qa.quantity(thedir[1][0],'rad'))
237 # convert AZEL to J2000 coordinates
238 _me.doframe(_me.epoch(rf='UTC', v0=_qa.quantity(thetime,'s')))
240 tbt.open(vis+'/ANTENNA')
241 thepos = tbt.getcell('POSITION',antid)
242 theposref = tbt.getcolkeyword('POSITION', 'MEASINFO')['Ref']
243 theposunits = tbt.getcolkeyword('POSITION', 'QuantumUnits')
244 tbt.close()
245 casalog.post( "Ref. antenna position is "+str(thepos)+' ('+theposunits[0]+', '+theposunits[1]+', '
246 +theposunits[2]+')('+theposref+')', 'NORMAL')
247 _me.doframe(_me.position(theposref,
248 _qa.quantity(thepos[0],theposunits[0]),
249 _qa.quantity(thepos[1],theposunits[1]),
250 _qa.quantity(thepos[2],theposunits[2]))
251 )
252 thedirmemod = _me.measure(v=thedirme, rf='J2000')
253 # casalog.post(thedirmemod)
254 thenewra_rad = thedirmemod['m0']['value']
255 thenewdec_rad = thedirmemod['m1']['value']
256 _me.done()
257 shutil.rmtree(fixplanetstemp2, ignore_errors=True)
259 else: # direction is not an empty string, use this instead of the pointing table information
260 if(type(direction)==str):
261 dirstr = direction.split(' ')
262 if len(dirstr)==1: # an ephemeris table was given
263 if(os.path.exists(dirstr[0])):
264 if os.path.isfile(dirstr[0]): # it is a file, i.e. not a CASA table
265 msg = "*** Error when interpreting parameter \'direction\':\n A file is given. Use of the JPL-Horizons "+\
266 "MIME format file is deprecated."
267 raise RuntimeError(msg)
268 else: # not a file, assume it is a CASA table
269 theephemeris = dirstr[0]
270 # add a check if it is going to replace the existing table in the MS
271 existingephemtab = _checkinternalephemtab(vis, field)
272 if existingephemtab != []:
273 casalog.post(f'Will replace existing ephemeris table {existingephemtab} in the MS with {direction}. '+\
274 'Ephemeris tables attached in the MS are assumed to be the ones used by the correlator and attaching' +\
275 ' a different ephemeris table may lead to scientifically wrong results.','WARN')
277 casalog.post('Will use ephemeris table '+theephemeris+' with offset (0,0)', 'NORMAL')
279 thenewra_rad = 0.
280 thenewdec_rad = 0.
282 else:
283 msg = "*** Error when interpreting parameter \'direction\':\n string is neither a direction nor a table."
284 raise RuntimeError(msg)
285 else:
286 if len(dirstr)==2: # a direction without ref frame was given
287 dirstr = ['J2000', dirstr[0], dirstr[1]]
288 casalog.post('Assuming reference frame J2000 for parameter \'direction\'', 'NORMAL')
289 elif not len(dirstr)==3:
290 msg = 'Incorrect format in parameter \'direction\''
291 raise RuntimeError(msg)
293 # process direction and refframe
294 try:
295 thedir = _me.direction(dirstr[0], dirstr[1], dirstr[2])
296 thenewra_rad = thedir['m0']['value']
297 thenewdec_rad = thedir['m1']['value']
298 except Exception as instance:
299 msg = "*** Error \'%s\' when interpreting parameter \'direction\': " % (instance)
300 raise RuntimeError(msg)
302 try:
303 tbt.open(vis+"/FIELD")
304 thenewrefindex = tbt.getcolkeyword('PHASE_DIR', 'MEASINFO')['TabRefTypes'].tolist().index(dirstr[0])
305 thenewref = tbt.getcolkeyword('PHASE_DIR', 'MEASINFO')['TabRefCodes'][thenewrefindex]
306 thenewrefstr = dirstr[0]
307 tbt.close()
308 except Exception as instance:
309 casalog.post("PHASE_DIR is not a variable reference column. Will leave reference as is.", 'WARN')
310 thenewref = -1
311 if(0<thenewref and thenewref<32):
312 msg = "*** Error when interpreting parameter \'direction\':\n presently only J2000 and solar system objects are supported."
313 raise RuntimeError(msg)
314 #endif
316 # modify FIELD table
317 tbt.open(vis+'/FIELD', nomodify=False)
318 numfields = tbt.nrows()
319 theolddir = tbt.getcell('PHASE_DIR',fld)
320 planetname = tbt.getcell('NAME',fld)
322 casalog.post( 'object: '+planetname, 'NORMAL')
323 casalog.post( 'old RA, DEC (rad) '+str(theolddir[0])+" " +str(theolddir[1]), 'NORMAL')
324 if(theephemeris==''):
325 casalog.post( 'new RA, DEC (rad) '+str(thenewra_rad)+" "+ str(thenewdec_rad), 'NORMAL')
326 else:
327 casalog.post( 'new RA, DEC to be taken from ephemeris table '+theephemeris
328 +' with offsets (rad) '+str(thenewra_rad)+" "+ str(thenewdec_rad), 'NORMAL')
330 pcol = tbt.getcol('PHASE_DIR')
331 pcol[0][0][fld] = thenewra_rad
332 pcol[1][0][fld] = thenewdec_rad
333 tbt.putcol('PHASE_DIR',pcol)
335 pcol = tbt.getcol('DELAY_DIR')
336 pcol[0][0][fld] = thenewra_rad
337 pcol[1][0][fld] = thenewdec_rad
338 tbt.putcol('DELAY_DIR',pcol)
340 pcol = tbt.getcol('REFERENCE_DIR')
341 pcol[0][0][fld] = thenewra_rad
342 pcol[1][0][fld] = thenewdec_rad
343 tbt.putcol('REFERENCE_DIR',pcol)
345 casalog.post("FIELD table PHASE_DIR, DELAY_DIR, and REFERENCE_DIR columns changed for field "+str(fld)+".", 'NORMAL')
347 if(theephemeris==''):
349 if ('EPHEMERIS_ID' in tbt.colnames()) and (tbt.getcell('EPHEMERIS_ID',fld)>=0): # originally an ephemeris was used
350 eidc = tbt.getcol('EPHEMERIS_ID')
351 eidc[fld] = -1
352 tbt.putcol('EPHEMERIS_ID', eidc) # remove the reference to it
353 casalog.post("FIELD table EPHEMERIS_ID column reset to -1 for field "+str(fld)+".", 'NORMAL')
355 if(thenewref!=-1):
356 # modify reference of the direction columns permanently
357 theoldref = -1
358 theoldref2 = -1
359 theoldref3 = -1
360 try:
361 pcol = tbt.getcol('PhaseDir_Ref')
362 theoldref = pcol[fld]
363 pcol[fld] = thenewref
364 oldrefcol = pcol
366 pcol2 = tbt.getcol('DelayDir_Ref')
367 theoldref2 = pcol2[fld]
368 pcol2[fld] = thenewref
370 pcol3 = tbt.getcol('RefDir_Ref')
371 theoldref3 = pcol3[fld]
372 pcol3[fld] = thenewref
374 if not (theoldref==theoldref2 and theoldref==theoldref3):
375 casalog.post("The three FIELD table direction reference frame entries for field "+str(fld)
376 +" are not identical in the input data: "
377 +str(theoldref)+", "+str(theoldref2)+", "+str(theoldref3)
378 +". Will try to continue ...", 'WARN')
379 theoldref=-1
380 else:
381 oldrefstrindex = tbt.getcolkeyword('PHASE_DIR', 'MEASINFO')['TabRefCodes'].tolist().index(theoldref)
382 oldrefstr = tbt.getcolkeyword('PHASE_DIR', 'MEASINFO')['TabRefTypes'][oldrefstrindex]
383 casalog.post("Original FIELD table direction reference frame entries for field "+str(fld)
384 +" are all "+str(theoldref)+" ("+oldrefstr+")", 'NORMAL')
385 tbt.putcol('PhaseDir_Ref', pcol)
386 tbt.putcol('DelayDir_Ref', pcol2)
387 tbt.putcol('RefDir_Ref', pcol3)
389 except Exception as instance:
390 msg = "*** Error \'%s\' when writing reference frames in FIELD table: " % (instance)
391 raise RuntimeError(msg)
393 if(theoldref != thenewref):
394 casalog.post("FIELD table direction reference frame entries for field "+str(fld)
395 +" set to "+str(thenewref)+" ("+thenewrefstr+")", 'NORMAL')
396 else:
397 casalog.post("FIELD table direction reference frame entries for field "+str(fld)
398 +" unchanged.", 'NORMAL')
399 else: # we are working with an ephemeris
400 try:
401 mst.open(vis, nomodify=False)
402 # sanitize the internal ephemeris table name
403 valid_chars = "-+_.%s%s" % (string.ascii_letters, string.digits)
404 planetname_for_table = ''.join(c if c in valid_chars else '_' for c in planetname)
405 # make sure theemepheris is given as absolute path to avoid MeasIERS warnings
406 mst.addephemeris(-1, os.path.abspath(theephemeris), planetname_for_table, fld) # -1 = take the next free ID
407 except Exception as instance:
408 raise RuntimeError("*** Error \'%s\' when attaching ephemeris: " %
409 (instance))
410 finally:
411 mst.close()
414 if(fixuvw and (len(oldrefcol)!=0) and (thenewref>0)):
415 # modify reference of phase dir for fixuvw
416 pcol = tbt.getcol('PhaseDir_Ref')
417 pcol[fld] = 0 # J2000
418 tbt.putcol('PhaseDir_Ref', pcol) # J2000
420 tbt.close()
422 #modify SOURCE table
423 newsra_rad = thenewra_rad
424 newsdec_rad = thenewdec_rad
425 if(theephemeris!=''): # get the nominal position from the ephemeris
426 mst.open(vis)
427 trec = mst.getfielddirmeas('PHASE_DIR',fld)
428 newsra_rad = trec['m0']['value']
429 newsdec_rad = trec['m1']['value']
430 mst.close()
431 tbt.open(vis+'/SOURCE', nomodify=False)
432 sdir = tbt.getcol('DIRECTION')
433 newsdir = sdir
434 sname = tbt.getcol('NAME')
436 for i in range(0,tbt.nrows()):
437 if(sname[i]==planetname):
438 #casalog.post('i old dir ' + i + " " + sdir[0][i] + sdir[1][i])
439 newsdir[0][i] = newsra_rad
440 newsdir[1][i] = newsdec_rad
441 #casalog.post(' new dir ' + newsdir[0][i] + newsdir[1][i])
442 tbt.putcol('DIRECTION', newsdir)
443 tbt.close()
444 casalog.post("SOURCE table DIRECTION column changed.", 'NORMAL')
446 # end for
448 if(fixuvw):
449 casalog.post("Fixing the UVW coordinates ...", 'NORMAL')
451 # similar to fixvis
453 fldids = []
454 for i in range(numfields):
455 if (i in fields):
456 fldids.append(i)
458 imt = imtool()
459 imt.open(vis, usescratch=False)
460 imt.calcuvw(fldids, refcode='J2000', reuse=False)
461 imt.close()
462 imt = None
464 if((len(oldrefcol)!=0) and (thenewref>0)):
465 tbt.open(vis+'/FIELD', nomodify=False)
466 tbt.putcol('PhaseDir_Ref', oldrefcol)
467 tbt.close()
469 else:
470 casalog.post("UVW coordinates not changed.", 'NORMAL')
472 if (ParallelTaskHelper.isParallelMS(vis)):
473 casalog.post("Tidying up the MMS subtables ...", 'NORMAL')
474 ParallelTaskHelper.restoreSubtableAgreement(vis)
476 finally:
477 mst = None
478 tbt = None
479 imt = None
481 # Write history to MS
482 try:
483 param_names = fixplanets.__code__.co_varnames[:fixplanets.__code__.co_argcount]
484 vars = locals()
485 param_vals = [vars[p] for p in param_names]
486 write_history(mstool(), vis, 'fixplanets', param_names,
487 param_vals, casalog)
488 except Exception as instance:
489 casalog.post("*** Error \'%s\' updating HISTORY" % (instance),
490 'WARN')
492 mst = None
493 tbt = None