Refactoring of vrfy - transferred bulk of plotting to sub-script

scripts/exgdas_global_marine_analysis_vrfy.py

@@ -22,189 +22,147 @@
 import gen_eva_obs_yaml
 import marine_eva_post
 import diag_statistics
-from soca_vrfy import plot_config, plot_horizontal_slice, plot_zonal_slice
+from soca_vrfy import statePlotter, plotConfig
 import subprocess
 from datetime import datetime, timedelta
 
 comout = os.getenv('COM_OCEAN_ANALYSIS')
 com_ice_history = os.getenv('COM_ICE_HISTORY_PREV')
 com_ocean_history = os.getenv('COM_OCEAN_HISTORY_PREV')
-data = os.getenv('DATA')
-pdy = os.getenv('PDY')
 cyc = os.getenv('cyc')
 RUN = os.getenv('RUN')
-bcyc = str((int(cyc) - 3) % 24).zfill(2)
 gcyc = str((int(cyc) - 6) % 24).zfill(2)
-gcdate = datetime.strptime(os.getenv('PDY')+os.getenv('cyc'), '%Y%m%d%H') - timedelta(hours=int(os.getenv('assim_freq')))
 
+bcyc = str((int(cyc) - 3) % 24).zfill(2)
 grid_file = os.path.join(comout, f'{RUN}.t'+bcyc+'z.ocngrid.nc')
 
 # for eva
 diagdir = os.path.join(comout, 'diags')
 HOMEgfs = os.getenv('HOMEgfs')
 
-
 #######################################
-# INCREMENT
+# ocean increment
 #######################################
-incr_cmap = 'RdBu'
+
 data_file = os.path.join(comout, f'{RUN}.t'+cyc+'z.ocninc.nc')
-config = plot_config(grid_file=grid_file,
-                     data_file=data_file,
-                     colormap=incr_cmap,
-                     comout=os.path.join(comout, 'vrfy', 'incr'))
+config = plotConfig(grid_file=grid_file,
+                    data_file=data_file,
+                    lats=np.arange(-60, 60, 10),
+                    variables_zonal=['Temp', 'Salt'],
+                    variables_horiz=['Temp', 'Salt', 'ave_ssh'],
+                    allbounds={'Temp': [-0.5, 0.5],
+                               'Salt': [-0.1, 0.1],
+                               'ave_ssh': [-0.1, 0.1]},
+                    colormap='RdBu',
+                    comout=os.path.join(comout, 'vrfy', 'incr'))
+ocnIncPlotter = statePlotter(config)
+ocnIncPlotter.plot()
 
 #######################################
-# zonal slices
-
-for lat in np.arange(-60, 60, 10):
-
-    for max_depth in [700.0, 5000.0]:
-        config['lats'] = [lat]
-        config['max depth'] = max_depth
-
-        # Temperature
-        config.update({'variable': 'Temp', 'levels': [1], 'bounds': [-.5, .5]})
-        plot_zonal_slice(config)
-
-        # Salinity
-        config.update({'variable': 'Salt', 'levels': [1], 'bounds': [-.1, .1]})
-        plot_zonal_slice(config)
-
+# sea ice increment
 #######################################
-# Horizontal slices
-
-# Temperature
-config.update({'variable': 'Temp', 'levels': [1], 'bounds': [-1, 1]})
-plot_horizontal_slice(config)
 
-# Salinity
-config.update({'variable': 'Salt', 'bounds': [-0.1, 0.1]})
-plot_horizontal_slice(config)
-
-# Sea surface height
-config.update({'variable': 'ave_ssh', 'bounds': [-0.1, 0.1]})
-plot_horizontal_slice(config)
-
-#######################################
-# Sea ice
 data_file = os.path.join(comout, f'{RUN}.t'+cyc+'z.ice.incr.nc')
-config = plot_config(grid_file=grid_file,
-                     data_file=data_file,
-                     colormap=incr_cmap,
-                     comout=os.path.join(comout, 'vrfy', 'incr'))
-
-for proj in ['North', 'South']:
-    # concentration
-    config.update({'variable': 'aicen', 'bounds': [-0.2, 0.2], 'proj': proj})
-    plot_horizontal_slice(config)
-
-    # thickness
-    config.update({'variable': 'hicen', 'bounds': [-0.5, 0.5], 'proj': proj})
-    plot_horizontal_slice(config)
+config = plotConfig(grid_file=grid_file,
+                    data_file=data_file,
+                    lats=np.arange(-60, 60, 10),
+                    variables_horiz=['aicen', 'hicen', 'hsnon'],
+                    allbounds={'aicen': [-0.2, 0.2],
+                               'hicen': [-0.5, 0.5],
+                               'hsnon': [-0.1, 0.1]},
+                    colormap='RdBu',
+                    projs=['North', 'South'],
+                    comout=os.path.join(comout, 'vrfy', 'incr'))
+iceIncPlotter = statePlotter(config)
+iceIncPlotter.plot()
 
-    # snow depth
-    config.update({'variable': 'hsnon', 'bounds': [-0.1, 0.1], 'proj': proj})
-    plot_horizontal_slice(config)
-
-#######################################
-# Analysis/Background
 #######################################
-
+# sea ice analysis
 #######################################
-# Sea ice
-data_files = [os.path.join(comout, f'{RUN}.t{cyc}z.iceana.nc'),
-              os.path.join(com_ice_history, f'{RUN}.t{gcyc}z.icef006.nc')]
-dirs_out = ['ana', 'bkg']
-ice_vars = {'bkg': ['aice_h', 'hs_h', 'hi_h'], 'ana': ['aicen', 'hicen', 'hsnon']}
-for data_file, dir_out in zip(data_files, dirs_out):
-    config = plot_config(grid_file=grid_file,
-                         data_file=data_file,
-                         colormap='jet',
-                         comout=os.path.join(comout, 'vrfy', dir_out))
-
-    for proj in ['North', 'South', 'Global']:
-        # concentration
-        var = ice_vars[dir_out]
-        config.update({'variable': var[0], 'bounds': [0.0, 1.0], 'proj': proj})
-        plot_horizontal_slice(config)
 
-        # thickness
-        config.update({'variable': var[1], 'bounds': [0.0, 4.0], 'proj': proj})
-        plot_horizontal_slice(config)
-
-        # snow depth
-        config.update({'variable': var[2], 'bounds': [0.0, 0.5], 'proj': proj})
-        plot_horizontal_slice(config)
+data_file = os.path.join(comout, f'{RUN}.t'+cyc+'z.iceana.nc')
+config = plotConfig(grid_file=grid_file,
+                    data_file=data_file,
+                    variables_horiz=['aicen', 'hicen', 'hsnon'],
+                    allbounds={'aicen': [0.0, 1.0],
+                               'hicen': [0.0, 4.0],
+                               'hsnon': [0.0, 0.5]},
+                    colormap='jet',
+                    projs=['North', 'South', 'Global'],
+                    comout=os.path.join(comout, 'vrfy', 'ana'))
+iceAnaPlotter = statePlotter(config)
+iceAnaPlotter.plot()
 
 #######################################
-# Ocean surface
-data_files = [os.path.join(comout, f'{RUN}.t'+cyc+'z.ocnana.nc'),
-              os.path.join(com_ocean_history, f'{RUN}.t{gcyc}z.ocnf006.nc')]
-dirs_out = ['ana', 'bkg']
-ocn_vars = ['ave_ssh', 'Temp', 'Salt']
-for data_file, dir_out in zip(data_files, dirs_out):
-    config = plot_config(grid_file=grid_file,
-                         data_file=data_file,
-                         colormap='jet',
-                         comout=os.path.join(comout, 'vrfy', dir_out))
-
-    # ssh
-    config.update({'variable': 'ave_ssh', 'bounds': [-1.8, 1.3], 'proj': proj, 'levels': [1]})
-    plot_horizontal_slice(config)
-
-    # sst
-    config.update({'variable': 'Temp', 'bounds': [-1.8, 34.0], 'proj': proj, 'levels': [1]})
-    plot_horizontal_slice(config)
+# sea ice background
+#######################################
 
-    # sss
-    config.update({'variable': 'Salt', 'bounds': [30, 38], 'proj': proj, 'levels': [1]})
-    plot_horizontal_slice(config)
+data_file = os.path.join(com_ice_history, f'{RUN}.t{gcyc}z.icef006.nc')
+config = plotConfig(grid_file=grid_file,
+                    data_file=data_file,
+                    variables_horiz=['aice_h', 'hs_h', 'hi_h'],
+                    allbounds={'aice_h': [0.0, 1.0],
+                               'hs_h': [0.0, 4.0],
+                               'hi_h': [0.0, 0.5]},
+                    colormap='jet',
+                    projs=['North', 'South', 'Global'],
+                    comout=os.path.join(comout, 'vrfy', 'bkg'))
+iceBkgPlotter = statePlotter(config)
+iceBkgPlotter.plot()
 
 #######################################
-# Std Bkg. Error
+# ocean surface analysis
 #######################################
-bmat_cmap = 'jet'
-data_file = os.path.join(comout, f'{RUN}.t'+cyc+'z.ocn.bkgerr_stddev.nc')
-config = plot_config(grid_file=grid_file,
-                     data_file=data_file,
-                     colormap=bmat_cmap,
-                     comout=os.path.join(comout, 'vrfy', 'bkgerr'))
-
-#######################################
-# zonal slices
-
-for lat in np.arange(-60, 60, 10):
-
-    for max_depth in [700.0, 5000.0]:
-        config['lats'] = [lat]
-        config['max depth'] = max_depth
 
-        # Temperature
-        config.update({'variable': 'Temp', 'levels': [1], 'bounds': [0, 1.5]})
-        plot_zonal_slice(config)
-
-        # Salinity
-        config.update({'variable': 'Salt', 'levels': [1], 'bounds': [0, .2]})
-        plot_zonal_slice(config)
+data_file = os.path.join(comout, f'{RUN}.t'+cyc+'z.ocnana.nc')
+config = plotConfig(grid_file=grid_file,
+                    data_file=data_file,
+                    variables_horiz=['ave_ssh', 'Temp', 'Salt'],
+                    allbounds={'ave_ssh': [-1.8, 1.3],
+                               'Temp': [-1.8, 34.0],
+                               'Salt': [30, 38]},
+                    colormap='jet',
+                    comout=os.path.join(comout, 'vrfy', 'ana'))
+ocnAnaPlotter = statePlotter(config)
+ocnAnaPlotter.plot()
 
 #######################################
-# Horizontal slices
+# ocean surface background
+#######################################
 
-# Temperature
-config.update({'variable': 'Temp', 'levels': [1], 'bounds': [0, 2]})
-plot_horizontal_slice(config)
+data_file = os.path.join(com_ocean_history, f'{RUN}.t{gcyc}z.ocnf006.nc')
+config = plotConfig(grid_file=grid_file,
+                    data_file=data_file,
+                    variables_horiz=['ave_ssh', 'Temp', 'Salt'],
+                    allbounds={'ave_ssh': [-1.8, 1.3],
+                               'Temp': [-1.8, 34.0],
+                               'Salt': [30, 38]},
+                    colormap='jet',
+                    comout=os.path.join(comout, 'vrfy', 'bkg'))
+ocnBkgPlotter = statePlotter(config)
+ocnBkgPlotter.plot()
 
-# Salinity
-config.update({'variable': 'Salt', 'bounds': [0, 0.2]})
-plot_horizontal_slice(config)
+#######################################
+# background error
+#######################################
 
-# Sea surface height
-config.update({'variable': 'ave_ssh', 'bounds': [0, 0.1]})
-plot_horizontal_slice(config)
+data_file = os.path.join(comout, f'{RUN}.t'+cyc+'z.ocn.bkgerr_stddev.nc')
+config = plotConfig(grid_file=grid_file,
+                    data_file=data_file,
+                    lats=np.arange(-60, 60, 10),
+                    variables_zonal=['Temp', 'Salt'],
+                    variables_horiz=['Temp', 'Salt', 'ave_ssh'],
+                    allbounds={'Temp': [0, 2],
+                               'Salt': [0, 0.2],
+                               'ave_ssh': [0, 0.1]},
+                    colormap='jet',
+                    comout=os.path.join(comout, 'vrfy', 'bkgerr'))
+bkgErrPlotter = statePlotter(config)
+bkgErrPlotter.plot()
 
 #######################################
 # eva plots
+#######################################
 
 evadir = os.path.join(HOMEgfs, 'sorc', f'{RUN}.cd', 'ush', 'eva')
 marinetemplate = os.path.join(evadir, 'marine_gdas_plots.yaml')
@@ -231,4 +189,8 @@
     print('running eva on', infile)
     subprocess.run(['eva', infile], check=True)
 
+#######################################
+# calculate diag statistics
+#######################################
+
 diag_statistics.get_diag_stats()