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

1#!/usr/bin/env python 

2 

3## ============================================================================= 

4## ============================================================================= 

5## 

6## tec_maps.py 

7## 

8## Created by J. E. Kooi 2014/04/25 v2.0 

9## Modified by gmoellen 2014/07/21 v2.1 (minor mods, mostly semantics) 

10## Modified by J. E. Kooi 2014/09/17 v2.2 (minor mods, formating and 

11## ~ Test for network connection 

12## ~ Heirarchy of IGS-related  

13## data products added 

14## ~ Added IGS RMS TEC maps 

15## ~ Made IGS default data product) 

16## Modified by gmoellen 2014/09/30 v2.3 (added simplified create method 

17## that hides all of the MAPGPS-related 

18## options, and thereby streamlines 

19## approved general usage for 

20## CASA 4.3) 

21## Modified by J. E. Kooi 2015/01/29 v2.4 (Changed get_IGS_TEC to accept 

22## IONEX format data with ANY grid  

23## spacing in deg. and time res in min) 

24## Modified by J. E. Kooi 2015/03/2 v2.5 (Changed ztec_value for doplot =  

25## True to plot a red bar over the VTEC 

26## to show the observing session) 

27## Modified by gmoellen 2017/05/30 v2.6 Generate plot disk file 

28## Modified by bkent 2017/10-31 v2.7 Fixed VisibleDeprecationWarning:  

29## using a non-integer  

30## number instead of an integer will  

31## result in an error in the future 

32## Modified by gmoellen 2020/08/28 v2.8 Updated TEC retrieval to use more-secure 

33## CDDIS server (ftp-ssl), as old 

34## service going offline 2020Oct31; 

35## other python cleanup 

36## 2020/10/13 python3-ification, for CASA 6. 

37##  

38## Tested in CASA 4.3.0 and 4.2.1 on RHEL release 6.4 (Santiago) 

39## Tested in CASA 4.3.0, 4.2.2, and 4.2.1 on Mac OS 10.8.5 

40## 

41## Modified by gmoellen 2023/12/12 v2.9 Updated to new CDDIS filename 

42## convention; use gunzip instead 

43## of uncompress; shrink labels 

44## in VLA-centric plot; misc. 

45## (incomplete) python cleanup 

46## 

47## 

48## The purpose of this python module is to retrieve vertical TEC/DTEC maps from 

49## either IGS (housed on the CDDIS servers) or MAPGPS (housed on the Madrigal  

50## servers). 

51## 

52## **** Currently, the MAPGPS is "turned off" and modification of line 306 **** 

53## **** is required to turn this functionality back "on." **** 

54## 

55## The IGS Final Product maps represent a combined TEC map: 

56## The different IGS Ionosphere Associate Analysis Centers IAAC TEC maps have  

57## been computed with different approaches but with a common formal resolution 

58## of 2 hours, 5 degrees, and 2.5 degrees in UT, longitude, and latitude  

59## (details can be found in e.g. Schaer, 1999; Feltens, 1998;  

60## Mannucci et al., 1998; Gao et al.; Hernandez-Pajares et al.,1999). The four  

61## IAACs TEC maps have been combined in an IGS combination product using weights 

62## obtained by two IGS Ionosphere Associate Validation Centers (IAVCs) from the 

63## corresponding performances in reproducing STEC and differences of STEC (IAVCs 

64## NRCAN and UPC respectively, see details in Feltens, 2002). 

65## ------ Krankowski & Sieradzki, 2013 (references can be found in this paper) 

66## 

67## The code has a hierarchy for IGS data products. It will initially search for 

68## the IGS Final product data (available ~ 10-14 days after an observation) and, 

69## if unavailable, then the IGS Rapid Product (available in ~ 1-2 days) and, 

70## if unavailable, then the JPL Rapid Product.  

71## 

72## The MAPGPS maps are computed with no pre-applied ionosphere model and, in  

73## this sense, represent 'raw TEC data.' The benefit is that they have a formal 

74## resolution of 5 minutes, 1.0 degree, and 1.0 degree in UT, longitude, and  

75## latitude (details can be found in Rideout & Coster, 2006). The drawback is 

76## that the data is sparsely gridded in many locations (e.g., there are no GPS 

77## ground stations in the middle of the ocean) and there are sporadic gaps in 

78## the data at some locations (e.g., even over the United States). 

79## 

80## For now, users are encouraged to use the IGS Final Product maps instead of 

81## MAPGPS products because there is data at all grid points and times,  

82## making it easier to deal with in the code. The code for MAPGPS is also set 

83## to make a 'patch' over North America and therefore is not a global map. 

84## However, the code is included and need only be uncommented. Feel free to 

85## experiment and compare the IGS and MAPGPS maps! 

86## 

87## This module calls several methods depending on the TEC/DTEC server and 

88## 

89## (1) Produces a CASA image of the TEC data for 

90## a) The whole world (IGS) (additional DTEC data image is output) 

91## b) A 15 deg. patch centered on the VLA (MAPGPS) 

92## (2) Optionally produces a vertical TEC/DTEC time series for the VLA 

93## 

94## !!!!! Important Note: In order to access the MAPGPS data, you must: 

95## (1) Have "madrigalWeb.py" and "globalDownload.py" and  

96## "madrigal.head" in the working directory 

97## (2) Provide your name, e-mail, and institution in the .create0 method to 

98## access the Madrigal server 

99## 

100## References: 

101## 

102## Krankowski, A., & Sieradzki, R. 2013, in International GNSS Service  

103## Technical Report, ed. R. Dach & Y. Jean, IGS Central Bureau  

104## (Astronomical Institute:University of Bern), 157 

105## Rideout, W., & Coster, A. 2006, GPS Solutions, 10, 219 

106## 

107## ============================================================================= 

108## 

109## Example Use: 

110## 

111## (0) The default use provides the IGS Product: 

112## > import tec_maps 

113## > msname = 'visibilities.ms' 

114## > CASA_image,CASA_RMS_image = tec_maps.create0(msname) 

115## 

116## (1) IGS Product (w/out the vertical TEC/DTEC time series at the VLA) 

117## > import tec_maps 

118## > msname = 'visibilities.ms' 

119## > CASA_image,CASA_RMS_image = tec_maps.create0(msname,'IGS') 

120## > viewer(CASA_image) 

121## > viewer(CASA_RMS_image) 

122## 

123## (2) IGS Product (WITH the vertical TEC/DTEC time series at the VLA) 

124## > import tec_maps 

125## > msname = 'visibilities.ms' 

126## > CASA_image,CASA_RMS_image = tec_maps.create0(msname,'IGS',True) 

127## > viewer(CASA_image) 

128## > viewer(CASA_RMS_image) 

129## 

130## (3) IGS Product (WITH the vertical TEC/DTEC time series at the VLA and for 

131## which the user specifies the output image name) 

132## > import tec_maps 

133## > msname = 'visibilities.ms' 

134## > imagename = 'my_image' 

135## > CASA_image,CASA_RMS_image = tec_maps.create0(msname,'IGS',True,imagename) 

136## > viewer(CASA_image) 

137## > viewer(CASA_RMS_image) 

138## 

139## (4) MAPGPS (WITH the vertical TEC/DTEC time series for the VLA) 

140## > import tec_maps 

141## > msname = 'visibilities.ms' 

142## > myname = 'John Doe' 

143## > myemail = 'john-doe@place.edu' 

144## > myplace = 'NRAO' 

145## > CASA_image,CASA_RMS_image =  

146## tec_maps.create0(msname,'MAPGPS',True,'',myname,myemail,myplace) 

147## > viewer(CASA_image) 

148## > viewer(CASA_RMS_image) 

149## 

150## 

151## ============================================================================= 

152## ============================================================================= 

153 

154 

155import glob, pylab, os, datetime 

156import numpy as np 

157from matplotlib import rc 

158import matplotlib.pyplot as plt 

159import ftplib 

160 

161from casatools import table, quanta, coordsys, image, measures 

162 

163tb = table() 

164qa = quanta() 

165cs = coordsys() 

166ia = image() 

167me = measures() 

168 

169workDir = os.getcwd()+'/' 

170 

171def create(vis,doplot=False,imname=''): 

172 """ 

173## ============================================================================= 

174## 

175## This method opens the .ms to determine the days of observation and calls the 

176## necessary functions to acquire the desired set of TEC/DTEC data. The method 

177## finally calls ztec_value and make_image once it has the data. 

178## 

179## ============================================================================= 

180## 

181## Inputs: 

182## vis type = string Name of the measurement set for which to 

183## acquire TEC/DTEC data 

184## doplot type = boolean When True, this will open a plot of the  

185## interpolated TEC/DTEC at the VLA. 

186## imname type = string Name of the output TEC Map optionally 

187## specified by the user 

188## 

189## Returns: 

190## Opens a plot showing the zenith TEC/DTEC for the VLA  

191## (if doplot = True) and the name of the CASA image file containing 

192## the TEC map. 

193## 

194## ============================================================================= 

195  

196 Usage:  

197 The default use provides the IGS Product: 

198 > vis = 'visibilities.ms' 

199 > CASA_image,CASA_RMS_image = tec_maps.create(vis) 

200 

201 IGS Product (WITH the vertical TEC/RMS TEC time series at the VLA) 

202 > vis = 'visibilities.ms' 

203 > CASA_image,CASA_RMS_image = tec_maps.create(msname,plot_vla_tec = True) 

204 

205 The TEC and RMS TEC images can then be examined using viewer: 

206 > viewer(CASA_image) 

207 > viewer(CASA_RMS_image) 

208 

209 The TEC image should then be used in gencal (caltype='tecim') to 

210 generate a sampled caltable that will nominally correct for 

211 ionospheric effects 

212 """ 

213 # call the more general method 

214 return create0(ms_name=vis,plot_vla_tec=doplot,im_name=imname) 

215 

216 

217def create0(ms_name,tec_server='IGS',plot_vla_tec=False,im_name='',username='',user_email='',affiliation=''): 

218 """ 

219## ============================================================================= 

220## 

221## This opens the .ms to determine the days of observation and calls the 

222## necessary functions to acquire the desired set of TEC/DTEC data. The method 

223## finally calls ztec_value and make_image once it has the data. 

224## 

225## ============================================================================= 

226## 

227## Inputs: 

228## ms_name type = string Name of the measurement set for which to 

229## acquire TEC/DTEC data 

230## tec_server type = string Server from which to retrieve TEC/DTEC 

231## plot_vla_tec type = boolean When True, this will open a plot of the  

232## interpolated TEC/DTEC at the VLA. 

233## im_name type = string Name of the output TEC Map optionally 

234## specified by the user 

235## username type = string MAPGPS only: full name of user accessing 

236## the site 

237## ex: '<First> <Last>' 

238## user_email type = string MAPGPS only: e-mail address at which the  

239## user can be reached 

240## ex: '<name>@<place>.com' 

241## affiliation type = string MAPGPS only: user's affiliated institute 

242## ex: 'NRAO' 

243## 

244## Returns: 

245## Opens a plot showing the zenith TEC/DTEC for the VLA  

246## (if plot_vla_tec = True) and the name of the CASA image file containing 

247## the TEC map. 

248## 

249## ============================================================================= 

250  

251 Usage:  

252 The default use provides the IGS Product: 

253 > msname = 'visibilities.ms' 

254 > CASA_image,CASA_RMS_image = tec_maps.create(msname) 

255 

256 IGS Product (WITH the vertical TEC/RMS TEC time series at the VLA) 

257 > msname = 'visibilities.ms' 

258 > CASA_image,CASA_RMS_image = tec_maps.create(msname,plot_vla_tec = True) 

259 

260 The TEC and RMS TEC images can then be examined using viewer: 

261 > viewer(CASA_image) 

262 > viewer(CASA_RMS_image) 

263 """ 

264 

265 ## Open the .ms to get the range in observation times 

266 tb.open(ms_name+'/OBSERVATION') 

267 obs_times = tb.getcol('TIME_RANGE') 

268 tb.close() 

269 

270 t_min = min(obs_times) 

271 t_max = max(obs_times) 

272 

273 ## Calculate the reference time for the TEC map to be generated. 

274 ref_time = 86400.*np.floor(t_min[0]/86400) 

275 ref_start = t_min[0]-ref_time 

276 ref_end = t_max[0]-ref_time 

277 

278 ## Gets the day string and the integer number of days of observation (only tested for two continuous days) 

279 begin_day = qa.time(str(t_min[0])+'s',form='ymd')[0][:10] 

280 end_day = qa.time(str(t_max[0])+'s',form='ymd')[0][:10] 

281 num_of_days = int(np.ceil((t_max[0]-ref_time)/86400.)) # must be an int as used below! 

282 

283 ## Set up the number of times we need to go get TEC files 

284 # (this is probably superfluous.... 2020Oct13 gmoellen) 

285 if begin_day == end_day: 

286 call_num = 1 

287 elif num_of_days == 0: 

288 call_num = 2 

289 else: 

290 call_num = num_of_days 

291 

292 ## Set up the days for which we need to go get TEC files 

293 day_list = [] 

294 next_day = begin_day 

295 for iter in range(call_num): 

296 day_list.append(next_day) 

297 next_day = qa.time(str(t_min[0]+86400.*(iter+1))+'s',form='ymd')[0][:10] 

298 

299 print('IGS files required for: '+str(day_list)) 

300 

301 # plot name declared here to avoid potential reference before assignment error 

302 plot_name = '' 

303 ## Runs the IGS methods 

304 if tec_server == 'IGS': 

305 ymd_date_num = 0 

306 #array = [] 

307 lo=0 

308 hi=0 

309 for ymd_date in day_list: 

310 points_long,points_lat,ref_long,ref_lat,incr_long,incr_lat,incr_time,num_maps,tec_array,tec_type = get_IGS_TEC(ymd_date) 

311 

312 ## Fill a new array with all the full set of TEC/DTEC values for all days in the observation set. 

313 if tec_type != '': 

314 if ymd_date_num == 0: 

315 # first pass make full output tec array 

316 # (this assumes all files have same number of timestamps (num_maps)) 

317 full_tec_array = np.zeros((2,int(points_long),int(points_lat),(int(num_maps)-1)*int(call_num)+1)) 

318 lo=0 

319 hi=int(num_maps) 

320 else: 

321 # First time of each subsequent file duplicates last of previous one, so overwrite it 

322 lo+=(int(num_maps)-1) 

323 hi=lo+int(num_maps) 

324 

325 full_tec_array[:,:,:,lo:hi] = tec_array[:,:,:,:] 

326 

327 ymd_date_num +=1 

328 

329 if tec_type != '': 

330 

331 if im_name == '': 

332 prefix = ms_name 

333 else: 

334 prefix = im_name 

335 

336 plot_name='' 

337 if plot_vla_tec: 

338 plot_name=prefix+'.IGS_TEC_at_site.png' 

339 

340 ztec_value(-107.6184,34.0790,points_long,points_lat,ref_long,ref_lat,incr_long,\ 

341 incr_lat,incr_time,ref_start,ref_end,int((num_maps-1)*call_num+1),\ 

342 full_tec_array,plot_vla_tec,plot_name) 

343 

344 CASA_image = make_image(prefix,ref_long,ref_lat,ref_time,incr_long,incr_lat,\ 

345 incr_time*60,full_tec_array[0],tec_type,appendix = '.IGS_TEC') 

346 CASA_RMS_image = make_image(prefix,ref_long,ref_lat,ref_time,incr_long,incr_lat,\ 

347 incr_time*60,full_tec_array[1],tec_type,appendix = '.IGS_RMS_TEC') 

348 else: 

349 CASA_image = '' 

350 CASA_RMS_image = '' 

351 

352 ## Runs the Madrigal methods; requires "madrigalWeb.py" and "globalDownload.py" and "madrigal.head" 

353 if tec_server == 'MAPGPS': 

354 print('\nCurrently, MAPGPS has been set to "False" \nand you must change tec_maps.py at line 306\n') 

355 CASA_image = '' 

356 CASA_RMS_image = '' 

357 

358 if (tec_server == 'MAPGPS' and False): ## replace with: if tec_server == 'MAPGPS': 

359 mad_file_front = 'gps'+str(begin_day.split('/')[0])[2:]+str(begin_day.split('/')[1])+str(begin_day.split('/')[2]) 

360 mad_data_file = ms_name[:-3]+'_MAPGPS_Data' 

361 

362 mad_file = '' 

363 mad_file = check_existence(mad_data_file) 

364 if username!='' and user_email!='' and affiliation!='': 

365 if mad_file == '': 

366 print('Retrieving the following MAPGPS files: '+begin_day+' to '+end_day) 

367 begin_mdy = str(begin_day.split('/')[1])+'/'+str(begin_day.split('/')[2])+'/'+str(begin_day.split('/')[0]) 

368 end_mdy = str(end_day.split('/')[1])+'/'+str(end_day.split('/')[2])+'/'+str(end_day.split('/')[0]) 

369 os.system('python globalDownload.py --url=http://madrigal.haystack.mit.edu/madrigal --outputDir='+workDir+\ 

370 ' --user_fullname="'+username+'" --user_email='+user_email+' --user_affiliation='+affiliation+\ 

371 ' --format=ascii --startDate='+str(begin_mdy)+' --endDate='+str(end_mdy)+' --inst=8000') 

372 

373 non_day = qa.time(str(t_min[0]-86400.)+'s',form='ymd')[0][:10] 

374 non_file_prefix = 'gps'+str(non_day.split('/')[0])[2:]+str(non_day.split('/')[1])+str(non_day.split('/')[2]) 

375 unwanted_file = check_existence(non_file_prefix) 

376 os.system('rm -rf '+unwanted_file+'.txt') 

377 

378 files = glob.glob(r''+workDir+'*.txt') 

379 outfile = open(workDir+mad_data_file+'.txt','a') 

380 for y in files: 

381 newfile = open(y,'r+') 

382 data = newfile.read() 

383 newfile.close() 

384 outfile.write(data) 

385 outfile.close() 

386 

387 

388 ## Making the CASA table is expensive, so we do it only once. 

389 print('Making a CASA table of: '+str(mad_data_file)) 

390 make_CASA_table(mad_data_file) 

391 

392 try: 

393 CASA_image = convert_MAPGPS_TEC(ms_name,mad_data_file,ref_time,ref_start,\ 

394 ref_end,plot_vla_tec,im_name) 

395 except: 

396 print('An error was encountered retrieving/interpreting the MAPGPS files.') 

397 CASA_image = '' 

398 CASA_RMS_image = '' 

399 else: 

400 print('You need to supply your username, e-mail, and affiliation to access the Madrigal server.') 

401 CASA_image = '' 

402 CASA_RMS_image = '' 

403 

404 ## Returns the name of the TEC image generated 

405 print('The following TEC map was generated: '+CASA_image+' & '+CASA_RMS_image) 

406 if len(plot_name)>0: 

407 print('The following TEC zenith plot was generated: '+plot_name) 

408 else: 

409 plot_name='none' 

410 return CASA_image,CASA_RMS_image,plot_name 

411 

412 

413 

414 

415 

416def get_IGS_TEC(ymd_date): 

417 """ 

418## ============================================================================= 

419## 

420## Retrieves the IGS data and, specifically, the IGS Final Product: this is 

421## the combined TEC map from the four IAACs TEC maps. 

422## 

423## ============================================================================= 

424## 

425## Inputs: 

426## ymd_date type = string The 'yyyy/mm/dd' date of observation for  

427## which this will retrieve data 

428## 

429## Returns: 

430## points_long type = integer Total number of points in longitude (deg)  

431## in the TEC map 

432## points_lat type = integer Total number of points in latitude (deg) 

433## in the TEC map 

434## start_long type = float Initial value for the longitude (deg) 

435## end_lat type = float Final value for the latitude (deg) 

436## incr_long type = float Increment by which longitude increases 

437## IGS data: 5 deg 

438## MAPGPS data: 1 deg 

439## incr_lat type = float Absolute value of the increment by which  

440## latitude increases 

441## IGS data: 2.5 deg 

442## MAPGPS data: 1 deg 

443## incr_time type = float Increment by which time increases 

444## IGS data: 120 min 

445## MAPGPS data: 5 min 

446## num_maps type = integer Number of maps (or time samples) of TEC 

447## tec_array type = array 4D array with axes consisting of  

448## [TEC_type,long.,lat.,time] that gives 

449## the TEC/DTEC in TECU 

450## tec_type type = string Specifies the origin of TEC data as a 

451## CASA table keyword 

452## 

453## ============================================================================= 

454 """ 

455 

456 # GPS weeks measured from 1980Jan06 (==MJD 44244.0) 

457 # (used to set filenames on remote server) 

458 gpsweek=(me.epoch('UTC',ymd_date)['m0']['value']-44244.0)/7 

459 

460 year = int(ymd_date.split('/')[0]) 

461 month = int(ymd_date.split('/')[1]) 

462 day = int(ymd_date.split('/')[2]) 

463 

464 ## Gives the day of the year of any given year 

465 dayofyear = datetime.datetime(year,month,day).timetuple().tm_yday 

466 

467 ## Convert dayofyear to 3-digit string, padded with zeros for IONEX filename format 

468 dayofyear = str(dayofyear).zfill(3) 

469 

470 ## ========================================================================= 

471 ## 

472 ## This goes to the CDDIS website, and downloads and uncompresses the IGS  

473 ## file. The preference is for the IGS Final product (IGSG), available 

474 ## ~ 12-14 days after data is collected (i.e. data for September 1 should 

475 ## be available by September 14). If this file is not already in the 

476 ## current working directory, it will retrieve it. If the file is  

477 ## unavailable, it will try to retrieve the IGS Rapid Product (IGRG), 

478 ## released ~ 2-3 days after data is collected. If this file is 

479 ## unavailable, it will try to retrieve the JPL Rapid Product (JPRG),  

480 ## released ~ 1 day after data is collected. While the 'uncompress' command 

481 ## is not necessary, it is the most straightforward on Linux. 

482 ## 

483 ## ========================================================================= 

484 

485 #CDDIS = 'ftp://cddis.gsfc.nasa.gov/gnss/products/ionex/' # pre-2020Nov01 version 

486 CDDIS = 'ftp://gdc.cddis.eosdis.nasa.gov/gps/products/ionex/' # new, more secure ftp-ssl server (2020Nov01) 

487 file_location = CDDIS+str(year)+'/'+str(dayofyear)+'/' 

488 #curlcmd='curl -u anonymous:casa-feedback@nrao.edu --ftp-ssl-reqd ' 

489 

490 

491 ## The name of the IONEX file you require.  

492 igs_file='igsg'+str(dayofyear)+'0.'+str(year)[2:4]+'i' if ( gpsweek<2238) else 'IGS0OPSFIN_'+str(year)+str(dayofyear)+'0000_01D_02H_GIM.INX' 

493 get_file=igs_file + ('.Z' if gpsweek<2238 else '.gz') 

494 

495 print('\nFor '+ymd_date+', the required IGS file is called: '+igs_file) 

496 

497 if len(glob.glob(igs_file))<1: # file does not yet exist locally 

498 #ftps login and navigation 

499 try: 

500 ftps = ftplib.FTP_TLS(host = 'gdc.cddis.eosdis.nasa.gov') # ftp-ssl version 

501 ftps.login(user='anonymous', passwd='casa-feedback@nrao.edu') 

502 ftps.prot_p() 

503 ftps.cwd('gps/products/ionex/'+str(year)+'/'+str(dayofyear)+'/') 

504 

505 if len(glob.glob(igs_file))<1: # file does not yet exist locally 

506 print('Attempting retrieval of IGS Final product file: '+str(get_file)) 

507 get_path = file_location+get_file 

508 if test_IONEX_connection(get_path): 

509 with open(get_file, 'wb') as fp: 

510 ftps.retrbinary("RETR " + get_file, fp.write) 

511 os.system('gunzip '+get_file) 

512 print('FILENAME: ', get_file) 

513 else: 

514 # IGS final product file does not exist; try rapid product file 

515 igs_file='igrg'+str(dayofyear)+'0.'+str(year)[2:4]+'i' if ( gpsweek<2238) else 'IGS0OPSRAP_'+str(year)+str(dayofyear)+'0000_01D_02H_GIM.INX' 

516 get_file=igs_file + ('.Z' if gpsweek<2238 else '.gz') 

517 

518 if len(glob.glob(igs_file))<1: # file does not yet exist 

519 print('Attempting retrieval of IGS Rapid product file: '+str(get_file)) 

520 get_path = file_location+get_file 

521 if test_IONEX_connection(get_path): 

522 with open(get_file, 'wb') as fp: 

523 ftps.retrbinary("RETR " + get_file, fp.write) 

524 os.system('gunzip '+get_file) 

525 

526 else: 

527 #igs_file = igs_file.replace('igr','jpr') 

528 igs_file= 'jprg'+str(dayofyear)+'0.'+str(year)[2:4]+'i' if (gpsweek<2274.5) else 'JPL0OPSRAP_'+str(year)+str(dayofyear)+'0000_01D_02H_GIM.INX' 

529 get_file=igs_file + ('.Z' if gpsweek<2274.5 else '.gz') 

530 

531 if len(glob.glob(igs_file))<1: # file does not yet exist 

532 print('Attempting retrieval of JPL Rapid product file: '+str(igs_file)) 

533 get_path = file_location+get_file 

534 if test_IONEX_connection(get_path): 

535 with open(get_file, 'wb') as fp: 

536 ftps.retrbinary("RETR " + get_file, fp.write) 

537 os.system('gunzip '+get_file) 

538 else: 

539 print('\nNo data products available. You may try to manually'+\ 

540 ' download the products at:\n'+\ 

541 'ftp://gdc.cddis.eosdis.nasa.gov/gps/products/ionex\n') 

542 return 0,0,0,0,0,0,0,0,[0],'' 

543 else: 

544 print('JPL Rapid product file: '+igs_file+' already available in current working directory.') 

545 else: 

546 print('IGS Rapid product file: '+igs_file+' already available in current working directory.') 

547 else: 

548 print('IGS Final product file: '+igs_file+' already available in current working directory.') 

549 

550 ftps.quit() 

551 except ftplib.all_errors as e: 

552 print('Failed to connect to ftp://gdc.cddis.eosdis.nasa.gov/gps/products/ionex/ with error: ', e) 

553 else: 

554 print('IGS Final product file: '+igs_file+' already available in current working directory.') 

555 

556 if igs_file.startswith('igs') or igs_file.startswith('IGS0OPSFIN'): 

557 tec_type = 'IGS_Final_Product' 

558 elif igs_file.startswith('igr') or igs_file.startswith('IGS0OPSRAP'): 

559 tec_type = 'IGS_Rapid_Product' 

560 elif igs_file.startswith('jpr') or igs_file.startswith('JPL0OPSRAP'): 

561 tec_type = 'JPL_Rapid_Product' 

562 

563 ## ========================================================================= 

564 ## 

565 ## The following section reads the lines of the ionex file for 1 day  

566 ## (13 maps total) into an array a[]. It also retrieves the thin-shell 

567 ## ionosphere height used by IGS, the lat./long. spacing, etc. for use 

568 ## later in this script. 

569 ## 

570 ## ========================================================================= 

571 print('Transforming IONEX data format to a TEC/DTEC array for '+str(ymd_date)) 

572 

573 ## Opening and reading the IONEX file into memory as a list 

574 linestring = open(igs_file, 'r').read() 

575 LongList = linestring.split('\n') 

576 

577 ## Create two lists without the header and only with the TEC and DTEC maps (based on code from ionFR.py) 

578 AddToList = 0 

579 TECLongList = [] 

580 DTECLongList = [] 

581 for i in range(len(LongList)-1): 

582 ## Once LongList[i] gets to the DTEC maps, append DTECLongList 

583 if LongList[i].split()[-1] == 'MAP': 

584 if LongList[i].split()[-2] == 'RMS': 

585 AddToList = 2 

586 if AddToList == 1: 

587 TECLongList.append(LongList[i]) 

588 if AddToList == 2: 

589 DTECLongList.append(LongList[i]) 

590 ## Determine the number of TEC/DTEC maps 

591 if LongList[i].split()[-1] == 'FILE': 

592 if LongList[i].split()[-3:-1] == ['MAPS','IN']: 

593 num_maps = float(LongList[i].split()[0]) 

594 ## Determine the shell ionosphere height (usually 450 km for IGS IONEX files) 

595 if LongList[i].split()[-1] == 'DHGT': 

596 ion_H = float(LongList[i].split()[0]) 

597 ## Determine the range in lat. and long. in the ionex file 

598 if LongList[i].split()[-1] == 'DLAT': 

599 start_lat = float(LongList[i].split()[0]) 

600 end_lat = float(LongList[i].split()[1]) 

601 incr_lat = float(LongList[i].split()[2]) 

602 if LongList[i].split()[-1] == 'DLON': 

603 start_long = float(LongList[i].split()[0]) 

604 end_long = float(LongList[i].split()[1]) 

605 incr_long = float(LongList[i].split()[2]) 

606 ## Find the end of the header so TECLongList can be appended 

607 if LongList[i].split()[0] == 'END': 

608 if LongList[i].split()[2] == 'HEADER': 

609 AddToList = 1 

610 

611 ## Variables that indicate the number of points in Lat. and Lon. 

612 points_long = ((end_long - start_long)/incr_long) + 1 

613 points_lat = ((end_lat - start_lat)/incr_lat) + 1 ## Note that incr_lat is defined as '-' here 

614 number_of_rows = int(np.ceil(points_long/16)) ## Note there are 16 columns of data in IONEX format 

615 

616 ## 4-D array that will contain TEC & DTEC (a[0] and a[1], respectively) values 

617 a = np.zeros((2,int(points_long),int(points_lat),int(num_maps))) 

618 

619 ## Selecting only the TEC/DTEC values to store in the 4-D array. 

620 for Titer in range(2): 

621 counterMaps = 1 

622 UseList = [] 

623 if Titer == 0: 

624 UseList = TECLongList 

625 elif Titer == 1: 

626 UseList = DTECLongList 

627 for i in range(len(UseList)): 

628 ## Pointing to first map (out of 13 maps) then by changing 'counterMaps' the other maps are selected 

629 if UseList[i].split()[0] == ''+str(counterMaps)+'': 

630 if UseList[i].split()[-4] == 'START': 

631 ## Pointing to the starting Latitude then by changing 'counterLat' we select TEC data 

632 ## at other latitudes within the selected map 

633 counterLat = 0 

634 newstartLat = float(str(start_lat)) 

635 for iLat in range(int(points_lat)): 

636 if UseList[i+2+counterLat].split()[0].split('-')[0] == ''+str(newstartLat)+'': 

637 ## Adding to array a[] a line of Latitude TEC data 

638 counterLon = 0 

639 for row_iter in range(number_of_rows): 

640 for item in range(len(UseList[i+3+row_iter+counterLat].split())): 

641 a[Titer,counterLon,iLat,counterMaps-1] = UseList[i+3+row_iter+counterLat].split()[item] 

642 counterLon = counterLon + 1 

643 if '-'+UseList[i+2+counterLat].split()[0].split('-')[1] == ''+str(newstartLat)+'': 

644 ## Adding to array a[] a line of Latitude TEC data. Same chunk as above but  

645 ## in this case we account for the TEC values at negative latitudes 

646 counterLon = 0 

647 for row_iter in range(number_of_rows): 

648 for item in range(len(UseList[i+3+row_iter+counterLat].split())): 

649 a[Titer,counterLon,iLat,counterMaps-1] = UseList[i+3+row_iter+counterLat].split()[item] 

650 counterLon = counterLon + 1 

651 counterLat = counterLat + row_iter + 2 

652 newstartLat = newstartLat + incr_lat 

653 counterMaps = counterMaps + 1 

654 

655 ## ========================================================================= 

656 ## 

657 ## The section creates a new array that is a copy of a[], but with the lower 

658 ## left-hand corner defined as the initial element (whereas a[] has the  

659 ## upper left-hand corner defined as the initial element). This also 

660 ## accounts for the fact that IONEX data is in 0.1*TECU. 

661 ## 

662 ## ========================================================================= 

663 

664 ## The native sampling of the IGS maps minutes 

665 incr_time = 24*60/(int(num_maps)-1) 

666 tec_array = np.zeros((2,int(points_long),int(points_lat),int(num_maps))) 

667 

668 for Titer in range(2): 

669 incr = 0 

670 for ilat in range(int(points_lat)): 

671 tec_array[Titer,:,ilat,:] = 0.1*a[Titer,:,int(points_lat)-1-ilat,:] 

672 

673 return points_long,points_lat,start_long,end_lat,incr_long,np.absolute(incr_lat),incr_time,num_maps,tec_array,tec_type 

674 

675 

676 

677 

678 

679def check_existence(file_prefix): 

680 """ 

681## ============================================================================= 

682## 

683## Checks to see if the IGS or MAPGPS file for a given day already exists and 

684## returns the file name prefix. Primarily used to ensure we only call the data 

685## and make the CASA table once. 

686## 

687## ============================================================================= 

688##  

689## Inputs: 

690## file_prefix type = string This is the prefix for the file name to 

691## search for in the current directorty 

692## IGS form: 'igsg'+'ddd'+'0.'+'yy'+'i' 

693## ex: August 6, 2011 is 'igsg2180.11i' 

694## MAPGPS form: 'gps'+'yy'+'mm'+'dd' 

695## ex: August 6, 2011 is 'gps110806' 

696## 

697## Returns: 

698## The name of the file it located. 

699## 

700## ============================================================================= 

701 """ 

702 file_name = '' 

703 return_file = '' 

704 for file_name in [doc for doc in os.listdir(workDir)]: 

705 if (file_prefix.startswith('igs') and file_name.startswith('igs')): 

706 if file_name.startswith(file_prefix): 

707 return_file = file_prefix 

708 if (file_prefix.startswith('igr') and file_name.startswith('igr')): 

709 if file_name.startswith(file_prefix): 

710 return_file = file_prefix 

711 if (file_prefix.startswith('jpr') and file_name.startswith('jpr')): 

712 if file_name.startswith(file_prefix): 

713 return_file = file_prefix 

714 if (file_prefix.startswith('gps') and file_name.startswith('gps')): 

715 if (file_name.startswith(file_prefix) and file_name.endswith('.txt')): 

716 return_file = file_name[:-4] 

717 if (file_prefix.endswith('_MAPGPS_Data') and file_name.endswith('_MAPGPS_Data.tab')): 

718 if (file_name.startswith(file_prefix) and file_name.endswith('.tab')): 

719 return_file = file_prefix 

720 return return_file 

721 

722 

723 

724 

725 

726def make_image(prefix,ref_long,ref_lat,ref_time,incr_long,incr_lat,incr_time,tec_array,tec_type,appendix=''): 

727 """ 

728## ============================================================================= 

729## 

730## Creates a new image file with the TEC data and then returns the image name.  

731## This also sets up the reference frame for use at the C++ level. 

732## 

733## ============================================================================= 

734## 

735## Inputs: 

736## prefix type = string Full file name for use in naming the image 

737## IGS form: 'igsg'+'ddd'+'0.'+'yy'+'i' 

738## ex: August 6, 2011 is 'igsg2180.11i' 

739## MAPGPS form: 'gps'+'yy'+'mm'+'dd'+'.###' 

740## ex: August 17, 2011 is 'gps110817g.001' 

741## ref_long type = float Reference long. (deg) for setting coordinates 

742## ref_lat type = float Reference lat. (deg) for setting coordinates 

743## ref_time type = float Reference time (s) for setting coordinates,  

744## UT 0 on the first day 

745## incr_long type = float Increment by which longitude increases 

746## IGS data: 5 deg 

747## MAPGPS data: 1 deg 

748## incr_lat type = float Increment by which latitude increases 

749## IGS data: 2.5 deg 

750## MAPGPS data: 1 deg 

751## incr_time type = float Increment by which time increases 

752## IGS data: 120 min 

753## MAPGPS data: 5 min 

754## tec_array type = array 3D array with axes consisting of  

755## [long.,lat.,time] giving the TEC in TECU 

756## tec_type type = string Specifies the origin of TEC data as a 

757## CASA table keyword 

758## appendix type = string Appendix to add to the end of the image name 

759## 

760## Returns: 

761## The name of the TEC map image 

762## 

763## ============================================================================= 

764 """ 

765 print('Making image: '+str(prefix+appendix)+'.im') 

766 

767 ## Set the coordinate system for the TEC image 

768 cs0=cs.newcoordsys(linear=3) 

769 cs0.setnames(value='Long Lat Time') 

770 cs0.setunits(type='linear',value='deg deg s',overwrite=True) 

771 cs0.setreferencevalue([ref_long,ref_lat,ref_time]) 

772 cs0.setincrement(type='linear',value=[incr_long,incr_lat,incr_time]) 

773 

774 ## Make and view the TEC image 

775 imname=prefix+appendix+'.im' 

776 ia.fromarray(outfile=imname,pixels=tec_array,csys=cs0.torecord(),overwrite=True) 

777 ia.summary() 

778 ## ia.statistics() 

779 ia.close() 

780 

781 ## Specify in the image where the TEC data came from 

782 tb.open(imname,nomodify=False) 

783 tb.putkeyword('TYPE',tec_type) 

784 tb.close() 

785 

786 return imname 

787 

788 

789 

790 

791def ztec_value(my_long,my_lat,points_long,points_lat,ref_long,ref_lat,incr_long,incr_lat,incr_time,ref_start,ref_end,num_maps,tec_array,PLOT=False,PLOTNAME=''): 

792 """ 

793## ============================================================================= 

794## 

795## Determine the TEC value for the coordinates given at every time sampling of 

796## the TEC. This locates the 4 points in the IONEX grid map which surround the 

797## coordinate for which you want to calculate the TEC value. 

798## 

799## ============================================================================= 

800## 

801## Inputs: 

802## my_long type = float Long. (deg) at which this interpolates TEC 

803## my_lat type = float Lat. (deg) at which this interpolates TEC 

804## points_long type = integer Total number of points in long. (deg) 

805## in the TEC map 

806## points_lat type = integer Total number of points in lat. (deg)  

807## in the TEC map 

808## ref_long type = float Initial value for the longitude (deg) 

809## ref_lat type = float Initial value for the latitude (deg) 

810## incr_long type = float Increment by which longitude increases 

811## IGS data: 5 deg 

812## MAPGPS data: 1 deg 

813## incr_lat type = float Increment by which latitude increases 

814## IGS data: 2.5 deg 

815## MAPGPS data: 1 deg 

816## incr_time type = float Increment by which time increases 

817## IGS data: 120 min 

818## MAPGPS data: 5 min 

819## ref_start type = float Beginning of observations (in seconds) 

820## ref_end type = float End of observations (in seconds) 

821## num_maps type = integer Number of maps (or time samples) of TEC 

822## tec_array type = array 3D array with axes consisting of  

823## [long.,lat.,time] giving TEC in TECU 

824## PLOT type = boolean Determines whether to plot or return the 

825## TEC time series of the local long./lat. 

826## 

827## Returns: 

828## site_tec type = array 2D array containing the TEC/DTEC values  

829## for the local long./lat. 

830## 

831## ============================================================================= 

832 """ 

833 indexLat = 0 

834 indexLon = 0 

835 n = 0 

836 m = 0 

837 ## Find the corners of the grid that surrounds the local long./lat. 

838 for lon in range(int(points_long)): 

839 if (my_long > (ref_long + (n+1)*incr_long) and my_long <= (ref_long + (n+2)*incr_long)) : 

840 lowerIndexLon = n + 1 

841 higherIndexLon = n + 2 

842 n = n + 1 

843 for lat in range(int(points_lat)): 

844 if (my_lat > (ref_lat + (m+1)*incr_lat) and my_lat <= (ref_lat + (m+2)*incr_lat)) : 

845 lowerIndexLat = m + 1 

846 higherIndexLat = m + 2 

847 m = m + 1 

848 ## Using the 4-point formula indicated in the IONEX manual, find the TEC value at the local coordinates 

849 diffLon = my_long - (ref_long + lowerIndexLon*incr_long) 

850 WLON = diffLon/incr_long 

851 diffLat = my_lat - (ref_lat + lowerIndexLat*incr_lat) 

852 WLAT = diffLat/incr_lat 

853 

854 site_tec = np.zeros((2,num_maps))