In [2]:
%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 [3]:
# 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 [10]:
load_file['dates_for_plotting']
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'][i,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 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 [5]:
Cross_Section_Title = 'Thiruvan_to_Bombay'
station_list_cs=[43371, 43353, 43285, 43192, 43003]
first_station=station_list_cs[0]
maths_bearing_to_rotate_winds=110.598427774 # 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}
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 [7]:
#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)
# 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 [23]:
# Weekly Climatology theta_e_minus_theta_e_sat with pbl and lcl, ROTATED winds
max_contour=315
min_contour=215
tick_interval=20
#variable='theta_e_minus_theta_e_sat'
variable='pot_temp'
cbar_number_format = '$%d$'
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)
cmap=plt.cm.jet
# 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
LinePlots(variable_list_line, station_list_cs, distances)
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.show()
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 [ ]:
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)
variable='pbl_pressure'
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='PBL')
variable='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], fmt="-", label='LCL')
variable='surface_pressure'
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='SP')
variable='T_eq_0'
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='T_eq_0')
plt.legend()
plt.title('%s' % station_name)
plt.gca().invert_yaxis()
plt.show()
In [28]:
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=315
min_contour=215
tick_interval=20
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(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 [ ]:
# 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()