import netCDF4 from netCDF4 import Dataset import numpy as np 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 import warnings import argparse from msd_regions import * import pygrib #Used in get_data() to handle files. 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 get_data(IVT_infile, MSL_infile ): #--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 IVT( fig,m, fig_name, IVT_infile, time, date, mslp_infile ): #---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------------------------------------------------------------------- #------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-- 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}$]', fontsize = 15) """ |------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("ECMWF-Analysis IVT (kg $m^{-1}$ $s^{-1}$; shaded), IVT Vector and MSLP(hPa; contours) \n " \ # +"Valid: "+time+" UTC, "+date_format, fontsize=11) plt.title( date_format + time+" UTC", fontsize=15, position=(0.17, 0.93), bbox=dict(facecolor='white', alpha=0.7), zorder=1000) def uptake(fig, m,infile): print("test") lat=np.arange(-90,90.5,0.5) lon=np.arange(-180,180.5,0.5) obj_up = Dataset(infile) msd=np.array(obj_up.variables['N'][0,0]*10**5) areas = return_cell_areas(lon, lat) #Normalisation of uptakes uptake_norm = np.array(normalized_uptake(msd, areas)) print(obj_up) flagmask=(msd<0.001) uptake_norm[flagmask]=np.nan masked_array=np.ma.array(uptake_norm,mask=np.isnan(uptake_norm)) longrid, latgrid = np.meshgrid(lon,lat) cmap=plt.cm.get_cmap('winter') levels=np.array([0.0001,0.005,0.01,0.02,0.04,0.06,0.08,0.1,0.12,0.14,0.16,0.18,0.2]) CS=m.contour(longrid, latgrid,uptake_norm, latlon=True,extend='max',levels=levels, cmap=cmap)#, colors = "g", linewidths = 2) #h1=m.plot(longrid, latgrid,uptake_norm, latlon=True,extend='max',levels=levels) cbaxes = fig.add_axes([0.19, 0.08, 0.67, 0.02]) #left, bottom, width, height], fractions that go from 0 to 1 ofthe plotting area. cbar = plt.colorbar(CS, fraction=0.046, cax = cbaxes, orientation="horizontal") cbar.ax.tick_params(labelsize=15) #NEW lines1=cbar.ax.get_children() print (lines1) numl = len(lines1[0].get_linewidths()) #lines1[0].set_color(['r']*11) lines1[0].set_linewidths([5]*numl) #plt.setp(CS.collections , linewidth=2) #cbar = plt.colorbar(CS, fraction=0.046 , orientation="vertical") cbar.set_label('% $(10^5$km$^2)^{-1}$ lower limit <0.0001 % $(10^5$km$^2)^{-1}$', labelpad=-1, size=15) def user_interface(): #----Sets up user interface from Terminal---------------------------- parser = argparse.ArgumentParser(description='Process some integers.') parser.add_argument( "--date", type = int, choices= [ 20171204,20171222,20171223,20171224,20160927, 20160928, 20160929, 20160930], 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------------------------------------------- event = "Walpurga" if event == "Walpurga": path_figures="/Users/ainajohannessen/Documents/Aina/skole/master/trajectoryruns/uptake/IVT_with_uptake/" #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 = path_norestore_files uptake_path= "/Users/ainajohannessen/Documents/Aina/skole/master/trajectoryruns/ETH_visit/Franziska/walpurga_large/" date_time = str(args.date) + "_" + args.time #example: 20160930_1800 analysis_date =datetime.strptime(date_time, '%Y%m%d_%H00') fig_outfile = path_figures+ "IVT_UPT_"+date_time+".png" IVT_infile = path_grib_files + "AR_2016/IVT/analysis/" + date_time + ".grib" mslp_infile = path_grib_files+"SCA/sfc/msl/" + "param_msl_"+date_time + ".grib" uptake_infile = uptake_path + "UTOT"+date_time[:-2]+".cdf" #uptake_infile = uptake_path + "UBL"+date_time[:-2]+".cdf" #uptake_infile = uptake_path + "UABL"+date_time[:-2]+".cdf" # -----start plots----- elif event =="Birk": path_figures="/Users/ainajohannessen/Documents/Aina/skole/master/trajectoryruns/uptake/IVT_with_uptake/" #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 = path_norestore_files uptake_path= "/Users/ainajohannessen/Documents/Aina/skole/master/trajectoryruns/ETH_visit/Franziska/birk/" date_time = str(args.date) + "_" + args.time #example: 20160930_1800 analysis_date =datetime.strptime(date_time, '%Y%m%d_%H00') fig_outfile = path_figures+ "IVT_UPT_"+date_time+".png" IVT_infile = path_grib_files + "Birk_2017/syn_IVT/" + "IVT_"+date_time + ".grib" mslp_infile = path_grib_files+ "Birk_2017/mslp/param_msl_"+date_time + ".grib" uptake_infile = uptake_path + "UTOT"+date_time[:-2]+".cdf" #uptake_infile = uptake_path + "UBL"+date_time[:-2]+".cdf" #uptake_infile = uptake_path + "UABL"+date_time[:-2]+".cdf" elif event =="Aina": path_figures="/Users/ainajohannessen/Documents/Aina/skole/master/trajectoryruns/uptake/IVT_with_uptake/" #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 = path_norestore_files uptake_path= "/Users/ainajohannessen/Documents/Aina/skole/master/trajectoryruns/ETH_visit/Franziska/aina/" date_time = str(args.date) + "_" + args.time #example: 20160930_1800 analysis_date =datetime.strptime(date_time, '%Y%m%d_%H00') fig_outfile = path_figures+ "IVT_UPT_"+date_time+".png" IVT_infile = path_grib_files + "Aina_2017/syn_IVT/" + "IVT_"+date_time + ".grib" mslp_infile = path_grib_files+ "Aina_2017/mslp/param_msl_"+date_time + ".grib" uptake_infile = uptake_path + "UTOT"+date_time[:-2]+".cdf" #uptake_infile = uptake_path + "UBL"+date_time[:-2]+".cdf" #uptake_infile = uptake_path + "UABL"+date_time[:-2]+".cdf" fig = plt.figure() m = set_background() #-------Calls main function------------------------------------- IVT(fig,m,fig_outfile, IVT_infile, args.time, args.date, mslp_infile) print ("IVT DONE") uptake(fig, m, uptake_infile) print ("DQ DONE") #traj(fig,m, analysis_date, date_time) #print ("taj DONE") print ("and noooooow") #plt.title(date_time) fig.set_size_inches( 12.80, 7.15 ) #plt.show() #-----SAVING FIG------------------------ fig_outfile = "./plotsss/grad.png" plt.savefig( fig_outfile, dpi = 200 ) print ("is it saved?") #plt.close( ) warnings.filterwarnings("ignore",category=matplotlib.mplDeprecation) user_interface()