Coverage for /wheeldirectory/casa-CAS-14623-1-py3.12.el8/lib/py/lib/python3.12/site-packages/casatasks/private/tec_maps.py: 61%

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 if len(CASA_image)>0: 

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

407 if len(plot_name)>0: 

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

409 else: 

410 plot_name='none' 

411 return CASA_image,CASA_RMS_image,plot_name 

412 

413 

414 

415 

416 

417def get_IGS_TEC(ymd_date): 

418 """ 

419## ============================================================================= 

420## 

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

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

423## 

424## ============================================================================= 

425## 

426## Inputs: 

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

428## which this will retrieve data 

429## 

430## Returns: 

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

432## in the TEC map 

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

434## in the TEC map 

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

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

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

438## IGS data: 5 deg 

439## MAPGPS data: 1 deg 

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

441## latitude increases 

442## IGS data: 2.5 deg 

443## MAPGPS data: 1 deg 

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

445## IGS data: 120 min 

446## MAPGPS data: 5 min 

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

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

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

450## the TEC/DTEC in TECU 

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

452## CASA table keyword 

453## 

454## ============================================================================= 

455 """ 

456 

457 # Ensure retrieved files are stored locally 

458 localDir=os.getcwd()+'/' 

459 

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

461 # (used to set filenames on remote server) 

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

463 

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

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

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

467 

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

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

470 

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

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

473 

474 ## ========================================================================= 

475 ## 

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

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

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

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

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

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

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

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

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

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

486 ## 

487 ## ========================================================================= 

488 

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

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

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

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

493 

494 

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

496 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' 

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

498 

499 

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

501 

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

503 #ftps login and navigation 

504 try: 

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

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

507 ftps.prot_p() 

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

509 

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

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

512 get_path = file_location+get_file 

513 if test_IONEX_connection(get_path): 

514 local_get_file=localDir+get_file 

515 with open(local_get_file, 'wb') as fp: 

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

517 os.system('gunzip '+local_get_file) 

518 else: 

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

520 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' 

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

522 

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

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

525 get_path = file_location+get_file 

526 if test_IONEX_connection(get_path): 

527 local_get_file=localDir+get_file 

528 with open(local_get_file, 'wb') as fp: 

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

530 os.system('gunzip '+local_get_file) 

531 

532 else: 

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

534 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' 

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

536 

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

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

539 get_path = file_location+get_file 

540 if test_IONEX_connection(get_path): 

541 local_get_file=localDir+get_file 

542 with open(local_get_file, 'wb') as fp: 

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

544 os.system('gunzip '+local_get_file) 

545 else: 

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

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

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

549 return 0,0,0,0,0,0,0,0,[0],'' 

550 else: 

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

552 else: 

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

554 else: 

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

556 

557 ftps.quit() 

558 

559 except ftplib.all_errors as e: 

560 if 'ftps' in locals(): 

561 # only quit ftps connection if it was assigned/constructed successfully above 

562 ftps.quit() 

563 print('Failed to successfully connect to and/or retrieve from ftp://gdc.cddis.eosdis.nasa.gov/gps/products/ionex/ with error: ', e) 

564 print('NB: A working internet connection is required!') 

565 print('NB: Ionex files may not be available for dates prior to circa June 1998.') 

566 return 0,0,0,0,0,0,0,0,[0],'' 

567 else: 

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

569 

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

571 tec_type = 'IGS_Final_Product' 

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

573 tec_type = 'IGS_Rapid_Product' 

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

575 tec_type = 'JPL_Rapid_Product' 

576 

577 ## ========================================================================= 

578 ## 

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

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

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

582 ## later in this script. 

583 ## 

584 ## ========================================================================= 

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

586 

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

588 igsf=open(igs_file, 'r') 

589 linestring = igsf.read() 

590 igsf.close() 

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

592 

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

594 AddToList = 0 

595 TECLongList = [] 

596 DTECLongList = [] 

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

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

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

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

601 AddToList = 2 

602 if AddToList == 1: 

603 TECLongList.append(LongList[i]) 

604 if AddToList == 2: 

605 DTECLongList.append(LongList[i]) 

606 ## Determine the number of TEC/DTEC maps 

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

625 AddToList = 1 

626 

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

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

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

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

631 

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

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

634 

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

636 for Titer in range(2): 

637 counterMaps = 1 

638 UseList = [] 

639 if Titer == 0: 

640 UseList = TECLongList 

641 elif Titer == 1: 

642 UseList = DTECLongList 

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

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

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

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

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

648 ## at other latitudes within the selected map 

649 counterLat = 0 

650 newstartLat = float(str(start_lat)) 

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

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

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

654 counterLon = 0 

655 for row_iter in range(number_of_rows): 

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

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

658 counterLon = counterLon + 1 

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

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

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

662 counterLon = 0 

663 for row_iter in range(number_of_rows): 

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

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

666 counterLon = counterLon + 1 

667 counterLat = counterLat + row_iter + 2 

668 newstartLat = newstartLat + incr_lat 

669 counterMaps = counterMaps + 1 

670 

671 ## ========================================================================= 

672 ## 

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

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

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

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

677 ## 

678 ## ========================================================================= 

679 

680 ## The native sampling of the IGS maps minutes 

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

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

683 

684 for Titer in range(2): 

685 incr = 0 

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

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

688 

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

690 

691 

692 

693 

694 

695def check_existence(file_prefix): 

696 """ 

697## ============================================================================= 

698## 

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

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

701## and make the CASA table once. 

702## 

703## ============================================================================= 

704##  

705## Inputs: 

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

707## search for in the current directorty 

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

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

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

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

712## 

713## Returns: 

714## The name of the file it located. 

715## 

716## ============================================================================= 

717 """ 

718 file_name = '' 

719 return_file = '' 

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

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

722 if file_name.startswith(file_prefix): 

723 return_file = file_prefix 

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

725 if file_name.startswith(file_prefix): 

726 return_file = file_prefix 

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

728 if file_name.startswith(file_prefix): 

729 return_file = file_prefix 

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

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

732 return_file = file_name[:-4] 

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

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

735 return_file = file_prefix 

736 return return_file 

737 

738 

739 

740 

741 

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

743 """ 

744## ============================================================================= 

745## 

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

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

748## 

749## ============================================================================= 

750## 

751## Inputs: 

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

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

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

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

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

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

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

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

760## UT 0 on the first day 

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

762## IGS data: 5 deg 

763## MAPGPS data: 1 deg 

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

765## IGS data: 2.5 deg 

766## MAPGPS data: 1 deg 

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

768## IGS data: 120 min 

769## MAPGPS data: 5 min 

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

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

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

773## CASA table keyword 

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

775## 

776## Returns: 

777## The name of the TEC map image 

778## 

779## ============================================================================= 

780 """ 

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

782 

783 ## Set the coordinate system for the TEC image 

784 cs0=cs.newcoordsys(linear=3) 

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

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

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

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

789 

790 ## Make and view the TEC image 

791 imname=prefix+appendix+'.im' 

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

793 ia.summary() 

794 ## ia.statistics() 

795 ia.close() 

796 

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

798 tb.open(imname,nomodify=False) 

799 tb.putkeyword('TYPE',tec_type) 

800 tb.close() 

801 

802 return imname 

803 

804 

805 

806 

807def 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=''): 

808 """ 

809## ============================================================================= 

810## 

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

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

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

814## 

815## ============================================================================= 

816## 

817## Inputs: 

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

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

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

821## in the TEC map 

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

823## in the TEC map 

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

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

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

827## IGS data: 5 deg 

828## MAPGPS data: 1 deg 

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

830## IGS data: 2.5 deg 

831## MAPGPS data: 1 deg 

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

833## IGS data: 120 min 

834## MAPGPS data: 5 min 

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

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

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

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

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

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

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

842## 

843## Returns: 

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

845## for the local long./lat. 

846## 

847## ============================================================================= 

848 """ 

849 indexLat = 0 

850 indexLon = 0 

851 n = 0 

852 m = 0 

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

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

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

856 lowerIndexLon = n + 1 

857 higherIndexLon = n + 2 

858 n = n + 1