diff --git a/calc_ref_state.py b/calc_ref_state.py index b393eab..0246953 100644 --- a/calc_ref_state.py +++ b/calc_ref_state.py @@ -12,6 +12,7 @@ import time # to measure elapsed time from get_data import get_data_woa2009 as read_data # plotting modules / functions from myplot import myplot_2dmap # plot of lat-lon fields +from myplot import myplot_2d # plot of crossections fields from myplot_inputs import myplot_create_pn # create plotname # ############################################################################################ @@ -23,11 +24,17 @@ from myplot_inputs import myplot_create_pn # create plotname # variables : variables necessary for calculations # MANUAL DEFINED OUTPUTS MISSING # ############################################################################################ +# information about input data print 'MANUAL DEFINED OUTPUTS MISSING' -dim = ('lon','lat','z','time',) -dir_plot = '/glusterfs/inspect/users/xg911182/Code/Python/plots_test/' -variables = {'s':{},'temp':{},'grd':[],'valid':[]} +dim = ('lon','lat','z','time',) +variables = {'s':{},'temp':{},'grd':[],'valid':[]} +# plotting input +saveplot = input('Do you want to save the plots? Enter 1 for yes or 0 for no') +dir_plot = '/glusterfs/inspect/users/xg911182/Code/Python/plots_test/' lonlat_reg = [25.0, 25.0+360.0, -90.0, 90.0] # plotting range for 2dmap +plot_2d = 1 # plot 2d fields (lat-lon) +plot_cs = 1 # plot crosssections +cs_plot = {'lon':[-47.5,172.5],'lat':[0]} # crosssections # Input 1 #################################################################################### # Physical constants # input : None @@ -53,13 +60,13 @@ lonlat_reg = [25.0, 25.0+360.0, -90.0, 90.0] # plotting range for 2dmap # ############################################################################################ beg_time = time.time() print '' -print ' ----------------------------------------------' +print '----------------------------------------------' print ' Read data' grd, data = read_data(variables,dim) -print ' ----------------------------------------------' +print '----------------------------------------------' end_time = time.time() print 'Elapsed time to read data: '+ '{:.2f}'.format(end_time-beg_time)+'s' -print ' ----------------------------------------------' +print '----------------------------------------------' # step b ##################################################################################### # Define grid for data from step a @@ -69,12 +76,12 @@ print ' ----------------------------------------------' # ############################################################################################ beg_time = time.time() print '' -print ' ----------------------------------------------' +print '----------------------------------------------' print 'Insert code to define grid of data' -print ' ----------------------------------------------' +print '----------------------------------------------' end_time = time.time() print 'Elapsed time to define grid: '+ '{:.2f}'.format(end_time-beg_time)+'s' -print ' ----------------------------------------------' +print '----------------------------------------------' # step c ##################################################################################### # Calculation of density @@ -87,12 +94,12 @@ print ' ----------------------------------------------' # ############################################################################################ beg_time = time.time() print '' -print ' ----------------------------------------------' +print '----------------------------------------------' print 'Insert code to calculate density' -print ' ----------------------------------------------' +print '----------------------------------------------' end_time = time.time() print 'Elapsed time to calculate density: '+ '{:.2f}'.format(end_time-beg_time)+'s' -print ' ----------------------------------------------' +print '----------------------------------------------' # step d ##################################################################################### # Calculate z-dependent ocean area @@ -105,12 +112,12 @@ print ' ----------------------------------------------' # ############################################################################################ beg_time = time.time() print '' -print ' ----------------------------------------------' +print '----------------------------------------------' print 'Insert code to calculate ocean area' -print ' ----------------------------------------------' +print '----------------------------------------------' end_time = time.time() print 'Elapsed time to calc. ocean area: '+ '{:.2f}'.format(end_time-beg_time)+'s' -print ' ----------------------------------------------' +print '----------------------------------------------' # step e ##################################################################################### # Calculate reference state @@ -123,12 +130,12 @@ print ' ----------------------------------------------' # ############################################################################################ beg_time = time.time() print '' -print ' ----------------------------------------------' +print '----------------------------------------------' print 'Insert code to calculate reference state' -print ' ----------------------------------------------' +print '----------------------------------------------' end_time = time.time() print 'Elapsed time to calculate reference state: '+ '{:.2f}'.format(end_time-beg_time)+'s' -print ' ----------------------------------------------' +print '----------------------------------------------' # step f ##################################################################################### # plot reference state @@ -139,31 +146,33 @@ print ' ----------------------------------------------' # ############################################################################################ beg_time = time.time() print '' -print ' ----------------------------------------------' +print '----------------------------------------------' print 'Insert code to plot reference state' -print ' ----------------------------------------------' +print '----------------------------------------------' end_time = time.time() print 'Elapsed time to plot and save reference state: '+ '{:.2f}'.format(end_time-beg_time)+'s' -print ' ----------------------------------------------' +print '----------------------------------------------' # step g ##################################################################################### # plot of full 2d fields -# input : None +# input : date_str, dir_plot, i_dep, i_t, lonlat_reg +# date_str : MISSING - DONE LOCALLY HERE +# dir_plot : output of Input 0 +# i_dep : MISSING - DONE LOCALLY HERE +# i_t : MISSING - DONE LOCALLY HERE +# lonlat_reg: output of Input 0 # output : None # saved/created : plots saved in dir_plot -# dir_plot : output of Input 0 # ############################################################################################ beg_time = time.time() print '' -print ' ----------------------------------------------' -print 'Insert code to plot full 2d fields' -# define plotname -list_allkeys = data.keys() +print '----------------------------------------------' +print 'Plot 2d fields (lat-lon)' +list_allkeys = data.keys() # all variables i_dep = 0 i_t = 0 date_str = '2009' for i_key in range(0,len(list_allkeys)): - print 'put this plotname stuff into a function' # create plotname plotname = myplot_create_pn({'dir_plot':dir_plot,'beg_str':'Map2d',\ 'varname':list_allkeys[i_key],'zlev':[grd.z[i_dep],grd.Uz],\ @@ -175,8 +184,66 @@ for i_key in range(0,len(list_allkeys)): myplot_2dmap(data_in,grd,lonlat_range=lonlat_reg,saveplot=1,\ title_c=title_in,plotname=plotname,unit_data=data[list_allkeys[i_key]]['units']) del data_in, title_in, plotname -print ' ----------------------------------------------' +print '----------------------------------------------' end_time = time.time() print 'Elapsed time to plot 2d fields (lat-lon): '+ '{:.2f}'.format(end_time-beg_time)+'s' -print ' ----------------------------------------------' +print '----------------------------------------------' +# step h ##################################################################################### +# plot of crossection fields +# input : date_str, dir_plot, i_t, plot_cs, saveplot +# date_str : MISSING - DONE LOCALLY HERE +# dir_plot : output of Input 0 +# i_t : MISSING - DONE LOCALLY HERE +# plot_cs : output of Input 0 +# saveplot : output of Input 0 +# output : None +# saved/created : plots saved in dir_plot +# ############################################################################################ +if plot_cs: + beg_time = time.time() + print '' + print '----------------------------------------------' + print 'Insert code to plot crossections' + list_allkeys = data.keys() # all variables + for i_key in range(0,len(list_allkeys)): + for i_cs in range(0,len(cs_plot['lat'])): + # load fixed latitude and associated index for grd.lat + lat_fix = cs_plot['lat'][i_cs] + i_lat = abs(grd.lat-cs_plot['lat'][i_cs]).argmin() + # create plotname + plotname = myplot_create_pn({'dir_plot':dir_plot,'beg_str':'CS',\ + 'varname':list_allkeys[i_key],'lonlat':[0,360,lat_fix,lat_fix],\ + 'date':date_str,'end_str':'.png'}) + # define plot input + title_in = 'Depth='+'{:.1f}'.format(grd.z[i_dep])+grd.Uz + cbarlabel_in = data[list_allkeys[i_key]]['standard_name'] + \ + ' ['+data[list_allkeys[i_key]]['units']+']' + data_in = np.transpose(np.squeeze(data[list_allkeys[i_key]]['val'][:,i_lat,:,i_t]),[1,0]) + data_in[data_in==data[list_allkeys[i_key]]['fill_value']] = np.nan + myplot_2d(grd.lon,-grd.z,data_in,saveplot=saveplot,FS=14,d_xtick=45,\ + xlabel_c='Longitude ['+grd.Ulon+']',ylabel_c='Depth ['+grd.Uz+']',\ + title_c=title_in,plotname=plotname,cbarlabel_c=cbarlabel_in) + del data_in, lat_fix, i_lat, title_in, plotname, cbarlabel_in + for i_cs in range(0,len(cs_plot['lon'])): + # load fixed longitude and associated index for grd.lon + lon_fix = cs_plot['lon'][i_cs] + i_lon = abs(grd.lon-cs_plot['lon'][i_cs]).argmin() + # create plotname + plotname = myplot_create_pn({'dir_plot':dir_plot,'beg_str':'CS',\ + 'varname':list_allkeys[i_key],'lonlat':[lon_fix,lon_fix,-90,90],\ + 'date':date_str,'end_str':'.png'}) + # define plot input + title_in = 'Depth='+'{:.1f}'.format(grd.z[i_dep])+grd.Uz + cbarlabel_in = data[list_allkeys[i_key]]['standard_name'] + \ + ' ['+data[list_allkeys[i_key]]['units']+']' + data_in = np.transpose(np.squeeze(data[list_allkeys[i_key]]['val'][i_lon,:,:,i_t]),[1,0]) + data_in[data_in==data[list_allkeys[i_key]]['fill_value']] = np.nan + myplot_2d(grd.lat,-grd.z,data_in,saveplot=saveplot,FS=14,d_xtick=45,\ + xlabel_c='Latitude ['+grd.Ulat+']',ylabel_c='Depth ['+grd.Uz+']',\ + title_c=title_in,plotname=plotname,cbarlabel_c=cbarlabel_in) + del data_in, lon_fix, i_lon, title_in, plotname, cbarlabel_in + print '----------------------------------------------' + end_time = time.time() + print 'Elapsed time to plot crossections: '+ '{:.2f}'.format(end_time-beg_time)+'s' + print '----------------------------------------------'