MIGHTI Level 2 Update

22 Apr 2015

from IPython.display import Image
Image(url='http://remote2.ece.illinois.edu/~bhardin2/MIGHTI/MIGHTI_blockdiag.png')

What to show? Check my notebook. Stuff below is from last week.

Copied from MIGHTI_L21_to_L22:

%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_20150420/L21/ICON_MIGHTI_A*.npz')
Bfns = glob.glob('/home/bhardin2/MIGHTI/Englandfiles/AB_20150420/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)


# Only use the first few files
Afns = Afns[:40]
Bfns = Bfns[:40]
# Only use the last few files
#Afns = Afns[-25:]
#Bfns = Bfns[-26:]

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
        npzfile.close()    
    
    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]])
            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 = x[0]
            v = x[1]
            
            # Remove co-rotation
            # (1) Convert location to ecef and determine rho (axial distance)
            lla_A = [lat_A[altbin,k_A], lon_A[altbin,k_A], alt_A[altbin,k_A]]
            ecef_xyz_A = ICON.wgs84_to_ecef(lla_A)
            rho = sqrt(ecef_xyz_A[0]**2 + ecef_xyz_A[1]**2)
            corot_vel = 2*pi*rho/(24*60*60)*1e3 # m/s
            # TEMP: THIS IS INTENTIONALLY WRONG (CM/S INSTEAD OF M/S)
            u = u - corot_vel*1e-2 # TODO: FIX THIS, REMOVE 1e-2
            
            U[altbin,k_A] = u
            V[altbin,k_A] = v

def plot_result(X,az=-70,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_20150420/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)

# 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))
plot_result(U)

True Zonal Wind from GCM

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

Difference (Recovered-True) Zonal Wind

plot_result(U-u_interp,vminmax=(-30,30))

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)

A vertical slice of the previous plot

i = 22

figure(figsize=(4,4))

subplot(1,2,1)
Urec = U[:,i]
Ugcm = u_interp[:,i]
alt = alt_A[:,i]
plot(Urec,alt,'r.-',label='Recovered')
plot(Ugcm,alt,'k-',label='From GCM file')
legend(loc='lower left',prop={'size':8});
xlabel('Zonal Wind [m/s]')
ylabel('Altitude [km]')

subplot(1,2,2)
Vrec = V[:,i]
Vgcm = v_interp[:,i]
alt = alt_A[:,i]
plot(Vrec,alt,'r.-',label='Recovered')
plot(Vgcm,alt,'k-',label='From GCM file')
#legend(loc='best',prop={'size':8})
xlabel('Meridional Wind [m/s]')

<matplotlib.text.Text at 0x7f1887406350>

Magnitude of wind

i = 18

figure(figsize=(2,4))

Urec = sqrt(U[:,i]**2 + V[:,i]**2)
Ugcm = sqrt(u_interp[:,i]**2 + v_interp[:,i]**2)
alt = alt_A[:,i]
plot(Urec,alt,'r.-',label='Recovered')
plot(Ugcm,alt,'k-',label='From GCM file')
xlim((0,nanmax(Ugcm)*1.3))
ylim((120,320))

legend(loc='best',prop={'size':8});
xlabel('Wind Speed [m/s]');
ylabel('Altitude [km]');

############# Look in Scott's file to see what u,v he actually used ###############

# Find source of the L2.1 file
npzfn = Afns[i]
npzfile = load(npzfn)
UIUCL1_src = npzfile['source_files'][0]
npzfile.close()

# Find source of that source, which is the original England file
npzfile = load(UIUCL1_src)
ncfn_src = npzfile['source_files'][0]
npzfile.close()

ncfn = ncfn_src

ncfile = Dataset(ncfn,'r')
vs = ncfile.variables.keys()
vs.sort()
#vs
print 'Looking at %s' % ncfn

Looking at /home/bhardin2/MIGHTI/Englandfiles/AB_20150420/England_raw/ICON_MIGHTI_A_ray_UT_555.000.nc

j = 1 # Use middle column only, so as not to care about horizontal resolution.

mighti_zonal_wind = ncfile.variables['MIGHTI_ZONAL_WIND'][0,:,:,j]
mighti_merid_wind = ncfile.variables['MIGHTI_MERID_WIND'][0,:,:,j]
mighti_tang_alts = ncfile.variables['MIGHTI_TANGENT_ALTITUDES_START'][0,:,j]


tang_point_along_ray_start = ncfile.variables['MIGHTI_TANGENT_POINT_ALONG_RAY_START'][0,:,j] # integer index

ny = shape(mighti_zonal_wind)[1]
nc_u = zeros(ny)
nc_v = zeros(ny)
for ii in range(ny):
    kstar = tang_point_along_ray_start[ii]
    nc_u[ii] = mighti_zonal_wind[kstar,ii]
    nc_v[ii] = mighti_merid_wind[kstar,ii]

figure(figsize=(4,4))

subplot(1,2,1)
Urec = U[:,i]
Ugcm = u_interp[:,i]
alt = alt_A[:,i]
plot(Urec,alt,'r.-',label='Recovered')
plot(Ugcm,alt,'k-',label='From GCM file')
plot(nc_u,mighti_tang_alts,'k--',label='From Scott\'s file') # TODO: is the altitude right???
legend(loc='best',prop={'size':8})
xlabel('Zonal Wind [m/s]')
ylabel('Altitude [km]')


subplot(1,2,2)
Vrec = V[:,i]
Vgcm = v_interp[:,i]
alt = alt_A[:,i]
plot(Vrec,alt,'r.-',label='Recovered')
plot(Vgcm,alt,'k-',label='From GCM file')
plot(nc_v,mighti_tang_alts,'k--',label='From Scott\'s file') # TODO: is the altitude right???
legend(loc='best',prop={'size':8});
xlabel('Meridional Wind [m/s]')

<matplotlib.text.Text at 0x7f18855f12d0>

What is the problem? Possibilities:

• Earth corotation. Correction is hundreds of m/s in zonal wind.
• Problem reading TIME-GCM file
• Problem with L2.2 algorithm

My extraction from the GCM matches Scott's recorded values (except for unexplained offset in zonal wind).

figure(figsize=(4,4))

subplot(1,2,1)
Urec = U[:,i]
Ugcm = u_interp[:,i]
alt = alt_A[:,i]
#plot(Urec,alt,'r.-',label='Recovered')
plot(Ugcm,alt,'k-',linewidth=2,label='From GCM file')
plot(nc_u,mighti_tang_alts,'b-',label='From Scott\'s file') # TODO: is the altitude right???
legend(loc='best',prop={'size':8});
xlabel('Zonal Wind [m/s]')
ylabel('Altitude [km]')


subplot(1,2,2)
Vrec = V[:,i]
Vgcm = v_interp[:,i]
alt = alt_A[:,i]
#plot(Vrec,alt,'r.-',label='Recovered')
plot(Vgcm,alt,'k-',linewidth=2,label='From GCM file')
plot(nc_v,mighti_tang_alts,'b-',label='From Scott\'s file') # TODO: is the altitude right???
legend(loc='best',prop={'size':8});
xlabel('Meridional Wind [m/s]')

<matplotlib.text.Text at 0x7f1889219590>

Backup slide: L2.1 LoS winds match OK

fnlist = concatenate((Afns,Bfns))

nsx = 6
nsy = len(fnlist)/nsx + 1
figure(figsize=(nsx*2,4*nsy))

for i in range(len(fnlist)):
    
    fn = fnlist[i]
    subplot(nsy,nsx,i+1)

    npzfile = load(fn)
    tang_alt = npzfile['tang_alt']
    tang_lat = npzfile['tang_lat']
    tang_lon = npzfile['tang_lon']
    los_wind = npzfile['los_wind']
    time = npzfile['time']
    datetime_created = npzfile['datetime_created']
    source_files = npzfile['source_files']
    npzfile.close()

    dirlist = fn.split('/')
    dirlist[-2] = 'Truth_raw'
    truthfn = '/'.join(dirlist)
    npzfile = load(truthfn)
    truth_tang_alt = npzfile['tang_alts']
    truth_winds = npzfile['truth_winds']
    npzfile.close()

    fns = fn.split('/')[-1][12:-4]
    title('%s\n(%.1f N,%.1f E)' % (fns,tang_lat[50],tang_lon[50]))
    plot(los_wind,tang_alt,'r.-',label='Recovered')
    plot(truth_winds,truth_tang_alt,'k-',label='Truth') # TODO: is the altitude right???
    legend(loc='best',prop={'size':8})
    xlabel('LoS Wind [m/s]')
    ylabel('Altitude [km]')
    #xlim((-80,80))
    
tight_layout()
savefig('/home/bhardin2/public_html/MIGHTI/MIGHTI_LoS_compare.png',dpi=80)

import subprocess
import shutil
nb_name = 'MIGHTI_L2_2015Apr22.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/')