Coverage for /wheeldirectory/casa-6.7.0-12-py3.10.el8/lib/py/lib/python3.10/site-packages/casatasks/private/tec_maps.py: 62%
335 statements
« prev ^ index » next coverage.py v7.6.4, created at 2024-11-01 07:19 +0000
« prev ^ index » next coverage.py v7.6.4, created at 2024-11-01 07:19 +0000
1#!/usr/bin/env python
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## =============================================================================
155import glob, pylab, os, datetime
156import numpy as np
157from matplotlib import rc
158import matplotlib.pyplot as plt
159import ftplib
161from casatools import table, quanta, coordsys, image, measures
163tb = table()
164qa = quanta()
165cs = coordsys()
166ia = image()
167me = measures()
169workDir = os.getcwd()+'/'
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## =============================================================================
196 Usage:
197 The default use provides the IGS Product:
198 > vis = 'visibilities.ms'
199 > CASA_image,CASA_RMS_image = tec_maps.create(vis)
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)
205 The TEC and RMS TEC images can then be examined using viewer:
206 > viewer(CASA_image)
207 > viewer(CASA_RMS_image)
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)
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## =============================================================================
251 Usage:
252 The default use provides the IGS Product:
253 > msname = 'visibilities.ms'
254 > CASA_image,CASA_RMS_image = tec_maps.create(msname)
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)
260 The TEC and RMS TEC images can then be examined using viewer:
261 > viewer(CASA_image)
262 > viewer(CASA_RMS_image)
263 """
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()
270 t_min = min(obs_times)
271 t_max = max(obs_times)
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
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!
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
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]
299 print('IGS files required for: '+str(day_list))
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)
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)
325 full_tec_array[:,:,:,lo:hi] = tec_array[:,:,:,:]
327 ymd_date_num +=1
329 if tec_type != '':
331 if im_name == '':
332 prefix = ms_name
333 else:
334 prefix = im_name
336 plot_name=''
337 if plot_vla_tec:
338 plot_name=prefix+'.IGS_TEC_at_site.png'
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)
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 = ''
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 = ''
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'
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')
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')
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()
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)
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 = ''
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
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 """
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
460 year = int(ymd_date.split('/')[0])
461 month = int(ymd_date.split('/')[1])
462 day = int(ymd_date.split('/')[2])
464 ## Gives the day of the year of any given year
465 dayofyear = datetime.datetime(year,month,day).timetuple().tm_yday
467 ## Convert dayofyear to 3-digit string, padded with zeros for IONEX filename format
468 dayofyear = str(dayofyear).zfill(3)
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 ## =========================================================================
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 '
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')
495 print('\nFor '+ymd_date+', the required IGS file is called: '+igs_file)
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)+'/')
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')
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)
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')
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.')
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.')
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'
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))
573 ## Opening and reading the IONEX file into memory as a list
574 linestring = open(igs_file, 'r').read()
575 LongList = linestring.split('\n')
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
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
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)))
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
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 ## =========================================================================
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)))
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,:]
673 return points_long,points_lat,start_long,end_lat,incr_long,np.absolute(incr_lat),incr_time,num_maps,tec_array,tec_type
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
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')
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])
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()
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()
786 return imname
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
854 site_tec = np.zeros((2,num_maps))
855 for Titer in range(2):
856 for m in range(num_maps):
857 site_tec[Titer,m] = (1.0-WLAT)*(1.0-WLON)*tec_array[Titer,lowerIndexLon,lowerIndexLat,m] +\
858 WLON*(1.0-WLAT)*tec_array[Titer,higherIndexLon,lowerIndexLat,m] +\
859 (1.0-WLON)*WLAT*tec_array[Titer,lowerIndexLon,higherIndexLat,m] +\
860 WLON*WLAT*tec_array[Titer,higherIndexLon,higherIndexLat,m]
862 if PLOT == True:
863 ## Set axis label size for the plots
864 rc('xtick', labelsize=8)
865 rc('ytick', labelsize=8)
866 plottimes = [x*incr_time for x in range(num_maps)]
867 plt.interactive(False)
868 plt.clf()
869 plt.errorbar(plottimes,site_tec[0],site_tec[1])
870 plt.axvspan(ref_start/60.0, ref_end/60.0, facecolor='r', alpha=0.5)
871 plt.xlabel(r'$\mathrm{Time}$ $\mathrm{(minutes)}$', fontsize=10)
872 plt.ylabel(r'$\mathrm{TEC}$ $\mathrm{(TECU)}$', fontsize=10)
873 plt.title(r'$\mathrm{TEC}$ $\mathrm{values}$ $\mathrm{for}$ $\mathrm{Long.}$ $\mathrm{=}$ '+\
874 '$\mathrm{'+str(my_long)+'}$ / $\mathrm{Lat.}$ $\mathrm{=}$ $\mathrm{'+str(my_lat)+'}$',\
875 fontsize=10)
876 plt.axis([min(plottimes),max(plottimes),0,1.1*max(site_tec[0])])
878 if len(PLOTNAME)>0:
879 plt.savefig( PLOTNAME )
881 if PLOT == False:
882 return site_tec
888def test_IONEX_connection(location):
889 """
890## =============================================================================
891##
892## Determines whether the machine is connected to the internet or not. Also
893## used to determine whether a given IONEX file exists.
894##
895## =============================================================================
896##
897## Inputs:
898## location type = string website/online file to attempt to access
899##
900## Returns:
901## True if the website/online file is accessible and False if not
902##
903## =============================================================================
904 """
905 try:
906 ok = os.system('curl -u anonymous:casa-feedback@nrao.edu --ftp-ssl-reqd -Isf -o/dev/null '+location)
907 if ok==0:
908 return True
909 else:
910 print(location+' unavailable! (or no internet?)')
911 return False
912 except:
913 print(location+' unavailable! (or no internet?)')
914 return False
920def convert_MAPGPS_TEC(ms_name,mad_data_file,ref_time,ref_start,ref_end,plot_vla_tec,im_name):
921 """
922## =============================================================================
923##
924## This opens the MAPGPS Data table and selects a subset of TEC/DTEC values
925## within a 15 deg square of the VLA. This then plots the zenith TEC/DTEC at the
926## VLA site and makes the TEC map for use at the C++ level. We chose to deal
927## with the MAPGPS data in this separate fashion because there are large
928## 'gaps' in the data where no TEC/DTEC values exist. Consequently, we use the
929## filled in CASA table to produce a TEC map and can not simply
930## concatenate arrays.
931##
932## =============================================================================
933##
934## Inputs:
935## ms_name type = string Name of the measurement set for which to
936## acquire TEC/DTEC data
937## mad_data_file type = string Name of the MAPGPS TEC/DTEC data table
938## ref_time type = float Reference time (s) for setting the
939## coordinates, UT 0 on the first day
940## plot_vla_tec type = boolean When True, this will open a plot of the
941## interpolated TEC/DTEC at the VLA.
942## im_name type = string Name of the output TEC Map optionally
943## specified by the user
944##
945## Returns:
946## Opens a plot showing the zenith TEC/DTEC at the VLA (if plot_vla_tec=True)
947## and the name of the CASA image file containing the TEC map.
948##
949## =============================================================================
950 """
951 ## Only retrieve data in a 15x15 deg. patch centered (more or less) at the VLA
952 tb.open(mad_data_file+'.tab')
953 st0=tb.query('GDLAT>19 && GDLAT<49 && GLON>-122 && GLON<-92',
954 ## If you want ALL the data to make a global map, use the line below:
955 #st0=tb.query('GDLAT>-90. && GDLAT<90. && GLON>-180. && GLON<180',
956 name='tecwindow')
958 utimes=pylab.unique(st0.getcol('UT1_UNIX'))
959 ulat=pylab.unique(st0.getcol('GDLAT'))
960 ulong=pylab.unique(st0.getcol('GLON'))
962 points_lat=len(ulat)
963 points_long=len(ulong)
964 num_maps=len(utimes)
966 ## Initialize the array which will be used to make the image
967 tec_array=pylab.zeros((2,points_long,points_lat,num_maps),dtype=float)
969 minlat=min(ulat)
970 minlong=min(ulong)
972 print('rows'+str(len(utimes)))
974 itime=0
975 for t in utimes:
976 st1=st0.query('UT1_UNIX=='+str(t),name='bytime')
977 n=st1.nrows()
978 if itime%100==0:
979 print(str(itime)+' '+str(n))
980 ilong=st1.getcol('GLON')-minlong
981 ilat=st1.getcol('GDLAT')-minlat
982 itec=st1.getcol('TEC')
983 idtec=st1.getcol('DTEC')
984 for i in range(n):
985 tec_array[0,int(ilong[i]),int(ilat[i]),itime]=itec[i]
986 tec_array[1,int(ilong[i]),int(ilat[i]),itime]=idtec[i]
987 st1.close()
989 ## Simply interpolate to cull as many zeros as possible
990 ## (median of good neighbors, if at least four of them)
991 thistec_array=tec_array[:,:,:,itime].copy()
992 thisgood=thistec_array[0]>0.0
993 for i in range(1,points_long-1):
994 for j in range(1,points_lat-1):
995 if not thisgood[i,j]:
996 mask=thisgood[(i-1):(i+2),(j-1):(j+2)]
997 if pylab.sum(mask)>4:
998 #print(str(itime)+' '+str(i)+' '+str(j)+str(pylab.sum(mask)))
999 tec_array[0,i,j,itime]=pylab.median(thistec_array[0,(i-1):(i+2),(j-1):(j+2)][mask])
1000 tec_array[1,i,j,itime]=pylab.median(thistec_array[1,(i-1):(i+2),(j-1):(j+2)][mask])
1001 itime+=1
1003 st0.close()
1004 tb.close()
1006 ztec_value(-107.6184,34.0790,points_long,points_lat,minlong,minlat,1,\
1007 1,5,ref_start,ref_end,int(num_maps),tec_array,plot_vla_tec)
1008 ## ref_time + 150 accounts for the fact that the MAPGPS map starts at 00:02:30 UT, not 00:00:00 UT
1009 if im_name == '':
1010 prefix = ms_name
1011 else:
1012 prefix = im_name
1013 CASA_image = make_image(prefix,minlong,minlat,ref_time+150.0,1,1,5*60,tec_array[0],'MAPGPS',appendix = '.MAPGPS_TEC')
1015 return CASA_image
1021def make_CASA_table(file_name):
1022 """
1023## =============================================================================
1024##
1025## Makes a CASA table for the MAPGPS file for a given day.
1026## It requires the 'madrigal.head' file be in the working directory.
1027##
1028## =============================================================================
1029##
1030## Inputs:
1031## file_name type = string This is the full MAPGPS file name
1032## MAPGPS form: 'gps'+'yy'+'mm'+'dd'+'.###'
1033## ex: August 17, 2011 is 'gps110817g.001'
1034##
1035## =============================================================================
1036 """
1037 os.system('rm -rf '+file_name+'.tab')
1038 tb.fromascii(tablename=file_name+'.tab',asciifile=file_name+'.txt',headerfile='madrigal.head',sep=" ",
1039 commentmarker='^YEAR MONT')
1041 # this will work when tb.fromascii is fixed...
1042 # (making the headerfile unnecessary)
1043 #tb.fromascii('madtest.tab','madtest.txt',sep=" ",
1044 # columnnames=['YEAR','MONTH','DAY','HOUR','MIN','SEC',
1045 # 'UT1_UNIX','UT2_UNIX','RECNO',
1046 # 'GDLAT','GLON','TEC','DTEC'],
1047 # datatypes=['I','I','I','I','I','I','I','I','I','R','R','R','R'],
1048 # autoheader=F,commentmarker='^YEAR MONT')