MIGHTI Update

• Level 2.1 - LoS winds
  • Improved accuracy
• Level 2.2 - Cardinal winds
  • Updates on algorithm development and validation

24 March 2015

Two updates to show: local horz projection calculation was wrong, so results got better. Results of comparison to truth cardinal winds in level 2.2

%pylab inline
from Scientific.IO.NetCDF import NetCDFFile as Dataset
from datetime import datetime, timedelta
from pyglow import pyglow
from mpl_toolkits.basemap import Basemap
import os
import glob
from scipy.interpolate import griddata
import subprocess
import MIGHTI
import ICON

# Make default figure size bigger
matplotlib.rcParams['savefig.dpi'] = 150
matplotlib.rcParams['figure.figsize'] = (4,3)
matplotlib.rcParams['font.size'] = 8

Populating the interactive namespace from numpy and matplotlib

WARNING: pylab import has clobbered these variables: ['datetime', 'griddata']
`%matplotlib` prevents importing * from pylab and numpy

Level 2.1 - LoS wind

• ~3 m/s discrepancy between true and recovered profiles.
• Scott and Brian narrowed it down to a bug in the forward model
• What we had before:

from IPython.display import Image
Image(url='http://remote2.ece.illinois.edu/~bhardin2/MIGHTI/ICON_MIGHTI_A_ray_UT_15sec_Spherical_Sym_L21_red.png',
     width=300,height=300)

• After fixing the forward model, the results seem better

from IPython.display import Image
Image(url='http://remote2.ece.illinois.edu/~bhardin2/MIGHTI/ICON_MIGHTI_A_ray_UT_Spherical_Symmetry_15sec_L21_red.png',
     width=300,height=300)

Copied from MIGHTI_L21_to_L22:

Just so I can make one slideshow

%pylab inline
from Scientific.IO.NetCDF import NetCDFFile as Dataset
from datetime import datetime, timedelta
from pyglow import pyglow
from mpl_toolkits.basemap import Basemap
import os
import glob
from scipy.interpolate import griddata
import subprocess
import MIGHTI
import ICON
from mpl_toolkits.mplot3d import Axes3D

# Make default figure size bigger
matplotlib.rcParams['savefig.dpi'] = 150
matplotlib.rcParams['figure.figsize'] = (4,3)
matplotlib.rcParams['font.size'] = 8

Populating the interactive namespace from numpy and matplotlib

WARNING: pylab import has clobbered these variables: ['datetime', 'griddata']
`%matplotlib` prevents importing * from pylab and numpy

Afns = glob.glob('/home/bhardin2/MIGHTI/Englandfiles/AB_20150324/L21/ICON_MIGHTI_A*.npz')
Bfns = glob.glob('/home/bhardin2/MIGHTI/Englandfiles/AB_20150324/L21/ICON_MIGHTI_B*.npz')

# Sort the files by time
tvec = [load(fn)['time'] for fn in Afns]
keydict = dict(zip(Afns, tvec))
Afns.sort(key=keydict.get)

tvec = [load(fn)['time'] for fn in Bfns]
keydict = dict(zip(Bfns, tvec))
Bfns.sort(key=keydict.get)

def load_l21(fns):
    '''
    Crawl the list of files given, and return N x K matrices of:
        latitude
        longitude
        altitude
        los_wind
        local_az
    where N = number of altitude bins 
          K = number of time samples (i.e., len(fns))
    '''
    N_altbins = len(load(fns[0])['tang_alt'])
    lat_A = zeros((N_altbins, len(fns)))
    lon_A = zeros((N_altbins, len(fns)))
    alt_A = zeros((N_altbins, len(fns)))
    los_wind_A = zeros((N_altbins, len(fns)))
    local_az_A = zeros((N_altbins, len(fns)))
    time_A = zeros(len(fns))

    for i in range(len(fns)):
        fn = fns[i]
        npzfile = load(fn)

        # Extract all information from the L2.1 file
        emission_color = str(npzfile['emission_color'])
        icon_alt = npzfile['icon_alt']
        icon_lat = npzfile['icon_lat']
        icon_lon = npzfile['icon_lon']
        tang_alt = npzfile['tang_alt']
        tang_lat = npzfile['tang_lat']
        tang_lon = npzfile['tang_lon']
        local_az = npzfile['local_az']
        los_wind = npzfile['los_wind']
        time = npzfile['time']
        datetime_created = npzfile['datetime_created']
        source_files = npzfile['source_files']

        lat_A[:,i] = tang_lat
        lon_A[:,i] = tang_lon
        alt_A[:,i] = tang_alt
        los_wind_A[:,i] = los_wind
        local_az_A[:,i] = local_az
        #time_A[i] = time[0] # something wrong here
    return lat_A, lon_A, alt_A, los_wind_A, local_az_A

lat_A, lon_A, alt_A, los_wind_A, local_az_A = load_l21(Afns)
lat_B, lon_B, alt_B, los_wind_B, local_az_B = load_l21(Bfns)

def plot_FOVs(az=-50,ze=30,rstride=1,cstride=1):
    
    lw = 2
    vmin = -100
    vmax = 100
    alpha = 1.0

    fig = plt.figure()
    ax = fig.add_subplot(111, projection='3d')

    ZA = (los_wind_A-vmin)/(vmax-vmin)
    ZA = ma.masked_array(ZA, isnan(ZA))
    col = cm.jet(ZA)
    a = col[:,:,3]
    a[a!=0] = alpha
    col[:,:,3] = a
    h = ax.plot_surface(lon_A, lat_A, alt_A, cstride=cstride,rstride=rstride, \
                    cmap=cm.jet, facecolors=col, linewidth=1,
                    antialiased=False)
    ax.plot(lon_A[0,:],lat_A[0,:],alt_A[0,:],'k-',linewidth=lw)
    ax.plot(lon_A[-1,:],lat_A[-1,:],alt_A[-1,:],'k-',linewidth=lw)
    ax.plot(lon_A[:,0],lat_A[:,0],alt_A[:,0],'k-',linewidth=lw)
    ax.plot(lon_A[:,-1],lat_A[:,-1],alt_A[:,-1],'k-',linewidth=lw)


    ZB = (los_wind_B-vmin)/(vmax-vmin)
    ZB = ma.masked_array(ZB, isnan(ZB))
    col = cm.jet(ZB)
    a = col[:,:,3]
    a[a!=0] = alpha
    col[:,:,3] = a
    #print a
    h = ax.plot_surface(lon_B, lat_B, alt_B, cstride=cstride,rstride=rstride, \
                    cmap=cm.jet, facecolors=col, linewidth=1,
                    antialiased=False)
    ax.plot(lon_B[0,:],lat_B[0,:],alt_B[0,:],'g-',linewidth=lw)
    ax.plot(lon_B[-1,:],lat_B[-1,:],alt_B[-1,:],'g-',linewidth=lw)
    ax.plot(lon_B[:,0],lat_B[:,0],alt_B[:,0],'g-',linewidth=lw)
    ax.plot(lon_B[:,-1],lat_B[:,-1],alt_B[:,-1],'g-',linewidth=lw)

    m = cm.ScalarMappable(cmap=cm.jet)
    m.set_array(ZA)
    m.set_clim((vmin,vmax))
    cb = plt.colorbar(m, shrink=0.5,aspect=5)
    cb.set_label('LoS velocity [m/s]')
    ax.view_init(ze, az)
    
    ax.set_xlabel('Longitude [deg]')
    ax.set_ylabel('Latitude [deg]')
    ax.set_zlabel('Altitude [km]')
    tight_layout()

    plt.show()    

# Simple idea: for each MIGHTI A measurement, match it with the closest MIGHTI B measurement (in longitude)
# Output variables:
U = nan*zeros(shape(lon_A))
V = nan*zeros(shape(lon_A))

sind = lambda(a): sin(a*np.pi/180.0)
cosd = lambda(a): cos(a*np.pi/180.0)
N_alts_A, N_times_A = shape(lon_A)
N_alts_B, N_times_B = shape(lon_B)
for altbin in range(N_alts_A):
    for k_A in range(N_times_A):
        # find time index of the B sample to the "left" (in longitude)
        k0_B = sum(lon_B[altbin,:] < lon_A[altbin,k_A])-1
        # find time index of the B sample to the "right" (in longitude)
        k1_B = N_times_B-sum(lon_B[altbin,:] > lon_A[altbin,k_A])
        if k0_B > -1 and k1_B < N_times_B:
            # Find closest MIGHTI B measurement
            k_B = argmin(abs(lon_B[altbin,:] - lon_A[altbin,k_A]))
            # Coordinate transformation
            y = array([los_wind_A[altbin,k_A], -los_wind_B[altbin,k_B]]) # negative so that y is los away
            az_A = local_az_A[altbin,k_A]
            az_B = local_az_B[altbin,k_B]
            A = array([[sind(az_A), cosd(az_A)],[sind(az_B), cosd(az_B)]])
            x = np.linalg.solve(A,y)
            U[altbin,k_A] = x[0]
            V[altbin,k_A] = x[1]

def plot_result(X,az=-50,ze=30,rstride=1,cstride=1,vminmax=None):
    '''
    X is either U or V
    '''
    
    lw = 2
    absmax = nanmax(abs(X))
    if vminmax is None:
        vmin = -absmax
        vmax = absmax
    else:
        vmin = vminmax[0]
        vmax = vminmax[1]
    alpha = 1.0
    cmap = cm.RdBu

    fig = plt.figure()
    ax = fig.add_subplot(111, projection='3d')

    ZA = (X-vmin)/(vmax-vmin)
    ZA = ma.masked_array(ZA, isnan(ZA))
    col = cmap(ZA)
    a = col[:,:,3]
    a[a!=0] = alpha
    col[:,:,3] = a
    h = ax.plot_surface(lon_A, lat_A, alt_A, cstride=cstride,rstride=rstride, \
                    cmap=cmap, facecolors=col, linewidth=1,
                    antialiased=False)
    ax.plot(lon_A[0,:],lat_A[0,:],alt_A[0,:],'k-',linewidth=lw)
    ax.plot(lon_A[-1,:],lat_A[-1,:],alt_A[-1,:],'k-',linewidth=lw)
    ax.plot(lon_A[:,0],lat_A[:,0],alt_A[:,0],'k-',linewidth=lw)
    ax.plot(lon_A[:,-1],lat_A[:,-1],alt_A[:,-1],'k-',linewidth=lw)
    ax.plot(lon_B[0,:],lat_B[0,:],alt_B[0,:],'g-',linewidth=lw)
    ax.plot(lon_B[-1,:],lat_B[-1,:],alt_B[-1,:],'g-',linewidth=lw)
    ax.plot(lon_B[:,0],lat_B[:,0],alt_B[:,0],'g-',linewidth=lw)
    ax.plot(lon_B[:,-1],lat_B[:,-1],alt_B[:,-1],'g-',linewidth=lw)

    m = cm.ScalarMappable(cmap=cmap)
    m.set_array(ZA)
    m.set_clim((vmin,vmax))
    cb = plt.colorbar(m, shrink=0.5,aspect=5)
    cb.set_label('Velocity [m/s]')
    ax.view_init(ze, az)
    
    ax.set_xlabel('Longitude [deg]')
    ax.set_ylabel('Latitude [deg]')
    ax.set_zlabel('Altitude [km]')
    tight_layout()

    plt.show()    

The simplest 2.1 to 2.2 conversion:

Assume that altitude bins of MIGHTI A match MIGHTI B (?). For each sample from MIGHTI A, find the closest sample in MIGHTI B. Then perform coordinate transformation (2x2 matrix multiplication).

$$ \begin{bmatrix} LoS~A \\ LoS~B \\ \end{bmatrix} = \begin{bmatrix} \sin(az_A) & \cos(az_A) \\ \sin(az_B) & \cos(az_B) \\ \end{bmatrix} \begin{bmatrix} u \\ v \\ \end{bmatrix} $$

Level 2.2 - Cardinal Winds

Compare with true values from GCM run

gcm_fn = '/home/bhardin2/MIGHTI/Englandfiles/AB_20150324/England_raw/sso076sept_mignmigtidegswm09.nc'
ncfile = Dataset(gcm_fn,'r')
vs = ncfile.variables.keys()
vs.sort()
#vs
# don't try to read "ne" says Scott

latvar = ncfile.variables['lat']
lonvar = ncfile.variables['lon']
altvar = ncfile.variables['ZG']
uvar = ncfile.variables['UN']
vvar = ncfile.variables['VN']
wvar = ncfile.variables['WN']
tvar = ncfile.variables['ut']

ti = 0 # decides which time in the GCM output to use
gcmlat = latvar.getValue() # deg north
gcmlon = lonvar.getValue() # deg east
gcmalt = altvar.getValue()[ti,:,:,:]*1e-5 # km
gcmu   = uvar.getValue()[ti,:,:,:]*1e-2 # m/s 
gcmv   = vvar.getValue()[ti,:,:,:]*1e-2 # m/s 
gcmw   = wvar.getValue()[ti,:,:,:]*1e-2 # m/s 
gcmt   = tvar.getValue()

print shape(gcmlat)
print shape(gcmlon)
print shape(gcmalt)
print shape(gcmu) # (97 alts x 72 lats x 144 lons)

(72,)
(144,)
(97, 72, 144)
(97, 72, 144)

# create meshed grid
gcmLAT = zeros(shape(gcmu))
for i in range(shape(gcmu)[0]):
    for k in range(shape(gcmu)[2]):
        gcmLAT[i,:,k] = gcmlat
gcmLON = zeros(shape(gcmu))
for i in range(shape(gcmu)[0]):
    for j in range(shape(gcmu)[1]):
        gcmLON[i,j,:] = gcmlon
gcmALT = gcmalt

lon_interp = lon_A.copy()
lat_interp = lat_A.copy()
alt_interp = alt_A.copy()
idx = lon_interp > 180.0
lon_interp[idx] -= 360.0

u_interp = zeros(shape(alt_interp))
v_interp = zeros(shape(alt_interp))
w_interp = zeros(shape(alt_interp))

gcm_pts = (gcmALT.flatten(),gcmLAT.flatten(),gcmLON.flatten())

for ii in range(shape(lon_interp)[0]):
    for jj in range(shape(lon_interp)[1]):
        alt = alt_interp[ii,jj]
        lat = lat_interp[ii,jj]
        lon = lon_interp[ii,jj]

        # first find lat index
        j = argmin(abs(gcmlat - lat))

        # find lon index
        k = argmin(abs(gcmlon - lon))

        # find alt index
        this_gcmalt = gcmalt[:,j,k]
        i = argmin(abs(this_gcmalt - alt))

        # read off values
        u_interp[ii,jj] = gcmu[i,j,k]
        v_interp[ii,jj] = gcmv[i,j,k]
        w_interp[ii,jj] = gcmw[i,j,k]

idxnan = isnan(U)
u_interp[idxnan] = nan
v_interp[idxnan] = nan
w_interp[idxnan] = nan

Recovered Zonal Wind

plot_result(U,vminmax=(-70,70))

True Zonal Wind from GCM

plot_result(u_interp,vminmax=(-70,70))

Difference (Recovered-True) Zonal Wind

plot_result(U-u_interp)

Recovered Meridional Wind

plot_result(V,vminmax=(-50,50))

True Meridional Wind from GCM

plot_result(v_interp,vminmax=(-50,50))

Difference (Recovered-True) Meridional Wind

plot_result(V-v_interp)

What is the problem? Possibilities:

• I haven't accounted for Earth corotation. (But shouldn't this effect be larger?)
• Problem reading TIME-GCM file
• Problem with L2.2 algorithm

import subprocess
import shutil
nb_name = 'MIGHTI_L2_2015Mar24.ipynb' # There's probably a way to get this programmatically
cmd = ['ipython', 'nbconvert', '--to', 'slides', nb_name, \
       '--reveal-prefix', '"http://cdn.jsdelivr.net/reveal.js/2.6.2"',\
      '--template','output_toggle']
# The --template command points to a custom template which opens and closes the input code cells with a click
subprocess.call(cmd)
html_name = '%s.slides.html' % (nb_name.split('.')[0])
shutil.copy(html_name, '/home/bhardin2/public_html/slideshows/')