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

1import glob 

2import os 

3 

4from casatools import table, measures, quanta, ms 

5from casatasks import casalog 

6 

7 

8# Conversion of TOPO ephemerides to GEO (ICRS) 

9# 

10# Example: 

11# import recipes.ephemerides.convertephem as ce 

12# ce.convert2geo('titan.ms', 'Titan') 

13# 

14# will find the ephemeris attached to field "Titan" in the MS "titan.ms" 

15# and convert RA, Dec, and RadVel to the geocentric ref frame (ICRS) 

16 

17 

18def converttopoephem2geo(tablename='', outtablename='', overwrite=True): 

19 """ 

20 converttopoephem2geo 

21 

22 Convert the given topo ephemeris table to the geocentric ref frame. 

23 Converted are the RA, DEC, and RadVel columns only 

24 

25 tablename -- name of the TOPO frame ephemeris table (in canonical CASA format) 

26 

27 outtablename -- name of the output GEO frame ephemeris table 

28 

29 """ 

30 

31 if type(tablename)!=str or tablename=='': 

32 casalog.post('Invalid parameter tablename', 'WARN') 

33 return False 

34 

35 if type(outtablename)!=str or outtablename=='': 

36 casalog.post('Invalid parameter outtablename', 'WARN') 

37 return False 

38 

39 #convert RA, DEC, and RadVel 

40 tbt = table() 

41 met = measures() 

42 qat = quanta() 

43 

44 tbt.open(tablename) 

45 ra = tbt.getcol('RA') 

46 dec = tbt.getcol('DEC') 

47 radvel = tbt.getcol('RadVel') 

48 radvelunit = 'km/s' 

49 tmpkw = tbt.getcolkeywords('RadVel') 

50 if 'UNIT' in tmpkw: 

51 radvelunit = tmpkw['UNIT'] 

52 elif 'QuantumUnits' in tmpkw: 

53 radvelunit = tmpkw['QuantumUnits'][0] 

54 else: 

55 casalog.post('Cannot determine units of radial velocity column. Assuming km/s.', 'WARN') 

56 mjd = tbt.getcol('MJD') 

57 kw = tbt.getkeywords() 

58 tbt.close() 

59 

60 geodist = kw['GeoDist'] # (km above reference ellipsoid) 

61 geolat = kw['GeoLat'] # (deg) 

62 geolong = kw['GeoLong'] # (deg) 

63 

64 if 'obsloc' in kw: 

65 obsloc = kw['obsloc'] 

66 else: 

67 casalog.post('Ephemeris does not have the obsloc keyword.', 'INFO') 

68 if (geodist==geolat==geolong==0.): 

69 casalog.post(' Assuming obsloc == GEOCENTRIC since lat, long, and dist are zero.', 'INFO') 

70 obsloc='GEOCENTRIC' 

71 else: 

72 obsloc='unknown' 

73 

74 oldref = 'J2000' 

75 newref = 'ICRS' 

76 

77 if 'posrefsys' in kw: 

78 posref = kw['posrefsys'] 

79 else: 

80 casalog.post('Ephemeris does not have the posrefsys keyword. Assuming ICRF/J2000.0', 'WARN') 

81 posref = 'ICRF/J2000.0' 

82 

83 if not ('J2000' in posref): 

84 if 'ICRS' in posref: 

85 oldref = 'ICRS' 

86 else: 

87 casalog.post('Position reference is '+posref+' is not supported, yet.', 'WARN') 

88 return False 

89 

90 newposref = 'ICRF/'+newref 

91 

92 if oldref!='ICRS': 

93 casalog.post('Position reference is '+oldref+'. Will convert positions to '+newref, 'INFO') 

94 

95 if obsloc=='GEOCENTRIC': 

96 casalog.post('Obsloc is already GEOCENTRIC. Nothing to be done.', 'INFO') 

97 return True 

98 

99 mepos = {'m0': {'value': geolong, 'unit': 'deg'}, 

100 'm1': {'value': geolat, 'unit': 'deg'}, # latitude 

101 'm2': {'value': geodist, 'unit': 'km'}, # alt above ref ellipsoid 

102 'refer': 'WGS84', 

103 'type': 'position'} 

104 met.doframe(mepos) 

105 

106 newra=[] 

107 newdec=[] 

108 newrho=[] 

109 newradvel=[] 

110 

111 for i in range(0, len(mjd)): 

112 memjd={'m0': {'value': mjd[i], 'unit': 'd'}, 

113 'refer': 'UTC', 

114 'type': 'epoch'} 

115 met.doframe(memjd) 

116 

117 olddir={'m0': {'value': ra[i], 'unit': 'deg'}, 

118 'm1': {'value': dec[i], 'unit': 'deg'}, 

119 'refer': oldref, 

120 'type': 'direction'} 

121 met.doframe(olddir) 

122 newdir=met.measure(olddir, newref) 

123 

124 tmpnewra = qat.convert(newdir['m0'],'deg')['value'] 

125 if tmpnewra<0: 

126 tmpnewra += 360. 

127 newra.append(tmpnewra) 

128 newdec.append(qat.convert(newdir['m1'],'deg')['value']) 

129 

130 oldradvel={'m0': {'value': radvel[i], 'unit': radvelunit}, 

131 'refer': 'TOPO', 

132 'type': 'radialvelocity'} 

133 

134 newradvelme = met.measure(oldradvel, 'GEO') 

135 newradvel.append(qat.convert(newradvelme['m0'], radvelunit)['value']) 

136 

137 # create the converted table 

138 safetycopyname=tablename+'.orig' 

139 if overwrite: 

140 if outtablename==tablename: 

141 os.system('cp -R '+tablename+' '+safetycopyname) 

142 else: 

143 os.system('rm -rf '+outtablename) 

144 os.system('cp -R '+tablename+' '+outtablename) 

145 else: 

146 if os.path.exists(outtablename): 

147 casalog.post('Output table '+outtablename+' already exists.', 'WARN') 

148 return False 

149 os.system('cp -R '+tablename+' '+outtablename) 

150 

151 try: 

152 tbt.open(outtablename, nomodify=False) 

153 tbt.putcol('RA', newra) 

154 tbt.putcol('DEC', newdec) 

155 tbt.putcol('RadVel', newradvel) 

156 tbt.putkeyword('GeoDist', 0.) 

157 tbt.putkeyword('GeoLat', 0.) 

158 tbt.putkeyword('GeoLong', 0.) 

159 tbt.putkeyword('obsloc', 'GEOCENTRIC') 

160 tbt.putkeyword('posrefsys', newposref) 

161 tbt.close() 

162 except Exception as instance: 

163 casalog.post("*** Error \'%s\' " % (instance), 'SEVERE') 

164 if overwrite and outtablename==tablename: 

165 casalog.post('Conversion in situ was not possible. Restoring original ephemeris ...', 'INFO') 

166 os.system('rm -rf '+tablename) 

167 os.system('mv '+safetycopyname+' '+tablename) 

168 return False 

169 

170 if overwrite and outtablename==tablename: 

171 os.system('rm -rf '+safetycopyname) 

172 

173 return True 

174 

175 

176def findattachedephemfields(vis='',field='*'): 

177 mst = ms() 

178 tbt = table() 

179 

180 tbt.open(vis+'/FIELD') 

181 fields = mst.msseltoindex(vis=vis,field=field)['field'] 

182 if(len(fields) == 0): 

183 casalog.post( "Field selection returned zero results.", 'WARN') 

184 return [] 

185 

186 thefields = [] 

187 

188 if ('EPHEMERIS_ID' in tbt.colnames()): 

189 for fld in fields: 

190 theid = tbt.getcell('EPHEMERIS_ID',fld) 

191 if theid>=0: # there is an ephemeris attached 

192 thefields.append(fld) 

193 

194 tbt.close() 

195 return thefields 

196 

197def convert2geo(vis='', field=''): 

198 

199 """ 

200 Convert the ephemeris attached to the given field of the given MS to GEO. 

201 Only RA, Dec, and RadVel are converted. 

202 If there is attached ephemeris or if the ephemeris is already in GEO, nothing is done. 

203 """ 

204 

205 mst = ms() 

206 tbt = table() 

207 

208 rval = True 

209 

210 try: 

211 fields = mst.msseltoindex(vis=vis,field=field)['field'] 

212 if(len(fields) == 0): 

213 casalog.post( "Field selection returned zero results.", 'WARN') 

214 return True 

215 

216 theephstoconvert = set([]) 

217 

218 tbt.open(vis+'/FIELD') 

219 

220 if ('EPHEMERIS_ID' in tbt.colnames()): 

221 for fld in fields: 

222 theid = tbt.getcell('EPHEMERIS_ID',fld) 

223 if theid>=0: # there is an ephemeris attached 

224 ephemname = glob.glob(vis+'/FIELD/EPHEM'+str(theid)+'_*.tab')[0] 

225 theephstoconvert.add(ephemname) 

226 casalog.post('Found ephemeris '+ephemname+' for field '+str(fld), 'INFO') 

227 

228 tbt.close() 

229 

230 if len(theephstoconvert)==0: 

231 casalog.post('No ephemerides attached.', 'INFO') 

232 else: 

233 for theeph in theephstoconvert: 

234 if converttopoephem2geo(theeph, theeph, overwrite=True): 

235 casalog.post('Converted '+theeph+' to GEO.', 'INFO') 

236 else: 

237 casalog.post('Error converting '+theeph+' to GEO.', 'WARN') 

238 rval = False 

239 

240 return rval 

241 

242 except Exception as instance: 

243 casalog.post("*** Error \'%s\' " % (instance), 'SEVERE') 

244 return False 

245 

246 

247