""" Aina Johannessen. 2017/2018 - Masterthesis """ from scipy.stats.kde import gaussian_kde #test import scipy.ndimage as ndimage import argparse import pygrib #To read grib files. import matplotlib import matplotlib.pyplot as plt import matplotlib.colors as colors from mpl_toolkits.basemap import Basemap #To make map from mpl_toolkits.basemap import shiftgrid #To make map from scipy.ndimage.filters import gaussian_filter # import numpy as np import math import warnings def get_data( obj, prm, lev, date, timelevel=0 ): """ ------------------------------------------------------------------ fatching the data and finds it latitude and longditude. parameters: - obj = object of the open file. - prm = parameter name that we want to find - lev = level. States which level we want to take or data from, feks 1000 hpa level. - date = which date to get - timelevel = which time we want. 0: time 0000UTC, 1: time 1200UTC, 2: time 1800UTC uses: Pygrib for reading grib file. returns: - lat, lon = lat,lon - data = values of the parameter. issues: Dont understand why I need to do lon-180 for it to use the whole map. Maybe this can be fixed when defining the drawing of the map or m object..? For now lon-180 works fine, but think it can be avoided. ------------------------------------------------------------------- """ parameter = obj( name = prm, level = lev, dataDate = date )[ timelevel ] print (parameter) #-----Checking grib type---------------------------------------------- if parameter.gridType == "sh": lat, lon, data = sh( parameter.values ) elif parameter.gridType == "reduced_gg": lat, lon = parameter.latlons() #very easy implementastion with a gg data = parameter.values elif parameter.gridType == "regular_gg": lat, lon = parameter.latlons() #very easy implementastion with a gg data = parameter.values else: print ( parameter.gridType ) x_dir = np.arange(0,lon.shape[0],1 ) y_dir = np.arange(0,lon.shape[1],1 ) for i in x_dir: data[i,:], lon[i,:] = shiftgrid(180, data[i,:], lon[i,:],start=False) #why is lon a 2dimensjonal thing..?? #data, lon = shiftgrid(180, data, lon, start=False) #for i in x_dir: # for k in y_dir: # if lon[i,k] > 180.: # lon[i,k] = lon[i,k] - 359. return lat, lon, data def plot_contourf( m, lat, lon, data, C, contur_val=None ): """ ------------------------------------------------------------------ Plot the data in a filed contour with color C Plots the colorbar parameters: - m = mab "object" from basemap. - lat = latitude - lon = longditude - data = the data values of our parameter - C = which Color to use. -contur_val = int:How many contours to plot. Alternatice colors: -#CS = m.contourf(x,y,data,15,cmap=plt.cm.jet) - See version 2.5 ------------------------------------------------------------------- """ x, y = m( lon,lat ) min = data.min() if contur_val.any(): CS = m.contourf( x, y, data, contur_val, colors = C, vmin = 264 , vmax=384 ) else: CS = m.contourf( x, y, data, colors = C, vmin = 264 , vmax=384 ) #cbar = plt.colorbar( drawedges = True ) # draw colorbar cbar = plt.colorbar(CS, fraction=0.046, pad=0.01) cbar.ax.tick_params(labelsize=15) #NEW cbar.set_label('Potential Temperature[ K ]', fontsize = 15) #cbar.set_label('[K]', fontsize=10, rotation=0, ) def plot_contour( m, lat, lon, data, contour, clr = 'k' ): """ ------------------------------------------------------------------ Plot the data in a contour with color xlr parameters: - m = mab "object" from basemap. - lat = latitude - lon = longditude - data = the data values of our parameter - contour - clr = Color to use. Issues: Had problems finding which contour values I should pot to make it nice. I did try different contour values for different levels, but didnt get any better Ended up with taking the "average" over the vertical levels and then filter it horizontally(see def DT), but its very time consuming, but looks better. Still not great though.. ------------------------------------------------------------------- """ x, y = m( lon,lat ) CS = m.contour( x, y, data, contour , linewidths = 0.8, colors = clr ) def map_area( m ): """ ------------------------------------------------------------------ Draws the background map based on m parameters: - m = mab "object" from basemap. It describes witch projection and latitudes to use for everything ------------------------------------------------------------------- """ m.drawcoastlines( linewidth = 1.5, linestyle = 'solid', color = [ 75./255., 75/255., 75/255. ] ) # ------draw parallels---------------- circles = np.arange( -90., 90. + 30, 20. ) #delat = 10. #circles = np.arange( -90., 90. + 30, 30. ) #delat = 30. #m.drawparallels( circles, labels = [ 1, 0, 0, 0 ] ) 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. #meridians = np.arange( 0., 360, 60. ) #delon = 60. #meridians = np.arange( -180., 180., 60. ) #delon = 60. #m.drawmeridians( meridians, labels = [ 0, 0, 0, 1 ] ) m.drawmeridians( meridians, color = [ 55./255., 55./255., 55./255. ], labels = [ 0, 0, 0, 1 ], linewidth = 0.1, fontsize=12 ) def plot_wind_bar( m, lat, lon, u, v ): """ ------------------------------------------------------------------ Plot wind bars for u and v parameters: - m = mab "object" from basemap. It describes witch projection and latitudes to use for everything - lat - lon - u - v source: #http://basemaptutorial.readthedocs.io/en/latest/plotting_data.html ------------------------------------------------------------------- """ x, y = m( lon, lat ) #---Have to limit number of windbarbs plotted----- yy = np.arange( 0, y.shape[ 0 ], 30 ) #skips over every 100th value xx = np.arange( 0, x.shape[ 1 ], 30 ) #skips over every 100th value points = np.meshgrid( yy, xx ) m.barbs( x[points], y[points], u[points], v[points], length = 5.5, pivot='middle', linewidth = 1., barbcolor = '#333333' ) def DT(time_lvl, date, DT_infile, path_figures): """ ------------------------------------------------------------------ Generates the Dynamical Tropopause(DT) map Parameter: -C: Found from measuring the color bar at the source. Uses: Functions: - map_area(m) - get_data( gfile, prm, lev, timelevel=0 ) - plot_contourf(m, lat, lon, data, C) - plot_contour(m, lat, lon, data, C) Issues: Tried using another basemap, one that is more curved, which is what I want, but the vorticity gets really weird in the North. I dont understand why, so I cant use it yet.. Source: http://www.atmos.albany.edu/student/abentley/realtime.html ------------------------------------------------------------------- """ #[ 250. , 257.33333333, 264.66666667, 272. , # 279.33333333, 286.66666667, 294. , 301.33333333, # 308.66666667, 316. , 323.33333333, 330.66666667, # 338. , 345.33333333, 352.66666667, 360. , # 367.33333333, 374.66666667, 382. , 389.33333333, # 396.66666667, 404. ]) #-------Customised color in RGB ------------ C = [[232,232,230],#grey 250. [203,203,203], #grey 257.33333333 [161,161,161], #grey 264.666 [130,130,130], #grey 272. [149,53,229], #lillac, 279.33333333, [39,64,197], #blue dark, 286.6666 [15,110,229], #blue 294. [80,149,240], #blue 301. [74,192,243], #blue 308.66666667, [152,219,248], #blue 316. [183,237,247], #blue 323.333 [251,217,198], #redish 330 [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], #389.3333 [217,155,210], # 396.66666667, [216,181,211]] # 404. C = np.array( C ) C = np.divide( C, 255. ) # RGB has to be between 0 and 1 in python #----------------------------------------------------------- fig = plt.figure() #-----Setting our map area and projection of interest------- #m = Basemap(projection='ortho', lat_0=-20, lon_0=-30, resolution='c') #m = Basemap( llcrnrlon = -170., llcrnrlat = 10., urcrnrlon = 50., urcrnrlat=80.,\ # resolution = 'h', area_thresh = 10000., projection = 'merc' ) #m = Basemap(projection='ortho',lat_0=90,lon_0=-50,resolution='l') m = Basemap( llcrnrlon = -90., llcrnrlat = 10., urcrnrlon = 50., urcrnrlat=70.,\ resolution = 'h', area_thresh = 10000., projection = 'merc' ) #m = Basemap(width=11500000,height=8500000,resolution='l',projection='lcc',\ # lat_1=07.,lat_2=40,lat_0=44,lon_0=-30. ) #m = Basemap(width=12000000,height=8000000, # resolution='l',projection='laea',\ # lat_ts=50,lat_0=50,lon_0=-30.) #projection='laea' """ #m=Basemap(llcrnrlon=-12.0, \ llcrnrlat=32.0, \ urcrnrlon=54.4, \ urcrnrlat=65.3, \ llcrnrx=None, \ llcrnry=None, \ urcrnrx=None, \ urcrnry=None, \ width=None, \ height=None, \ projection=projection, \ resolution='h', \ area_thresh=1000, \ rsphere=6370997.0, \ ellps=None, \ lat_ts=None, \ lat_1=None, \ lat_2=None, \ lat_0=60.0, \ lon_0=13.0, \ lon_1=None, \ lon_2=None, \ o_lon_p=None, \ o_lat_p=None, \ k_0=None, \ no_rot=True, \ suppress_ticks=True, \ satellite_height=35786000, \ boundinglat=None, \ fix_aspect=True, \ anchor='C', \ celestial=False, \ round=False, \ epsg=None, \ ax=None) """ #m = Basemap(width=190000,height=2200000,resolution='l', projection='tmerc',lon_0=-30,lat_0=44) map_area( m ) # ploting background #file = path_gribs +"DT_var.grib" obj = pygrib.open( DT_infile ) #-FETCHING ALL THE VALUES---------------------------------------- #-----Potential temperature--------------------------------------- lat, lon, data = get_data( obj,'Potential temperature', 2000, date, timelevel = time_lvl ) #contour_val = np.linspace( 264, 384, 22 ) #contours for potential tempeature contour_val = np.linspace( 250, 404, 22 ) #contours for potential tempeature contour_val = np.linspace( 250, 404, 22 ) #contours for potential tempeature #contour_val = np.linspace( 228, 404, 22 ) #contours for potential tempeature #[ 250. , 257.33333333, 264.66666667, 272. , # 279.33333333, 286.66666667, 294. , 301.33333333, # 308.66666667, 316. , 323.33333333, 330.66666667, # 338. , 345.33333333, 352.66666667, 360. , # 367.33333333, 374.66666667, 382. , 389.33333333, # 396.66666667, 404. ]) plot_contourf( m, lat, lon, data, C, contour_val ) #-----Relative vorticity, diff level------------------------------0.25 10 4 s 1 intervals above 1 10 4 s 1 in black contours) #contour=[ 1.E-4, 1.5E-4, 2.E-4, 2.5E-4, 2.8E-4, 3.5E-4, 4.5E-4, 6.5E-4, 7.E-4, 7.5E-4, 8.E-4 ] #1.5E-4,2.5E-4]# #contour = np.arange(1.E-4, 8.E-4, 0.25E-4 ) #contour = np.arange(1.E-4, 9.E-4, 1.E-4 ) contour = np.arange(0.5E-4, 7.5E-4, 0.5E-4 ) #contour = [1.E-4, 1.5E-4, 2.E-4, 2.5E-4] #contour=[ 1.E-4, 1.5E-4, 2.E-4, 2.5E-4, 2.8E-4, 3.5E-4, 4.5E-4, 6.5E-4, 7.E-4, 7.5E-4, 8.E-4 ] #1.5E-4,2.5E-4]# lat, lon, data925 = get_data( obj, 'Vorticity (relative)', 925, date, timelevel = time_lvl ) lat, lon, data900 = get_data( obj, 'Vorticity (relative)', 900, date, timelevel = time_lvl ) lat, lon, data850 = get_data( obj, 'Vorticity (relative)', 850, date, timelevel = time_lvl ) print ("size of lat ",lat.size) print ("shape of lat",lat.shape) #->--->---->--mean value over height and filtering---------------- #data = np.sqrt( data900**2 + 2*data850**2 + data925**2 ) #Vertical "average", weightet values at 850hpa double. data = ( data900 + data850 + data925 )/3 data[np.where(data < 1.E-5 )] = 0 data[np.where(data > 10 )] = 0 #footprint = np.array([[0,0,0,1,1,1,1,0,0,0], #footprint=np.ones((3,10)) # [0,0,1,1,1,2,1,1,0,0], # [1,1,1,2,2,1,2,1,1,1], # [0,1,1,1,1,2,1,1,1,0], # [0,0,1,1,1,1,1,1,0,0]]) #data = ndimage.generic_filter( data, np.mean, footprint = footprint, mode='wrap' ) data =gaussian_filter( data, sigma=1 ) plot_contour( m, lat,lon, data,contour, clr = 'k') #linewidth #-----Wind barbs---------------------------------------------------- lat, lon, data_u = get_data( obj , 'U component of wind', 2000, date, timelevel = time_lvl ) lat, lon, data_v = get_data( obj , 'V component of wind', 2000, date, timelevel = time_lvl ) plot_wind_bar( m, lat, lon, data_u, data_v ) #----------------------------------------------- #----------------------------------------------- #-SAVE AND CLOSE---------------------------------------------------- #------------------------------------------------------------------ obj.close() if time_lvl == 0: t = "0000" elif time_lvl == 1: t = "0600" elif time_lvl == 2: t = "1200" elif time_lvl == 3: t = "1800" else: t = "t_not_set" fig_name = path_figures+"DT_" + str( date ) + "_" + str( t )#+ ".TIFF" #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 DT; Potential temperature at pv level = 2 (K; shaded) \n \ # Wind at (m $s^{-1}$; barbs), Low level relative vorticity($s^{-1}$; contours) \n " \ # +"Valid: "+t+" UTC, "+date_format, fontsize=10) #plt.title( t+" UTC, "+date_format, fontsize=15, position=(0.15, 0.9), bbox=dict(facecolor='white', alpha=0.5)) plt.title( date_format + t+" UTC", fontsize=15, position=(0.17, 0.93), bbox=dict(facecolor='white', alpha=0.5)) ax = plt.gca( ) plt.rc( 'font', size = 6 ) fig.set_size_inches( 12.80, 7.15 ) fig.savefig( fig_name, dpi = 200 ) plt.close( ) #plt.show() #-------------------------- #---------------------------- def user_interface(): parser = argparse.ArgumentParser(description='Process some integers.') t = [] #timelevel parser.add_argument( "--date", type = int, choices= [ 20160920,20160921,20160922,20160923,20160924,20160925,20160926,20160927, 20160928, 20160929, 20160930, 20161001], help = "the date you want. " ) parser.add_argument( "time", type = str, help = "the times you want. ex 160924, 160926, 160922, 160930", default = "all" ) args = parser.parse_args( ) date = [ 20160920, 20160921, 20160922, 20160923, 20160924, 20160925, 20160926, 20160927, 20160928, 20160929, 20160930, 20161001 ] if args.date: date = [args.date] if args.time == "all": t = [ 0, 1, 2, 3 ] elif args.time == "0000" or args.time == "00" or args.time == "0": t = [ 0 ] elif args.time == "0600" or args.time == "06" or args.time == "6" or args.time == "1": t = [ 1 ] elif args.time == "1200" or args.time == "12" or args.time == "2": t = [ 2 ] elif args.time == "1800" or args.time == "18" or args.time == "3": t = [ 3 ] print ("\n-------------------------------------------------------") print ("you chosed, time: "+args.time+"corresponding to timelevel: "+str(t)) print ("you chosed, date: "+str(args.date) ) print ("---------------------------------------------------------") print("Creating a dynamic tropopause....\n ...") #---PATHS----------------------------------------------------------------------------------- path_gribs="../gribs/" path_figures="../figs/" DT_infile = path_gribs +"DT_var.grib" #------------------------------------------------------------------------------------------- #------Setting time and date for plots---------------------- for d in date: print ("one day made") for n in t: DT(n, d, DT_infile, path_figures) print("one time made") #------------------------------------------------------------ #matplotlib warns me that I am using a function(hold = on) that is not supported in updatet version of matplotlib. #Though I am not using it, so probably a package that uses it. Might be a problem if I dont update all the packages correctly. warnings.filterwarnings("ignore",category=matplotlib.mplDeprecation) user_interface()