In [5]:

    

%matplotlib inline

import linecache
import sys
import gc
import pdb

import matplotlib
import matplotlib.pyplot as plt
import matplotlib.mlab as ml
from matplotlib.colors import from_levels_and_colors
import datetime
from dateutil.relativedelta import relativedelta
import re
import numpy as np
from math import *
import scipy.interpolate


    

In [6]:

    

# Exception handling, with line number and stuff

def PrintException():
    exc_type, exc_obj, tb = sys.exc_info()
    f = tb.tb_frame
    lineno = tb.tb_lineno
    filename = f.f_code.co_filename
    linecache.checkcache(filename)
    line = linecache.getline(filename, lineno, f.f_globals)
    print 'EXCEPTION IN ({}, LINE {} "{}"): {}'.format(filename, lineno, line.strip(), exc_obj)


    

In [7]:

    

def rotate_u_v_wind_array(u_wind, v_wind, brng):
    '''
    Rotate winds with given angle (which is a maths angle - positive counter-clockwise from positive x to positive y)
    Returns roated u and v masked array 
    '''
  
  #nan_mask_u_wind = np.ma.masked_array(np.array(u_wind[:,i,0:-1:10], dtype=float).flatten(), np.isnan(np.array(u_wind[:,i,0:-1:10], dtype=float).flatten()))
  #nan_mask_v_wind = np.ma.masked_array(np.array(v_wind[:,i,0:-1:10], dtype=float).flatten(), np.isnan(np.array(v_wind[:,i,0:-1:10], dtype=float).flatten()))
    nan_mask_u_wind = np.ma.masked_array(np.array(u_wind, dtype=float).flatten(), np.isnan(np.array(u_wind, dtype=float).flatten()))
    nan_mask_v_wind = np.ma.masked_array(np.array(v_wind, dtype=float).flatten(), np.isnan(np.array(v_wind, dtype=float).flatten()))
    
    #print nan_mask
 
    u_rot = nan_mask_u_wind * cos(brng) + nan_mask_v_wind * sin(brng)
    v_rot = -nan_mask_u_wind * sin(brng) + nan_mask_v_wind * cos(brng)
    
    return u_rot,v_rot

def u_v_winds(wind_direction, wind_speed):
    wind_rad = np.radians(wind_direction)
    u_wind=-((wind_speed)*np.sin(wind_rad))
    v_wind=-((wind_speed)*np.cos(wind_rad))
    return u_wind,v_wind

def grid_data_cs(pressure, distance, param):
    xi=np.linspace(0, max(distance), 200) 
 #yi=np.linspace(np.nanmin(pressure), np.nanmax(pressure), 500) 
    yi=np.linspace(5000, 100000, 50) # Points for pressure interpolation
    
    try:
        zi = ml.griddata(distance, pressure,param,xi, yi, interp='linear')
        #zi = scipy.interpolate.griddata((distance, pressure),  param, (xi[None,:],yi[:,None]), method='linear')
    except Exception, e:
        PrintException()
    return xi,yi,zi 
    

    

def station_info_search(stat):
    for line in station_metadata:
         if "%s" % stat in line: 
             st = line[2].lower().title().replace('_',' ')
             lo = float(line[3])
             la = float(line[4])
    return st,la,lo

def calculate_distance_from_first_station(stat, first_station_lon, first_station_lat, station_lat, station_lon):

    fslat_rad = radians(first_station_lat)
    fslon_rad = radians(first_station_lon)
    lat_rad = radians(station_lat)
    lon_rad = radians(station_lon)

 #Haversine Formula
    
    a = sin((lat_rad-fslat_rad)/2)**2 + cos(lat_rad) * cos(fslat_rad) * sin((lon_rad-fslon_rad)/2)**2
    c = 2 * atan2(sqrt(a), sqrt(1-a))
    d = 6371 * c

    return d

def variable_name_index_match(variable, variable_list):
    for key, value in variable_list.iteritems():   # iter on both keys and values
        if key.startswith('%s' % variable) and key.endswith('%s' % variable):
            arr_index_var=value 
    return arr_index_var

def variable_cat(var_index, station_list_cs):
    '''
    Load date binned radiosonde data, saved as output for each station from previous script, and gather
    info of stations in line, calculate dist from first station and concatente station data
    '''
    
    var_cat=[]
    distances=[]
    date_min_max=[]
    
    for stat in station_list_cs:
        load_file = np.load('/nfs/a90/eepdw/Data/Observations/Radiosonde_Numpy/Radiosonde_Cross_Section_'
                           'IND_SOUNDING_INTERP_MEAN_Climat_%s_%s_%s_%s.npz' 
                           % (date_min.strftime('%Y%m%d'), date_max.strftime('%Y%m%d'), delta, stat))
  
         #print load_file['date_bin_mean_all_dates_one_station'].shape

        if date_min_max ==[]:
            date_min_max=np.empty(load_file['min_max_date_bin'].shape)
        
        station_title, station_lon, station_lat = station_info_search(stat)
    
        dist_from_first_station = calculate_distance_from_first_station(stat, first_station_lon, first_station_lat, station_lat, station_lon)
    
        var_cat.append(load_file['date_bin_mean_all_dates_one_station'][:,var_index,:])
        distances.append(dist_from_first_station)
    
        date_min_max = np.where((load_file['min_max_date_bin']>date_min) & (load_file['min_max_date_bin']<date_max), load_file['min_max_date_bin'], date_min_max )
        
    dates_savefig_name=load_file['dates_for_plotting']
  
    return np.array(var_cat), np.array(distances, dtype=float), date_min_max, dates_savefig_name

def variable_cat_line(var_index, station_list_cs):
    var_cat=[]
    distances=[]
    date_min_max=[]
    
    for stat in station_list_cs:
     load_file = np.load('/nfs/a90/eepdw/Data/Observations/Radiosonde_Numpy/Radiosonde_Cross_Section_'
                        'IND_SOUNDING_INTERP_MEAN_Climat_%s_%s_%s_%s.npz' 
                        % (date_min.strftime('%Y%m%d'), date_max.strftime('%Y%m%d'), delta, stat))
  
     #print load_file['date_bin_mean_all_dates_one_station'].shape

     if date_min_max ==[]:
      date_min_max=np.empty(load_file['min_max_date_bin'].shape)
        
     station_title, station_lon, station_lat = station_info_search(stat)
    
     dist_from_first_station = calculate_distance_from_first_station(stat, first_station_lon, first_station_lat, station_lat, station_lon)
    
     #print dist_from_first_station
     #print load_file['date_bin_mean_all_dates_one_station'][:,var_index,:].shape
     var_cat.append(load_file['date_bin_mean_all_dates_one_station_single'][:,0,var_index])
     distances.append(dist_from_first_station)
     #pdb.set_trace()
     #if load_file['min_max_date_bin'].any() != np.NAN:
     #date_min_max=np.ma.masked_outside(load_file['min_max_date_bin'], date_min, date_max ).data
     date_min_max = np.where((load_file['min_max_date_bin']>date_min) & (load_file['min_max_date_bin']<date_max), load_file['min_max_date_bin'], date_min_max )
    #print np.array(var_cat).shape
    #print date_min_max
    return np.array(var_cat), np.array(distances, dtype=float), date_min_max

def PlotRadSondeTransNoWindLoop(i, xi,yi,zi, min_contour, max_contour, cmap, cbar_number_format):
    
    num_cols=32
   
    plt.figure(figsize=(14,8.5))
       
    clevs = np.linspace(min_contour, max_contour, num_cols)
    ticks = (np.arange(min_contour, max_contour+tick_interval,tick_interval))
    
 #cmap=plt.cm.RdBu_r
    
    cont = plt.contourf(xi,yi/100, zi, clevs, cmap=cmap, extend='both')

 #plt.clabel(cont_wind, fontsize=9, inline=1, fmt='%1d')
  
    cbar = plt.colorbar(cont, orientation='vertical', pad=0.05, extend='both')
 #cbar.set_label('$W m^{-2}$') 
    cbar.set_ticks(np.arange(min_contour, max_contour+tick_interval,tick_interval))
    cbar.set_ticklabels([cbar_number_format % i for i in ticks])
 
    plt.xlim(np.min(xi)-200,np.max(xi)+200)
    plt.gca().invert_yaxis()
    plt.gca().invert_xaxis()
    plt.ylabel('Pressure (hPa)')
    plt.xlabel('km from first station')
    
    return cont,cbar

def plot_rad_cs_winds(i, xi,yi,zi, min_contour, max_contour, xiw, yiw, uiw, viw, cmap):
    
    num_cols=32
   
    plt.figure(figsize=(14,8.5))
       
    clevs = np.linspace(min_contour, max_contour, num_cols)
    ticks = (np.arange(min_contour, max_contour+tick_interval,tick_interval))
    
 #cmap=plt.cm.RdBu_r
    
    cont = plt.contourf(xi,yi/100, zi, clevs, cmap=cmap, extend='both')
 
    cont_wind = plt.barbs(xiw,yiw/100, uiw, viw)
 #plt.clabel(cont_wind, fontsize=9, inline=1, fmt='%1d')
  
    cbar = plt.colorbar(cont, orientation='vertical', pad=0.05, extend='both')
 #cbar.set_label('$W m^{-2}$') 
    cbar.set_ticks(np.arange(min_contour, max_contour+tick_interval,tick_interval))
    cbar.set_ticklabels([cbar_number_format % i for i in ticks])
 
    plt.xlim(np.min(xi)-200,np.max(xi)+200)
    plt.gca().invert_yaxis()
    plt.gca().invert_xaxis()
    plt.ylabel('Pressure (hPa)')
    plt.xlabel('km from first station')
    
    return cont,cbar

def plot_rad_cs_winds_normalize(i, xi,yi,zi, min_contour, max_contour, xiw, yiw, uiw, viw, cmap):
   
    num_cols=32
    
    plt.figure(figsize=(14,8))
   
    clevs = np.linspace(min_contour, max_contour, num_cols)
    ticks = (np.arange(min_contour, max_contour+tick_interval,tick_interval))
    
 #cmap=plt.cm.RdBu_r
    
    midpoint=0
    midp = np.mean(np.c_[clevs[:-1], clevs[1:]], axis=1)
    vals = np.interp(midp, [min_contour, midpoint, max_contour], [0, 0.5, 1])
    cols = cmap(vals)
    clevs_extend = np.linspace(min_contour, max_contour,num_cols-2)
    cmap, norm = from_levels_and_colors(clevs_extend, cols, extend='both')

    cont = plt.contourf(xi,yi/100, zi, clevs, cmap=cmap, extend='both')
 
    cont_wind = plt.barbs(xiw,yiw/100, uiw, viw)
 #plt.clabel(cont_wind, fontsize=9, inline=1, fmt='%1d')
  
    cbar = plt.colorbar(cont, orientation='vertical', pad=0.05, extend='both')
 #cbar.set_label('$W m^{-2}$') 
    cbar.set_ticks(np.arange(min_contour, max_contour+tick_interval,tick_interval))
    cbar.set_ticklabels([cbar_number_format % i for i in ticks])
 
    plt.xlim(np.min(xi)-200,np.max(xi)+200)
    plt.gca().invert_yaxis()
    plt.gca().invert_xaxis()
    plt.ylabel('Pressure (hPa)')
    plt.xlabel('km from first station')
    
    return cont,cbar


def plot_rad_cs_winds_rot_and_unrot(i, xi,yi,zi, min_contour, max_contour, xiw, yiw, uiw, viw, uiw_rot, viw_rot, cmap, cbar_number_format):
    
    num_cols=32
   
    plt.figure(figsize=(14,8.5))
       
    clevs = np.linspace(min_contour, max_contour, num_cols)
    ticks = (np.arange(min_contour, max_contour+tick_interval,tick_interval))
    
 #cmap=plt.cm.RdBu_r
    
    cont = plt.contourf(xi,yi/100, zi, clevs, cmap=cmap, extend='both')
 
    cont_wind = plt.barbs(xiw,yiw/100, uiw, viw)
    
    cont_wind_rot = plt.quiver(xiw,yiw/100, uiw_rot, viw_rot,
               scale_units='inches',
               angles='uv',
               scale=40,
               headlength=1,
               headaxislength=1,
               width=0.004,
               alpha=0.8)
 #plt.clabel(cont_wind, fontsize=9, inline=1, fmt='%1d')
  
    cbar = plt.colorbar(cont, orientation='vertical', pad=0.05, extend='both')
 #cbar.set_label('$W m^{-2}$') 
    cbar.set_ticks(np.arange(min_contour, max_contour+tick_interval,tick_interval))
    cbar.set_ticklabels([cbar_number_format % i for i in ticks])
 
    plt.xlim(np.min(xi)-200,np.max(xi)+200)
    plt.gca().invert_yaxis()
    plt.gca().invert_xaxis()
    plt.ylabel('Pressure (hPa)')
    plt.xlabel('km from first station')
    
    return cont,cbar

def plot_rad_cs_winds_rot_and_unrot_normalize(i, xi,yi,zi, min_contour, max_contour, xiw, yiw, uiw, viw, uiw_rot, viw_rot, cmap, cbar_number_format):
    
    num_cols=32
   
    plt.figure(figsize=(14,8.5))
       
    clevs = np.linspace(min_contour, max_contour, num_cols)
    ticks = (np.arange(min_contour, max_contour+tick_interval,tick_interval))
    
 #cmap=plt.cm.RdBu_r
    
    midpoint=0
    midp = np.mean(np.c_[clevs[:-1], clevs[1:]], axis=1)
    vals = np.interp(midp, [min_contour, midpoint, max_contour], [0, 0.5, 1])
    cols = cmap(vals)
    clevs_extend = np.linspace(min_contour, max_contour,num_cols-2)
    cmap, norm = from_levels_and_colors(clevs_extend, cols, extend='both')
 #cmap=plt.cm.RdBu_r
    
    cont = plt.contourf(xi,yi/100, zi, clevs, cmap=cmap, extend='both')
 
    cont_wind = plt.barbs(xiw,yiw/100, uiw, viw)
    
    cont_wind_rot = plt.quiver(xiw,yiw/100, uiw_rot, viw_rot,
               scale_units='inches',
               angles='uv',
               scale=40,
               headlength=1,
               headaxislength=1,
               width=0.004,
               alpha=0.8)
 #plt.clabel(cont_wind, fontsize=9, inline=1, fmt='%1d')
  
    cbar = plt.colorbar(cont, orientation='vertical', pad=0.05, extend='both')
 #cbar.set_label('$W m^{-2}$') 
    cbar.set_ticks(np.arange(min_contour, max_contour+tick_interval,tick_interval))
    cbar.set_ticklabels([cbar_number_format % i for i in ticks])
 
    plt.xlim(np.min(xi)-200,np.max(xi)+200)
    plt.gca().invert_yaxis()
    plt.gca().invert_xaxis()
    plt.ylabel('Pressure (hPa)')
    plt.xlabel('km from first station')
    
    return cont,cbar


import re

station_list_search='/nfs/a90/eepdw/Data/Observations/Radiosonde_downloaded_from_NOAA_GUAN/igra-stations.txt'
station_metadata=[]
f = open(station_list_search,'r')
for line in f:
     line = line.strip()
     line=re.sub(r'([A-Z])\s([A-Z])', r'\1_\2',line)
     line=re.sub(r'([A-Z])\s\s([A-Z])', r'\1_\2',line)
     station_metadata.append(line.split())
f.close()

def LinePlots(variable_list_line, station_list_cs, distances):
        legendEntries=[]
        legendtext=[]
    
        variable='surface_pressure'
        var_index = variable_name_index_match(variable, variable_list_line)
        surface_pressure, distances, date_min_max = variable_cat_line(var_index, station_list_cs)
        l, = plt.plot(distances, surface_pressure/100, color='k', linewidth=2.)
        
        legendEntries.append(l)
        legendtext.append('Surface Pressure')
        
        variable='T_eq_0'
        var_index = variable_name_index_match(variable, variable_list_line)
        T_eq_0, distances_line, date_min_max_line = variable_cat_line(var_index, station_list_cs)
        l, = plt.plot(distances, T_eq_0/100, color='red', linewidth=2., linestyle='-')
  
        legendEntries.append(l)
        legendtext.append('T=0')   
  
        variable='lcl_vpt'

        var_index = variable_name_index_match(variable, variable_list_line)
        lcl_vpt, distances_line, date_min_max_line = variable_cat_line(var_index, station_list_cs)
        l, = plt.plot(distances, lcl_vpt/100, color='k', linewidth=2., linestyle=':')
  
        legendEntries.append(l)
        legendtext.append('LCL')   
  
        variable='pbl_pressure'

        var_index = variable_name_index_match(variable, variable_list_line)
        pbl_pressure, distances_line, date_min_max_line = variable_cat_line(var_index, station_list_cs)
        l,= plt.plot(distances, pbl_pressure/100, color='k', linewidth=2., linestyle='--')
   
        legendEntries.append(l)
        legendtext.append('PBL')

        l0=plt.legend(legendEntries, legendtext,title='', loc='upper left', bbox_to_anchor=(1.2, 1.05),
                      ncol=3, fancybox=True, shadow=True )


    

In [8]:

    

load_file['date_bin_mean_all_dates_one_station_single'].shape


    

    




---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-8-3bc5af8753df> in <module>()
----> 1 load_file['date_bin_mean_all_dates_one_station_single'].shape

NameError: name 'load_file' is not defined

In [ ]:

    

Cross_Section_Title = 'Vizag_to_Kandahar'
station_list_cs=[43150, 42867, 43014, 42339, 40990]
first_station=station_list_cs[0]

maths_bearing_to_rotate_winds=-42.469648282 # See 'Radiosonde Location Map of Picked Sites' notebook for calculation

date_min=datetime.datetime(1960,5,1,0,0,0)
date_max=datetime.datetime(2014,10,1,0,0,0)

delta = relativedelta(weeks=+1)

variable_list={'pressures': 0, 'temps':1, 'dewpoints':2, 'winddirs':3, 'windspeeds':4, 'pot_temp':5, 
               'sat_vap_pres':6, 'vap_press':7, 'rel_hum':8, 'wvmr':9, 'sp_hum':10, 'sat_temp':11, 
               'theta_e':12, 'theta_e_sat':13, 'theta_e_minus_theta_e_sat':14}
variable_list_line={'lcl_temp': 0, 'lcl_vpt':1, 'pbl_pressure':2, 'surface_pressure':3, 'T_eq_0':4, 
                    'CAPE':5, 'CIN':6, 'lclp':7, 'lclt':8, 'lfcp':9, 'equil_level':10}

first_station_title, first_station_lon, first_station_lat = station_info_search(first_station)

# Create plot arrays for pressure, winds etc

variable='pressures'

var_index = variable_name_index_match(variable, variable_list)
pressures, distances, date_min_max,dates_savefig = variable_cat(var_index, station_list_cs)

variable='windspeeds'

var_index = variable_name_index_match(variable, variable_list)
wind_speed, distances, date_min_max,dates_savefig = variable_cat(var_index, station_list_cs)

variable='winddirs'

var_index = variable_name_index_match(variable, variable_list)
wind_direction, distances, date_min_max,dates_savefig = variable_cat(var_index, station_list_cs)

u_wind,v_wind = u_v_winds(wind_direction, wind_speed)


    

In [ ]:

    

#Weekly Climatology with pbl and lcl, ROTATED winds

max_contour=100
min_contour=0
tick_interval=10

cmap=plt.cm.jet


variable_to_plot='rel_hum' # From variable_list

var_index = variable_name_index_match(variable_to_plot, variable_list)
concat_plot_variable, distances, date_min_max,dates_savefig = variable_cat(var_index, station_list_cs)

cbar_number_format = '$%d$'
# For each time bin

for i, date_bin in enumerate(dates_savefig):

# Rotate winds
            
    try:    
   
# Grid data
        u_rot,v_rot = rotate_u_v_wind_array(u_wind[:,i,0:-1:10], v_wind[:,i,0:-1:10], radians(maths_bearing_to_rotate_winds))
        nan_mask = np.ma.masked_array(np.array(concat_plot_variable[:,i,:], dtype=float).flatten(), np.isnan(np.array(concat_plot_variable[:,i,:], dtype=float).flatten()))
        if nan_mask.mask.all() == False:
            xi,yi,zi = grid_data_cs(np.array(pressures[:,i,:], dtype=float).flatten(), np.repeat(distances, concat_plot_variable[:,i,:].shape[1]).flatten(), nan_mask) 
            
# Contour plot
            
            cont,cbar = plot_rad_cs_winds_rot_and_unrot(i, xi, yi, zi, min_contour, max_contour, np.repeat(distances, concat_plot_variable[:,i,0:-1:10].shape[1]).flatten(),
                                  np.array(pressures[:,i,0:-1:10], dtype=float).flatten(), u_wind[:,i,0:-1:10], v_wind[:,i,0:-1:10], u_rot, v_rot, cmap, cbar_number_format)
           # station_name_plot(station_list_cs, first_station, yi)


# Line plots
        legendEntries=[]
        legendtext=[]
    
        variable='surface_pressure'
        var_index = variable_name_index_match(variable, variable_list_line)
        surface_pressure, distances, date_min_max = variable_cat_line(var_index, station_list_cs)
        l, = plt.plot(distances, surface_pressure/100, color='k', linewidth=2.)
        
        legendEntries.append(l)
        legendtext.append('Surface Pressure')
        
        variable='T_eq_0'
        var_index = variable_name_index_match(variable, variable_list_line)
        T_eq_0, distances_line, date_min_max_line = variable_cat_line(var_index, station_list_cs)
        l, = plt.plot(distances, T_eq_0/100, color='red', linewidth=2., linestyle='-')
  
        legendEntries.append(l)
        legendtext.append('T=0')   
  
        variable='lcl_vpt'

        var_index = variable_name_index_match(variable, variable_list_line)
        lcl_vpt, distances_line, date_min_max_line = variable_cat_line(var_index, station_list_cs)
        l, = plt.plot(distances, lcl_vpt/100, color='k', linewidth=2., linestyle=':')
  
        legendEntries.append(l)
        legendtext.append('LCL')   
  
        variable='pbl_pressure'

        var_index = variable_name_index_match(variable, variable_list_line)
        pbl_pressure, distances_line, date_min_max_line = variable_cat_line(var_index, station_list_cs)
        l,= plt.plot(distances, pbl_pressure/100, color='k', linewidth=2., linestyle='--')
   
        legendEntries.append(l)
        legendtext.append('PBL')

        l0=plt.legend(legendEntries, legendtext,title='', loc='upper left', bbox_to_anchor=(1.2, 1.05),
                      ncol=3, fancybox=True, shadow=True )  
    
        cbar.set_label('\%', rotation=90)
    
# Plot station names and vertical lines at their location

        y_offset_text=0
        first_station_title, first_station_lon, first_station_lat = station_info_search(first_station)
    
        for stat in station_list_cs: 
     
            station_title, station_lon, station_lat = station_info_search(stat)
    
            dist_from_first_station = calculate_distance_from_first_station(stat, first_station_lon, first_station_lat, station_lat, station_lon)
            
            plt.axvline(x=dist_from_first_station, ymin=0, ymax=1, label=station_title, color='k')
            plt.text(dist_from_first_station+0.1,max(yi)/100+20,station_title,rotation=-45)

            y_offset_text=+1 
        
  
        #plt.title('Week Beg.%s %s Cross-Section of Relative Humidity and Winds along Cross-Section Plane from Radiosonde Soundings' % (date_bin.strftime("%d %B"), Cross_Section_Title.replace('_',' ') ))
        plt.title('Week Beg.%s %s Relative Humidity' % (date_bin.strftime("%d %B"), Cross_Section_Title.replace('_',' ') ))

        #plt.show()
        plt.savefig('/nfs/a90/eepdw/Figures/Radiosonde/Cross_Sections/%s_%s_%s_Relative_Humidity_No_Title.png' % (Cross_Section_Title, date_bin.strftime("%y"), date_bin.strftime("%d_%B")),  format='png', bbox_inches='tight')

    except Exception, e:
        PrintException()
  
 # plt.savefig('/nfs/a90/eepdw/Figures/Radiosonde/Cross_Sections/%s_%s_%s_Relative_Humidity_V_Wind.png' % (Cross_Section_Title, date_bin.strftime("%y"), date_bin.strftime("%d_%B")),  format='png', bbox_inches='tight')
 # plt.close()
 # plt.clf()
 # gc.collect()


    

In [ ]:

    

# Weekly Climatology POT_TEMP  with pbl and lcl, ROTATED winds

max_contour=350
min_contour=300
tick_interval=10

#variable='theta_e_minus_theta_e_sat'
def station_name_vertical_line(station_list_cs, first_station_lon, first_station_lat, date_bin, Cross_Section_Title):
    
    for stat in station_list_cs: 
     
        station_title, station_lon, station_lat = station_info_search(stat)
    
        dist_from_first_station = calculate_distance_from_first_station(stat, first_station_lon, first_station_lat, station_lat, station_lon)
            
        plt.axvline(x=dist_from_first_station, ymin=0, ymax=1, label=station_title, color='k')
        plt.text(dist_from_first_station+0.1,max(yi)/100+20,station_title,rotation=-45)

        y_offset_text=+1 
        
        
variable='pressures'

var_index = variable_name_index_match(variable, variable_list)
concat_plot_variable, distances, date_min_max,dates_savefig = variable_cat(var_index, station_list_cs)

#Potential Temperature

variable='pot_temp'

var_index = variable_name_index_match(variable, variable_list)
concat_plot_variable, distances, date_min_max,dates_savefig = variable_cat(var_index, station_list_cs)
print var_index
print variable_list
    
cbar_number_format = '$%d$'
cmap=plt.cm.jet

for i, date_bin in enumerate(dates_savefig):
        
 try:
  u_rot,v_rot = rotate_u_v_wind_array(u_wind[:,i,0:-1:10], v_wind[:,i,0:-1:10], radians(maths_bearing_to_rotate_winds))
  nan_mask = np.ma.masked_array(np.array(concat_plot_variable[:,i,:], dtype=float).flatten(), np.isnan(np.array(concat_plot_variable[:,i,:], dtype=float).flatten()))
  #if nan_mask.mask.all() == False:
  xi,yi,zi = grid_data_cs(np.array(pressures[:,i,:], dtype=float).flatten(), np.repeat(distances, concat_plot_variable[:,i,:].shape[1]).flatten(), nan_mask) 
            
# Contour plot
            
  cont,cbar = plot_rad_cs_winds_rot_and_unrot(i, xi, yi, zi, min_contour, max_contour, np.repeat(distances, concat_plot_variable[:,i,0:-1:10].shape[1]).flatten(),
                                  np.array(pressures[:,i,0:-1:10], dtype=float).flatten(), u_wind[:,i,0:-1:10], v_wind[:,i,0:-1:10], u_rot, v_rot, cmap, cbar_number_format)
        
  LinePlots(variable_list_line, station_list_cs, distances) 
        
  cbar.set_label('\%', rotation=90)
  plt.title('Week beg. %s %s Cross-Section of Potential Temperature from Radiosonde Soundings' % (date_bin.strftime("%d %B"), Cross_Section_Title.replace('_',' ') ))

  station_name_vertical_line(station_list_cs, first_station_lon, first_station_lat, date_bin, Cross_Section_Title)
        
  #plt.show()
  plt.savefig('/nfs/a90/eepdw/Figures/Radiosonde/Cross_Sections/%s_%s_%s_%s_pottemp.png' % (Cross_Section_Title,  date_min.strftime('%Y%m%d'), date_max.strftime('%Y%m%d'), date_bin.strftime("%d_%B")),  format='png', bbox_inches='tight')
  plt.close()
  plt.clf()
  gc.collect()
 except Exception, e:
   PrintException()


    

In [ ]:

    

grid=[]

date_min_month_bin=datetime.datetime(2011,5,1,0,0,0)
date_max_month_bin=datetime.datetime(2011,10,1,0,0,0)
delta = relativedelta(weeks=+1)

d = date_min_month_bin
while (d <= date_max_month_bin): 
 grid.append(d.timetuple().tm_yday)
 d += delta
    

for stat in station_list_cs:

  load_file = np.load('/nfs/a90/eepdw/Data/Observations/Radiosonde_Numpy/Radiosonde_Cross_Section_'
                        'IND_SOUNDING_INTERP_MEAN_Climat_%s_%s_%s_%s.npz' 
                        % (date_min.strftime('%Y%m%d'), date_max.strftime('%Y%m%d'), delta, stat))
  load_file['date_bin_mean_all_dates_one_station_single'][:,0,1]

  date_plot = [datetime.datetime.fromordinal(d).replace(year=2011) for d in grid]
    
  station_name,la,lo=station_info_search(stat)
  
  fig=plt.figure()
    
  variable='CAPE'
  var_index = variable_name_index_match(variable, variable_list_line)
  plt.plot_date(date_plot[1:], load_file['date_bin_mean_all_dates_one_station_single'][1:,0,var_index], fmt="-", label=variable)  
  
  variable='CIN'
  var_index = variable_name_index_match(variable, variable_list_line)
  plt.plot_date(date_plot[1:], load_file['date_bin_mean_all_dates_one_station_single'][1:,0,var_index], fmt="-", label=variable)  
  
  plt.legend()
  plt.title('%s' % station_name)
    
  fig.autofmt_xdate()

  #plt.gca().invert_yaxis()
  plt.show()


    

In [ ]:

    

grid=[]

date_min_month_bin=datetime.datetime(2011,5,1,0,0,0)
date_max_month_bin=datetime.datetime(2011,10,1,0,0,0)
delta = relativedelta(weeks=+1)

d = date_min_month_bin
while (d <= date_max_month_bin): 
 grid.append(d.timetuple().tm_yday)
 d += delta
    

for stat in station_list_cs:

  load_file = np.load('/nfs/a90/eepdw/Data/Observations/Radiosonde_Numpy/Radiosonde_Cross_Section_'
                        'IND_SOUNDING_INTERP_MEAN_Climat_%s_%s_%s_%s.npz' 
                        % (date_min.strftime('%Y%m%d'), date_max.strftime('%Y%m%d'), delta, stat))
  load_file['date_bin_mean_all_dates_one_station_single'][:,0,1]

  date_plot = [datetime.datetime.fromordinal(d).replace(year=2011) for d in grid]
    
  station_name,la,lo=station_info_search(stat)
  
  fig=plt.figure()
    
  for variable in ('pbl_pressure', 'surface_pressure', 'T_eq_0', 'lclp', 'lfcp', 'equil_level', 'lcl_vpt'):
      if variable in ('equil_level', 'lclp', 'lfcp'):                                  
          var_index = variable_name_index_match(variable, variable_list_line)
          plt.plot_date(date_plot[1:], load_file['date_bin_mean_all_dates_one_station_single'][1:,0,var_index], fmt="-", label=variable)
      elif variable in ('pbl_pressure', 'surface_pressure', 'T_eq_0', 'lcl_vpt'):                                  
          var_index = variable_name_index_match(variable, variable_list_line)
          plt.plot_date(date_plot[1:], load_file['date_bin_mean_all_dates_one_station_single'][1:,0,var_index]/100, fmt="-", label=variable)
      
  plt.legend()
  plt.title('%s' % station_name)
  plt.gca().invert_yaxis()
                
  fig.autofmt_xdate()
                
  plt.show()


    

In [ ]:

    

def station_name_vertical_line(station_list_cs, first_station_lon, first_station_lat, date_bin, Cross_Section_Title):
    
    for stat in station_list_cs: 
     
        station_title, station_lon, station_lat = station_info_search(stat)
    
        dist_from_first_station = calculate_distance_from_first_station(stat, first_station_lon, first_station_lat, station_lat, station_lon)
            
        plt.axvline(x=dist_from_first_station, ymin=0, ymax=1, label=station_title, color='k')
        plt.text(dist_from_first_station+0.1,max(yi)/100+20,station_title,rotation=-45)

        y_offset_text=+1 
        
        
variable='pressures'

var_index = variable_name_index_match(variable, variable_list)
concat_plot_variable, distances, date_min_max,dates_savefig = variable_cat(var_index, station_list_cs)

variable='rel_hum'

var_index = variable_name_index_match(variable, variable_list)
concat_plot_variable, distances, date_min_max,dates_savefig = variable_cat(var_index, station_list_cs)

max_contour=100
min_contour=0
tick_interval=10

cbar_number_format = '$%d$'
cmap=plt.cm.jet
    
for i, date_bin in enumerate(dates_savefig):
        
 try:
  # Grid data
  u_rot,v_rot = rotate_u_v_wind_array(u_wind[:,i,0:-1:10], v_wind[:,i,0:-1:10], radians(maths_bearing_to_rotate_winds))
  nan_mask = np.ma.masked_array(np.array(concat_plot_variable[:,i,:], dtype=float).flatten(), np.isnan(np.array(concat_plot_variable[:,i,:], dtype=float).flatten()))
  if nan_mask.mask.all() == False:
      xi,yi,zi = grid_data_cs(np.array(pressures[:,i,:], dtype=float).flatten(), np.repeat(distances, concat_plot_variable[:,i,:].shape[1]).flatten(), nan_mask) 
            
# Contour plot
            
  cont,cbar = plot_rad_cs_winds_rot_and_unrot(i, xi, yi, zi, min_contour, max_contour, np.repeat(distances, concat_plot_variable[:,i,0:-1:10].shape[1]).flatten(),
                                  np.array(pressures[:,i,0:-1:10], dtype=float).flatten(), u_wind[:,i,0:-1:10], v_wind[:,i,0:-1:10], u_rot, v_rot, cmap, cbar_number_format)
           # station_name_plot(station_list_cs, first_station, yi)
        
  LinePlots(variable_list_line, station_list_cs, distances) 
        
  cbar.set_label('\%', rotation=90)

  #print date_bin
  
  y_offset_text=0
  first_station_title, first_station_lon, first_station_lat = station_info_search(first_station)
  
  station_name_vertical_line(station_list_cs, first_station_lon, first_station_lat, date_bin, Cross_Section_Title)
    
  plt.title('%s %s Cross-Section of Relative Humidity from Radiosonde Soundings' % (date_bin.strftime("%d %B"), Cross_Section_Title.replace('_',' ') ))
    
  #plt.show()
  
  plt.savefig('/nfs/a90/eepdw/Figures/Radiosonde/Cross_Sections/%s_%s_%s_%s_Relative_Humidity.png' % (Cross_Section_Title,  date_min.strftime('%Y%m%d'), date_max.strftime('%Y%m%d'), date_bin.strftime("%d_%B")),  format='png', bbox_inches='tight')
  plt.close()
  plt.clf()
  gc.collect()
 except Exception, e:
   PrintException()
        
#Potential Temperature

variable='pot_temp'

var_index = variable_name_index_match(variable, variable_list)
concat_plot_variable, distances, date_min_max,dates_savefig = variable_cat(var_index, station_list_cs)

max_contour=350
min_contour=300
tick_interval=10
    
cbar_number_format = '$%d$'
cmap=plt.cm.jet

for i, date_bin in enumerate(dates_savefig):
        
 try:
  u_rot,v_rot = rotate_u_v_wind_array(u_wind[:,i,0:-1:10], v_wind[:,i,0:-1:10], radians(maths_bearing_to_rotate_winds))
  nan_mask = np.ma.masked_array(np.array(concat_plot_variable[:,i,:], dtype=float).flatten(), np.isnan(np.array(concat_plot_variable[:,i,:], dtype=float).flatten()))
  if nan_mask.mask.all() == False:
      xi,yi,zi = grid_data_cs(np.array(pressures[:,i,:], dtype=float).flatten(), np.repeat(distances, concat_plot_variable[:,i,:].shape[1]).flatten(), nan_mask) 
            
# Contour plot
            
  cont,cbar = plot_rad_cs_winds_rot_and_unrot(i, xi, yi, zi, min_contour, max_contour, np.repeat(distances, concat_plot_variable[:,i,0:-1:10].shape[1]).flatten(),
                                  np.array(pressures[:,i,0:-1:10], dtype=float).flatten(), u_wind[:,i,0:-1:10], v_wind[:,i,0:-1:10], u_rot, v_rot, cmap, cbar_number_format)
        
  LinePlots(variable_list_line, station_list_cs, distances) 
        
  cbar.set_label('\%', rotation=90)
  plt.title('Week beg. %s %s Cross-Section of Potential Temperature from Radiosonde Soundings' % (date_bin.strftime("%d %B"), Cross_Section_Title.replace('_',' ') ))

  station_name_vertical_line(station_list_cs, first_station_lon, first_station_lat, date_bin, Cross_Section_Title)
        
  #plt.show()
  plt.savefig('/nfs/a90/eepdw/Figures/Radiosonde/Cross_Sections/%s_%s_%s_%s_pottemp.png' % (Cross_Section_Title,  date_min.strftime('%Y%m%d'), date_max.strftime('%Y%m%d'), date_bin.strftime("%d_%B")),  format='png', bbox_inches='tight')
  plt.close()
  plt.clf()
  gc.collect()
 except Exception, e:
   PrintException()


#Theta E

variable='theta_e'

var_index = variable_name_index_match(variable, variable_list)
concat_plot_variable, distances, date_min_max,dates_savefig = variable_cat(var_index, station_list_cs)

max_contour=380
min_contour=320
tick_interval=10
    
cbar_number_format = '$%d$'
cmap=plt.cm.jet

for i, date_bin in enumerate(dates_savefig):
        
 try:
  u_rot,v_rot = rotate_u_v_wind_array(u_wind[:,i,0:-1:10], v_wind[:,i,0:-1:10], radians(maths_bearing_to_rotate_winds))
  nan_mask = np.ma.masked_array(np.array(concat_plot_variable[:,i,:], dtype=float).flatten(), np.isnan(np.array(concat_plot_variable[:,i,:], dtype=float).flatten()))
  if nan_mask.mask.all() == False:
      xi,yi,zi = grid_data_cs(np.array(pressures[:,i,:], dtype=float).flatten(), np.repeat(distances, concat_plot_variable[:,i,:].shape[1]).flatten(), nan_mask) 
            
# Contour plot
            
  cont,cbar = plot_rad_cs_winds_rot_and_unrot(i, xi, yi, zi, min_contour, max_contour, np.repeat(distances, concat_plot_variable[:,i,0:-1:10].shape[1]).flatten(),
                                  np.array(pressures[:,i,0:-1:10], dtype=float).flatten(), u_wind[:,i,0:-1:10], v_wind[:,i,0:-1:10], u_rot, v_rot, cmap, cbar_number_format)
        
  LinePlots(variable_list_line, station_list_cs, distances) 
        
  cbar.set_label('K', rotation=90)
  plt.title('Week beg. %s %s Cross-Section of Equivalent Potential Temperature from Radiosonde Soundings' % (date_bin.strftime("%d %B"), Cross_Section_Title.replace('_',' ') ))

  station_name_vertical_line(station_list_cs, first_station_lon, first_station_lat, date_bin, Cross_Section_Title)
    
  #plt.show()
  plt.savefig('/nfs/a90/eepdw/Figures/Radiosonde/Cross_Sections/%s_%s_%s_%s_thetae.png' % (Cross_Section_Title,  date_min.strftime('%Y%m%d'), date_max.strftime('%Y%m%d'), date_bin.strftime("%d_%B")),  format='png', bbox_inches='tight')
  plt.close()
  plt.clf()
  gc.collect()
 except Exception, e:
   PrintException()

#Theta Es

variable='theta_e_sat'

var_index = variable_name_index_match(variable, variable_list)
concat_plot_variable, distances, date_min_max,dates_savefig = variable_cat(var_index, station_list_cs)

max_contour=420
min_contour=330
tick_interval=10
    
cbar_number_format = '$%d$'
cmap=plt.cm.jet

for i, date_bin in enumerate(dates_savefig):
        
 try:
  u_rot,v_rot = rotate_u_v_wind_array(u_wind[:,i,0:-1:10], v_wind[:,i,0:-1:10], radians(maths_bearing_to_rotate_winds))
  nan_mask = np.ma.masked_array(np.array(concat_plot_variable[:,i,:], dtype=float).flatten(), np.isnan(np.array(concat_plot_variable[:,i,:], dtype=float).flatten()))
  if nan_mask.mask.all() == False:
      xi,yi,zi = grid_data_cs(np.array(pressures[:,i,:], dtype=float).flatten(), np.repeat(distances, concat_plot_variable[:,i,:].shape[1]).flatten(), nan_mask) 
            
# Contour plot
            
  cont,cbar = plot_rad_cs_winds_rot_and_unrot(i, xi, yi, zi, min_contour, max_contour, np.repeat(distances, concat_plot_variable[:,i,0:-1:10].shape[1]).flatten(),
                                  np.array(pressures[:,i,0:-1:10], dtype=float).flatten(), u_wind[:,i,0:-1:10], v_wind[:,i,0:-1:10], u_rot, v_rot, cmap, cbar_number_format)
        
  LinePlots(variable_list_line, station_list_cs, distances) 
        
  cbar.set_label('K', rotation=90)
  plt.title('Week beg. %s %s Cross-Section of Saturated Equivalent Potential Temperature from Radiosonde Soundings' % (date_bin.strftime("%d %B"), Cross_Section_Title.replace('_',' ') ))

  station_name_vertical_line(station_list_cs, first_station_lon, first_station_lat, date_bin, Cross_Section_Title)
    
  #plt.show()
  plt.savefig('/nfs/a90/eepdw/Figures/Radiosonde/Cross_Sections/%s_%s_%s_%s_thetaesat.png' % (Cross_Section_Title,  date_min.strftime('%Y%m%d'), date_max.strftime('%Y%m%d') , date_bin.strftime("%d_%B")),  format='png', bbox_inches='tight')
  plt.close()
  plt.clf()
  gc.collect()
 except Exception, e:
   PrintException()
        
        
# WVMR

variable = 'wvmr'

var_index = variable_name_index_match(variable, variable_list)
concat_plot_variable, distances, date_min_max,dates_savefig = variable_cat(var_index, station_list_cs)

max_contour=0.02
min_contour=0.
tick_interval=0.001

cbar_number_format = '$%.3f$'
cmap=plt.cm.jet

for i, date_bin in enumerate(dates_savefig):
        
 try:
  u_rot,v_rot = rotate_u_v_wind_array(u_wind[:,i,0:-1:10], v_wind[:,i,0:-1:10], radians(maths_bearing_to_rotate_winds))
  nan_mask = np.ma.masked_array(np.array(concat_plot_variable[:,i,:], dtype=float).flatten(), np.isnan(np.array(concat_plot_variable[:,i,:], dtype=float).flatten()))
  if nan_mask.mask.all() == False:
      xi,yi,zi = grid_data_cs(np.array(pressures[:,i,:], dtype=float).flatten(), np.repeat(distances, concat_plot_variable[:,i,:].shape[1]).flatten(), nan_mask) 
            
# Contour plot
            
  cont,cbar = plot_rad_cs_winds_rot_and_unrot(i, xi, yi, zi, min_contour, max_contour, np.repeat(distances, concat_plot_variable[:,i,0:-1:10].shape[1]).flatten(),
                                  np.array(pressures[:,i,0:-1:10], dtype=float).flatten(), u_wind[:,i,0:-1:10], v_wind[:,i,0:-1:10], u_rot, v_rot, cmap, cbar_number_format)
        
  LinePlots(variable_list_line, station_list_cs, distances) 
        
  cbar.set_label('', rotation=90)
  plt.title('Week beg. %s %s Cross-Section of WVMR from Radiosonde Soundings' % (date_bin.strftime("%d %B"), Cross_Section_Title.replace('_',' ') ))

  station_name_vertical_line(station_list_cs, first_station_lon, first_station_lat, date_bin, Cross_Section_Title)
    
    #plt.show()
  plt.savefig('/nfs/a90/eepdw/Figures/Radiosonde/Cross_Sections/%s_%s_%s_%s_wvmr.png' % (Cross_Section_Title,  date_min.strftime('%Y%m%d'), date_max.strftime('%Y%m%d'), date_bin.strftime("%d_%B")),  format='png', bbox_inches='tight')
  plt.close()
  plt.clf()
  gc.collect()
 except Exception, e:
   PrintException()



#Theta Es - Theta E - Buoyancy

variable='theta_e_minus_theta_e_sat'

var_index = variable_name_index_match(variable, variable_list)
concat_plot_variable, distances, date_min_max,dates_savefig = variable_cat(var_index, station_list_cs)

cbar_number_format = '$%d$'

cmap=plt.cm.RdBu_r

max_contour=10
min_contour=-50
tick_interval=10
for i, date_bin in enumerate(dates_savefig):
        
 try:
  u_rot,v_rot = rotate_u_v_wind_array(u_wind[:,i,0:-1:10], v_wind[:,i,0:-1:10], radians(maths_bearing_to_rotate_winds))
  nan_mask = np.ma.masked_array(np.array(concat_plot_variable[:,i,:], dtype=float).flatten(), np.isnan(np.array(concat_plot_variable[:,i,:], dtype=float).flatten()))
  if nan_mask.mask.all() == False:
      xi,yi,zi = grid_data_cs(np.array(pressures[:,i,:], dtype=float).flatten(), np.repeat(distances, concat_plot_variable[:,i,:].shape[1]).flatten(), nan_mask) 
            
# Contour plot
            
  cont,cbar = plot_rad_cs_winds_rot_and_unrot_normalize(i, xi, yi, zi, min_contour, max_contour, np.repeat(distances, concat_plot_variable[:,i,0:-1:10].shape[1]).flatten(),
                                  np.array(pressures[:,i,0:-1:10], dtype=float).flatten(), u_wind[:,i,0:-1:10], v_wind[:,i,0:-1:10], u_rot, v_rot, cmap, cbar_number_format)
        
  LinePlots(i, variable_list_line, station_list_cs, distances) 
        
        
  cbar.set_label('K', rotation=90)
  plt.title('Week beg. %s %s Cross-Section of Theta-e BL minus theta-es local from Radiosonde Soundings' % (date_bin.strftime("%d %B"), Cross_Section_Title.replace('_',' ') ))

  station_name_vertical_line(station_list_cs, first_station_lon, first_station_lat, date_bin, Cross_Section_Title)
    
  #plt.show()
  plt.savefig('/nfs/a90/eepdw/Figures/Radiosonde/Cross_Sections/%s_%s_%s_%s_thetaeminthetaesat_buoy.png' % (Cross_Section_Title, date_min.strftime('%Y%m%d'), date_max.strftime('%Y%m%d'), date_bin.strftime("%d_%B")),  format='png', bbox_inches='tight')
  plt.close()
  plt.clf()
  gc.collect()
 except Exception, e:
   PrintException()


    

In [ ]:

    

def station_name_vertical_line(station_list_cs, first_station_lon, first_station_lat, date_bin, Cross_Section_Title):
    
    for stat in station_list_cs: 
     
        station_title, station_lon, station_lat = station_info_search(stat)
    
        dist_from_first_station = calculate_distance_from_first_station(stat, first_station_lon, first_station_lat, station_lat, station_lon)
            
        plt.axvline(x=dist_from_first_station, ymin=0, ymax=1, label=station_title, color='k')
        plt.text(dist_from_first_station+0.1,max(yi)/100+20,station_title,rotation=-45)

        y_offset_text=+1 
        
def LinePlots(i, variable_list_line, station_list_cs, distances):
        legendEntries=[]
        legendtext=[]
    
        variable='surface_pressure'
        var_index = variable_name_index_match(variable, variable_list_line)
        surface_pressure, distances, date_min_max = variable_cat_line(var_index, station_list_cs)
        l, = plt.plot(distances, surface_pressure[:,i]/100, color='k', linewidth=2.)
        
        legendEntries.append(l)
        legendtext.append('Surface Pressure')
        
        variable='T_eq_0'
        var_index = variable_name_index_match(variable, variable_list_line)
        T_eq_0, distances_line, date_min_max_line = variable_cat_line(var_index, station_list_cs)
        l, = plt.plot(distances, T_eq_0[:,i]/100, color='red', linewidth=2., linestyle='-')
  
        legendEntries.append(l)
        legendtext.append('T=0')   
  
        variable='lcl_vpt'

        var_index = variable_name_index_match(variable, variable_list_line)
        lcl_vpt, distances_line, date_min_max_line = variable_cat_line(var_index, station_list_cs)
        l, = plt.plot(distances, lcl_vpt[:,i]/100, color='k', linewidth=2., linestyle=':')
  
        legendEntries.append(l)
        legendtext.append('LCL')   
  
        variable='pbl_pressure'

        var_index = variable_name_index_match(variable, variable_list_line)
        pbl_pressure, distances_line, date_min_max_line = variable_cat_line(var_index, station_list_cs)
        l,= plt.plot(distances, pbl_pressure[:,i]/100, color='k', linewidth=2., linestyle='--')
   
        legendEntries.append(l)
        legendtext.append('PBL')

        l0=plt.legend(legendEntries, legendtext,title='', loc='upper left', bbox_to_anchor=(1.2, 1.05),
                      ncol=3, fancybox=True, shadow=True ) 
        
def PlotRadSondeTransNoWindLoop(i, xi,yi,zi, min_contour, max_contour, cmap, cbar_number_format):
    
    num_cols=32
   
    plt.figure(figsize=(14,8.5))
       
    clevs = np.linspace(min_contour, max_contour, num_cols)
    ticks = (np.arange(min_contour, max_contour+tick_interval,tick_interval))
    
 #cmap=plt.cm.RdBu_r
    
    cont = plt.contourf(xi,yi/100, zi, clevs, cmap=cmap, extend='both')

 #plt.clabel(cont_wind, fontsize=9, inline=1, fmt='%1d')
  
    cbar = plt.colorbar(cont, orientation='vertical', pad=0.05, extend='both')
 #cbar.set_label('$W m^{-2}$') 
    cbar.set_ticks(np.arange(min_contour, max_contour+tick_interval,tick_interval))
    cbar.set_ticklabels([cbar_number_format % i for i in ticks])
 
    plt.xlim(np.min(xi)-200,np.max(xi)+200)
    plt.gca().invert_yaxis()
    plt.gca().invert_xaxis()
    plt.ylabel('Pressure (hPa)')
    plt.xlabel('km from first station')
    
    return cont,cbar

def station_info_search(stat):
    for line in station_metadata:
         if "%s" % stat in line: 
             st = line[2].lower().title().replace('_',' ')
             lo = float(line[3])
             la = float(line[4])
    return st,la,lo

variable='pressures'

var_index = variable_name_index_match(variable, variable_list)
pressures, distances, date_min_max,dates_savefig = variable_cat(var_index, station_list_cs)

#Theta Es

variable='theta_e_sat'

var_index = variable_name_index_match(variable, variable_list)
theta_e_sat, distances, date_min_max,dates_savefig = variable_cat(var_index, station_list_cs)

#Theta E

variable='theta_e'

var_index = variable_name_index_match(variable, variable_list)
theta_e, distances, date_min_max,dates_savefig = variable_cat(var_index, station_list_cs)

#Theta e boundary minus Theta Es

variable='theta_e_minus_theta_e_sat'

var_index = variable_name_index_match(variable, variable_list)
theta_e_minus_theta_e_sat, distances, date_min_max,dates_savefig = variable_cat(var_index, station_list_cs)

#Theta E

variable='theta_e'

var_index = variable_name_index_match(variable, variable_list)
theta_e, distances, date_min_max,dates_savefig = variable_cat(var_index, station_list_cs)

# Surface pressure

variable='surface_pressure'

var_index = variable_name_index_match(variable, variable_list_line)
surface_pressure, distances, date_min_max = variable_cat_line(var_index, station_list_cs)

# Get array of theta_e 50 hPa above the surface

i_idx_50hPa = np.argmin(np.abs(np.rollaxis(pressures,axis=-1)-(surface_pressure-5000.)),axis=0)  # Closest point to 50 hPa less than surface

first_axis_length = len(i_idx_50hPa[:,0])
second_axis_length = len(i_idx_50hPa[0,:])

first_axis = np.repeat(np.arange(0, first_axis_length), second_axis_length).T.reshape(i_idx_50hPa.shape)
second_axis= np.repeat(np.arange(0, second_axis_length),first_axis_length).reshape(i_idx_50hPa.shape[::-1]).T

theta_e = theta_e[first_axis,second_axis, i_idx_50hPa]
         
theta_e_time_axis_diff = np.diff(theta_e, axis=1)
theta_e_sat_time_axis_diff = np.diff(theta_e_sat, axis=1)
theta_e_minus_theta_e_sat_time_axis_diff = np.diff(theta_e_minus_theta_e_sat, axis=1)

cmap=plt.cm.jet
cbar_number_format = '$%d$'

max_contour=10
min_contour=0
tick_interval=1

fig = plt.figure(figsize=(14,8.5))

for s in range(len(station_list_cs)):
    
    st,la, lo = station_info_search(station_list_cs[s])

    plt.plot_date(dates_savefig[1:], theta_e_time_axis_diff[s], fmt="-", label=st)
    #
    #plt.yticks(arange(0,24,4))
    #plt.ylim(0,24)
    
#plt.ylim(360,380)    
#plt.title('%s to %s hPa' % (hpa_level_lower, hpa_level_higher))
plt.xlabel('Week')
#plt.ylabel('Theta-e at %s hPa' % hpa_level)    
fig.autofmt_xdate()    
plt.legend()
#plt.show()
plt.savefig('/nfs/a90/eepdw/Figures/Radiosonde/Cross_Sections/%s_%s_%s_Theta_e_time_axis_diff.png' % (Cross_Section_Title, date_bin.strftime("%y"), date_bin.strftime("%d_%B")),  format='png', bbox_inches='tight')


for i, date_bin in enumerate(dates_savefig[1:]):
    
    nan_mask_theta_e_sat = np.abs(np.rollaxis(
                            theta_e_sat_time_axis_diff[:,i,:], axis=-1))
    nan_mask_theta_e = np.abs(theta_e_time_axis_diff[:,i])
    
    #plot_diff = np.rollaxis(
     #           (nan_mask_theta_e_sat)
      #          *
       #         (nan_mask_theta_e_sat
        #        /
         #       (nan_mask_theta_e_sat + nan_mask_theta_e))
          #      ,axis=0, start=2)
    plot_diff=np.abs(theta_e_sat_time_axis_diff[:,i,:])
    #plot_diff2 = theta_e_minus_theta_e_sat_time_axis_diff[:,i,:]/theta_e_sat_time_axis_diff[:,i,:]
    nan_mask = np.ma.masked_array(plot_diff, 
                                  np.isnan(plot_diff))
    
    
    if nan_mask.mask.all() == False:
        xi,yi,zi = grid_data_cs(pressures[:,i,:].flatten(), np.repeat(distances, nan_mask[0,:].flatten().shape[0]), nan_mask.flatten()) 

        cont, cbar = PlotRadSondeTransNoWindLoop(i, xi,yi,zi, min_contour, max_contour, cmap, cbar_number_format)
    
        cbar.set_label('K', rotation=90)
  

        plt.title('Week beg. %s %s Cross-Section of Theta-es time axis diff from Radiosonde Soundings' % (date_bin.strftime("%d %B"), Cross_Section_Title.replace('_',' ') ))

        station_name_vertical_line(station_list_cs, first_station_lon, first_station_lat, date_bin, Cross_Section_Title)

        LinePlots(i, variable_list_line, station_list_cs, distances) 
    
        #plt.show()
        plt.savefig('/nfs/a90/eepdw/Figures/Radiosonde/Cross_Sections/%s_%s_%s_Theta_e_sat_time_axis_diff.png' % (Cross_Section_Title, date_bin.strftime("%y"), date_bin.strftime("%d_%B")),  format='png', bbox_inches='tight')