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

1import os 

2import sys 

3 

4 

5from casatools import image, quanta 

6from casatasks import casalog 

7from .ialib import write_image_history 

8 

9_qa = quanta() 

10 

11def importfits(fitsimage,imagename,whichrep,whichhdu,zeroblanks,overwrite,defaultaxes,defaultaxesvalues,beam): 

12 """Convert an image FITS file into a CASA image: 

13 

14 Keyword arguments: 

15 fitsimage -- Name of input image FITS file 

16 default: none; example='3C273XC1.fits' 

17 imagename -- Name of output CASA image 

18 default: none; example: imagename='3C273XC1.image' 

19 whichrep -- If fits image has multiple coordinate reps, choose one. 

20 default: 0 means first; example: whichrep=1 

21 whichhdu -- If its file contains multiple images, choose one (0 = first HDU, -1 = first valid image) 

22 default=-1 ; example: whichhdu=1 

23 zeroblanks -- Set blanked pixels to zero (not NaN) 

24 default=True; example: zeroblanks=True 

25 overwrite -- Overwrite pre-existing imagename 

26 default=False; example: overwrite=True 

27 defaultaxes -- Add the default 4D coordinate axes where they are missing 

28 default=False, example: defaultaxes=True 

29 defaultaxesvalues -- List of values to assign to added degenerate axes when defaultaxes==True (ra,dec,freq,stokes) 

30 default = [], example: defaultaxesvalues=['13.5h', '-2.5deg', '88.5GHz', 'Q']  

31 beam -- List of values to be used to define the synthesized beam [BMAJ,BMIN,BPA] (as in the FITS keywords) 

32 default = [] (i.e. take from FITS file), example: beam=['0.35arcsec', '0.24arcsec', '25deg'] 

33 

34 """ 

35 

36 #Python script 

37 casalog.origin('importfits') 

38 _myia = image() 

39 tmpname = imagename 

40 reorder = False 

41 addaxes = False 

42 adddir = False 

43 addstokes = False 

44 addfreq=False 

45 defaultorder = ['Right Ascension', 'Declination', 'Stokes', 'Frequency'] 

46 addbeam = False 

47 

48 try: 

49 _myia.dohistory(False) 

50 if os.path.exists(imagename): 

51 if not overwrite: 

52 raise RuntimeError('Output image exists already and you did not set overwrite to True.') 

53 else: 

54 os.system('rm -rf '+imagename) 

55 

56 if defaultaxes: 

57 if len(defaultaxesvalues)!=4: 

58 raise TypeError('When defaultaxes==True, parameter defaultaxesvalues must be provided as a list of 4 values: RA, Dec, Freq, Stokes,\n e.g. [\'13.5h\', \'-2.5deg\', \'88.5GHz\', \'I\']\nFor existing axes, empty strings can be given as values.') 

59 _myia.open(fitsimage) 

60 _mycs=_myia.coordsys() 

61 acts = _mycs.axiscoordinatetypes() 

62 cnames = _mycs.names() 

63 _myia.close() 

64 if ('Direction' in acts and not ('Right Ascension' in cnames and 'Declination' in cnames)): 

65 raise RuntimeError('Non-standard direction axes. Cannot add default axes.') 

66 if ('Spectral' in acts and not 'Frequency' in cnames): 

67 raise RuntimeError('Non-standard spectral axis. Cannot add default axes.') 

68 if ('Stokes' in acts and not 'Stokes' in cnames): 

69 raise RuntimeError('Non-standard Stokes axis. Cannot add default axes.') 

70 

71 if not ('Right Ascension' in cnames and 'Declination' in cnames): 

72 addaxes = True 

73 adddir = True 

74 if not ('Frequency' in cnames): 

75 addaxes = True 

76 addfreq = True 

77 if not ('Stokes' in cnames): 

78 addaxes = True 

79 addstokes = True 

80 if not addaxes and cnames!=defaultorder: 

81 reorder = True 

82 if addaxes or reorder: 

83 tmpname = imagename+'.tmp' 

84 os.system('rm -rf '+tmpname) 

85 

86 if beam!=[]: # user has set beam 

87 if type(beam)!=list or len(beam)!=3 or not (type(beam[0])==str and type(beam[1])==str and type(beam[2])): 

88 raise TypeError("Parameter beam is invalid (should be list of 3 strings or empty): %s" % beam) 

89 qabeam = [] 

90 for i in range(0,3): 

91 try: 

92 qabeam.append(_qa.quantity(beam[i])) 

93 tmp = _qa.toangle(qabeam[i]) 

94 except: 

95 raise TypeError("Parameter beam[%s] is invalid (should be an angle): %s" % (i,beam[i])) 

96 if (_qa.convert(qabeam[0], 'arcsec')['value'] >= _qa.convert(qabeam[1], 'arcsec')['value']): 

97 addbeam = True 

98 else: 

99 raise TypeError("Parameter beam[%s] is invalid (major axis must be >= minor axis): %s" % (i,beam)) 

100 

101 _myia.fromfits(tmpname,fitsimage,whichrep,whichhdu,zeroblanks) 

102 _myia.close() 

103 

104 if addaxes: 

105 casalog.post('Adding missing coodinate axes ...', 'INFO') 

106 tmpname2 = imagename+'.tmp2' 

107 os.system('rm -rf '+tmpname2) 

108 

109 _myia.open(tmpname) 

110 ia2 = _myia.adddegaxes(outfile=tmpname2, direction=True, spectral=True, stokes=defaultaxesvalues[3], 

111 silent=True) 

112 _myia.close() 

113 os.system('rm -rf '+tmpname) 

114 ia2.close() 

115 ia2.open(tmpname2) 

116 

117 # set the right reference values in the added axes 

118 _mynewcs=ia2.coordsys() 

119 raval = 0. 

120 decval = 0. 

121 freqval = 0. 

122 if adddir: 

123 ra = defaultaxesvalues[0] 

124 if type(ra)==int or type(ra)==float: 

125 raval = ra 

126 else: 

127 qara = _qa.quantity(_qa.angle(ra)[0]) 

128 if qara['unit'].find('deg') < 0: 

129 raise TypeError("RA default value is not a valid angle quantity: %s" % ra) 

130 raval = qara['value'] 

131 

132 dec = defaultaxesvalues[1] 

133 if type(dec)==int or type(dec)==float: 

134 decval = dec 

135 else: 

136 qadec = _qa.quantity(_qa.angle(dec)[0]) 

137 if qadec['unit'].find('deg') < 0: 

138 raise TypeError("DEC default value is not a valid angle quantity: %s" % dec) 

139 decval = qadec['value'] 

140 

141 _mynewcs.setunits(value='deg deg', type='direction') 

142 _mynewcs.setreferencevalue(type='direction', value=[raval,decval]) 

143 

144 if addfreq: 

145 freq = defaultaxesvalues[2] 

146 if type(freq)==int or type(freq)==float: 

147 freqval = freq 

148 else: 

149 qafreq = _qa.quantity(freq) 

150 if qafreq['unit'].find('Hz') < 0: 

151 raise TypeError("Freq default value is not a valid frequency quantity: %s" % freq) 

152 freqval = _qa.convertfreq(qafreq,'Hz')['value'] 

153 _mynewcs.setunits(value='Hz', type='spectral') 

154 _mynewcs.setreferencevalue(type='spectral', value=freqval) 

155 _mynewcs.setrestfrequency(freqval) 

156 

157 # Note: stokes default value was already set in adddegaxes 

158 

159 if adddir or addfreq: 

160 ia2.setcoordsys(_mynewcs.torecord()) 

161 

162 ia2.close() 

163 

164 cnames = _mynewcs.names() 

165 

166 if len(cnames)==4 and not (cnames == defaultorder): 

167 # need to reorder 

168 reorder = True 

169 tmpname = tmpname2 

170 else: 

171 os.system('mv '+tmpname2+' '+imagename) 

172 

173 if reorder: 

174 casalog.post('Transposing coodinate axes ...', 'INFO') 

175 _myia.open(tmpname) 

176 ia2 = _myia.transpose(outfile=imagename, order=defaultorder) 

177 _myia.close() 

178 ia2.close() 

179 os.system('rm -rf '+tmpname) 

180 

181 _myia.open(imagename) 

182 if addbeam: 

183 _myia.setrestoringbeam(beam[0], beam[1], beam[2]) 

184 else: 

185 mybeam = _myia.restoringbeam() 

186 if mybeam =={}: # the fits image had no beam 

187 casalog.post("This image has no beam or angular resolution provided, so you will not receive warnings from\n" 

188 "tasks such as imregrid if your image pixels do not sample the the angular resolution well.\n" 

189 "(This only affects warnings, not any functionality).\n" 

190 "Providing a beam and brightness units in an image can also be useful for flux calculations.\n" 

191 "If you wish to add a beam or brightness units to your image, please use\n" 

192 "the \"beam\" parameter or ia.setrestoringbeam() and ia.setbrightnessunit()", 'WARN') 

193 try: 

194 param_names = importfits.__code__.co_varnames[:importfits.__code__.co_argcount] 

195 vars = locals( ) 

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

197 write_image_history(_myia, sys._getframe().f_code.co_name, param_names, param_vals, casalog) 

198 except Exception as instance: 

199 casalog.post("*** Error \'%s\' updating HISTORY" % instance, 'WARN') 

200 

201 finally: 

202 _myia.done()