""" Aina Johannessen 2017/2018 Plot trajectories in the vertical """ import pygrib #Used in get_data() to handle files. import matplotlib #Used for further calling import matplotlib.pyplot as plt #Uses in IVT() to set up a plot from mpl_toolkits.basemap import Basemap #Used in set_background() To make map import numpy as np #Used in IVT() for array operations from dypy.plotting import Mapfigure from dypy.small_tools import * #from dypy.lagranto import Tra from lagranto import LagrantoRun from lagranto import Tra import pandas as pd import warnings from matplotlib.collections import LineCollection import matplotlib.dates as mdates import matplotlib.cm as cm import matplotlib.pyplot as plt import pandas as pd from matplotlib.collections import LineCollection from matplotlib.colors import ListedColormap, BoundaryNorm from dypy.small_tools import CrossSection def select_more(trajs, time, path): new_trajs = Tra() #below 30 in atlantic if time =="20160927_18": if path == "all": new_trajs = trajs if path == "A": new_trajs_idx = np.where( (np.max(trajs["lat"], axis=1) > 70 ) & (np.min(trajs["lon"], axis=1) > -40) | # & (np.max(trajs["lon"], axis=1) < -20)) pickingsec = new_trajs_idx[0] new_trajs.set_array(trajs[pickingsec, :]) if path == "B": print("o") new_trajs_idx = np.where( (trajs["lat"] < 50 ) & ( trajs["lon"] > -55)& (trajs["lon"] < -53) ) # & pickingsec = new_trajs_idx[0] new_trajs.set_array(trajs[pickingsec, :]) new_trajs_idx = np.where( (np.min(new_trajs["lat"], axis=1) > 40 ) & (np.max(new_trajs["lat"], axis=1) <67)) pickingsec = new_trajs_idx[0] new_trajs.set_array(new_trajs[pickingsec, :]) #new_trajs_idx = np.where(((new_trajs["Q"][:, :1] - new_trajs["Q"] ) > 5.0)) #other? of the low levels ##new_trajs_idx = np.where( new_trajs["lat"][:,117] > 55 ) #& (np.max(new_trajs["p"], axis=1) >700)) #pickingsec = new_trajs_idx[0] #new_trajs.set_array(new_trajs[pickingsec, :]) #two of the low levels #new_trajs_idx = np.where( new_trajs["lat"][:,117] < 50 ) #& (np.max(new_trajs["p"], axis=1) >700)) #pickingsec = new_trajs_idx[0] #new_trajs.set_array(new_trajs[pickingsec, :]) #1 of the low levels #new_trajs_idx = np.where( new_trajs["p"][:,0] < 700 ) #under or above gives the two diff #pickingsec = new_trajs_idx[0] #new_trajs.set_array(new_trajs[pickingsec, :]) #seperate upper and lower path-- see thesis/heinis comment #new_trajs_idx = np.where( (np.min(new_trajs["Q"], axis=1) > 2) )#& (np.max(new_trajs["p"], axis=1) >700)) #new_trajs_idx = np.where( (np.max(new_trajs["Q"], axis=1) < 6) )#& (np.max(new_trajs["p"], axis=1) >700)) if time =="20160929_00": print ("WAAA") if path =="C": new_trajs_idx = np.where( (np.min(trajs["lat"], axis=1) < 20) & (np.min(trajs["lon"], axis=1) > -90) ) #& #(np.max(trajs["lon"], axis=1) < -50)) pickingsec = new_trajs_idx[0] new_trajs.set_array(trajs[pickingsec, :]) #if path =="C1": # new_trajs_idx = np.where( (np.min(trajs["lat"], axis=1) < 20) & # (np.min(trajs["lon"], axis=1) > -90) & # (np.max(trajs["lon"], axis=1) < 50)) # #(np.mean( trajs["lon"][:,:]**2 ) + trajs["lat"][:,:])) # pickingsec = new_trajs_idx[0] # new_trajs.set_array(trajs[pickingsec, :]) return new_trajs def horizontal_traj(trajs, fig_h, param, t,clr_para): if clr_para=="TH": print ("YEEE") C = [[232,232,230],#grey [203,203,203], #grey [161,161,161], #grey [130,130,130], #grey #[149,53,229], #lillac, 39 64 197 149,53,229 [39,64,197], #blue dark,7,67,194 [15,110,229], #blue [80,149,240], #blue [74,192,243], #blue [152,219,248], #blue [183,237,247], #blue [251,217,198], #redish [255,197,166], #redish [255,172,164], #redish [253,139,142], #redish [253,101,105], #redish [255,66,74], #redish [238,13,28], #red [214,78,166], #pink [214,102,201], [217,155,210], [216,181,211]] C = np.array( C ) C = np.divide( C, 255. ) # RGB has to be between 0 and 1 in python cmap = ListedColormap(C) norm = plt.Normalize(vmin=250,vmax=404) #norm = plt.Normalize(vmin=280,vmax=300) clr_label = "'potential temperature [$^\circ$C]" if clr_para =="Q": #cmap = plt.get_cmap('coolwarm_r') cmap = plt.get_cmap('RdYlGn') #PiYG Spectral RdYlGn norm = plt.Normalize(vmin=0,vmax=12) clr_label = "Specific humidity [g kg$^{-1}$]" elif clr_para =="T": cmap = plt.get_cmap('coolwarm') norm = plt.Normalize(vmin=-15,vmax=15) clr_label = "Temperature $^\circ$C] " elif clr_para =="PV": cmap = ListedColormap(['g','b', 'c', 'y','r', 'r']) norm = BoundaryNorm([ -0.5 ,0.5, 1, 1.5, 2, 3], cmap.N) clr_label = "Potential Vorticity" elif clr_para=="RH": cmap = plt.get_cmap('coolwarm_r') norm = plt.Normalize(vmin=0,vmax=100) clr_label = "Relative humidity" elif clr_para=="RWC": cmap = plt.get_cmap('Blues') norm = plt.Normalize(vmin=0,vmax=50) clr_label = "rain water content [ mg/kg ]" #fig.set_facecolor('grey') else: cmap = plt.get_cmap('coolwarm_r') norm = None clr_label = "??" m = Mapfigure(basemap = Basemap( llcrnrlon = -90., llcrnrlat = 10., urcrnrlon = 50., urcrnrlat=70., resolution = 'l', area_thresh = 10000., projection = 'merc' )) #m = Mapfigure() m.drawcoastlines( linewidth = 0.5, linestyle = 'solid', color = "k", zorder=3)#[ 75./255., 75/255., 75/255. ] ) m.drawmapboundary()#fill_color='aqua') m.fillcontinents(color=[ 75./255., 75/255., 75/255. ], zorder=1)#'coral',lake_color='aqua') # --------draw parallels------ circles = np.arange( -90., 90. + 30, 10. ) #delat = 10. m.drawparallels( circles, color = [ 55./255., 55/255., 55/255. ], labels = [ 1, 0, 0, 0 ], linewidth = 0.1, fontsize=7 ) # --------draw meridians----- meridians = np.arange( 0., 360, 10. ) #delon = 10. m.drawmeridians( meridians, color = [ 55./255., 55./255., 55./255. ], labels = [ 0, 0, 0, 1 ], linewidth = 0.1, fontsize=7 ) values = [-7, -2, -1, -0.5, 0, 0.5, 1, 2, 3] CS = m.plot_traj(trajs , param, linewidths=0.6 , cmap=cmap, norm = norm)# plt.title( " "+t +" - 10 days: UTC \n \ DQ(0-6h)>0.2: RH(0-6h)>80: P(traj)>700=BL \n", fontsize=9 ) cbaxes = fig_h.add_axes([0.19, 0.09, 0.67, 0.03]) #left, bottom, width, height], fractions that go from 0 to 1 ofthe plotting cbar = plt.colorbar(CS, fraction=0.046, cax = cbaxes, orientation="horizontal") cbar.set_label(clr_label, labelpad=-1) #cbar.set_label('PV:Potential vorticity', labelpad=-1) #fig_h.savefig( outfile, dpi = 600 ) #plt.show() return fig_h def traj(time, infile, y_para, clr_para, fig, ax, traj_inter, path, trajs): for traj in trajs: ###Handle datetime time### x = range(-len(traj['time']),0)[::-1] time = traj["time"] x = np.cumsum(np.random.normal(size=len(time))) s = pd.Series(x, index=time) inxval = mdates.date2num(s.index.to_pydatetime()) ######################### y = traj[y_para] clr = traj[clr_para] points = np.array([inxval, y ]).T.reshape(-1,1,2) segments = np.concatenate([points[:-1], points[1:]], axis=1) if clr_para =="Q": #cmap = plt.get_cmap('coolwarm_r') cmap = plt.get_cmap('RdYlGn') #PiYG Spectral RdYlGn norm = plt.Normalize(vmin=0,vmax=12) clr_label = "Specific humidity [g kg$^{-1}$]" elif clr_para =="T": cmap = plt.get_cmap('coolwarm') norm = plt.Normalize(vmin=-15,vmax=15) clr_label = "Temperature $^\circ$C] " elif clr_para =="PV": cmap = ListedColormap(['g','b', 'c', 'y','r', 'r']) norm = BoundaryNorm([ -0.5 ,0.5, 1, 1.5, 2, 3], cmap.N) clr_label = "Potential Vorticity" elif clr_para=="TH": C = [[232,232,230],#grey [203,203,203], #grey [161,161,161], #grey [130,130,130], #grey #[149,53,229], #lillac, 39 64 197 149,53,229 [39,64,197], #blue dark,7,67,194 [15,110,229], #blue [80,149,240], #blue [74,192,243], #blue [152,219,248], #blue [183,237,247], #blue [251,217,198], #redish [255,197,166], #redish [255,172,164], #redish [253,139,142], #redish [253,101,105], #redish [255,66,74], #redish [238,13,28], #red [214,78,166], #pink [214,102,201], [217,155,210], [216,181,211]] C = np.array( C ) C = np.divide( C, 255. ) # RGB has to be between 0 and 1 in python cmap = ListedColormap(C) #norm = plt.Normalize(vmin=250,vmax=404) #norm = plt.Normalize(vmin=250,vmax=404) #norm = plt.Normalize(vmin=280,vmax=300) norm=None clr_label = "'potential temperature [$^\circ$C]" elif clr_para=="RH": cmap = plt.get_cmap('coolwarm_r') norm = plt.Normalize(vmin=0,vmax=100) clr_label = "Relative humidity [%]" elif clr_para=="RWC": cmap = plt.get_cmap('Blues') norm = plt.Normalize(vmin=0,vmax=50) clr_label = "rain water content [ mg/kg ]" ax.set_facecolor('xkcd:grey') else: cmap = plt.get_cmap('coolwarm_r') norm = None clr_label = "??" lc = LineCollection(segments, cmap=cmap, norm = norm) lc.set_array(clr) lc.set_linewidth(1) plt.gca().add_collection(lc) if y_para =="p": y_label="Pressure [hpa]" plt.gca().invert_yaxis() elif y_para =="PV": y_label="Potential vorticity [pvu]" elif y_para=="TH": y_label = "'potential temperature [$^\circ$K]" elif y_para=="RWC": y_label = "Rain water content" else: y_label="???" ###SET date time to make x axis have nice dates#### hours = mdates.HourLocator(interval = 24) h_fmt = mdates.DateFormatter('%d Sept %H:00') ax.xaxis.set_major_locator( mdates.DayLocator() ) ax.xaxis.set_minor_locator( mdates.HourLocator() ) #interval = 1 ax.xaxis.set_major_locator(hours) ax.xaxis.set_major_formatter(h_fmt) ################################################### axcb = fig.colorbar(lc, fraction=0.046, pad=0.001) axcb.set_label('\n' + clr_label, fontsize=20, labelpad=0.01) #axcb.legend(handletextpad=0.1) axcb.ax.tick_params(labelsize=20) ax.autoscale_view() fig.autofmt_xdate() ax.set_ylabel(y_label,fontsize=20, labelpad=-1) ax.yaxis.set_tick_params(labelsize=20) ax.xaxis.set_tick_params(labelsize=20) plt.grid() return fig, ax def main(): """ time time in h relative to start of trajectories lon longitude in deg lat latitude in deg p pressure in hPa Q specific humidity in g/kg LWC liquid water content in mg/kg IWC ice water content in mg/kg RWC rain water content in mg/kg SWC snow water content in mg/kg RH relative humidity in % T temperature in degC TH potential temperature in K PV potential vorticity in pvu """ event = "Walpurga" #######CHANGE ONE, CHANGE BOTH!! time = "20160929_00" path = "C" ######### y_para = "p" clr_para = "RH" #Available fields: time/lon/lat/p/Q/LWC/IWC/RWC/SWC/RH/T/TH/PV total duration: -14400.0 minutes traj_inter = 2 #traj_inter = 25 #A every 2nd = 115 #B was every 15th = 193 path_nc="/Users/ainajohannessen/Documents/Aina/skole/master/trajectoryruns/sel_traj/" #infile = path_nc + "rainout_sources/1_all/netcdf/moist_"+time+".nc" infile = "./trajs_in_netcdf/moist_20160929_00.nc" outfile = path_nc+"rainout_sources/1_all/figs/moist_"+y_para+ "_VS_" +clr_para+"_"+time+".png" outfile = path_nc+"ME_moist_"+y_para+ "_VS_" +clr_para+"_"+time+".png" trajs = Tra() trajs.load_netcdf(infile) print("NUMBER OF TRAJS") print(np.shape(trajs)) trajs = select_more(trajs, time, path) trajs = trajs[::traj_inter,:] #118(perfect) #62 perf), 63, 64 65,66,67,68, 69,70,71,73,74,75,76, 77, 78, 79, 80,81, 82 83 84 85 86 87 88 89 90 91 92 93 94 95 print("NUMBER OF TRAJS after sel!") print(np.shape(trajs)) filename= clr_para + time + "path_"+path+".png" fig_h = plt.figure() fig_h = horizontal_traj(trajs, fig_h, clr_para, time, clr_para) plt.savefig("MAP"+filename, bbox_inches='tight', dpi = 200) fig, ax = plt.subplots(figsize=(10,4)) fig, ax = traj(time, infile, y_para, clr_para, fig, ax,traj_inter, path, trajs) filename= clr_para + time + "path_"+path+".png" #fig.set_size_inches( 11, 3 ) plt.savefig(filename, bbox_inches='tight', dpi = 200) #plt.savefig(outfile, bbox_inches='tight', dpi = 200) plt.show() warnings.filterwarnings("ignore",category=matplotlib.mplDeprecation) main()