""" Aina Johannessen. 2017/2018 - Masterthesis """ import argparse #Used in user_interface() to do terminal input 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 from mpl_toolkits.basemap import shiftgrid #Used in IVT() to shift grid data from scipy.ndimage.filters import gaussian_filter #Used in IVT() to evenout bumps in data from scipy import ndimage, misc import numpy as np #Used in IVT() for array operations import warnings #Used before user_interface() call to remove uneccasary warnings from matplotlib.patches import Rectangle from mpl_toolkits.axes_grid1 import make_axes_locatable #import matplotlib.patches as patches #import matplotlib.colors as colors # If you want to spesify aditional colors from matplitlib #import spharm # Conda has it. Converts from spectral grid to gaussian grid, http://code.google.com/p/pyspharm """ ------PROGRAM USAGE------------------------------------------------------------------------------------------------------------- -------------------------------------------------------------------------------------------------------------------------------- This python program is independent of bash skript. It assumes all IVT files has been generated already(from bash program). It plots IVT and saves file. Path is set in user_interface() function Run from terminal, example: $python IVT_independent.py --date 20160920 0000 -------------------------------------------------------------------------------------------------------------------------------- -------------------------------------------------------------------------------------------------------------------------------- """ def get_data(IVT_infile, MSL_infile ): """ ---------------------------------------------------------------------------------------------------------------------------- Mission of function: - Extracts data from files and makes them suitable for further usage - Called from my custom IVT() function PARAMETERS inputs: - IVT_infile: Contains integrated water vapor transport generated from bash file: integrate_expr.sh. - MSL_infile: Contains mean sea level pressure. Gets directly from data_norestore server returns: - uq_data: IVT values in u direction [kg/(ms)]. Got from IVT_infile: - vq_data: IVT values in v direction [kg/(ms)]. Got from IVT_infile: - msl_data: mslp values [PA]=[kg/..]. Got from MSL_infile: - lat, lon Coordinates from IVT_infile other: - None USES installed functions: - pygrib - array.values() - array.latlons() custom functions: - None source: - None ------------------------------------------------------------------------------------------------------------------------------ """ #--INITIATING-------------------------------------------------------------------------------------------------- obj_uq = pygrib.open( IVT_infile ) #File containing IVT obj_msl= pygrib.open( MSL_infile ) print (IVT_infile) uq_obj = obj_uq(name = "U component of wind")[0] vq_obj = obj_uq(name = "V component of wind")[0] msl_obj = obj_msl(name = "Mean sea level pressure")[0] print (uq_obj) #---Extract data------------------------------------------------------------------------------------------------ msl_data = msl_obj.values/100. #divide by 100 to go from PA to hPa uq_data = uq_obj.values vq_data = vq_obj.values plat,plon= msl_obj.latlons() #lat and lon from mslp file:MSL_infile lat,lon = uq_obj.latlons() #lat and lon from u component of wind latv,lonv = vq_obj.latlons() ##lat and lon from v component of wind """ |->-----------NOTE for expression above------------------------------------------------------------------------| | Had to split up lat and lon to use shiftgrid on each parameter without shifting lon and lat multiple times. | | Probably a more efficient way were I can extract just one lat and lon, but for now it is what it is. | |--------------------------------------------------------------------------------------------------------------| """ #---Shifts grid in lon dir to get values where I need--------------------------- x_dir = np.arange(0,lon.shape[0],1 ) #y_dir = np.arange(0,lon.shape[1],1 ) for i in x_dir: uq_data[i,:], lon[i,:] = shiftgrid(180, uq_data[i,:], lon[i,:],start=False) vq_data[i,:], lonv[i,:] = shiftgrid(180, vq_data[i,:], lonv[i,:],start=False) msl_data[i,:], plon[i,:] = shiftgrid(180, msl_data[i,:], plon[i,:],start=False) """ |->-----------NOTE on for loop above---------------------------------------------------------------------------| | Not all basemap needs this, ut the one I have does. This is because the data starts at lon=0 (over england) | | and continues to plot with increasing lon. But that means I do not get data at negative vanlues of lon | \ which it is west of england. | |--------------------------------------------------------------------------------------------------------------| """ obj_uq.close() # obj_msl.close() #return uq_data, vq_data, msl_data, lat, lon return uq_data, vq_data, lat, lon, msl_data def set_background(): """ ------------------------------------------------------------------ Mission of function: - Generates the background map - Called from my custom IVT() function PARAMETERS inputs: - None returns: - m: other: - None USES installed functions: - Basemap - custom functions: - None source: - None ------------------------------------------------------------------- """ # ------------Atlantic----------------- m = Basemap( llcrnrlon = -90., llcrnrlat = 10., urcrnrlon = 50., urcrnrlat=70.,\ resolution = 'l', area_thresh = 10000., projection = 'merc' ) # ------------Norway------------------- #m = Basemap( llcrnrlon = -30., llcrnrlat = 40., urcrnrlon = 40., urcrnrlat=70.,\ # resolution = 'h', area_thresh = 10000., projection = 'merc' ) m.drawcoastlines( linewidth = 0.5, linestyle = 'solid', color = "k", zorder=5)#[ 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, 20. ) #delat = 10. m.drawparallels( circles, color = [ 55./255., 55/255., 55/255. ], labels = [ 1, 0, 0, 0 ], linewidth = 0.1, fontsize=12 ) # --------draw meridians----- meridians = np.arange( 0., 360, 20. ) #delon = 10. m.drawmeridians( meridians, color = [ 55./255., 55./255., 55./255. ], labels = [ 0, 0, 0, 1 ], linewidth = 0.1, fontsize=12 ) return m def IVT( fig_name, IVT_infile, time, date, mslp_infile ): """ ------------------------------------------------------------------ Mission of function: - The Brain of getting the IVT plots - Called from my custom user_interface() function - Can also be called indivually by: IVT(fig_name, IVT_infile, MSL_infile, time=0000, date=20160920) PARAMETERS inputs: - fig_name: - IVT_infile: Contains integrated water vapor transport generated from bash file: integrate_expr.sh. - MSL_infile: Contains mean sea level pressure. Gets directly from data_norestore server - time: - date: returns: - file: one figure in path "fig_name" other: - None USES installed functions: - numpy as np - np.arange() - np.divide() - np.meshgrid() - np.srt() - matplotlib.pyplot as plt - plt.clabel() - plt.colorbar() - quiver() - contour() - contourf() - scipy.ndimage.filters import gaussian_filter custom functions: - get_data() - set_background() as m source: - None ------------------------------------------------------------------- """ #---Extract data->custom function get_data-------------------------------------------------- #uq_data, vq_data, msl_data, lat, lon = get_data(IVT_infile, MSL_infile) uq_data, vq_data, lat, lon, msl_data = get_data(IVT_infile, mslp_infile) print ("get data done") #-------Customised color in RGB ------------ C = [[255,255,60], #grey #[255,255,255],#gre [247,228,53], #grey [239,199,47], #grey [232,172,40], #lillac, 39 64 197 149,53,229 [224,129,30], #blue dark,7,67,194 [218,81,14], #blue [215,35,6], #blue [202,20,18], #blue [158,9,59], #blue [114,1,97], #blue [76,0,133]] C = np.array( C ) C = np.divide( C, 255. ) # RGB has to be between 0 and 1 in python #------BAKGROUND MAP------------------------------------------------------------------- fig = plt.figure() m = set_background() #------PLOT DATA------------------------------------------------------------------- x, y = m( lon,lat ) IVT_size = np.sqrt( vq_data*vq_data+uq_data*uq_data ) #blurred = gaussian_filter( IVT_size, sigma=2 ) blurred = IVT_size #no blurr #blurred = ndimage.uniform_filter(IVT_size, size=7) pblurred = msl_data #pblurred=gaussian_filter( msl_data, sigma=7 ) """ |------Notes to expressions above------------------------------------------------------------------------------------------| | Should read about what is actually does..But it makes output smooth and nice | |--------------------------------------------------------------------------------------------------------------------------| """ pcont_val = np.arange( 800,2000,5 ) #defines which contourlines to plot contour_val = [250,300,400,500,600,700,800,1000,1200,1400,1600,1900] #defines which contourlines to plot #--->PLOTS PRESSURE ISOBARS----------------------------------------- CS = m.contour( x,y,pblurred,pcont_val, colors = "k", linewidths = 0.5, zorder=4 ) plt.clabel( CS, fontsize=7, inline=1,fmt = '%1.0f')#,ticks=Contourrange ) #--->PLOT IVT------------------------------------------------------ #----->Generates thinlines at contour interface-- print (blurred) CS = m.contour( x,y,blurred,contour_val, colors = "k", vmin = 0 , vmax=1600,linewidths = 0.2, zorder=3 ) #CS = m.contour( x,y,blurred, colors = "k", vmin = 0 , vmax=1600,linewidths = 0.2, zorder=3 ) #----->Generates the filled contours for IVT----- CS = m.contourf( x,y,blurred,contour_val, colors = C, vmin = 0 , vmax=1600, zorder=2 ) #CS = m.contourf( x,y,blurred, colors = C, vmin = 0 , vmax=1600, zorder=2 ) """ |------Notes to zorder------------------------------------------------------------------------------------------------------| | zorder makes it possible to tell which layer of plots should be above or under another layer | |---------------------------------------------------------------------------------------------------------------------------| """ cbar = plt.colorbar(CS, fraction=0.046, pad=0.01) cbar.ax.tick_params(labelsize=15) #NEW cbar.set_label('Integrated water vapor transport(IVT): [kg $m^{-1}$ $s^{-1}$]') """ |------Notes to colorbar-----------------------------------------------------------------------------------------------------| | Goal was to make colorbare constant on each plot so that it is easier to compare plots. To do this I had to specify which | | color to use(I made a customized color similar to what is made with GFS: http://cw3e.ucsd.edu/iwv-and-ivt-forecasts/) | | I had to define specific contours(which should match approx the same number as different colors). | | So parameters contour_val was made for this reason. | | I also specified vmin and vmax for the same reason of keeping colorbar constant, but dont know if it actually does | | anything | |----------------------------------------------------------------------------------------------------------------------------| """ #----VECTOR PLOT---------------------- yy = np.arange( 0, y.shape[ 0 ], 30 ) #skips over every 100th value 100 = inbetween xx = np.arange( 0, x.shape[ 1 ], 30 ) #skips over every 100th value points = np.meshgrid( yy, xx ) """ |------Notes to expressions above------------------------------------------------------------------------------------------| | - yy, xx spesifies which IVT vector to plot. We can't plot every vector(to messy), so specifies at which intervall we | | should plot. | | - points: |--------------------------------------------------------------------------------------------------------------------------| """ for i in range(0,IVT_size.shape[0]): for j in range(0,IVT_size.shape[1]): if IVT_size[i,j] <= 250: uq_data[i,j] = None vq_data[i,j] = None """ |-------Notes to for loop above-------------------------------------------------------------------------------------------| | Here I remove vectors with size less than 250, which is defines as the threshold value of atmospheric River. | | This will let us only plot IVT in those possible Atmospheric River regimes. | |-------------------------------------------------------------------------------------------------------------------------| """ scale = 3000. i = 1000. legend_size = i/scale + i qv = m.quiver( x[points], y[points], uq_data[points], vq_data[points], scale = scale,scale_units='inches',zorder=6) #-----#box--- ax = plt.gca() #f = fig.patches.extend([plt.Rectangle((0.767,0.89),0.12,0.042, # fill=False, color='k', alpha=0.8, zorder=9, # transform=fig.transFigure, figure=fig)]) #p = plt.quiverkey(qv, 0.96, 1.05 , legend_size , str(i) + " [kg $m^{-1}$$s^{-1}$]", coordinates='axes', zorder = 1, # color='k', labelcolor='k', labelpos = "S", labelsep= 0.02) f = fig.patches.extend([plt.Rectangle((0.863,0.89),0.092,0.042, fill=False, color='k', alpha=0.8, zorder=9, transform=fig.transFigure, figure=fig)]) p = plt.quiverkey(qv, 1.085, 1.05 , legend_size , "$10^{3}$kg $m^{-1}$$s^{-1}$", coordinates='axes', zorder = 1, color='k', labelcolor='k', labelpos = "S", labelsep= 0.02) #--Making Title nice formated------ date_format = "{6}{7} / {4}{5} - {0}{1}{2}{3}".format(*str(date)) date_format = "{0}{1}{2}{3}-{4}{5}-{6}{7} ".format(*str(date)) plt.title( date_format + time+" UTC", fontsize=15, position=(0.17, 0.93), bbox=dict(facecolor='white', alpha=0.7), zorder=1000) #plt.title("ECMWF-Analysis IVT (kg $m^{-1}$ $s^{-1}$; shaded), IVT Vector and MSLP(hPa; contours) \n " \ # +"Valid: "+time+" UTC, "+date_format, fontsize=11) fig.set_size_inches( 12.80, 7.15 ) fig.savefig( fig_name, dpi = 200 ) plt.close( ) def user_interface(): """ ------------------------------------------------------------------ Mission of function: - Handles Terminal input from user. - Makes it possible to run for a user defined date and time EXAMPLE OF TERMINAL RUN: - $ python IVT_independent.py --date 20160930 1800 PARAMETERS inputs: - None returns: - None other: - None USES installed functions: - arparse() custom functions: - IVT() source: - None ------------------------------------------------------------------- """ #----Sets up user interface from Terminal---------------------------- parser = argparse.ArgumentParser(description='Process some integers.') parser.add_argument( "--date", type = int, choices=[20171215,20171216, 20171217, 20171218, 20171219, 20171220, 20171221, 20171222, 20171223, 20171224, 20171225, 20171226], help = "the date you want. " ) parser.add_argument( "time", type = str, help = "the times you want. ex 0000, 0600, 1200, 1800", default = "all" ) args = parser.parse_args( ) #----------SET PATHS------------------------------------------- path_figures="../figs/" #path_grib_files="/Users/ainajohannessen/Documents/Aina/skole/master/Aina_2017storm/IVT/IVT_grib/" #path_norestore_files="/Users/ainajohannessen/Documents/Aina/skole/master/Masterthesis/data_norstore/" path_grib_files = "../gribs/" main_name = str(args.date) + "_" + args.time #example: 20160930_1800 fig_path_and_name = path_figures+ "IVT_"+main_name IVT_infile = path_grib_files + "IVT_"+main_name + ".grib" mslp_infile = path_grib_files + "param_msl_"+main_name + ".grib" #-------Calls main function------------------------------------- IVT(fig_path_and_name, IVT_infile, args.time, args.date, mslp_infile) warnings.filterwarnings("ignore",category=matplotlib.mplDeprecation) user_interface()