Update to establish running version with plotting of p-lat and p-lon crosssections

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':[]}
+# 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 '----------------------------------------------'
-print 'Insert code to plot full 2d fields'
+print 'Plot 2d fields (lat-lon)'
-# define plotname
+list_allkeys = data.keys() # all variables
-list_allkeys = data.keys()
 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 '----------------------------------------------'