MIGHTI A/B Simulation

From level 2.1 to level 2.2

25 Feb 2015

This notebook will read in Scott England's observations file taken from a long-term TIME-GCM run, and look at how to match up MIGHTI A and B

%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

# 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_20150221/England_raw/ICON_MIGHTI_A_*.nc')
Afns.sort()
Bfns = glob.glob('/home/bhardin2/MIGHTI/Englandfiles/AB_20150221/England_raw/ICON_MIGHTI_B_*.nc')
Bfns.sort()

def quick_read(ncfn, j=5):
    '''
    Return the time and tangent location of each ray in this exposure.
    (j indexes which of the 11 columns of Scott's simulation to use)
    '''    
    ncfile = Dataset(ncfn,'r')
    
    tang_alt_start = ncfile.variables['MIGHTI_TANGENT_ALTITUDES_START'][0,:,:]
    tang_alt_end = ncfile.variables['MIGHTI_TANGENT_ALTITUDES_END'][0,:,:]
    tang_lat_start = ncfile.variables['MIGHTI_TANGENT_LATITUDES_START'][0,:,:]
    tang_lat_end = ncfile.variables['MIGHTI_TANGENT_LATITUDES_END'][0,:,:]
    tang_lon_start = ncfile.variables['MIGHTI_TANGENT_LONGITUDES_START'][0,:,:]
    tang_lon_end = ncfile.variables['MIGHTI_TANGENT_LONGITUDES_END'][0,:,:]
    

    tang_alts = (tang_alt_start[:,j] + tang_alt_end[:,j])/2
    tang_lats = (tang_lat_start[:,j] + tang_lat_end[:,j])/2
    tang_lons = (tang_lon_start[:,j] + tang_lon_end[:,j])/2
    
    t = float(ncfn.split('/')[-1].split('_')[-1][:-3])

    return t, tang_lats, tang_lons, tang_alts
    

fn = Afns[10]
t, lats, lons, alts = quick_read(fn) 
plot(lats,lons,'k.')

[<matplotlib.lines.Line2D at 0x7fd082769c90>]

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

tvec = [quick_read(fn)[0] for fn in Bfns]
keydict = dict(zip(Bfns, tvec))
Bfns.sort(key=keydict.get)

t, lats, lons, alts = quick_read(Afns[0]) 

Atvec = zeros(len(Afns))
Aaltvec = zeros((len(Afns),len(alts))) # assume all exposures have the same tang alt bins
Alatvec = zeros((len(Afns),len(alts))) 
Alonvec = zeros((len(Afns),len(alts))) 
for i in range(len(Afns)):
    fn = Afns[i]
    t, lats, lons, alts = quick_read(fn) 
    Atvec[i] = t
    Alatvec[i,:] = lats
    Alonvec[i,:] = lons
    Aaltvec[i,:] = alts

t, lats, lons, alts = quick_read(Bfns[0]) 

Btvec = zeros(len(Bfns))
Baltvec = zeros((len(Bfns),len(alts))) # assume all exposures have the same tang alt bins
Blatvec = zeros((len(Bfns),len(alts))) 
Blonvec = zeros((len(Bfns),len(alts))) 
for i in range(len(Bfns)):
    fn = Bfns[i]
    t, lats, lons, alts = quick_read(fn) 
    Btvec[i] = t
    Blatvec[i,:] = lats
    Blonvec[i,:] = lons
    Baltvec[i,:] = alts

Scott:

• TIME-GCM
• ICON geometry and orbit
• 17 minutes at 30-sec cadence

Brian:

• Simulate L1 data products for each exposure (eventually Kenneth)
• Convert L1 to L2.1 (interferograms --> LoS winds)
• Convert L2.1 to L2.2 (LoS winds --> cardinal winds)

# Define the Basemap
fig = plt.figure(figsize=(4,4))
meanlon = mean((mean(Alonvec),mean(Blonvec)))
meanlat = mean((mean(Alatvec),mean(Blatvec)))
m = Basemap(projection='nsper',lon_0=meanlon,lat_0=meanlat,satellite_height=20000e3,resolution='l')
m.drawcoastlines(linewidth=0.5)
m.drawlsmask(land_color='none',ocean_color=array([187,235,251])/255.,lakes=True)
# draw parallels and meridians.
m.drawparallels(np.arange(-90.,120.,10.))
m.drawmeridians(np.arange(0.,420.,10.))
m.drawmapboundary(fill_color='coral')

altbin = 0
xa,ya = m(Alonvec[:,altbin], Alatvec[:,altbin])
xb,yb = m(Blonvec[:,altbin], Blatvec[:,altbin])
m.plot(xa,ya,'k.',label='A',markersize=5)
m.plot(xb,yb,'g.',label='B',markersize=5)
legend(loc='best',numpoints=1,prop={'size':7})
title('Tangent Points for lowest ray (~90 km)')

savefig('/home/bhardin2/public_html/MIGHTI/AB_tracks.png',dpi=200)

# Define the Basemap
fig = plt.figure(figsize=(4,4))
meanlon = mean((mean(Alonvec),mean(Blonvec)))
meanlat = mean((mean(Alatvec),mean(Blatvec)))
m = Basemap(projection='nsper',lon_0=meanlon,lat_0=meanlat,satellite_height=20000e3,resolution='l')
m.drawcoastlines(linewidth=0.5)
m.drawlsmask(land_color='none',ocean_color=array([187,235,251])/255.,lakes=True)
# draw parallels and meridians.
m.drawparallels(np.arange(-90.,120.,10.))
m.drawmeridians(np.arange(0.,420.,10.))
m.drawmapboundary(fill_color='coral')

altbin = -1
xa,ya = m(Alonvec[:,altbin], Alatvec[:,altbin])
xb,yb = m(Blonvec[:,altbin], Blatvec[:,altbin])
m.plot(xa,ya,'k.',label='A',markersize=5)
m.plot(xb,yb,'g.',label='B',markersize=5)
legend(loc='best',numpoints=1,prop={'size':7});
title('Tangent Points for highest ray (~300 km)');

#savefig('/home/bhardin2/public_html/MIGHTI/AB_tracks.png',dpi=200)

<matplotlib.text.Text at 0x7fd0760fd110>

# Define the Basemap
fig = plt.figure(figsize=(4,4))
meanlon = mean((mean(Alonvec),mean(Blonvec)))
meanlat = mean((mean(Alatvec),mean(Blatvec)))
m = Basemap(projection='nsper',lon_0=meanlon,lat_0=meanlat,satellite_height=20000e3,resolution='l')
m.drawcoastlines(linewidth=0.5)
m.drawlsmask(land_color='none',ocean_color=array([187,235,251])/255.,lakes=True)
# draw parallels and meridians.
m.drawparallels(np.arange(-90.,120.,10.))
m.drawmeridians(np.arange(0.,420.,10.))
m.drawmapboundary(fill_color='coral')

altbin = 0
xa,ya = m(Alonvec.flatten(), Alatvec.flatten())
xb,yb = m(Blonvec.flatten(), Blatvec.flatten())
m.plot(xa,ya,'k.',label='A',markersize=1)
m.plot(xb,yb,'g.',label='B',markersize=1)
legend(loc='best',numpoints=1,prop={'size':7})
title('Tangent Points of all rays')

#savefig('/home/bhardin2/public_html/MIGHTI/AB_tracks.png',dpi=200)

<matplotlib.text.Text at 0x7fd07b771a50>

t, lats, lons, alts = quick_read(Afns[0]) 

lAtvec = zeros(len(Afns))
lAaltvec = zeros((len(Afns),len(alts))) # assume all exposures have the same tang alt bins
lAlatvec = zeros((len(Afns),len(alts))) 
lAlonvec = zeros((len(Afns),len(alts))) 
rAtvec = zeros(len(Afns))
rAaltvec = zeros((len(Afns),len(alts))) # assume all exposures have the same tang alt bins
rAlatvec = zeros((len(Afns),len(alts))) 
rAlonvec = zeros((len(Afns),len(alts))) 
for i in range(len(Afns)):
    fn = Afns[i]
    t, lats, lons, alts = quick_read(fn,j=0) 
    lAtvec[i] = t
    lAlatvec[i,:] = lats
    lAlonvec[i,:] = lons
    lAaltvec[i,:] = alts
    t, lats, lons, alts = quick_read(fn,j=10) 
    rAtvec[i] = t
    rAlatvec[i,:] = lats
    rAlonvec[i,:] = lons
    rAaltvec[i,:] = alts

Fun Fact: Effect of Horizontal FOV $\approx$ Effect of smearing due to S/C velocity

# Define the Basemap
fig = plt.figure(figsize=(4,4))
meanlon = mean((mean(Alonvec),mean(Blonvec)))
meanlat = mean((mean(Alatvec),mean(Blatvec)))
m = Basemap(projection='nsper',lon_0=meanlon,lat_0=meanlat,satellite_height=1000e3,resolution='l')
m.drawcoastlines(linewidth=0.5)
m.drawlsmask(land_color='none',ocean_color=array([187,235,251])/255.,lakes=True)
# draw parallels and meridians.
m.drawparallels(np.arange(-90.,120.,10.))
m.drawmeridians(np.arange(0.,420.,10.))
m.drawmapboundary(fill_color='coral')

altbin = 0
xa,ya = m(rAlonvec[10,:].flatten(), rAlatvec[10,:].flatten())
xb,yb = m(lAlonvec[10,:].flatten(), lAlatvec[10,:].flatten())
m.plot(concatenate((xa,xb)),concatenate((ya,yb)),'g.',label='A Exp 1',markersize=2)
xa,ya = m(rAlonvec[11,:].flatten(), rAlatvec[11,:].flatten())
xb,yb = m(lAlonvec[11,:].flatten(), lAlatvec[11,:].flatten())
m.plot(concatenate((xa,xb)),concatenate((ya,yb)),'k.',label='A Exp 2',markersize=2)

legend(loc='upper right',numpoints=1,prop={'size':7})
title('Tangent Points and FOV of\nconsecutive MIGHTI A Exposures')

#savefig('/home/bhardin2/public_html/MIGHTI/AB_tracks.png',dpi=200)

<matplotlib.text.Text at 0x7fd07b158a10>

MIGHTI A/B Lineup

figure(figsize=(6,6))
plot(Alonvec[:-16,:].flatten(), Aaltvec[:-16,:].flatten(),'k.',markersize=2,label='A');
plot(Blonvec[8:,:].flatten(), Baltvec[8:,:].flatten(),'g.',markersize=2,label='B');
legend(loc='best',prop={'size':8},numpoints=1)
xlabel('Longitude [deg]')
ylabel('Altitude [km]')
grid()

figure(figsize=(6,6))
plot(Alatvec[:-16,:].flatten(), Aaltvec[:-16,:].flatten(),'k.',markersize=2,label='A');
plot(Blatvec[8:,:].flatten(), Baltvec[8:,:].flatten(),'g.',markersize=2,label='B');
legend(loc='best',prop={'size':8},numpoints=1)
xlabel('Latitude [deg]')
ylabel('Altitude [km]')
grid()

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

From L2.1 analysis: line of sight wind from A and B

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)

vmin = -100
vmax = 100
figure()
ZA = ma.masked_array(los_wind_A, mask = isnan(los_wind_A))
pcolor(lon_A,alt_A,ZA,vmin=vmin,vmax=vmax)
xlim((min((lon_A.min(),lon_B.min())),max((lon_A.max(),lon_B.max()))))
ylim((90,320))
xlabel('Longitude [deg]')
ylabel('Altitude [km]')
title('MIGHTI A LOS wind')
cb = colorbar()
cb.set_label('LOS velocity [m/s]')

figure()
ZB = ma.masked_array(los_wind_B, mask = isnan(los_wind_B))
pcolor(lon_B,alt_B,ZB,vmin=vmin,vmax=vmax)
xlim((min((lon_A.min(),lon_B.min())),max((lon_A.max(),lon_B.max()))))
ylim((90,320))
xlabel('Longitude [deg]')
ylabel('Altitude [km]')
title('MIGHTI B LOS wind')
cb = colorbar()
cb.set_label('LOS velocity [m/s]')

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()    

Same plots in 3D:

plot_FOVs()

# 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[altbin,k_A] = x[0]
            V[altbin,k_A] = x[1]

def plot_result(X,az=-50,ze=30,rstride=1,cstride=1):
    '''
    X is either U or V
    '''
    
    lw = 2
    vmin = nanmin(X)
    vmax = nanmax(X)
    alpha = 1.0

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

    ZA = (X-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)
    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('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} $$

L2.2 Zonal Wind

plot_result(U)

L2.2 Meridional Wind

plot_result(V)

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