Added meridional section in the marine vrfy (#583)

* dded meridional section

* Added meridional section

---------

Co-authored-by: Guillaume Vernieres <[email protected]>

scripts/exgdas_global_marine_analysis_vrfy.py

@@ -49,11 +49,14 @@
 config = plotConfig(grid_file=grid_file,
                     data_file=data_file,
                     lats=np.arange(-60, 60, 10),
+                    lons=np.arange(-280, 80, 10),
                     variables_zonal={'Temp': [-0.5, 0.5],
                                      'Salt': [-0.1, 0.1]},
                     variables_horiz={'Temp': [-0.5, 0.5],
                                      'Salt': [-0.1, 0.1],
                                      'ave_ssh': [-0.1, 0.1]},
+                    variables_meridional={'Temp': [-0.5, 0.5],
+                                          'Salt': [-0.1, 0.1]},
                     colormap='RdBu',
                     comout=os.path.join(comout, 'vrfy', 'incr'))
 ocnIncPlotter = statePlotter(config)
@@ -146,11 +149,14 @@
 config = plotConfig(grid_file=grid_file,
                     data_file=data_file,
                     lats=np.arange(-60, 60, 10),
+                    lons=np.arange(-280, 80, 10),
                     variables_zonal={'Temp': [0, 2],
                                      'Salt': [0, 0.2]},
                     variables_horiz={'Temp': [0, 2],
                                      'Salt': [0, 0.2],
                                      'ave_ssh': [0, 0.1]},
+                    variables_meridional={'Temp': [0, 2],
+                                          'Salt': [0, 0.2]},
                     colormap='jet',
                     comout=os.path.join(comout, 'vrfy', 'bkgerr'))
 bkgErrPlotter = statePlotter(config)

ush/soca/soca_vrfy.py

@@ -29,8 +29,11 @@ def plotConfig(grid_file=[],
                comout=[],
                variables_horiz={},
                variables_zonal={},
+               variables_meridional={},
                lat=np.nan,
                lats=np.arange(-60, 60, 10),
+               lon=np.nan,
+               lons=np.arange(-280, 80, 10),
                proj='set me',
                projs=['Global']):
 
@@ -49,10 +52,13 @@ def plotConfig(grid_file=[],
     config['bounds'] = bounds
     config['lats'] = lats  # all the lats to plot
     config['lat'] = lat  # the lat being currently plotted
+    config['lons'] = lons
+    config['lon'] = lon
     config['max depths'] = max_depths  # all the max depths to plot
     config['max depth'] = max_depth  # the max depth currently plotted
     config['horiz variables'] = variables_horiz  # all the vars for horiz plots
     config['zonal variables'] = variables_zonal  # all the vars for zonal plots
+    config['meridional variables'] = variables_meridional  # all the vars for meridional plots
     config['variable'] = variable  # the variable currently plotted
     config['projs'] = projs  # all the projections etc.
     config['proj'] = proj
@@ -106,6 +112,7 @@ def plotHorizontalSlice(config):
         ax.set_extent([-180, 180, 50, 90], ccrs.PlateCarree())
     # ax.add_feature(cartopy.feature.LAND)  # TODO: make this work on hpc
     plt.savefig(figname, bbox_inches='tight', dpi=600)
+    plt.close(fig)
 
 
 def plotZonalSlice(config):
@@ -139,6 +146,35 @@ def plotZonalSlice(config):
     figname = os.path.join(dirname, variable +
                            'zonal_lat_'+str(int(lat)) + '_' + str(int(config['max depth'])) + 'm')
     plt.savefig(figname, bbox_inches='tight', dpi=600)
+    plt.close(fig)
+
+
+def plotMeridionalSlice(config):
+    """
+    pcolormesh of a Meridional slice of an ocean field
+    """
+    lon = float(config['lon'])
+    grid = xr.open_dataset(config['grid file'])
+    data = xr.open_dataset(config['fields file'])
+    lon_index = np.argmin(np.array(np.abs(np.squeeze(grid.lon)[:, 0]-lon)))
+    slice_data = np.squeeze(np.array(data[config['variable']]))[:, lon_index, :]
+    depth = np.squeeze(np.array(grid['h']))[:, lon_index, :]
+    depth[np.where(np.abs(depth) > 10000.0)] = 0.0
+    depth = np.cumsum(depth, axis=0)
+    bounds = config['bounds']
+    y = np.tile(np.squeeze(grid.lon[:, lon_index]), (np.shape(depth)[0], 1))
+    fig, ax = plt.subplots(figsize=(8, 5))
+    plt.pcolormesh(y, -depth, slice_data,
+                   vmin=bounds[0], vmax=bounds[1],
+                   cmap=config['colormap'])
+    plt.colorbar(label=config['variable']+' Lon '+str(lon), shrink=0.5, orientation='horizontal')
+    ax.set_ylim(-config['max depth'], 0)
+    dirname = os.path.join(config['comout'], config['variable'])
+    os.makedirs(dirname, exist_ok=True)
+    figname = os.path.join(dirname, config['variable'] +
+                           'meridional_lon_'+str(int(lon)) + '_' + str(int(config['max depth'])) + 'm')
+    plt.savefig(figname, bbox_inches='tight', dpi=600)
+    plt.close(fig)
 
 
 class statePlotter:
@@ -164,6 +200,21 @@ def plot(self):
                     self.config.update({'variable': variable, 'bounds': bounds})
                     plotZonalSlice(self.config)
 
+        #######################################
+        # Meridional slices
+
+        for lon in self.config['lons']:
+            self.config['lon'] = lon
+
+            for max_depth in self.config['max depths']:
+                self.config['max depth'] = max_depth
+
+                variableBounds = self.config['meridional variables']
+                for variable in variableBounds.keys():
+                    bounds = variableBounds[variable]
+                    self.config.update({'variable': variable, 'bounds': bounds})
+                    plotMeridionalSlice(self.config)
+
         #######################################
         # Horizontal slices
         for proj in self.config['projs']: