ionex_download_and_display.py

"""
Donload and display
==================

Downloads IONEX data from IGS and displays a TEC map.

Uses the ionex library: https://github.com/gnss-lab/ionex

"""
import os
import sys
import ftputil
import numpy as np
import ionex

from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt


if len(sys.argv) < 4:
    print("Usage: ionex_download_and_display.py <year> <day> <hour>")
    sys.exit(0)    

year = int(sys.argv[1])
day = int(sys.argv[2])
hour = int(sys.argv[3])


year_str = '{:d}'.format(year)
day_str = '{:03d}'.format(day)
file_name = 'igsg{}0.{}i'.format(day_str,year_str[-2:])
zfile_name = file_name+'.Z'

if not file_name in os.listdir('.'):
    with ftputil.FTPHost('cddis.gsfc.nasa.gov', 'anonymous', '') as host:
        print("Conected ")
        host.download(host.path.join('gnss/products/ionex',year_str,day_str,zfile_name),
                      zfile_name)
        print("downloaded ", end=' ', flush=True)        
        os.system('uncompress '+  zfile_name)
        print('done')


found = False
with open(file_name) as file:
    inx = ionex.reader(file)
    for ionex_map in inx:
            h = ionex_map.epoch.hour
            if h == hour:
                la = ionex_map.grid.latitude;
                la_n = int((la.lat2 - la.lat1) // la.dlat + 1)
                lo = ionex_map.grid.longitude 
                lo_n = int((lo.lon2 - lo.lon1) // lo.dlon + 1)              
                tecmap = np.reshape(ionex_map.tec,(la_n, lo_n))
                found = True
                break
if not found:
    print("Map not found")
    sys.exit(0)

lats = np.linspace(la.lat1, la.lat2, la_n)
lons = np.linspace(lo.lon1, lo.lon2, lo_n)
lons, lats = np.meshgrid(lons,lats)
m = Basemap(
    projection= 'cyl', 
    lat_1= la.lat1, 
    lat_2= la.lat2, 
    lon_1= lo.lon1, 
    lon_2= lo.lon2,
    )
m.pcolormesh(lons, lats, tecmap,cmap='plasma',latlon=True)
m.drawparallels(np.linspace( -90, +90,  7))
m.drawmeridians(np.linspace(-180,+180, 13))
m.drawcoastlines()
plt.title("IGS VTEC Map at {}/{} {:2d}:00".format(year_str,day,h))
plt.tight_layout()
plt.show()