From 253ffcbeca3eb16b8340a3cc7fa14828dc3af523 Mon Sep 17 00:00:00 2001
From: sloosvel <saskia.loosveldt@bsc.es>
Date: Tue, 21 May 2019 14:04:29 +0200
Subject: [PATCH 01/30] Add diagonals functions

---
 earthdiagnostics/ocean/regionmean.py | 190 ++++++++++++++-------------
 1 file changed, 102 insertions(+), 88 deletions(-)

diff --git a/earthdiagnostics/ocean/regionmean.py b/earthdiagnostics/ocean/regionmean.py
index 75c283db..63367706 100644
--- a/earthdiagnostics/ocean/regionmean.py
+++ b/earthdiagnostics/ocean/regionmean.py
@@ -8,13 +8,19 @@ import iris.exceptions
 
 import numpy as np
 
+import netCDF4
+
 from earthdiagnostics.box import Box
 from earthdiagnostics.constants import Basins
 from earthdiagnostics.diagnostic import Diagnostic, DiagnosticOption, DiagnosticIntOption, DiagnosticDomainOption, \
-    DiagnosticBoolOption, DiagnosticBasinOption, DiagnosticVariableOption
+    DiagnosticBoolOption, DiagnosticBasinListOption, DiagnosticVariableOption
 from earthdiagnostics.modelingrealm import ModelingRealms
 from earthdiagnostics.utils import Utils, TempFile
 
+import diagonals.regmean as regmean
+from diagonals.mesh_helpers.nemo import Nemo
+
+
 
 class RegionMean(Diagnostic):
     """
@@ -45,7 +51,7 @@ class RegionMean(Diagnostic):
     "Diagnostic alias for the configuration file"
 
     def __init__(self, data_manager, startdate, member, chunk, domain, variable, box, save3d,
-                 variance, basin, grid_point):
+                 variance, basins, grid_point):
         Diagnostic.__init__(self, data_manager)
         self.startdate = startdate
         self.member = member
@@ -55,7 +61,7 @@ class RegionMean(Diagnostic):
         self.box = box
         self.save3d = save3d
         self.variance = variance
-        self.basin = basin
+        self.basins = basins
         self.grid_point = grid_point
 
         self.declared = {}
@@ -90,9 +96,13 @@ class RegionMean(Diagnostic):
         options_available = (DiagnosticDomainOption(),
                              DiagnosticVariableOption(diags.data_manager.config.var_manager),
                              DiagnosticOption('grid_point', 'T'),
-                             DiagnosticBasinOption('basin', Basins().Global),
+                             DiagnosticBasinListOption('basins', Basins().Global),
                              DiagnosticIntOption('min_depth', -1),
                              DiagnosticIntOption('max_depth', -1),
+                             DiagnosticIntOption('min_lat', None),
+                             DiagnosticIntOption('max_lat', None),
+                             DiagnosticIntOption('min_lon', None),
+                             DiagnosticIntOption('max_lon', None),
                              DiagnosticBoolOption('save3D', True),
                              DiagnosticBoolOption('variance', False),
                              DiagnosticOption('grid', ''))
@@ -101,12 +111,22 @@ class RegionMean(Diagnostic):
         box = Box()
         box.min_depth = options['min_depth']
         box.max_depth = options['max_depth']
+        box.min_lat = options['min_lat']
+        box.max_lat = options['max_lat']
+        box.min_lon = options['min_lon']
+        box.max_lon = options['max_lon']
+
+        basins = options['basins']
+        if not basins:
+            Log.error('Basins not recognized')
+            return()
+
 
         job_list = list()
         for startdate, member, chunk in diags.config.experiment.get_chunk_list():
             job = RegionMean(diags.data_manager, startdate, member, chunk,
                              options['domain'], options['variable'], box,
-                             options['save3D'], options['variance'], options['basin'],
+                             options['save3D'], options['variance'], options['basins'],
                              options['grid_point'].lower())
             job_list.append(job)
 
@@ -126,98 +146,48 @@ class RegionMean(Diagnostic):
         self._declare_var('mean', False, box_save)
         self._declare_var('mean', True, box_save)
 
-        if self.variance:
-            self._declare_var('var', False, box_save)
-            self._declare_var('var', True, box_save)
-
     def compute(self):
         """Run the diagnostic"""
         has_levels = self._fix_file_metadata()
         data = self._load_data()
+        masks = {}
+        self.basins.sort()
+        for basin in self.basins:
+            masks[basin] = Utils.get_mask(basin)
+
+        mesh = Nemo('mesh_hgr.nc', 'mask_regions.nc')
+        if self.box.min_lat is not None and self.box.max_lat is not None and self.box.min_lon is not None and self.box.max_lat is not None:
+            name = '{0}_{1}'.format(Box.get_lat_str(self.box), Box.get_lon_str(self.box))
+
+            masks[name] = mesh.get_region_mask(self.box.min_lat,
+                                               self.box.max_lat,
+                                               self.box.min_lon,
+                                               self.box.max_lon)
         if has_levels:
-            self._meand_3d_variable(data)
+            self._meand_3d_variable(data, mesh, masks)
         else:
-            self._mean_2d_var(data)
-
-    def _mean_2d_var(self, data):
-        mean = iris.cube.CubeList()
-        var = iris.cube.CubeList()
-        mask = np.squeeze(Utils.get_mask(self.basin, True))
-        if len(mask.shape) == 3:
-            mask = mask[0, ...]
-
-        e1 = self._try_load_cube(1)
-        e2 = self._try_load_cube(2)
-        weights = e1 * e2 * mask
-        for time_slice in data.slices_over('time'):
-            mean.append(time_slice.collapsed(['latitude', 'longitude'], iris.analysis.MEAN, weights=weights.data))
-            if self.variance:
-                var.append(time_slice.collapsed(['latitude', 'longitude'], iris.analysis.VARIANCE))
-        self._send_var('mean', False, mean)
-        if self.variance:
-            self._send_var('var', False, var)
-
-    def _meand_3d_variable(self, data):
-        mean = iris.cube.CubeList()
-        mean3d = iris.cube.CubeList()
-        var = iris.cube.CubeList()
-        var3d = iris.cube.CubeList()
-        mask = np.squeeze(Utils.get_mask(self.basin, True))
-        e1 = self._try_load_cube(1)
-        e2 = self._try_load_cube(2)
+            self._mean_2d_var(data, mesh, masks)
+
+    def _mean_2d_var(self, data, mesh, masks):
+        areacello = mesh.get_areacello()
+        mean = regmean.compute_regmean_2D(data.data, masks, areacello)
+        self._save_result_2D('mean', mean, data)
+
+
+    def _meand_3d_variable(self, data, mesh, masks):
+        areacello = mesh.get_areacello()
         e3 = self._try_load_cube(3)
         e3 = self._rename_depth(e3)
-        weights = e1 * e2
-        weights = e3 * weights.data * mask
         depth_constraint = iris.Constraint(depth=lambda c: self.box.min_depth <= c <= self.box.max_depth)
-        weights = weights.extract(depth_constraint).data
+        volcello = areacello*e3.extract(depth_constraint).data
         data = data.extract(depth_constraint)
+        mean = regmean.compute_regmean_3D(data.data, masks, volcello)
+        self._save_result_2D('mean', mean, data)
+        if self.save3d:
+            mean3d = regmean.compute_regmean_levels(data.data, masks, volcello)
+            self._save_result_3D('mean3d', mean3d, data)
+
 
-        for time_slice in data.slices_over('time'):
-            mean.append(time_slice.collapsed(
-                ['latitude', 'longitude', 'depth'],
-                iris.analysis.MEAN,
-                weights=weights
-            ))
-            if self.save3d:
-                mean3d.append(time_slice.collapsed(
-                    ['latitude', 'longitude'],
-                    iris.analysis.MEAN,
-                    weights=weights
-                ))
-            if self.variance:
-                var.append(time_slice.collapsed(
-                    ['latitude', 'longitude', 'depth'],
-                    iris.analysis.VARIANCE
-                ))
-                if self.save3d:
-                    var3d.append(time_slice.collapsed(
-                        ['latitude', 'longitude'],
-                        iris.analysis.VARIANCE
-                    ))
-        self._send_var('mean', True, mean3d)
-        self._send_var('mean', False, mean)
-        if self.variance:
-            self._send_var('var', True, var3d)
-
-            self._send_var('var', False, var)
-
-    def _try_load_cube(self, number):
-        try:
-            cube = iris.load_cube('mesh_hgr.nc', 'e{0}{1}'.format(number, self.grid_point))
-        except iris.exceptions.ConstraintMismatchError:
-            cube = iris.load_cube('mesh_hgr.nc', 'e{0}{1}_0'.format(number, self.grid_point))
-        cube = iris.util.squeeze(cube)
-        dims = len(cube.shape)
-        try:
-            cube.coord('i')
-        except iris.exceptions.CoordinateNotFoundError:
-            cube.add_dim_coord(iris.coords.DimCoord(np.arange(cube.shape[dims - 1]), var_name='i'), dims - 1)
-        try:
-            cube.coord('j')
-        except iris.exceptions.CoordinateNotFoundError:
-            cube.add_dim_coord(iris.coords.DimCoord(np.arange(cube.shape[dims - 2]), var_name='j'), dims - 2)
-        return cube
 
     def _load_data(self):
         coords = []
@@ -266,7 +236,7 @@ class RegionMean(Diagnostic):
             final_name = '{1}{0}'.format(var, self.variable)
 
         self.declared[final_name] = self.declare_chunk(ModelingRealms.ocean, final_name, self.startdate, self.member,
-                                                       self.chunk, box=box_save, region=self.basin)
+                                                       self.chunk, box=box_save, region=self.basins)
 
     def _send_var(self, var, threed, cube_list):
         if threed:
@@ -289,4 +259,48 @@ class RegionMean(Diagnostic):
             pass
         temp = TempFile.get()
         iris.save(cube, temp)
-        self.declared[final_name].set_local_file(temp, diagnostic=self, rename_var='result', region=self.basin)
+        self.declared[final_name].set_local_file(temp, diagnostic=self, rename_var='result', region=self.basins)
+
+
+    def _save_result_2D(self, var, result, data):
+        final_name = '{1}{0}'.format(var, self.variable)
+        temp = TempFile.get()
+        handler_source = Utils.open_cdf(self.variable_file.local_file)
+        handler_temp = Utils.open_cdf(temp, 'w')
+        Utils.copy_variable(handler_source, handler_temp, 'time', True, True)
+        handler_temp.createDimension('region', len(result))
+        handler_temp.createDimension('region_length', 50)
+        var_region = handler_temp.createVariable('region', 'S1',
+                                                 ('region', 'region_length'))
+        var = handler_temp.createVariable('{1}{0}'.format(var, self.variable),
+                                          float, ('time', 'region',),)
+        var.units = '{0}'.format(data.units)
+        for i, basin in enumerate(result):
+            var_region[i, ...] = netCDF4.stringtoarr(str(basin), 50)
+            var[..., i] = result[basin]
+        handler_temp.close()
+        self.declared[final_name].set_local_file(temp, diagnostic=self)
+
+
+
+    def _save_result_3D(self, var, result, data):
+        final_name = '{1}{0}'.format(var, self.variable)
+        temp = TempFile.get()
+        handler_source = Utils.open_cdf(self.variable_file.local_file)
+        handler_temp = Utils.open_cdf(temp, 'w')
+        Utils.copy_variable(handler_source, handler_temp, 'time', True, True)
+        Utils.copy_variable(handler_source, handler_temp, 'lev', True, True)
+        handler_temp.createDimension('region', len(result))
+        handler_temp.createDimension('region_length', 50)
+        var_region = handler_temp.createVariable('region', 'S1',
+                                                 ('region', 'region_length'))
+        var = handler_temp.createVariable('{1}{0}'.format(var, self.variable),
+                                          float, ('time', 'lev', 'region',),)
+        var.units = '{0}'.format(data.units)
+        for i, basin in enumerate(result):
+            var_region[i, ...] = netCDF4.stringtoarr(str(basin), 50)
+            var[..., i] = result[basin]
+        handler_temp.close()
+        self.declared[final_name].set_local_file(temp, diagnostic=self)
+
+
-- 
GitLab


From 34859c8e1f09a00797e82345b66806dfaccaafe4 Mon Sep 17 00:00:00 2001
From: sloosvel <saskia.loosveldt@bsc.es>
Date: Tue, 21 May 2019 16:44:54 +0200
Subject: [PATCH 02/30] Add functions for 3d vars

---
 earthdiagnostics/ocean/regionmean.py | 50 ++++++++++++++++++++--------
 1 file changed, 36 insertions(+), 14 deletions(-)

diff --git a/earthdiagnostics/ocean/regionmean.py b/earthdiagnostics/ocean/regionmean.py
index 63367706..8aafdf84 100644
--- a/earthdiagnostics/ocean/regionmean.py
+++ b/earthdiagnostics/ocean/regionmean.py
@@ -99,10 +99,10 @@ class RegionMean(Diagnostic):
                              DiagnosticBasinListOption('basins', Basins().Global),
                              DiagnosticIntOption('min_depth', -1),
                              DiagnosticIntOption('max_depth', -1),
-                             DiagnosticIntOption('min_lat', None),
-                             DiagnosticIntOption('max_lat', None),
-                             DiagnosticIntOption('min_lon', None),
-                             DiagnosticIntOption('max_lon', None),
+                             DiagnosticIntOption('min_lat', ''),
+                             DiagnosticIntOption('max_lat', ''),
+                             DiagnosticIntOption('min_lon', ''),
+                             DiagnosticIntOption('max_lon', ''),
                              DiagnosticBoolOption('save3D', True),
                              DiagnosticBoolOption('variance', False),
                              DiagnosticOption('grid', ''))
@@ -156,7 +156,7 @@ class RegionMean(Diagnostic):
             masks[basin] = Utils.get_mask(basin)
 
         mesh = Nemo('mesh_hgr.nc', 'mask_regions.nc')
-        if self.box.min_lat is not None and self.box.max_lat is not None and self.box.min_lon is not None and self.box.max_lat is not None:
+        if self.box.min_lat is not '' and self.box.max_lat is not '' and self.box.min_lon is not '' and self.box.max_lat is not '':
             name = '{0}_{1}'.format(Box.get_lat_str(self.box), Box.get_lon_str(self.box))
 
             masks[name] = mesh.get_region_mask(self.box.min_lat,
@@ -169,24 +169,41 @@ class RegionMean(Diagnostic):
             self._mean_2d_var(data, mesh, masks)
 
     def _mean_2d_var(self, data, mesh, masks):
-        areacello = mesh.get_areacello()
+        areacello = mesh.get_areacello(cell_point=self.grid_point)
         mean = regmean.compute_regmean_2D(data.data, masks, areacello)
         self._save_result_2D('mean', mean, data)
 
 
     def _meand_3d_variable(self, data, mesh, masks):
-        areacello = mesh.get_areacello()
+        areacello = mesh.get_areacello(cell_point=self.grid_point)
         e3 = self._try_load_cube(3)
         e3 = self._rename_depth(e3)
+        e3.coord('depth').bounds = data.coord('depth').bounds
         depth_constraint = iris.Constraint(depth=lambda c: self.box.min_depth <= c <= self.box.max_depth)
-        volcello = areacello*e3.extract(depth_constraint).data
+        volcello = areacello*e3.extract(depth_constraint).data.astype(np.float32)
         data = data.extract(depth_constraint)
         mean = regmean.compute_regmean_3D(data.data, masks, volcello)
         self._save_result_2D('mean', mean, data)
         if self.save3d:
             mean3d = regmean.compute_regmean_levels(data.data, masks, volcello)
-            self._save_result_3D('mean3d', mean3d, data)
+            self._save_result_3D('mean', mean3d, data)
 
+    def _try_load_cube(self, number):
+        try:
+            cube = iris.load_cube('mesh_hgr.nc', 'e{0}{1}'.format(number, self.grid_point))
+        except iris.exceptions.ConstraintMismatchError:
+            cube = iris.load_cube('mesh_hgr.nc', 'e{0}{1}_0'.format(number, self.grid_point))
+        cube = iris.util.squeeze(cube)
+        dims = len(cube.shape)
+        try:
+            cube.coord('i')
+        except iris.exceptions.CoordinateNotFoundError:
+            cube.add_dim_coord(iris.coords.DimCoord(np.arange(cube.shape[dims - 1]), var_name='i'), dims - 1)
+        try:
+            cube.coord('j')
+        except iris.exceptions.CoordinateNotFoundError:
+            cube.add_dim_coord(iris.coords.DimCoord(np.arange(cube.shape[dims - 2]), var_name='j'), dims - 2)
+        return cube
 
 
     def _load_data(self):
@@ -282,19 +299,25 @@ class RegionMean(Diagnostic):
         self.declared[final_name].set_local_file(temp, diagnostic=self)
 
 
-
     def _save_result_3D(self, var, result, data):
-        final_name = '{1}{0}'.format(var, self.variable)
+        final_name = '{1}3d{0}'.format(var, self.variable)
         temp = TempFile.get()
         handler_source = Utils.open_cdf(self.variable_file.local_file)
         handler_temp = Utils.open_cdf(temp, 'w')
         Utils.copy_variable(handler_source, handler_temp, 'time', True, True)
-        Utils.copy_variable(handler_source, handler_temp, 'lev', True, True)
         handler_temp.createDimension('region', len(result))
         handler_temp.createDimension('region_length', 50)
+        handler_temp.createDimension('lev', data.shape[1])
+        var_level = handler_temp.createVariable('lev', float, 'lev')
+        var_level[...] = data.coord('depth').points
+        var_level.units = 'm'
+        var_level.axis = 'Z'
+        var_level.positive = 'down'
+        var_level.long_name = 'ocean depth coordinate'
+        var_level.standard_name = 'depth'
         var_region = handler_temp.createVariable('region', 'S1',
                                                  ('region', 'region_length'))
-        var = handler_temp.createVariable('{1}{0}'.format(var, self.variable),
+        var = handler_temp.createVariable('{1}3d{0}'.format(var, self.variable),
                                           float, ('time', 'lev', 'region',),)
         var.units = '{0}'.format(data.units)
         for i, basin in enumerate(result):
@@ -303,4 +326,3 @@ class RegionMean(Diagnostic):
         handler_temp.close()
         self.declared[final_name].set_local_file(temp, diagnostic=self)
 
-
-- 
GitLab


From f13246e66c66e4bebeec292744444e57a06b40e8 Mon Sep 17 00:00:00 2001
From: sloosvel <saskia.loosveldt@bsc.es>
Date: Tue, 21 May 2019 18:10:07 +0200
Subject: [PATCH 03/30] Change default value for box limits

---
 earthdiagnostics/ocean/regionmean.py | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/earthdiagnostics/ocean/regionmean.py b/earthdiagnostics/ocean/regionmean.py
index 8aafdf84..ed24e759 100644
--- a/earthdiagnostics/ocean/regionmean.py
+++ b/earthdiagnostics/ocean/regionmean.py
@@ -99,10 +99,10 @@ class RegionMean(Diagnostic):
                              DiagnosticBasinListOption('basins', Basins().Global),
                              DiagnosticIntOption('min_depth', -1),
                              DiagnosticIntOption('max_depth', -1),
-                             DiagnosticIntOption('min_lat', ''),
-                             DiagnosticIntOption('max_lat', ''),
-                             DiagnosticIntOption('min_lon', ''),
-                             DiagnosticIntOption('max_lon', ''),
+                             DiagnosticIntOption('min_lat', -1),
+                             DiagnosticIntOption('max_lat', -1),
+                             DiagnosticIntOption('min_lon', -1),
+                             DiagnosticIntOption('max_lon', -1),
                              DiagnosticBoolOption('save3D', True),
                              DiagnosticBoolOption('variance', False),
                              DiagnosticOption('grid', ''))
@@ -156,7 +156,7 @@ class RegionMean(Diagnostic):
             masks[basin] = Utils.get_mask(basin)
 
         mesh = Nemo('mesh_hgr.nc', 'mask_regions.nc')
-        if self.box.min_lat is not '' and self.box.max_lat is not '' and self.box.min_lon is not '' and self.box.max_lat is not '':
+        if self.box.min_lat is not -1 and self.box.max_lat is not -1 and self.box.min_lon is not -1 and self.box.max_lat is not -1:
             name = '{0}_{1}'.format(Box.get_lat_str(self.box), Box.get_lon_str(self.box))
 
             masks[name] = mesh.get_region_mask(self.box.min_lat,
-- 
GitLab


From 9fb87985827fe82cf0c304b634302595148ced04 Mon Sep 17 00:00:00 2001
From: sloosvel <saskia.loosveldt@bsc.es>
Date: Tue, 21 May 2019 18:10:36 +0200
Subject: [PATCH 04/30] Add diagonals functions, test 3d pending

---
 earthdiagnostics/ocean/regionsum.py | 252 +++++++++++++++++++++++-----
 1 file changed, 210 insertions(+), 42 deletions(-)

diff --git a/earthdiagnostics/ocean/regionsum.py b/earthdiagnostics/ocean/regionsum.py
index 9e2c9b95..2f431ff0 100644
--- a/earthdiagnostics/ocean/regionsum.py
+++ b/earthdiagnostics/ocean/regionsum.py
@@ -2,21 +2,35 @@
 """Diagnostic to calculate a region total"""
 import os
 
+import iris
+import iris.util
+import iris.coords
+import iris.analysis
+import iris.exceptions
+
+import numpy as np
+
+import netCDF4
+
 from earthdiagnostics import cdftools
 from earthdiagnostics.box import Box
 from earthdiagnostics.constants import Basins
-from earthdiagnostics.diagnostic import Diagnostic, DiagnosticOption, DiagnosticIntOption, DiagnosticDomainOption, \
-    DiagnosticBoolOption, DiagnosticBasinOption, DiagnosticVariableOption
+from earthdiagnostics.diagnostic import Diagnostic, DiagnosticOption, \
+    DiagnosticIntOption, DiagnosticDomainOption, \
+    DiagnosticBoolOption, DiagnosticBasinListOption, DiagnosticVariableOption
 from earthdiagnostics.modelingrealm import ModelingRealms
 from earthdiagnostics.utils import Utils, TempFile
 
+import diagonals.regsum as regsum
+from diagonals.mesh_helpers.nemo import Nemo
 
 class RegionSum(Diagnostic):
     """
     Computes the sum of the field (3D, weighted).
 
-    For 3D fields, a horizontal mean for each level is also given. If a spatial window
-    is specified, the mean value is computed only in this window.
+    For 3D fields, a horizontal mean for each level is also given.
+    If a spatial window is specified, the mean value is computed only
+    in this window.
 
     :original author: Javier Vegas-Regidor <javier.vegas@bsc.es>
 
@@ -39,7 +53,8 @@ class RegionSum(Diagnostic):
     alias = 'regsum'
     "Diagnostic alias for the configuration file"
 
-    def __init__(self, data_manager, startdate, member, chunk, domain, variable, grid_point, box, save3d, basin,
+    def __init__(self, data_manager, startdate, member, chunk, domain,
+                 variable, grid_point, box, save3d, basins,
                  grid):
         Diagnostic.__init__(self, data_manager)
         self.startdate = startdate
@@ -50,7 +65,7 @@ class RegionSum(Diagnostic):
         self.grid_point = grid_point.upper()
         self.box = box
         self.save3d = save3d
-        self.basin = basin
+        self.basins = basins
         self.grid = grid
         self.declared = {}
 
@@ -61,14 +76,22 @@ class RegionSum(Diagnostic):
         if self._different_type(other):
             return False
 
-        return self.startdate == other.startdate and self.member == other.member and self.chunk == other.chunk and \
+        return self.startdate == other.startdate and \
+            self.member == other.member and self.chunk == other.chunk and \
             self.box == other.box and self.variable == other.variable and \
-            self.grid_point == other.grid_point and self.grid == other.grid and self.basin == other.basin
+            self.grid_point == other.grid_point and self.grid == other.grid \
+            and self.basin == other.basin
 
     def __str__(self):
-        return 'Region sum Startdate: {0.startdate} Member: {0.member} Chunk: {0.chunk} Variable: {0.variable} ' \
+        return 'Region sum Startdate: {0.startdate} Member: {0.member}' \
+               'Chunk: {0.chunk} Variable: {0.variable} ' \
                'Grid point: {0.grid_point} Box: {0.box} Save 3D: {0.save3d}' \
-               'Original grid: {0.grid} Basin: {0.basin}'.format(self)
+               'Original grid: {0.grid} Basin: {0.basins}'.format(self)
+
+
+    def __hash__(self):
+        return hash(str(self))
+
 
     @classmethod
     def generate_jobs(cls, diags, options):
@@ -84,9 +107,13 @@ class RegionSum(Diagnostic):
         options_available = (DiagnosticDomainOption(),
                              DiagnosticVariableOption(diags.data_manager.config.var_manager),
                              DiagnosticOption('grid_point', 'T'),
-                             DiagnosticBasinOption('basin', Basins().Global),
-                             DiagnosticIntOption('min_depth', 0),
-                             DiagnosticIntOption('max_depth', 0),
+                             DiagnosticBasinListOption('basins', Basins().Global),
+                             DiagnosticIntOption('min_depth', -1),
+                             DiagnosticIntOption('max_depth', -1),
+                             DiagnosticIntOption('min_lat', -1),
+                             DiagnosticIntOption('max_lat', -1),
+                             DiagnosticIntOption('min_lon', -1),
+                             DiagnosticIntOption('max_lon', -1),
                              DiagnosticBoolOption('save3D', True),
                              DiagnosticOption('grid', ''))
         options = cls.process_options(options, options_available)
@@ -94,23 +121,33 @@ class RegionSum(Diagnostic):
         box = Box()
         box.min_depth = options['min_depth']
         box.max_depth = options['max_depth']
+        box.min_lat = options['min_lat']
+        box.max_lat = options['max_lat']
+        box.min_lon = options['min_lon']
+        box.max_lon = options['max_lon']
+
+        basins = options['basins']
+        if not basins:
+            Log.error('Basins not recognized')
+            return()
 
         job_list = list()
         for startdate, member, chunk in diags.config.experiment.get_chunk_list():
             job_list.append(RegionSum(diags.data_manager, startdate, member, chunk,
                                       options['domain'], options['variable'], options['grid_point'], box,
-                                      options['save3D'], options['basin'], options['grid']))
+                                      options['save3D'], options['basins'], options['grid']))
         return job_list
 
     def request_data(self):
         """Request data required by the diagnostic"""
-        self.variable_file = self.request_chunk(self.domain, self.variable, self.startdate, self.member, self.chunk,
+        self.variable_file = self.request_chunk(self.domain, self.variable,
+                                                self.startdate,
+                                                self.member, self.chunk,
                                                 grid=self.grid)
 
     def declare_data_generated(self):
         """Declare data to be generated by the diagnostic"""
-        if self.box.min_depth == 0:
-            # To cdftools, this means all levels
+        if self.box.min_depth == -1 and self.box.max_depth == -1:
             box_save = None
         else:
             box_save = self.box
@@ -120,34 +157,117 @@ class RegionSum(Diagnostic):
 
     def compute(self):
         """Run the diagnostic"""
-        mean_file = TempFile.get()
-
-        variable_file = self.variable_file.local_file
-
-        handler = Utils.open_cdf(variable_file)
-        self.save3d &= 'lev' in handler.dimensions
-        if "latitude" in handler.variables:
-            self.lat_name = 'latitude'
-        if "longitude" in handler.variables:
-            self.lon_name = 'longitude'
-
+        has_levels = self._fix_file_metadata()
+        data = self._load_data()
+        masks = {}
+        self.basins.sort()
+        for basin in self.basins:
+            masks[basin] = Utils.get_mask(basin)
+
+        mesh = Nemo('mesh_hgr.nc', 'mask_regions.nc')
+        if self.box.min_lat is not '' and self.box.max_lat is not '' and self.box.min_lon is not '' and self.box.max_lat is not '':
+            name = '{0}_{1}'.format(Box.get_lat_str(self.box), Box.get_lon_str(self.box))
+
+            masks[name] = mesh.get_region_mask(self.box.min_lat,
+                                               self.box.max_lat,
+                                               self.box.min_lon,
+                                               self.box.max_lon)
+        if has_levels:
+            self._sum_3d_variable(data, mesh, masks)
+        else:
+            self._sum_2d_var(data, mesh, masks)
+
+
+    def _sum_2d_var(self, data, mesh, masks):
+        areacello = mesh.get_areacello(cell_point=self.grid_point)
+        varsum = regsum.compute_regsum_2D(data.data, masks, areacello)
+        self._save_result_2D('sum', varsum, data)
+
+
+    def _sum_3d_variable(self, data, mesh, masks):
+        areacello = mesh.get_areacello(cell_point=self.grid_point)
+        tmask = mesh.get_landsea_mask(cell_point=self.grid_point)
+        e3 = self._try_load_cube(3)
+        e3 = self._rename_depth(e3)
+        e3.coord('depth').bounds = data.coord('depth').bounds
+        depth_constraint = iris.Constraint(depth=lambda c: self.box.min_depth <= c <= self.box.max_depth)
+        volcello = areacello*e3.extract(depth_constraint).data.astype(np.float32)
+        data = data.extract(depth_constraint)
+        varsum = regsum.compute_regsum_3D(data.data, masks, volcello, tmask)
+        self._save_result_2D('sum', varsum, data)
+        if self.save3d:
+            varsum3d = regsum.compute_regsum_levels(data.data, masks,
+                                                    volcello, tmask)
+            self._save_result_3D('aum', varsum3d, data)
+
+
+    def _try_load_cube(self, number):
+        try:
+            cube = iris.load_cube('mesh_hgr.nc', 'e{0}{1}'.format(number, self.grid_point))
+        except iris.exceptions.ConstraintMismatchError:
+            cube = iris.load_cube('mesh_hgr.nc', 'e{0}{1}_0'.format(number, self.grid_point))
+        cube = iris.util.squeeze(cube)
+        dims = len(cube.shape)
+        try:
+            cube.coord('i')
+        except iris.exceptions.CoordinateNotFoundError:
+            cube.add_dim_coord(iris.coords.DimCoord(np.arange(cube.shape[dims - 1]), var_name='i'), dims - 1)
+        try:
+            cube.coord('j')
+        except iris.exceptions.CoordinateNotFoundError:
+            cube.add_dim_coord(iris.coords.DimCoord(np.arange(cube.shape[dims - 2]), var_name='j'), dims - 2)
+        return cube
+
+    def _load_data(self):
+        coords = []
+        handler = Utils.open_cdf(self.variable_file.local_file)
+        for variable in handler.variables:
+            if variable in ('time', 'lev', 'lat', 'lon', 'latitude', 'longitude', 'leadtime', 'time_centered'):
+                coords.append(variable)
+            if variable == 'time_centered':
+                handler.variables[variable].standard_name = ''
+
+        handler.variables[self.variable].coordinates = ' '.join(coords)
         handler.close()
 
-        cdfmean_options = ['-v', self.variable, '-p', self.grid_point,
-                           '-zoom', 0, 0, 0, 0, self.box.min_depth, self.box.max_depth]
-        if self.basin != Basins().Global:
-            cdfmean_options.append('-M')
-            cdfmean_options.append('mask_regions.3d.nc')
-            cdfmean_options.append(self.basin.name)
-
-        cdftools.run('cdfsum', input_file=variable_file, input_option='-f', output_file=mean_file,
-                     options=cdfmean_options)
-        Utils.rename_variables(mean_file, {'gdept': 'lev', 'gdepw': 'lev'}, must_exist=False)
+        data = iris.load_cube(self.variable_file.local_file)
+        return self._rename_depth(data)
+
+
+    def _rename_depth(self, data):
+        for coord_name in ('model_level_number', 'Vertical T levels', 'lev'):
+            if data.coords(coord_name):
+                coord = data.coord(coord_name)
+                coord.standard_name = 'depth'
+                coord.long_name = 'depth'
+                break
+        return data
+
+
+    def _fix_file_metadata(self):
+        handler = Utils.open_cdf(self.variable_file.local_file)
+        var = handler.variables[self.variable]
+        coordinates = ''
+        has_levels = False
+        for dimension in handler.variables.keys():
+            if dimension in ['time', 'lev', 'lat', 'latitude', 'lon', 'longitude', 'i', 'j']:
+                coordinates += ' {0}'.format(dimension)
+            if dimension == 'lev':
+                has_levels = True
+        var.coordinates = coordinates
+        handler.close()
+        return has_levels
 
-        self._send_var(False, mean_file)
-        self._send_var(True, mean_file)
+    def _declare_var(self, var, threed, box_save):
+        if threed:
+            if not self.save3d:
+                return False
+            final_name = '{1}3d{0}'.format(var, self.variable)
+        else:
+            final_name = '{1}{0}'.format(var, self.variable)
 
-        os.remove(mean_file)
+        self.declared[final_name] = self.declare_chunk(ModelingRealms.ocean, final_name, self.startdate, self.member,
+                                                       self.chunk, box=box_save, region=self.basins)
 
     def _send_var(self, threed, mean_file):
         var = 'sum'
@@ -186,4 +306,52 @@ class RegionSum(Diagnostic):
             final_name = '{1}{0}'.format(var, self.variable)
 
         self.declared[final_name] = self.declare_chunk(ModelingRealms.ocean, final_name, self.startdate, self.member,
-                                                       self.chunk, box=box_save, region=self.basin, grid=self.grid)
+                                                       self.chunk, box=box_save, region=self.basins, grid=self.grid)
+
+
+    def _save_result_2D(self, var, result, data):
+        final_name = '{1}{0}'.format(var, self.variable)
+        temp = TempFile.get()
+        handler_source = Utils.open_cdf(self.variable_file.local_file)
+        handler_temp = Utils.open_cdf(temp, 'w')
+        Utils.copy_variable(handler_source, handler_temp, 'time', True, True)
+        handler_temp.createDimension('region', len(result))
+        handler_temp.createDimension('region_length', 50)
+        var_region = handler_temp.createVariable('region', 'S1',
+                                                 ('region', 'region_length'))
+        var = handler_temp.createVariable('{1}{0}'.format(var, self.variable),
+                                          float, ('time', 'region',),)
+        var.units = '{0}'.format(data.units)
+        for i, basin in enumerate(result):
+            var_region[i, ...] = netCDF4.stringtoarr(str(basin), 50)
+            var[..., i] = result[basin]
+        handler_temp.close()
+        self.declared[final_name].set_local_file(temp, diagnostic=self)
+
+
+    def _save_result_3D(self, var, result, data):
+        final_name = '{1}3d{0}'.format(var, self.variable)
+        temp = TempFile.get()
+        handler_source = Utils.open_cdf(self.variable_file.local_file)
+        handler_temp = Utils.open_cdf(temp, 'w')
+        Utils.copy_variable(handler_source, handler_temp, 'time', True, True)
+        handler_temp.createDimension('region', len(result))
+        handler_temp.createDimension('region_length', 50)
+        handler_temp.createDimension('lev', data.shape[1])
+        var_level = handler_temp.createVariable('lev', float, 'lev')
+        var_level[...] = data.coord('depth').points
+        var_level.units = 'm'
+        var_level.axis = 'Z'
+        var_level.positive = 'down'
+        var_level.long_name = 'ocean depth coordinate'
+        var_level.standard_name = 'depth'
+        var_region = handler_temp.createVariable('region', 'S1',
+                                                 ('region', 'region_length'))
+        var = handler_temp.createVariable('{1}3d{0}'.format(var, self.variable),
+                                          float, ('time', 'lev', 'region',),)
+        var.units = '{0}'.format(data.units)
+        for i, basin in enumerate(result):
+            var_region[i, ...] = netCDF4.stringtoarr(str(basin), 50)
+            var[..., i] = result[basin]
+        handler_temp.close()
+        self.declared[final_name].set_local_file(temp, diagnostic=self)
-- 
GitLab


From 15a6d3fef2d2dca2aa1c0757149c2b9827ca51c6 Mon Sep 17 00:00:00 2001
From: sloosvel <saskia.loosveldt@bsc.es>
Date: Wed, 22 May 2019 09:49:47 +0200
Subject: [PATCH 05/30] Change default value for region box

---
 earthdiagnostics/ocean/regionsum.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/earthdiagnostics/ocean/regionsum.py b/earthdiagnostics/ocean/regionsum.py
index 2f431ff0..ab2de26e 100644
--- a/earthdiagnostics/ocean/regionsum.py
+++ b/earthdiagnostics/ocean/regionsum.py
@@ -165,7 +165,7 @@ class RegionSum(Diagnostic):
             masks[basin] = Utils.get_mask(basin)
 
         mesh = Nemo('mesh_hgr.nc', 'mask_regions.nc')
-        if self.box.min_lat is not '' and self.box.max_lat is not '' and self.box.min_lon is not '' and self.box.max_lat is not '':
+        if self.box.min_lat is not -1 and self.box.max_lat is not -1 and self.box.min_lon is not -1 and self.box.max_lat is not -1:
             name = '{0}_{1}'.format(Box.get_lat_str(self.box), Box.get_lon_str(self.box))
 
             masks[name] = mesh.get_region_mask(self.box.min_lat,
-- 
GitLab


From 0ef83cd087d47286f4982a5ebf048d5bd579e923 Mon Sep 17 00:00:00 2001
From: sloosvel <saskia.loosveldt@bsc.es>
Date: Wed, 22 May 2019 11:15:34 +0200
Subject: [PATCH 06/30] Fix depth constrain bug

---
 earthdiagnostics/ocean/regionmean.py |  8 +++++---
 earthdiagnostics/ocean/regionsum.py  | 18 +++++++++++-------
 2 files changed, 16 insertions(+), 10 deletions(-)

diff --git a/earthdiagnostics/ocean/regionmean.py b/earthdiagnostics/ocean/regionmean.py
index ed24e759..b6ced885 100644
--- a/earthdiagnostics/ocean/regionmean.py
+++ b/earthdiagnostics/ocean/regionmean.py
@@ -179,9 +179,11 @@ class RegionMean(Diagnostic):
         e3 = self._try_load_cube(3)
         e3 = self._rename_depth(e3)
         e3.coord('depth').bounds = data.coord('depth').bounds
-        depth_constraint = iris.Constraint(depth=lambda c: self.box.min_depth <= c <= self.box.max_depth)
-        volcello = areacello*e3.extract(depth_constraint).data.astype(np.float32)
-        data = data.extract(depth_constraint)
+        if self.box.min_depth is not -1 and self.box.max_depth is not -1:
+            depth_constraint = iris.Constraint(depth=lambda c: self.box.min_depth <= c <= self.box.max_depth)
+            e3 = e3.extract(depth_constraint)
+            data = data.extract(depth_constraint)
+        volcello = areacello*e3.data.astype(np.float32)
         mean = regmean.compute_regmean_3D(data.data, masks, volcello)
         self._save_result_2D('mean', mean, data)
         if self.save3d:
diff --git a/earthdiagnostics/ocean/regionsum.py b/earthdiagnostics/ocean/regionsum.py
index ab2de26e..bc58650a 100644
--- a/earthdiagnostics/ocean/regionsum.py
+++ b/earthdiagnostics/ocean/regionsum.py
@@ -62,7 +62,7 @@ class RegionSum(Diagnostic):
         self.chunk = chunk
         self.domain = domain
         self.variable = variable
-        self.grid_point = grid_point.upper()
+        self.grid_point = grid_point
         self.box = box
         self.save3d = save3d
         self.basins = basins
@@ -134,8 +134,10 @@ class RegionSum(Diagnostic):
         job_list = list()
         for startdate, member, chunk in diags.config.experiment.get_chunk_list():
             job_list.append(RegionSum(diags.data_manager, startdate, member, chunk,
-                                      options['domain'], options['variable'], options['grid_point'], box,
-                                      options['save3D'], options['basins'], options['grid']))
+                                      options['domain'], options['variable'],
+                                      options['grid_point'].lower(), box,
+                                      options['save3D'], options['basins'],
+                                      options['grid']))
         return job_list
 
     def request_data(self):
@@ -190,15 +192,17 @@ class RegionSum(Diagnostic):
         e3 = self._try_load_cube(3)
         e3 = self._rename_depth(e3)
         e3.coord('depth').bounds = data.coord('depth').bounds
-        depth_constraint = iris.Constraint(depth=lambda c: self.box.min_depth <= c <= self.box.max_depth)
-        volcello = areacello*e3.extract(depth_constraint).data.astype(np.float32)
-        data = data.extract(depth_constraint)
+        if self.box.min_depth is not -1 and self.box.max_depth is not -1:
+            depth_constraint = iris.Constraint(depth=lambda c: self.box.min_depth <= c <= self.box.max_depth)
+            e3 = e3.extract(depth_constraint)
+            data = data.extract(depth_constraint)
+        volcello = areacello*e3.data.astype(np.float32)
         varsum = regsum.compute_regsum_3D(data.data, masks, volcello, tmask)
         self._save_result_2D('sum', varsum, data)
         if self.save3d:
             varsum3d = regsum.compute_regsum_levels(data.data, masks,
                                                     volcello, tmask)
-            self._save_result_3D('aum', varsum3d, data)
+            self._save_result_3D('sum', varsum3d, data)
 
 
     def _try_load_cube(self, number):
-- 
GitLab


From 214d7bf0d51bc74a165495df40302cda6e0c521d Mon Sep 17 00:00:00 2001
From: sloosvel <saskia.loosveldt@bsc.es>
Date: Wed, 22 May 2019 11:47:17 +0200
Subject: [PATCH 07/30] Add constrain for tmask

---
 earthdiagnostics/ocean/regionsum.py | 8 +++++---
 1 file changed, 5 insertions(+), 3 deletions(-)

diff --git a/earthdiagnostics/ocean/regionsum.py b/earthdiagnostics/ocean/regionsum.py
index bc58650a..f9207e59 100644
--- a/earthdiagnostics/ocean/regionsum.py
+++ b/earthdiagnostics/ocean/regionsum.py
@@ -188,7 +188,7 @@ class RegionSum(Diagnostic):
 
     def _sum_3d_variable(self, data, mesh, masks):
         areacello = mesh.get_areacello(cell_point=self.grid_point)
-        tmask = mesh.get_landsea_mask(cell_point=self.grid_point)
+        tmask = iris.load_cube('mesh_hgr.nc', '{0}mask'.format(self.grid_point))
         e3 = self._try_load_cube(3)
         e3 = self._rename_depth(e3)
         e3.coord('depth').bounds = data.coord('depth').bounds
@@ -196,8 +196,10 @@ class RegionSum(Diagnostic):
             depth_constraint = iris.Constraint(depth=lambda c: self.box.min_depth <= c <= self.box.max_depth)
             e3 = e3.extract(depth_constraint)
             data = data.extract(depth_constraint)
+            tmask = tmask.extract(depth_constraint)
         volcello = areacello*e3.data.astype(np.float32)
-        varsum = regsum.compute_regsum_3D(data.data, masks, volcello, tmask)
+        varsum = regsum.compute_regsum_3D(data.data, masks, volcello,
+                                          tmask.data)
         self._save_result_2D('sum', varsum, data)
         if self.save3d:
             varsum3d = regsum.compute_regsum_levels(data.data, masks,
@@ -299,7 +301,7 @@ class RegionSum(Diagnostic):
         if hasattr(var_handler, 'valid_max'):
             del var_handler.valid_max
         handler.close()
-        self.declared[final_name].set_local_file(temp2, diagnostic=self, rename_var=original_name, region=self.basin)
+        self.declared[final_name].set_local_file(temp2, diagnostic=self, rename_var=original_name, region=self.basins)
 
     def _declare_var(self, var, threed, box_save):
         if threed:
-- 
GitLab


From 469973dd43740d5fd6cbfa89fe27596949242ab3 Mon Sep 17 00:00:00 2001
From: sloosvel <saskia.loosveldt@bsc.es>
Date: Thu, 23 May 2019 11:17:26 +0200
Subject: [PATCH 08/30] Call diagonals functions, issue with missing values

---
 earthdiagnostics/ocean/moc.py | 135 ++++++++++++++++++----------------
 1 file changed, 72 insertions(+), 63 deletions(-)

diff --git a/earthdiagnostics/ocean/moc.py b/earthdiagnostics/ocean/moc.py
index fefc8e14..44e1301d 100644
--- a/earthdiagnostics/ocean/moc.py
+++ b/earthdiagnostics/ocean/moc.py
@@ -9,11 +9,17 @@ from iris.coords import DimCoord, AuxCoord
 from iris.cube import CubeList
 import iris.analysis
 
+
+import netCDF4
+
 from earthdiagnostics.constants import Basins
 from earthdiagnostics.diagnostic import Diagnostic, DiagnosticBasinListOption
 from earthdiagnostics.modelingrealm import ModelingRealms
 from earthdiagnostics.utils import Utils, TempFile
 
+import diagonals.moc as moc
+from diagonals.mesh_helpers.nemo import Nemo
+
 
 class Moc(Diagnostic):
     """
@@ -47,9 +53,11 @@ class Moc(Diagnostic):
         self.chunk = chunk
         self.basins = basins
 
+        self.results = {}
+
     def __str__(self):
         basins = []
-        basins.extend(self.basins.keys())
+        basins.extend(self.basins)
         return 'MOC Startdate: {0.startdate} Member: {0.member} ' \
             'Chunk: {0.chunk} Basins: {1}'.format(self, basins)
 
@@ -86,31 +94,9 @@ class Moc(Diagnostic):
             Log.error('Basins not recognized')
             return ()
 
-        try:
-            e1v = iris.load_cube('mesh_hgr.nc', 'e1v')
-        except iris.exceptions.ConstraintMismatchError:
-            e1v = iris.load_cube('mesh_hgr.nc', 'e1v_0')
-        try:
-            e3v = iris.load_cube('mesh_hgr.nc', 'e3t_0')
-        except iris.exceptions.ConstraintMismatchError:
-            e3v = iris.load_cube('mesh_hgr.nc', 'e3t_0')
-        e1v = iris.util.squeeze(e1v).data
-        e3v = iris.util.squeeze(e3v).data
-        if len(e3v.shape) == 1:
-            e3v = np.expand_dims(e3v.data, 1)
-            e3v = np.expand_dims(e3v, 2)
-        else:
-            e3v = e3v.data
-        mesh = - e1v * e3v
-
-        masks = {}
-        basins.sort()
-        for basin in basins:
-            masks[basin.name] = Utils.get_mask(basin) * mesh / 1e6
-
         job_list = list()
         for startdate, member, chunk in diags.config.experiment.get_chunk_list():
-            job_list.append(Moc(diags.data_manager, startdate, member, chunk, masks))
+            job_list.append(Moc(diags.data_manager, startdate, member, chunk, basins))
         return job_list
 
     def request_data(self):
@@ -123,43 +109,66 @@ class Moc(Diagnostic):
 
     def compute(self):
         """Run the diagnostic"""
-        data = iris.load_cube(self.variable_file.local_file)
-        Log.debug(str(data))
-        try:
-            data.coord('i')
-        except iris.exceptions.CoordinateNotFoundError:
-            dims = len(data.shape)
-            data.add_dim_coord(iris.coords.DimCoord(np.arange(data.shape[dims - 1]), var_name='i'), dims - 1)
-        try:
-            data.coord('j')
-        except iris.exceptions.CoordinateNotFoundError:
-            dims = len(data.shape)
-            data.add_dim_coord(iris.coords.DimCoord(np.arange(data.shape[dims - 2]), var_name='j'), dims - 2)
-
-        for coord_name in ('model_level_number', 'Vertical V levels', 'lev'):
-            if data.coords(coord_name):
-                coord = data.coord(coord_name)
-                coord.standard_name = 'depth'
-                coord.long_name = 'depth'
-                break
-
-        moc_results = CubeList()
-        for map_slice in data.slices_over('time'):
-            # Force data loading
-            map_slice.data
-            Log.debug(str(map_slice))
-            for basin, mask in six.iteritems(self.basins):
-                moc = map_slice.collapsed(('i', 'depth'), iris.analysis.SUM, weights=mask)
-                moc.add_aux_coord(
-                    AuxCoord([basin], var_name='region', standard_name='region')
-                )
-                moc_results.append(moc)
-        results = moc_results.merge_cube()
+        vo = iris.load_cube(self.variable_file.local_file)
+        mesh = Nemo('mesh_hgr.nc', 'mask_regions.nc')
+        e1v = mesh.get_i_length(cell_point='V')
+        e3v = mesh.get_k_length(cell_point='V')
+
+        masks = {}
+        self.basins.sort()
+        for basin in self.basins:
+            if basin is 'Global':
+                global_mask = mesh.get_landsea_mask(cell_point='V')
+                global_mask[..., 0] = 0.0
+                global_mask[..., -1] = 0.0
+                masks[basin] = global_mask
+            else:
+                masks[basin] = Utils.get_mask(basin)
+
+        moc_results = moc.compute(masks, e1v, e3v, vo.data)
+        del vo, e1v, e3v
+        self._save_result(moc_results, mesh)
+
+
+    def _save_result(self, result, mesh):
         temp = TempFile.get()
-        results.var_name = Moc.vsftmyz
-        results.remove_coord('i')
-        results.remove_coord('depth')
-        results.remove_coord('longitude')
-        results.units = 'Sverdrup'
-        iris.save(results, temp)
-        self.results.set_local_file(temp)
+        handler_source = Utils.open_cdf(self.variable_file.local_file)
+        handler_temp = Utils.open_cdf(temp, 'w')
+        gphiv = np.squeeze(mesh.get_grid_latitude(cell_point='V'))
+        max_gphiv = np.unravel_index(np.argmax(gphiv), gphiv.shape)[1]
+
+        Utils.copy_variable(handler_source, handler_temp, 'time', True, True)
+        Utils.copy_variable(handler_source, handler_temp, 'lev', True, True)
+        handler_temp.createDimension('i', 1)
+        handler_temp.createDimension('j', gphiv.shape[0])
+        handler_temp.createDimension('region', len(result))
+        handler_temp.createDimension('region_length', 50)
+
+
+        var_region = handler_temp.createVariable('region', 'S1',
+                                                 ('region', 'region_length'))
+
+        lat = handler_temp.createVariable('lat', float, ('j', 'i'))
+        lat[...] = gphiv[:, max_gphiv]
+        lat.units = 'degrees_north'
+        lat.long_name = "Latitude"
+
+        lon = handler_temp.createVariable('lon', float, ('j', 'i'))
+        lon[...] = 0
+        lon.units = 'degrees_east'
+        lon.long_name = "Longitude"
+
+        var = handler_temp.createVariable('vsftmyz', float, ('time', 'lev',
+                                                             'i', 'j',
+                                                             'region'))
+        var.units = 'Sverdrup'
+        var.coordinates = 'lev time'
+        var.long_name = 'Ocean meridional overturning volume streamfunction'
+        var.missing_value = 1e20
+        var.fill_value = 1e20
+
+        for i, basin in enumerate(result):
+            var_region[i, ...] = netCDF4.stringtoarr(str(basin), 50)
+            var[..., i] = result[basin]
+        handler_temp.close()
+        self.results.set_local_file(temp, diagnostic=self)
-- 
GitLab


From e5a5039920f3da496993388e5002b4b0dbb8ba3d Mon Sep 17 00:00:00 2001
From: sloosvel <saskia.loosveldt@bsc.es>
Date: Fri, 24 May 2019 10:01:52 +0200
Subject: [PATCH 09/30] Add diagonals functions

---
 earthdiagnostics/ocean/siasiesiv.py | 77 ++++++++++-------------------
 1 file changed, 27 insertions(+), 50 deletions(-)

diff --git a/earthdiagnostics/ocean/siasiesiv.py b/earthdiagnostics/ocean/siasiesiv.py
index 763abbaf..c1102b02 100644
--- a/earthdiagnostics/ocean/siasiesiv.py
+++ b/earthdiagnostics/ocean/siasiesiv.py
@@ -1,10 +1,12 @@
 # coding=utf-8
 """Compute the sea ice extent , area and volume  in both hemispheres or a specified region"""
-import os
 import six
+import os
 
 import numpy as np
 
+import netCDF4
+
 import iris
 import iris.analysis
 import iris.coords
@@ -16,6 +18,9 @@ from earthdiagnostics.diagnostic import Diagnostic, DiagnosticBasinListOption, D
 from earthdiagnostics.modelingrealm import ModelingRealms
 from earthdiagnostics.utils import Utils, TempFile
 
+import diagonals.siasie as siasie
+from diagonals.mesh_helpers.nemo import Nemo
+
 
 # noinspection PyUnresolvedReferences
 
@@ -40,10 +45,6 @@ class Siasiesiv(Diagnostic):
     alias = 'siasiesiv'
     "Diagnostic alias for the configuration file"
 
-    e1t = None
-    e2t = None
-    gphit = None
-
     def __init__(self, data_manager, startdate, member, chunk, masks, var_manager, data_convention, omit_vol):
         Diagnostic.__init__(self, data_manager)
         self.startdate = startdate
@@ -58,7 +59,7 @@ class Siasiesiv(Diagnostic):
         self.data_convention = data_convention
 
         self.results = {}
-        for var in ('siarean', 'siareas', 'sivoln', 'sivols', 'siextentn', 'siextents'):
+        for var in ('siarean', 'siareas', 'siextentn', 'siextents'):
             self.results[var] = {}
 
     def __str__(self):
@@ -78,7 +79,7 @@ class Siasiesiv(Diagnostic):
         :return:
         """
         options_available = (DiagnosticBasinListOption('basins', [Basins().Global]),
-                             DiagnosticBoolOption('omit_volume', False))
+                             DiagnosticBoolOption('omit_volume', True))
         options = cls.process_options(options, options_available)
 
         basins = options['basins']
@@ -97,10 +98,6 @@ class Siasiesiv(Diagnostic):
                                       diags.config.var_manager, diags.config.data_convention,
                                       options['omit_volume']))
 
-        e1t = iris.load_cube('mesh_hgr.nc', 'e1t')
-        e2t = iris.load_cube('mesh_hgr.nc', 'e2t')
-        Siasiesiv.area = e1t * e2t
-
         return job_list
 
     def request_data(self):
@@ -137,26 +134,14 @@ class Siasiesiv(Diagnostic):
         if sic.units.origin == '%' and sic.data.max() < 2:
             sic.units = '1.0'
 
-        extent = sic.copy((sic.data >= 0.15).astype(np.int8)) * Siasiesiv.area.data
-        sic *= Siasiesiv.area.data
-
-        if not self.omit_volume:
-            self._fix_coordinates_attribute(
-                self.sit.local_file, self.sit_varname
-            )
-            sit = iris.load_cube(self.sit.local_file)
-
-        for basin, mask in six.iteritems(self.masks):
-            self.results['siarean'][basin] = self.sum(sic, mask, north=True)
-            self.results['siareas'][basin] = self.sum(sic, mask, north=False)
-
-            if not self.omit_volume:
-                volume = sic * sit.data
-                self.results['sivoln'][basin] = self.sum(volume, mask, north=True)
-                self.results['sivols'][basin] = self.sum(volume, mask, north=False)
-
-            self.results['siextentn'][basin] = self.sum(extent, mask, north=True)
-            self.results['siextents'][basin] = self.sum(extent, mask, north=False)
+        sic_slices = []
+        for sic_data in sic.slices_over('time'):
+            sic_data.data = np.ma.filled(sic_data.data, 0.0).astype(np.float32)
+            sic_slices.append(sic_data)
+        mesh = Nemo('mesh_hgr.nc', 'mask_regions.nc')
+        areacello = mesh.get_areacello(cell_point='T')
+        gphit = mesh.get_grid_latitude(cell_point='T')
+        self.results['siextentn'], self.results['siextents'], self.results['siarean'], self.results['siareas'] = siasie.compute(gphit, areacello, sic_slices, self.masks)
 
         self.save()
 
@@ -171,29 +156,21 @@ class Siasiesiv(Diagnostic):
             ' '.join(set(coordinates) | add_coordinates)
         handler.close()
 
-    def sum(self, data, mask, north=True):
-        if north:
-            condition = data.coord('latitude').points > 0
-        else:
-            condition = data.coord('latitude').points < 0
-        weights = iris.util.broadcast_to_shape(condition, data.shape, data.coord_dims('latitude')) * mask
-        return data.collapsed(('latitude', 'longitude'), iris.analysis.SUM, weights=weights)
 
     def save(self):
-        for var, basins in six.iteritems(self.results):
+        for var, time in six.iteritems(self.results):
             results = iris.cube.CubeList()
-            for basin, result in six.iteritems(basins):
-                result.var_name = var
-                if var.startswith('sivol'):
-                    result.units = 'm^3'
-                else:
-                    result.units = 'm^2'
-                result.add_aux_coord(iris.coords.AuxCoord(basin.name, var_name='region'))
-                results.append(result)
-            if not results:
-                continue
+            for time, basin in six.iteritems(time):
+                for basin, result in six.iteritems(basin):
+                    cube = iris.cube.Cube(result)
+                    cube.var_name = var
+                    cube.units = 'm^2'
+                    cube.add_aux_coord(iris.coords.AuxCoord(basin.name, var_name='region'))
+                    cube.add_aux_coord(iris.coords.AuxCoord(time, var_name='time'))
+                    results.append(cube)
             self._save_file(results.merge_cube(), var)
 
+
     def _save_file(self, data, var):
         generated_file = self.generated[var]
         temp = TempFile.get()
@@ -201,7 +178,7 @@ class Siasiesiv(Diagnostic):
         data.remove_coord('region')
         iris.save(data, temp, zlib=True)
         if len(region) > 1:
-            Utils.rename_variable(temp, 'dim0', 'region', False)
+            Utils.rename_variable(temp, 'dim1', 'region', False)
             handler = Utils.open_cdf(temp)
             var = handler.createVariable('region2', str, ('region',))
             var[...] = region
-- 
GitLab


From 2527f469f7ac44f656182079dff8b92c88c08c4a Mon Sep 17 00:00:00 2001
From: sloosvel <saskia.loosveldt@bsc.es>
Date: Fri, 24 May 2019 10:16:24 +0200
Subject: [PATCH 10/30] Remove unused functions

---
 earthdiagnostics/ocean/regionmean.py | 89 ++++++++++++++--------------
 1 file changed, 45 insertions(+), 44 deletions(-)

diff --git a/earthdiagnostics/ocean/regionmean.py b/earthdiagnostics/ocean/regionmean.py
index b6ced885..792793a5 100644
--- a/earthdiagnostics/ocean/regionmean.py
+++ b/earthdiagnostics/ocean/regionmean.py
@@ -12,7 +12,8 @@ import netCDF4
 
 from earthdiagnostics.box import Box
 from earthdiagnostics.constants import Basins
-from earthdiagnostics.diagnostic import Diagnostic, DiagnosticOption, DiagnosticIntOption, DiagnosticDomainOption, \
+from earthdiagnostics.diagnostic import Diagnostic, DiagnosticOption, \
+    DiagnosticIntOption, DiagnosticDomainOption, \
     DiagnosticBoolOption, DiagnosticBasinListOption, DiagnosticVariableOption
 from earthdiagnostics.modelingrealm import ModelingRealms
 from earthdiagnostics.utils import Utils, TempFile
@@ -26,8 +27,9 @@ class RegionMean(Diagnostic):
     """
     Computes the mean value of the field (3D, weighted).
 
-    For 3D fields, a horizontal mean for each level is also given. If a spatial window
-    is specified, the mean value is computed only in this window.
+    For 3D fields, a horizontal mean for each level is also given.
+    If a spatial window is specified, the mean value is computed
+    only in this window.
 
     :original author: Javier Vegas-Regidor <javier.vegas@bsc.es>
 
@@ -50,8 +52,8 @@ class RegionMean(Diagnostic):
     alias = 'regmean'
     "Diagnostic alias for the configuration file"
 
-    def __init__(self, data_manager, startdate, member, chunk, domain, variable, box, save3d,
-                 variance, basins, grid_point):
+    def __init__(self, data_manager, startdate, member, chunk, domain,
+                 variable, box, save3d, variance, basins, grid_point):
         Diagnostic.__init__(self, data_manager)
         self.startdate = startdate
         self.member = member
@@ -72,12 +74,15 @@ class RegionMean(Diagnostic):
     def __eq__(self, other):
         if self._different_type(other):
             return False
-        return self.startdate == other.startdate and self.member == other.member and self.chunk == other.chunk and \
+        return self.startdate == other.startdate and \
+            self.member == other.member and self.chunk == other.chunk and \
             self.box == other.box and self.variable == other.variable
 
     def __str__(self):
-        return 'Region mean Startdate: {0.startdate} Member: {0.member} Chunk: {0.chunk} Variable: {0.variable} ' \
-               'Box: {0.box} Save 3D: {0.save3d} Save variance: {0.variance} Grid point: {0.grid_point}'.format(self)
+        return 'Region mean Startdate: {0.startdate} Member: {0.member}' \
+               'Chunk: {0.chunk} Variable: {0.variable} ' \
+               'Box: {0.box} Save 3D: {0.save3d} Save variance: {0.variance}' \
+               'Grid point: {0.grid_point}'.format(self)
 
     def __hash__(self):
         return hash(str(self))
@@ -96,7 +101,8 @@ class RegionMean(Diagnostic):
         options_available = (DiagnosticDomainOption(),
                              DiagnosticVariableOption(diags.data_manager.config.var_manager),
                              DiagnosticOption('grid_point', 'T'),
-                             DiagnosticBasinListOption('basins', Basins().Global),
+                             DiagnosticBasinListOption('basins',
+                                                       Basins().Global),
                              DiagnosticIntOption('min_depth', -1),
                              DiagnosticIntOption('max_depth', -1),
                              DiagnosticIntOption('min_lat', -1),
@@ -126,7 +132,8 @@ class RegionMean(Diagnostic):
         for startdate, member, chunk in diags.config.experiment.get_chunk_list():
             job = RegionMean(diags.data_manager, startdate, member, chunk,
                              options['domain'], options['variable'], box,
-                             options['save3D'], options['variance'], options['basins'],
+                             options['save3D'], options['variance'],
+                             options['basins'],
                              options['grid_point'].lower())
             job_list.append(job)
 
@@ -134,7 +141,9 @@ class RegionMean(Diagnostic):
 
     def request_data(self):
         """Request data required by the diagnostic"""
-        self.variable_file = self.request_chunk(self.domain, self.variable, self.startdate, self.member, self.chunk)
+        self.variable_file = self.request_chunk(self.domain, self.variable,
+                                                self.startdate, self.member,
+                                                self.chunk)
 
     def declare_data_generated(self):
         """Declare data to be generated by the diagnostic"""
@@ -156,8 +165,10 @@ class RegionMean(Diagnostic):
             masks[basin] = Utils.get_mask(basin)
 
         mesh = Nemo('mesh_hgr.nc', 'mask_regions.nc')
-        if self.box.min_lat is not -1 and self.box.max_lat is not -1 and self.box.min_lon is not -1 and self.box.max_lat is not -1:
-            name = '{0}_{1}'.format(Box.get_lat_str(self.box), Box.get_lon_str(self.box))
+        if self.box.min_lat is not -1 and self.box.max_lat is not -1 and \
+            self.box.min_lon is not -1 and self.box.max_lat is not -1:
+            name = '{0}_{1}'.format(Box.get_lat_str(self.box),
+                                    Box.get_lon_str(self.box))
 
             masks[name] = mesh.get_region_mask(self.box.min_lat,
                                                self.box.max_lat,
@@ -192,19 +203,23 @@ class RegionMean(Diagnostic):
 
     def _try_load_cube(self, number):
         try:
-            cube = iris.load_cube('mesh_hgr.nc', 'e{0}{1}'.format(number, self.grid_point))
+            cube = iris.load_cube('mesh_hgr.nc', 'e{0}{1}'.format(number,
+                                                                  self.grid_point))
         except iris.exceptions.ConstraintMismatchError:
-            cube = iris.load_cube('mesh_hgr.nc', 'e{0}{1}_0'.format(number, self.grid_point))
+            cube = iris.load_cube('mesh_hgr.nc', 'e{0}{1}_0'.format(number,
+                                                                    self.grid_point))
         cube = iris.util.squeeze(cube)
         dims = len(cube.shape)
         try:
             cube.coord('i')
         except iris.exceptions.CoordinateNotFoundError:
-            cube.add_dim_coord(iris.coords.DimCoord(np.arange(cube.shape[dims - 1]), var_name='i'), dims - 1)
+            cube.add_dim_coord(iris.coords.DimCoord(np.arange(cube.shape[dims - 1]),
+                                                    var_name='i'), dims - 1)
         try:
             cube.coord('j')
         except iris.exceptions.CoordinateNotFoundError:
-            cube.add_dim_coord(iris.coords.DimCoord(np.arange(cube.shape[dims - 2]), var_name='j'), dims - 2)
+            cube.add_dim_coord(iris.coords.DimCoord(np.arange(cube.shape[dims - 2]),
+                                                    var_name='j'), dims - 2)
         return cube
 
 
@@ -212,7 +227,8 @@ class RegionMean(Diagnostic):
         coords = []
         handler = Utils.open_cdf(self.variable_file.local_file)
         for variable in handler.variables:
-            if variable in ('time', 'lev', 'lat', 'lon', 'latitude', 'longitude', 'leadtime', 'time_centered'):
+            if variable in ('time', 'lev', 'lat', 'lon', 'latitude',
+                            'longitude', 'leadtime', 'time_centered'):
                 coords.append(variable)
             if variable == 'time_centered':
                 handler.variables[variable].standard_name = ''
@@ -238,7 +254,8 @@ class RegionMean(Diagnostic):
         coordinates = ''
         has_levels = False
         for dimension in handler.variables.keys():
-            if dimension in ['time', 'lev', 'lat', 'latitude', 'lon', 'longitude', 'i', 'j']:
+            if dimension in ['time', 'lev', 'lat', 'latitude',
+                             'lon', 'longitude', 'i', 'j']:
                 coordinates += ' {0}'.format(dimension)
             if dimension == 'lev':
                 has_levels = True
@@ -254,31 +271,13 @@ class RegionMean(Diagnostic):
         else:
             final_name = '{1}{0}'.format(var, self.variable)
 
-        self.declared[final_name] = self.declare_chunk(ModelingRealms.ocean, final_name, self.startdate, self.member,
-                                                       self.chunk, box=box_save, region=self.basins)
-
-    def _send_var(self, var, threed, cube_list):
-        if threed:
-            if not self.save3d and threed:
-                return False
-            final_name = '{1}3d{0}'.format(var, self.variable)
-        else:
-            final_name = '{1}{0}'.format(var, self.variable)
-        cube = cube_list.merge_cube()
-        cube.var_name = 'result'
-        cube.remove_coord('latitude')
-        cube.remove_coord('longitude')
-        try:
-            cube.remove_coord('depth')
-        except iris.exceptions.CoordinateNotFoundError:
-            pass
-        try:
-            cube.remove_coord('lev')
-        except iris.exceptions.CoordinateNotFoundError:
-            pass
-        temp = TempFile.get()
-        iris.save(cube, temp)
-        self.declared[final_name].set_local_file(temp, diagnostic=self, rename_var='result', region=self.basins)
+        self.declared[final_name] = self.declare_chunk(ModelingRealms.ocean,
+                                                       final_name,
+                                                       self.startdate,
+                                                       self.member,
+                                                       self.chunk,
+                                                       box=box_save,
+                                                       region=self.basins)
 
 
     def _save_result_2D(self, var, result, data):
@@ -319,8 +318,10 @@ class RegionMean(Diagnostic):
         var_level.standard_name = 'depth'
         var_region = handler_temp.createVariable('region', 'S1',
                                                  ('region', 'region_length'))
+
         var = handler_temp.createVariable('{1}3d{0}'.format(var, self.variable),
                                           float, ('time', 'lev', 'region',),)
+                                          
         var.units = '{0}'.format(data.units)
         for i, basin in enumerate(result):
             var_region[i, ...] = netCDF4.stringtoarr(str(basin), 50)
-- 
GitLab


From d45911813d141b124eb4d0eed52e49f943caeac4 Mon Sep 17 00:00:00 2001
From: sloosvel <saskia.loosveldt@bsc.es>
Date: Fri, 24 May 2019 10:17:52 +0200
Subject: [PATCH 11/30] Fix some style issues

---
 earthdiagnostics/ocean/moc.py | 11 ++++++++---
 1 file changed, 8 insertions(+), 3 deletions(-)

diff --git a/earthdiagnostics/ocean/moc.py b/earthdiagnostics/ocean/moc.py
index 44e1301d..b20ea03e 100644
--- a/earthdiagnostics/ocean/moc.py
+++ b/earthdiagnostics/ocean/moc.py
@@ -67,7 +67,8 @@ class Moc(Diagnostic):
     def __eq__(self, other):
         if self._different_type(other):
             return False
-        return self.startdate == other.startdate and self.member == other.member and self.chunk == other.chunk
+        return self.startdate == other.startdate and \
+               self.member == other.member and self.chunk == other.chunk
 
     @classmethod
     def generate_jobs(cls, diags, options):
@@ -101,11 +102,15 @@ class Moc(Diagnostic):
 
     def request_data(self):
         """Request data required by the diagnostic"""
-        self.variable_file = self.request_chunk(ModelingRealms.ocean, 'vo', self.startdate, self.member, self.chunk)
+        self.variable_file = self.request_chunk(ModelingRealms.ocean, 'vo',
+                                                self.startdate, self.member,
+                                                self.chunk)
 
     def declare_data_generated(self):
         """Declare data to be generated by the diagnostic"""
-        self.results = self.declare_chunk(ModelingRealms.ocean, Moc.vsftmyz, self.startdate, self.member, self.chunk)
+        self.results = self.declare_chunk(ModelingRealms.ocean, Moc.vsftmyz,
+                                          self.startdate, self.member,
+                                          self.chunk)
 
     def compute(self):
         """Run the diagnostic"""
-- 
GitLab


From b838ba4eb3ae0e1d3c60273eefa829d72610c970 Mon Sep 17 00:00:00 2001
From: sloosvel <saskia.loosveldt@bsc.es>
Date: Fri, 24 May 2019 10:41:05 +0200
Subject: [PATCH 12/30] Load missing values as 0

---
 earthdiagnostics/ocean/moc.py | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/earthdiagnostics/ocean/moc.py b/earthdiagnostics/ocean/moc.py
index b20ea03e..5967b606 100644
--- a/earthdiagnostics/ocean/moc.py
+++ b/earthdiagnostics/ocean/moc.py
@@ -114,7 +114,8 @@ class Moc(Diagnostic):
 
     def compute(self):
         """Run the diagnostic"""
-        vo = iris.load_cube(self.variable_file.local_file)
+        vo_cube = iris.load_cube(self.variable_file.local_file)
+        vo = np.ma.filled(vo_cube.data, 0.0).astype(np.float32)
         mesh = Nemo('mesh_hgr.nc', 'mask_regions.nc')
         e1v = mesh.get_i_length(cell_point='V')
         e3v = mesh.get_k_length(cell_point='V')
@@ -130,7 +131,7 @@ class Moc(Diagnostic):
             else:
                 masks[basin] = Utils.get_mask(basin)
 
-        moc_results = moc.compute(masks, e1v, e3v, vo.data)
+        moc_results = moc.compute(masks, e1v, e3v, vo)
         del vo, e1v, e3v
         self._save_result(moc_results, mesh)
 
-- 
GitLab


From 8ca727c4284fadae963c4477f6b3a0a003a65607 Mon Sep 17 00:00:00 2001
From: sloosvel <saskia.loosveldt@bsc.es>
Date: Fri, 24 May 2019 10:43:12 +0200
Subject: [PATCH 13/30] Fix style issue

---
 earthdiagnostics/ocean/moc.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/earthdiagnostics/ocean/moc.py b/earthdiagnostics/ocean/moc.py
index 5967b606..7732726d 100644
--- a/earthdiagnostics/ocean/moc.py
+++ b/earthdiagnostics/ocean/moc.py
@@ -97,7 +97,8 @@ class Moc(Diagnostic):
 
         job_list = list()
         for startdate, member, chunk in diags.config.experiment.get_chunk_list():
-            job_list.append(Moc(diags.data_manager, startdate, member, chunk, basins))
+            job_list.append(Moc(diags.data_manager, startdate, member, chunk,
+                                basins))
         return job_list
 
     def request_data(self):
-- 
GitLab


From 7c893232350252f49b4a7afa3338f93b6ad9b7b1 Mon Sep 17 00:00:00 2001
From: sloosvel <saskia.loosveldt@bsc.es>
Date: Fri, 24 May 2019 10:46:40 +0200
Subject: [PATCH 14/30] Fix some style issues

---
 earthdiagnostics/ocean/siasiesiv.py | 43 +++++++++++++++++++----------
 1 file changed, 28 insertions(+), 15 deletions(-)

diff --git a/earthdiagnostics/ocean/siasiesiv.py b/earthdiagnostics/ocean/siasiesiv.py
index c1102b02..b79efd89 100644
--- a/earthdiagnostics/ocean/siasiesiv.py
+++ b/earthdiagnostics/ocean/siasiesiv.py
@@ -1,7 +1,7 @@
 # coding=utf-8
 """Compute the sea ice extent , area and volume  in both hemispheres or a specified region"""
-import six
 import os
+import six
 
 import numpy as np
 
@@ -14,7 +14,8 @@ import iris.util
 from bscearth.utils.log import Log
 
 from earthdiagnostics.constants import Basins
-from earthdiagnostics.diagnostic import Diagnostic, DiagnosticBasinListOption, DiagnosticBoolOption
+from earthdiagnostics.diagnostic import Diagnostic, \
+    DiagnosticBasinListOption, DiagnosticBoolOption
 from earthdiagnostics.modelingrealm import ModelingRealms
 from earthdiagnostics.utils import Utils, TempFile
 
@@ -27,7 +28,8 @@ from diagonals.mesh_helpers.nemo import Nemo
 
 class Siasiesiv(Diagnostic):
     """
-    Compute the sea ice extent , area and volume  in both hemispheres or a specified region.
+    Compute the sea ice extent , area and volume  in both hemispheres or a
+    specified region.
 
     Parameters
     ----------
@@ -45,7 +47,8 @@ class Siasiesiv(Diagnostic):
     alias = 'siasiesiv'
     "Diagnostic alias for the configuration file"
 
-    def __init__(self, data_manager, startdate, member, chunk, masks, var_manager, data_convention, omit_vol):
+    def __init__(self, data_manager, startdate, member, chunk, masks,
+                 var_manager, data_convention, omit_vol):
         Diagnostic.__init__(self, data_manager)
         self.startdate = startdate
         self.member = member
@@ -78,7 +81,8 @@ class Siasiesiv(Diagnostic):
         :type options: list[str]
         :return:
         """
-        options_available = (DiagnosticBasinListOption('basins', [Basins().Global]),
+        options_available = (DiagnosticBasinListOption('basins',
+                                                       [Basins().Global]),
                              DiagnosticBoolOption('omit_volume', True))
         options = cls.process_options(options, options_available)
 
@@ -94,8 +98,9 @@ class Siasiesiv(Diagnostic):
 
         job_list = list()
         for startdate, member, chunk in diags.config.experiment.get_chunk_list():
-            job_list.append(Siasiesiv(diags.data_manager, startdate, member, chunk, masks,
-                                      diags.config.var_manager, diags.config.data_convention,
+            job_list.append(Siasiesiv(diags.data_manager, startdate, member,
+                                      chunk, masks, diags.config.var_manager,
+                                      diags.config.data_convention,
                                       options['omit_volume']))
 
         return job_list
@@ -103,8 +108,11 @@ class Siasiesiv(Diagnostic):
     def request_data(self):
         """Request data required by the diagnostic"""
         if not self.omit_volume:
-            self.sit = self.request_chunk(ModelingRealms.seaIce, self.sit_varname,
-                                          self.startdate, self.member, self.chunk)
+            self.sit = self.request_chunk(ModelingRealms.seaIce,
+                                          self.sit_varname,
+                                          self.startdate,
+                                          self.member,
+                                          self.chunk)
         self.sic = self.request_chunk(ModelingRealms.seaIce, self.sic_varname,
                                       self.startdate, self.member, self.chunk)
 
@@ -121,8 +129,9 @@ class Siasiesiv(Diagnostic):
         self._declare_var('siextentn')
 
     def _declare_var(self, var_name):
-        self.generated[var_name] = self.declare_chunk(ModelingRealms.seaIce, var_name,
-                                                      self.startdate, self.member, self.chunk)
+        self.generated[var_name] = self.declare_chunk(ModelingRealms.seaIce,
+                                                      var_name, self.startdate,
+                                                      self.member, self.chunk)
 
     def compute(self):
         """Run the diagnostic"""
@@ -165,8 +174,10 @@ class Siasiesiv(Diagnostic):
                     cube = iris.cube.Cube(result)
                     cube.var_name = var
                     cube.units = 'm^2'
-                    cube.add_aux_coord(iris.coords.AuxCoord(basin.name, var_name='region'))
-                    cube.add_aux_coord(iris.coords.AuxCoord(time, var_name='time'))
+                    cube.add_aux_coord(iris.coords.AuxCoord(basin.name,
+                                                            var_name='region'))
+                    cube.add_aux_coord(iris.coords.AuxCoord(time,
+                                                            var_name='time'))
                     results.append(cube)
             self._save_file(results.merge_cube(), var)
 
@@ -199,12 +210,14 @@ class Siasiesiv(Diagnostic):
                         continue
                     Utils.copy_variable(handler, new_handler, variable)
                 old_var = handler.variables[var]
-                new_var = new_handler.createVariable(var, old_var.dtype, ('region',) + old_var.dimensions,
+                new_var = new_handler.createVariable(var, old_var.dtype,
+                                                     ('region',) + old_var.dimensions,
                                                      zlib=True, fill_value=1.0e20)
                 Utils.copy_attributes(new_var, old_var)
                 new_var[0, :] = old_var[:]
 
-                new_var = new_handler.createVariable('region', str, ('region',))
+                new_var = new_handler.createVariable('region', str,
+                                                     ('region',))
                 new_var[0] = region[0]
 
                 new_handler.close()
-- 
GitLab


From fed94ae1a9df0a07b41ccb4f46e004ddc1b9a2ad Mon Sep 17 00:00:00 2001
From: sloosvel <saskia.loosveldt@bsc.es>
Date: Fri, 24 May 2019 18:21:44 +0200
Subject: [PATCH 15/30] Save creating netcdf file manually

---
 earthdiagnostics/ocean/siasiesiv.py | 78 ++++++++---------------------
 1 file changed, 20 insertions(+), 58 deletions(-)

diff --git a/earthdiagnostics/ocean/siasiesiv.py b/earthdiagnostics/ocean/siasiesiv.py
index b79efd89..5cf4a905 100644
--- a/earthdiagnostics/ocean/siasiesiv.py
+++ b/earthdiagnostics/ocean/siasiesiv.py
@@ -167,61 +167,23 @@ class Siasiesiv(Diagnostic):
 
 
     def save(self):
-        for var, time in six.iteritems(self.results):
-            results = iris.cube.CubeList()
-            for time, basin in six.iteritems(time):
-                for basin, result in six.iteritems(basin):
-                    cube = iris.cube.Cube(result)
-                    cube.var_name = var
-                    cube.units = 'm^2'
-                    cube.add_aux_coord(iris.coords.AuxCoord(basin.name,
-                                                            var_name='region'))
-                    cube.add_aux_coord(iris.coords.AuxCoord(time,
-                                                            var_name='time'))
-                    results.append(cube)
-            self._save_file(results.merge_cube(), var)
-
-
-    def _save_file(self, data, var):
-        generated_file = self.generated[var]
-        temp = TempFile.get()
-        region = data.coord('region').points
-        data.remove_coord('region')
-        iris.save(data, temp, zlib=True)
-        if len(region) > 1:
-            Utils.rename_variable(temp, 'dim1', 'region', False)
-            handler = Utils.open_cdf(temp)
-            var = handler.createVariable('region2', str, ('region',))
-            var[...] = region
-            handler.close()
-            Utils.rename_variable(temp, 'region2', 'region', True)
-        else:
-            handler = Utils.open_cdf(temp)
-            if 'region' not in handler.dimensions:
-                new_file = TempFile.get()
-                new_handler = Utils.open_cdf(new_file, 'w')
-
-                new_handler.createDimension('region', 1)
-                for dimension in handler.dimensions:
-                    Utils.copy_dimension(handler, new_handler, dimension)
-
-                for variable in handler.variables.keys():
-                    if variable in (var, 'region'):
-                        continue
-                    Utils.copy_variable(handler, new_handler, variable)
-                old_var = handler.variables[var]
-                new_var = new_handler.createVariable(var, old_var.dtype,
-                                                     ('region',) + old_var.dimensions,
-                                                     zlib=True, fill_value=1.0e20)
-                Utils.copy_attributes(new_var, old_var)
-                new_var[0, :] = old_var[:]
-
-                new_var = new_handler.createVariable('region', str,
-                                                     ('region',))
-                new_var[0] = region[0]
-
-                new_handler.close()
-                os.remove(temp)
-                temp = new_file
-            handler.close()
-        generated_file.set_local_file(temp)
+        for var in self.results.keys():
+            res = self.results[var]
+            temp = TempFile.get()
+            handler_source = Utils.open_cdf(self.sic.local_file)
+            handler_temp = Utils.open_cdf(temp, 'w')
+            Utils.copy_variable(handler_source, handler_temp, 'time', True, True)
+            handler_temp.createDimension('region', len(self.masks))
+            handler_temp.createDimension('region_length', 50)
+            var_region = handler_temp.createVariable('region', 'S1',
+                                                 ('region', 'region_length'))
+            var_res = handler_temp.createVariable('{0}'.format(var), float,
+                                                  ('time', 'region',))
+            var_res.units = 'm^2'
+            for i, basin in enumerate(self.masks):
+                var_region[i, ...] = netCDF4.stringtoarr(str(basin), 50)
+                var_res[..., i] = res[i, ...]
+            handler_temp.close()
+            self.generated[var].set_local_file(temp, diagnostic=self)
+
+
-- 
GitLab


From 35d5165bc9cfd0b2aca04cf91f545f5bae6461c1 Mon Sep 17 00:00:00 2001
From: sloosvel <saskia.loosveldt@bsc.es>
Date: Mon, 27 May 2019 11:01:45 +0200
Subject: [PATCH 16/30] Add diagonals to list of packages

---
 setup.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/setup.py b/setup.py
index e8e0e31d..32640d3f 100644
--- a/setup.py
+++ b/setup.py
@@ -78,6 +78,7 @@ setup(
         'scitools-iris>=2.2',
         'six',
         'xxhash',
+        'diagonals'
     ],
     packages=find_packages(),
     include_package_data=True,
-- 
GitLab


From a23923a75ce9fc75b93a0c0981e3b4aabfa52d49 Mon Sep 17 00:00:00 2001
From: Javier Vegas-Regidor <javier.vegas@bsc.es>
Date: Mon, 27 May 2019 11:33:36 +0200
Subject: [PATCH 17/30] Update CI to install diagonals from repository

---
 .gitlab-ci.yml | 25 +++++++++++++++++--------
 1 file changed, 17 insertions(+), 8 deletions(-)

diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 813b67e7..ea47b884 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -1,6 +1,7 @@
 before_script:
   - export GIT_SSL_NO_VERIFY=1
   - export PATH="$HOME/miniconda2/bin:$PATH"
+  - export DIAGONALS_BRANCH="production"
 
 stages:
   - prepare
@@ -24,18 +25,26 @@ test_python2:
     - git submodule update --init --recursive
     - conda env update -f environment.yml -n earthdiagnostics2 python=2.7
     - source activate earthdiagnostics2
+    - cd ..
+    - git clone -b ${DIAGONALS_BRANCH} https://earth.bsc.es/gitlab/es/diagonals.git
+    - pip install ./diagonals
+    - cd earthdiagnostics
     - pip install -e .
     - python setup.py test
 
 test_python3:
- stage: test
- script:
-   - git submodule sync --recursive
-   - git submodule update --init --recursive
-   - conda env update -f environment.yml -n earthdiagnostics3 python=3.7
-   - source activate earthdiagnostics3
-   - pip install -e .
-   - python setup.py test
+  stage: test
+  script:
+    - git submodule sync --recursive
+    - git submodule update --init --recursive
+    - conda env update -f environment.yml -n earthdiagnostics3 python=3.7
+    - source activate earthdiagnostics3
+    - cd ..
+    - git clone -b ${DIAGONALS_BRANCH} https://earth.bsc.es/gitlab/es/diagonals.git
+    - pip install ./diagonals
+    - cd earthdiagnostics
+    - pip install -e .
+    - python setup.py test
 
 report_codacy:
   stage: report
-- 
GitLab


From ca16b49c01585868d0516a15f85d7451877c7455 Mon Sep 17 00:00:00 2001
From: Javier Vegas-Regidor <javier.vegas@bsc.es>
Date: Mon, 27 May 2019 11:47:14 +0200
Subject: [PATCH 18/30] Delete diagonals folder

---
 .gitlab-ci.yml | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index ea47b884..61d2b4f3 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -40,8 +40,10 @@ test_python3:
     - conda env update -f environment.yml -n earthdiagnostics3 python=3.7
     - source activate earthdiagnostics3
     - cd ..
+    - rm -rf ./diagonals
     - git clone -b ${DIAGONALS_BRANCH} https://earth.bsc.es/gitlab/es/diagonals.git
     - pip install ./diagonals
+    - rm -r ./diagonals
     - cd earthdiagnostics
     - pip install -e .
     - python setup.py test
-- 
GitLab


From d699d92e287ac6eadfc5a6e040a5ebe16c6c8758 Mon Sep 17 00:00:00 2001
From: Javier Vegas-Regidor <javier.vegas@bsc.es>
Date: Mon, 27 May 2019 11:53:16 +0200
Subject: [PATCH 19/30] Let CI fail until released in pip

---
 .gitlab-ci.yml | 15 ++-------------
 1 file changed, 2 insertions(+), 13 deletions(-)

diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 61d2b4f3..e7a4660d 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -1,7 +1,6 @@
 before_script:
   - export GIT_SSL_NO_VERIFY=1
   - export PATH="$HOME/miniconda2/bin:$PATH"
-  - export DIAGONALS_BRANCH="production"
 
 stages:
   - prepare
@@ -25,11 +24,7 @@ test_python2:
     - git submodule update --init --recursive
     - conda env update -f environment.yml -n earthdiagnostics2 python=2.7
     - source activate earthdiagnostics2
-    - cd ..
-    - git clone -b ${DIAGONALS_BRANCH} https://earth.bsc.es/gitlab/es/diagonals.git
-    - pip install ./diagonals
-    - cd earthdiagnostics
-    - pip install -e .
+    - pip install .
     - python setup.py test
 
 test_python3:
@@ -39,13 +34,7 @@ test_python3:
     - git submodule update --init --recursive
     - conda env update -f environment.yml -n earthdiagnostics3 python=3.7
     - source activate earthdiagnostics3
-    - cd ..
-    - rm -rf ./diagonals
-    - git clone -b ${DIAGONALS_BRANCH} https://earth.bsc.es/gitlab/es/diagonals.git
-    - pip install ./diagonals
-    - rm -r ./diagonals
-    - cd earthdiagnostics
-    - pip install -e .
+    - pip install .
     - python setup.py test
 
 report_codacy:
-- 
GitLab


From 7e830fed0582091eb731b0e5465ffc5d13fa7ea6 Mon Sep 17 00:00:00 2001
From: sloosvel <saskia.loosveldt@bsc.es>
Date: Fri, 21 Jun 2019 15:29:43 +0200
Subject: [PATCH 20/30] Add sivol computing functions from diagonals

---
 earthdiagnostics/ocean/siasiesiv.py | 10 +++++++++-
 1 file changed, 9 insertions(+), 1 deletion(-)

diff --git a/earthdiagnostics/ocean/siasiesiv.py b/earthdiagnostics/ocean/siasiesiv.py
index 5cf4a905..0bf84947 100644
--- a/earthdiagnostics/ocean/siasiesiv.py
+++ b/earthdiagnostics/ocean/siasiesiv.py
@@ -150,7 +150,13 @@ class Siasiesiv(Diagnostic):
         mesh = Nemo('mesh_hgr.nc', 'mask_regions.nc')
         areacello = mesh.get_areacello(cell_point='T')
         gphit = mesh.get_grid_latitude(cell_point='T')
-        self.results['siextentn'], self.results['siextents'], self.results['siarean'], self.results['siareas'] = siasie.compute(gphit, areacello, sic_slices, self.masks)
+
+        if not self.omit_volume:
+            sit = iris.load_cube(self.sit.local_file)
+            self.results['siextentn'], self.results['siextents'], self.results['siarean'], self.results['siareas'], self.results['sivoln'], self.results['sivols'] = siasie.compute(gphit, areacello, sic_slices, self.masks, sit)
+        else:
+            self.results['siextentn'], self.results['siextents'], self.results['siarean'], self.results['siareas'] = siasie.compute(gphit, areacello, sic_slices, self.masks, None)
+
 
         self.save()
 
@@ -180,6 +186,8 @@ class Siasiesiv(Diagnostic):
             var_res = handler_temp.createVariable('{0}'.format(var), float,
                                                   ('time', 'region',))
             var_res.units = 'm^2'
+            if var in ('voln', 'vols'):
+                var_res.units = 'm^3'
             for i, basin in enumerate(self.masks):
                 var_region[i, ...] = netCDF4.stringtoarr(str(basin), 50)
                 var_res[..., i] = res[i, ...]
-- 
GitLab


From 19ce7f19c660ce1d2d64434a653bc9c2761a7844 Mon Sep 17 00:00:00 2001
From: sloosvel <saskia.loosveldt@bsc.es>
Date: Tue, 25 Jun 2019 11:53:49 +0200
Subject: [PATCH 21/30] Remove ugly hack

---
 bin/earthdiags | 4 ----
 1 file changed, 4 deletions(-)

diff --git a/bin/earthdiags b/bin/earthdiags
index 200086ee..99a4a377 100644
--- a/bin/earthdiags
+++ b/bin/earthdiags
@@ -5,10 +5,6 @@
 import os
 import sys
 
-scriptdir = os.path.abspath(os.path.dirname(sys.argv[0]))
-assert sys.path[0] == scriptdir
-sys.path[0] = os.path.normpath(os.path.join(scriptdir, os.pardir))
-
 # noinspection PyUnresolvedReferences,PyPep8
 from earthdiagnostics.earthdiags import EarthDiags
 
-- 
GitLab


From 8ffc67d7c133fb680354cbd4899cbf596f4ce28d Mon Sep 17 00:00:00 2001
From: sloosvel <saskia.loosveldt@bsc.es>
Date: Fri, 28 Jun 2019 16:02:39 +0200
Subject: [PATCH 22/30] Remove things related to hash

---
 earthdiagnostics/datafile.py |  4 +-
 earthdiagnostics/utils.py    | 88 ++++++++++++++++++------------------
 2 files changed, 46 insertions(+), 46 deletions(-)

diff --git a/earthdiagnostics/datafile.py b/earthdiagnostics/datafile.py
index a3d0ad22..7c903e1e 100644
--- a/earthdiagnostics/datafile.py
+++ b/earthdiagnostics/datafile.py
@@ -648,11 +648,11 @@ class NetCDFFile(DataFile):
             Log.debug('Downloading file {0}...', self.remote_file)
             if not self.local_file:
                 self.local_file = TempFile.get()
-            Utils.get_file_hash(self.remote_file, use_stored=True, save=True)
+            #Utils.get_file_hash(self.remote_file, use_stored=True, save=True)
             try:
                 Utils.copy_file(self.remote_file, self.local_file, retrials=1)
             except Utils.CopyException:
-                Utils.get_file_hash(self.remote_file, use_stored=False, save=True)
+             #   Utils.get_file_hash(self.remote_file, use_stored=False, save=True)
                 Utils.copy_file(self.remote_file, self.local_file, retrials=2)
 
             if self.data_convention == 'meteofrance':
diff --git a/earthdiagnostics/utils.py b/earthdiagnostics/utils.py
index c08c1f08..0920a365 100644
--- a/earthdiagnostics/utils.py
+++ b/earthdiagnostics/utils.py
@@ -328,21 +328,21 @@ class Utils(object):
                 # This can be due to a race condition. If directory already exists, we don have to do nothing
                 if not os.path.exists(dirname_path):
                     raise ex
-        hash_destiny = None
-        Log.debug('Hashing original file... {0}', datetime.datetime.now())
-        hash_original = Utils.get_file_hash(source, use_stored=use_stored_hash)
-
-        if retrials < 1:
-            retrials = 1
-
-        while hash_original != hash_destiny:
-            if retrials == 0:
-                raise Utils.CopyException('Can not copy {0} to {1}'.format(source, destiny))
-            Log.debug('Copying... {0}', datetime.datetime.now())
-            shutil.copyfile(source, destiny)
-            Log.debug('Hashing copy ... {0}', datetime.datetime.now())
-            hash_destiny = Utils.get_file_hash(destiny, save=save_hash)
-            retrials -= 1
+        #hash_destiny = None
+        #Log.debug('Hashing original file... {0}', datetime.datetime.now())
+        #hash_original = Utils.get_file_hash(source, use_stored=use_stored_hash)
+
+       # if retrials < 1:
+       #     retrials = 1
+
+       # while hash_original != hash_destiny:
+       #     if retrials == 0:
+       #         raise Utils.CopyException('Can not copy {0} to {1}'.format(source, destiny))
+        Log.debug('Copying... {0}', datetime.datetime.now())
+        shutil.copyfile(source, destiny)
+       #     Log.debug('Hashing copy ... {0}', datetime.datetime.now())
+            #hash_destiny = Utils.get_file_hash(destiny, save=save_hash)
+        #    retrials -= 1
         Log.debug('Finished {0}', datetime.datetime.now())
 
     @staticmethod
@@ -426,8 +426,8 @@ class Utils(object):
         time.sleep(2)
         shutil.rmtree(source)
 
-    @staticmethod
-    def get_file_hash(filepath, use_stored=False, save=False):
+ #   @staticmethod
+   # def get_file_hash(filepath, use_stored=False, save=False):
         """
         Get the xxHash hash for a given file
 
@@ -439,33 +439,33 @@ class Utils(object):
         save: bool, optional
             If True, saves the hash to a file
         """
-        if use_stored:
-            hash_file = Utils._get_hash_filename(filepath)
-            if os.path.isfile(hash_file):
-                hash_value = open(hash_file, 'r').readline()
-                return hash_value
-
-        blocksize = 104857600
-        hasher = xxhash.xxh64()
-        with open(filepath, 'rb') as afile:
-            buf = afile.read(blocksize)
-            while len(buf) > 0:
-                hasher.update(buf)
-                buf = afile.read(blocksize)
-        hash_value = hasher.hexdigest()
-        if save:
-            hash_file = open(Utils._get_hash_filename(filepath), 'w')
-            hash_file.write(hash_value)
-            hash_file.close()
-
-        return hash_value
-
-    @staticmethod
-    def _get_hash_filename(filepath):
-        folder = os.path.dirname(filepath)
-        filename = os.path.basename(filepath)
-        hash_file = os.path.join(folder, '.{0}.xxhash64.hash'.format(filename))
-        return hash_file
+      #  if use_stored:
+      #      hash_file = Utils._get_hash_filename(filepath)
+      #      if os.path.isfile(hash_file):
+      #          hash_value = open(hash_file, 'r').readline()
+      #          return hash_value
+
+     #   blocksize = 104857600
+     #   hasher = xxhash.xxh64()
+     #   with open(filepath, 'rb') as afile:
+     #       buf = afile.read(blocksize)
+     #       while len(buf) > 0:
+     #           hasher.update(buf)
+     #           buf = afile.read(blocksize)
+     #   hash_value = hasher.hexdigest()
+     #   if save:
+     #       hash_file = open(Utils._get_hash_filename(filepath), 'w')
+     #       hash_file.write(hash_value)
+     #       hash_file.close()
+
+    #    return hash_value
+
+  #  @staticmethod
+  #  def _get_hash_filename(filepath):
+  #      folder = os.path.dirname(filepath)
+  #      filename = os.path.basename(filepath)
+  #      hash_file = os.path.join(folder, '.{0}.xxhash64.hash'.format(filename))
+  #      return hash_file
 
     @staticmethod
     def execute_shell_command(command, log_level=Log.DEBUG):
-- 
GitLab


From b30026b5fd0bfaf6a72f9bd3252b4464ffeae1dd Mon Sep 17 00:00:00 2001
From: sloosvel <saskia.loosveldt@bsc.es>
Date: Fri, 28 Jun 2019 16:24:40 +0200
Subject: [PATCH 23/30] Fix bug when saving volume variables

---
 earthdiagnostics/ocean/siasiesiv.py | 14 +++++++++++---
 1 file changed, 11 insertions(+), 3 deletions(-)

diff --git a/earthdiagnostics/ocean/siasiesiv.py b/earthdiagnostics/ocean/siasiesiv.py
index 0bf84947..cc65228e 100644
--- a/earthdiagnostics/ocean/siasiesiv.py
+++ b/earthdiagnostics/ocean/siasiesiv.py
@@ -70,6 +70,9 @@ class Siasiesiv(Diagnostic):
                'Basins: {1} Omit volume: {0.omit_volume}'.format(self,
                                                                  ','.join(str(basin) for basin in self.masks.keys()))
 
+    def __hash__(self):
+        return hash(str(self))
+
     @classmethod
     def generate_jobs(cls, diags, options):
         """
@@ -153,7 +156,11 @@ class Siasiesiv(Diagnostic):
 
         if not self.omit_volume:
             sit = iris.load_cube(self.sit.local_file)
-            self.results['siextentn'], self.results['siextents'], self.results['siarean'], self.results['siareas'], self.results['sivoln'], self.results['sivols'] = siasie.compute(gphit, areacello, sic_slices, self.masks, sit)
+            sit_slices = []
+            for sit_data in sit.slices_over('time'):
+                sit_data.data = np.ma.filled(sit_data.data, 0.0).astype(np.float32)
+                sit_slices.append(sit_data)
+            self.results['siextentn'], self.results['siextents'], self.results['siarean'], self.results['siareas'], self.results['sivoln'], self.results['sivols'] = siasie.compute(gphit, areacello, sic_slices, self.masks, sit_slices)
         else:
             self.results['siextentn'], self.results['siextents'], self.results['siarean'], self.results['siareas'] = siasie.compute(gphit, areacello, sic_slices, self.masks, None)
 
@@ -185,9 +192,10 @@ class Siasiesiv(Diagnostic):
                                                  ('region', 'region_length'))
             var_res = handler_temp.createVariable('{0}'.format(var), float,
                                                   ('time', 'region',))
-            var_res.units = 'm^2'
-            if var in ('voln', 'vols'):
+            if var in ('sivoln', 'sivols'):
                 var_res.units = 'm^3'
+            else:
+                var_res.units = 'm^2'
             for i, basin in enumerate(self.masks):
                 var_region[i, ...] = netCDF4.stringtoarr(str(basin), 50)
                 var_res[..., i] = res[i, ...]
-- 
GitLab


From 9c592aa539805b83db6838cbe81eb6234ea16c03 Mon Sep 17 00:00:00 2001
From: sloosvel <saskia.loosveldt@bsc.es>
Date: Tue, 2 Jul 2019 09:32:17 +0200
Subject: [PATCH 24/30] Integrate diagonals

---
 earthdiagnostics/ocean/zonalmean.py | 128 +++++++++++++---------------
 1 file changed, 59 insertions(+), 69 deletions(-)

diff --git a/earthdiagnostics/ocean/zonalmean.py b/earthdiagnostics/ocean/zonalmean.py
index 5e917790..b7b9908e 100644
--- a/earthdiagnostics/ocean/zonalmean.py
+++ b/earthdiagnostics/ocean/zonalmean.py
@@ -11,13 +11,18 @@ from iris.cube import Cube, CubeList
 import numpy as np
 import numba
 
+import netCDF4
+
 from earthdiagnostics.box import Box
 from earthdiagnostics.constants import Basins
 from earthdiagnostics.diagnostic import Diagnostic, DiagnosticOption, DiagnosticIntOption, DiagnosticDomainOption, \
-    DiagnosticBoolOption, DiagnosticBasinOption, DiagnosticVariableOption
+    DiagnosticBoolOption, DiagnosticBasinListOption, DiagnosticVariableOption
 from earthdiagnostics.modelingrealm import ModelingRealms
 from earthdiagnostics.utils import Utils, TempFile
 
+import diagonals.zonmean as zonmean
+from diagonals.mesh_helpers.nemo import Nemo
+
 
 class ZonalMean(Diagnostic):
     """
@@ -40,14 +45,14 @@ class ZonalMean(Diagnostic):
     alias = 'zonmean'
     "Diagnostic alias for the configuration file"
 
-    def __init__(self, data_manager, startdate, member, chunk, domain, variable, basin, grid_point):
+    def __init__(self, data_manager, startdate, member, chunk, domain, variable, basins, grid_point):
         Diagnostic.__init__(self, data_manager)
         self.startdate = startdate
         self.member = member
         self.chunk = chunk
         self.domain = domain
         self.variable = variable
-        self.basin = basin
+        self.basins = basins
         self.grid_point = grid_point
 
         self.declared = {}
@@ -83,14 +88,14 @@ class ZonalMean(Diagnostic):
             DiagnosticDomainOption(),
             DiagnosticVariableOption(diags.data_manager.config.var_manager),
             DiagnosticOption('grid_point', 'T'),
-            DiagnosticBasinOption('basin', Basins().Global),
+            DiagnosticBasinListOption('basins', Basins().Global),
         )
         options = cls.process_options(options, options_available)
 
         job_list = list()
         for startdate, member, chunk in diags.config.experiment.get_chunk_list():
             job = ZonalMean(diags.data_manager, startdate, member, chunk,
-                            options['domain'], options['variable'], options['basin'],
+                            options['domain'], options['variable'], options['basins'],
                             options['grid_point'].lower())
             job_list.append(job)
 
@@ -106,76 +111,24 @@ class ZonalMean(Diagnostic):
         """Declare data to be generated by the diagnostic"""
         self.zonal_mean = self.declare_chunk(
             ModelingRealms.ocean, self.variable + 'zonal',
-            self.startdate, self.member, self.chunk, region=self.basin
+            self.startdate, self.member, self.chunk, region=self.basins
         )
 
     def compute(self):
         """Run the diagnostic"""
         self._fix_file_metadata()
+        mesh = Nemo('mesh_hgr.nc', 'mask_regions.nc')
+        areacello = mesh.get_areacello(cell_point=self.grid_point)
+        lats = mesh.get_mesh_var('latitude', dtype=np.float32)
         data = self._load_data()
-        self._meand_3d_variable(data)
-
-    def _meand_3d_variable(self, data):
-        e1 = self._try_load_cube(1)
-        e2 = self._try_load_cube(2)
-        mask = np.squeeze(Utils.get_mask(self.basin, True))
-        mask = e1.data * e2.data * mask
-        if len(mask.shape) == 2:
-            data.add_aux_coord(
-                iris.coords.AuxCoord(mask.data, long_name='mask'),
-                data.coord_dims('latitude')
-            )
-        else:
-            data.add_aux_coord(
-                iris.coords.AuxCoord(mask.data, long_name='mask'),
-                data.coord_dims('depth') + data.coord_dims('latitude')
-            )
-
-        @numba.njit()
-        def get_zonal_mean(variable, weight, latitude):
-            total = np.zeros(180, np.float64)
-            weights = np.zeros(180, np.float64)
-            for i in range(variable.shape[0]):
-                for j in range(variable.shape[1]):
-                    if weight[i, j] == 0:
-                        continue
-                    bin_value = int(round(latitude[i, j]) + 90)
-                    weights[bin_value] += weight[i, j]
-                    total[bin_value] += variable[i, j] * weight[i, j]
-            return total / weights
-
-        mean = iris.cube.CubeList()
-        lat_coord = None
-        for map_slice in data.slices_over('time'):
-            # Force data loading
-            map_slice.data
-            surface_cubes = iris.cube.CubeList()
-            for surface_slice in map_slice.slices_over('depth'):
-                value = get_zonal_mean(
-                    surface_slice.data,
-                    surface_slice.coord('mask').points,
-                    surface_slice.coord('latitude').points,
-                )
-                cube = Cube(value)
-                cube.add_aux_coord(surface_slice.coord('depth'))
-                if lat_coord is None:
-                    lat_coord = surface_slice.coord('latitude')
-                    lat_coord = lat_coord.copy(
-                        np.arange(-90, 90, dtype=np.float32)
-                    )
-                    lat_coord = iris.coords.DimCoord.from_coord(lat_coord)
-                cube.add_dim_coord(lat_coord, 0)
-                surface_cubes.append(cube)
-            time_cube = surface_cubes.merge_cube()
-            time_cube.add_aux_coord(map_slice.coord('time'))
-            mean.append(time_cube)
-        cube = mean.merge_cube()
-        cube.var_name = 'result'
-        cube.units = data.units
-        cube.attributes = data.attributes
-        temp = TempFile.get()
-        iris.save(cube, temp)
-        self.zonal_mean.set_local_file(temp, rename_var='result', region=self.basin)
+        masks = {}
+        self.basins.sort()
+        for basin in self.basins:
+            masks[basin] = Utils.get_mask(basin)
+        area_basin = zonmean.get_basin_area(areacello, masks)
+        zonalmean = zonmean.compute_zonmean(data.data, lats, area_basin)
+        self._save_result(zonalmean, data)
+
 
     def _try_load_cube(self, number):
         try:
@@ -231,3 +184,40 @@ class ZonalMean(Diagnostic):
         var.coordinates = coordinates
         handler.close()
         return has_levels
+
+
+    def _save_result(self, result, data):
+        final_name = '{0}zonal'.format(self.variable)
+        temp = TempFile.get()
+        handler_source = Utils.open_cdf(self.variable_file.local_file)
+        handler_temp = Utils.open_cdf(temp, 'w')
+        Utils.copy_variable(handler_source, handler_temp, 'time', True, True)
+        handler_temp.createDimension('region', len(result))
+        handler_temp.createDimension('region_length', 50)
+        handler_temp.createDimension('lev', data.shape[1])
+        handler_temp.createDimension('lat', 180)
+
+        var_lat = handler_temp.createVariable('lat', float, 'lat')
+        var_lat[...] = np.arange(-90, 90, dtype=np.float32)
+        var_lat.units = 'degrees_north'
+        var_lat.long_name = 'latitude'
+        var_lat.standard_name = 'latitude'
+
+        var_level = handler_temp.createVariable('lev', float, 'lev')
+        var_level[...] = data.coord('depth').points
+        var_level.units = 'm'
+        var_level.axis = 'Z'
+        var_level.positive = 'down'
+        var_level.long_name = 'ocean depth coordinate'
+        var_level.standard_name = 'depth'
+        var_region = handler_temp.createVariable('region', 'S1',
+                                                 ('region', 'region_length'))
+        var = handler_temp.createVariable('{0}zonal'.format(self.variable),
+                                          float, ('time', 'lev', 'lat','region',),)
+        var.units = '{0}'.format(data.units)
+        for i, basin in enumerate(result):
+            var_region[i, ...] = netCDF4.stringtoarr(str(basin), 50)
+            var[..., i] = result[basin]
+        handler_temp.close()
+        self.declared[final_name].set_local_file(temp, diagnostic=self)
+
-- 
GitLab


From 85af8ce83518b728ae6b670ff95986a766502eb3 Mon Sep 17 00:00:00 2001
From: sloosvel <saskia.loosveldt@bsc.es>
Date: Fri, 12 Jul 2019 15:36:14 +0200
Subject: [PATCH 25/30] Add hfcorr and wfcorr to default.csv table

---
 earthdiagnostics/cmor_tables/default.csv | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/earthdiagnostics/cmor_tables/default.csv b/earthdiagnostics/cmor_tables/default.csv
index 7133c04f..67cc35ff 100644
--- a/earthdiagnostics/cmor_tables/default.csv
+++ b/earthdiagnostics/cmor_tables/default.csv
@@ -1,5 +1,5 @@
 Variable,Shortname,Name,Long name,Domain,Basin,Units,Valid min,Valid max,Grid,Tables
-hfxout,hfxout,,Ocean surface heat flux,ocean,,W m-2,,,,,
+hfxout,hfxout,,Ocean surface heat flux,ocean,,W m-2,,,,
 iiceages:siage:iice_otd,ageice,age_of_sea_ice,Age of sea ice,seaIce,,,,,,
 al,al,surface_albedo,Albedo,atmos,,,,,,
 bgfrcsal,bgfrcsal,change_over_time_in_heat_content_from_forcing,Change over time in salt content from forcing,ocean,,,,,,
@@ -347,4 +347,6 @@ zqlw,rlntds,surface_net_downward_longwave_flux,Surface Net Downward Longwave Rad
 var78,tclw,total_column_liquid_water,Total column liquid water,atmos,,kg m-2,,,,
 var79,tciw,total_column_ice_water,Total column ice water,atmos,,kg m-2,,,,
 rho,rhopoto,sea_water_potential_density,Sea Water Potential Density,ocean,,kg m-3,,,,
-hc300,hc300,,Heat Content to 300m depth,ocean,J m-2,,,,
\ No newline at end of file
+hc300,hc300,,Heat Content to 300m depth,ocean,,J m-2,,,,
+hfcorr,hfcorr,heat_flux_into_sea_water_due_to_newtonian_relaxation,heat_flux_correction,ocean,,W m-2,,,,
+wfcorr,wfcorr,water_flux_out_of_sea_water_due_to_newtonian_relaxation,water_flux_correction,ocean,,kg m-2 s-1,,,,
-- 
GitLab


From 1dd12e0689ccb80099b75a33f9b3507d9c1f603e Mon Sep 17 00:00:00 2001
From: sloosvel <saskia.loosveldt@bsc.es>
Date: Fri, 12 Jul 2019 15:56:44 +0200
Subject: [PATCH 26/30] Add units

---
 earthdiagnostics/cmor_tables/default.csv | 1 +
 1 file changed, 1 insertion(+)

diff --git a/earthdiagnostics/cmor_tables/default.csv b/earthdiagnostics/cmor_tables/default.csv
index 67cc35ff..eaa24064 100644
--- a/earthdiagnostics/cmor_tables/default.csv
+++ b/earthdiagnostics/cmor_tables/default.csv
@@ -350,3 +350,4 @@ rho,rhopoto,sea_water_potential_density,Sea Water Potential Density,ocean,,kg m-
 hc300,hc300,,Heat Content to 300m depth,ocean,,J m-2,,,,
 hfcorr,hfcorr,heat_flux_into_sea_water_due_to_newtonian_relaxation,heat_flux_correction,ocean,,W m-2,,,,
 wfcorr,wfcorr,water_flux_out_of_sea_water_due_to_newtonian_relaxation,water_flux_correction,ocean,,kg m-2 s-1,,,,
+,,,,ocean,,kg m-2 s-1,,,,
-- 
GitLab


From 8aae6828e9752792062df14fab487c08d770a54d Mon Sep 17 00:00:00 2001
From: sloosvel <saskia.loosveldt@bsc.es>
Date: Thu, 18 Jul 2019 13:17:23 +0200
Subject: [PATCH 27/30] Add hfxin to cmor table

---
 earthdiagnostics/cmor_tables/default.csv | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/earthdiagnostics/cmor_tables/default.csv b/earthdiagnostics/cmor_tables/default.csv
index eaa24064..2fe50114 100644
--- a/earthdiagnostics/cmor_tables/default.csv
+++ b/earthdiagnostics/cmor_tables/default.csv
@@ -350,4 +350,4 @@ rho,rhopoto,sea_water_potential_density,Sea Water Potential Density,ocean,,kg m-
 hc300,hc300,,Heat Content to 300m depth,ocean,,J m-2,,,,
 hfcorr,hfcorr,heat_flux_into_sea_water_due_to_newtonian_relaxation,heat_flux_correction,ocean,,W m-2,,,,
 wfcorr,wfcorr,water_flux_out_of_sea_water_due_to_newtonian_relaxation,water_flux_correction,ocean,,kg m-2 s-1,,,,
-,,,,ocean,,kg m-2 s-1,,,,
+hfxin,hfxin,,total heat fluxes at the ice/ocean surface,ocean,,W m-2,,,,
-- 
GitLab


From f16f3e4bc88b1b407d00d9ed211400e7b1e06cee Mon Sep 17 00:00:00 2001
From: Javier Vegas-Regidor <javier.vegas@bsc.es>
Date: Mon, 22 Jul 2019 14:13:38 +0200
Subject: [PATCH 28/30] Fix startdate appending for DCPP

---
 .gitlab-ci.yml                              |  1 +
 earthdiagnostics/config.py                  |  1 +
 earthdiagnostics/data_convention.py         |  2 +
 earthdiagnostics/ocean/moc.py               |  8 ++--
 earthdiagnostics/ocean/regionmean.py        | 41 ++++++++----------
 earthdiagnostics/ocean/regionsum.py         | 13 ++----
 earthdiagnostics/ocean/siasiesiv.py         | 16 +++----
 test/unit/data_convention/test_primavera.py | 14 ++++++
 test/unit/ocean/test_region_mean.py         | 48 +++++++++++++--------
 9 files changed, 81 insertions(+), 63 deletions(-)

diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index e7a4660d..f3e8d8c5 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -40,6 +40,7 @@ test_python3:
 report_codacy:
   stage: report
   script:
+    - '[ -z "$CODACY_PROJECT_TOKEN" ] && exit 0'
     - source activate earthdiagnostics3
     - pip install codacy-coverage --upgrade
     - python-codacy-coverage -r test/report/python3/coverage.xml
diff --git a/earthdiagnostics/config.py b/earthdiagnostics/config.py
index 6bdf281a..838bb223 100644
--- a/earthdiagnostics/config.py
+++ b/earthdiagnostics/config.py
@@ -250,6 +250,7 @@ class CMORConfig(object):
         self.activity = parser.get_option('CMOR', 'ACTIVITY', 'CMIP')
         self.min_cmorized_vars = parser.get_int_option('CMOR', 'MIN_CMORIZED_VARS', 10)
         self.append_startdate = parser.get_bool_option('CMOR', 'APPEND_STARTDATE', False)
+        self.append_startdate_year_only = parser.get_bool_option('CMOR', 'APPEND_STARTDATE_YEAR_ONLY', False)
 
         vars_string = parser.get_option('CMOR', 'VARIABLE_LIST', '')
         self.var_manager = var_manager
diff --git a/earthdiagnostics/data_convention.py b/earthdiagnostics/data_convention.py
index 9e855dd1..1c54d35c 100644
--- a/earthdiagnostics/data_convention.py
+++ b/earthdiagnostics/data_convention.py
@@ -585,6 +585,8 @@ class Cmor3Convention(DataConvention):
             else:
                 grid = self.config.cmor.default_atmos_grid
         if self.config.cmor.append_startdate:
+            if self.config.cmor.append_startdate_year_only:
+                startdate = startdate[0:4]
             subexp_id = 's{}-'.format(startdate)
         else:
             subexp_id = ''
diff --git a/earthdiagnostics/ocean/moc.py b/earthdiagnostics/ocean/moc.py
index 7732726d..a294b824 100644
--- a/earthdiagnostics/ocean/moc.py
+++ b/earthdiagnostics/ocean/moc.py
@@ -67,8 +67,10 @@ class Moc(Diagnostic):
     def __eq__(self, other):
         if self._different_type(other):
             return False
-        return self.startdate == other.startdate and \
-               self.member == other.member and self.chunk == other.chunk
+        return (
+            self.startdate == other.startdate and
+            self.member == other.member and self.chunk == other.chunk
+        )
 
     @classmethod
     def generate_jobs(cls, diags, options):
@@ -136,7 +138,6 @@ class Moc(Diagnostic):
         del vo, e1v, e3v
         self._save_result(moc_results, mesh)
 
-
     def _save_result(self, result, mesh):
         temp = TempFile.get()
         handler_source = Utils.open_cdf(self.variable_file.local_file)
@@ -151,7 +152,6 @@ class Moc(Diagnostic):
         handler_temp.createDimension('region', len(result))
         handler_temp.createDimension('region_length', 50)
 
-
         var_region = handler_temp.createVariable('region', 'S1',
                                                  ('region', 'region_length'))
 
diff --git a/earthdiagnostics/ocean/regionmean.py b/earthdiagnostics/ocean/regionmean.py
index 792793a5..c70c9500 100644
--- a/earthdiagnostics/ocean/regionmean.py
+++ b/earthdiagnostics/ocean/regionmean.py
@@ -22,7 +22,6 @@ import diagonals.regmean as regmean
 from diagonals.mesh_helpers.nemo import Nemo
 
 
-
 class RegionMean(Diagnostic):
     """
     Computes the mean value of the field (3D, weighted).
@@ -79,10 +78,10 @@ class RegionMean(Diagnostic):
             self.box == other.box and self.variable == other.variable
 
     def __str__(self):
-        return 'Region mean Startdate: {0.startdate} Member: {0.member}' \
-               'Chunk: {0.chunk} Variable: {0.variable} ' \
-               'Box: {0.box} Save 3D: {0.save3d} Save variance: {0.variance}' \
-               'Grid point: {0.grid_point}'.format(self)
+        return ('Region mean Startdate: {0.startdate} Member: {0.member} '
+                'Chunk: {0.chunk} Variable: {0.variable} '
+                'Box: {0.box} Save 3D: {0.save3d} Save variance: {0.variance} '
+                'Grid point: {0.grid_point}'.format(self))
 
     def __hash__(self):
         return hash(str(self))
@@ -101,8 +100,7 @@ class RegionMean(Diagnostic):
         options_available = (DiagnosticDomainOption(),
                              DiagnosticVariableOption(diags.data_manager.config.var_manager),
                              DiagnosticOption('grid_point', 'T'),
-                             DiagnosticBasinListOption('basins',
-                                                       Basins().Global),
+                             DiagnosticBasinListOption('basins', 'global'),
                              DiagnosticIntOption('min_depth', -1),
                              DiagnosticIntOption('max_depth', -1),
                              DiagnosticIntOption('min_lat', -1),
@@ -117,17 +115,18 @@ class RegionMean(Diagnostic):
         box = Box()
         box.min_depth = options['min_depth']
         box.max_depth = options['max_depth']
-        box.min_lat = options['min_lat']
-        box.max_lat = options['max_lat']
-        box.min_lon = options['min_lon']
-        box.max_lon = options['max_lon']
+        if options['min_lat'] != -1:
+            box.min_lat = options['min_lat']
+            box.max_lat = options['max_lat']
+        if options['min_lon'] != -1 or options['max_lon'] != -1:
+            box.min_lon = options['min_lon']
+            box.max_lon = options['max_lon']
 
         basins = options['basins']
         if not basins:
             Log.error('Basins not recognized')
             return()
 
-
         job_list = list()
         for startdate, member, chunk in diags.config.experiment.get_chunk_list():
             job = RegionMean(diags.data_manager, startdate, member, chunk,
@@ -165,10 +164,13 @@ class RegionMean(Diagnostic):
             masks[basin] = Utils.get_mask(basin)
 
         mesh = Nemo('mesh_hgr.nc', 'mask_regions.nc')
-        if self.box.min_lat is not -1 and self.box.max_lat is not -1 and \
-            self.box.min_lon is not -1 and self.box.max_lat is not -1:
-            name = '{0}_{1}'.format(Box.get_lat_str(self.box),
-                                    Box.get_lon_str(self.box))
+        if (
+            self.box.min_lat is not -1 and self.box.max_lat is not -1 and
+            self.box.min_lon is not -1 and self.box.max_lat is not -1
+        ):
+            name = '{0}_{1}'.format(
+                Box.get_lat_str(self.box), Box.get_lon_str(self.box)
+            )
 
             masks[name] = mesh.get_region_mask(self.box.min_lat,
                                                self.box.max_lat,
@@ -184,7 +186,6 @@ class RegionMean(Diagnostic):
         mean = regmean.compute_regmean_2D(data.data, masks, areacello)
         self._save_result_2D('mean', mean, data)
 
-
     def _meand_3d_variable(self, data, mesh, masks):
         areacello = mesh.get_areacello(cell_point=self.grid_point)
         e3 = self._try_load_cube(3)
@@ -222,7 +223,6 @@ class RegionMean(Diagnostic):
                                                     var_name='j'), dims - 2)
         return cube
 
-
     def _load_data(self):
         coords = []
         handler = Utils.open_cdf(self.variable_file.local_file)
@@ -279,7 +279,6 @@ class RegionMean(Diagnostic):
                                                        box=box_save,
                                                        region=self.basins)
 
-
     def _save_result_2D(self, var, result, data):
         final_name = '{1}{0}'.format(var, self.variable)
         temp = TempFile.get()
@@ -299,7 +298,6 @@ class RegionMean(Diagnostic):
         handler_temp.close()
         self.declared[final_name].set_local_file(temp, diagnostic=self)
 
-
     def _save_result_3D(self, var, result, data):
         final_name = '{1}3d{0}'.format(var, self.variable)
         temp = TempFile.get()
@@ -321,11 +319,10 @@ class RegionMean(Diagnostic):
 
         var = handler_temp.createVariable('{1}3d{0}'.format(var, self.variable),
                                           float, ('time', 'lev', 'region',),)
-                                          
+
         var.units = '{0}'.format(data.units)
         for i, basin in enumerate(result):
             var_region[i, ...] = netCDF4.stringtoarr(str(basin), 50)
             var[..., i] = result[basin]
         handler_temp.close()
         self.declared[final_name].set_local_file(temp, diagnostic=self)
-
diff --git a/earthdiagnostics/ocean/regionsum.py b/earthdiagnostics/ocean/regionsum.py
index f9207e59..99b46ce1 100644
--- a/earthdiagnostics/ocean/regionsum.py
+++ b/earthdiagnostics/ocean/regionsum.py
@@ -24,6 +24,7 @@ from earthdiagnostics.utils import Utils, TempFile
 import diagonals.regsum as regsum
 from diagonals.mesh_helpers.nemo import Nemo
 
+
 class RegionSum(Diagnostic):
     """
     Computes the sum of the field (3D, weighted).
@@ -88,11 +89,9 @@ class RegionSum(Diagnostic):
                'Grid point: {0.grid_point} Box: {0.box} Save 3D: {0.save3d}' \
                'Original grid: {0.grid} Basin: {0.basins}'.format(self)
 
-
     def __hash__(self):
         return hash(str(self))
 
-
     @classmethod
     def generate_jobs(cls, diags, options):
         """
@@ -167,7 +166,8 @@ class RegionSum(Diagnostic):
             masks[basin] = Utils.get_mask(basin)
 
         mesh = Nemo('mesh_hgr.nc', 'mask_regions.nc')
-        if self.box.min_lat is not -1 and self.box.max_lat is not -1 and self.box.min_lon is not -1 and self.box.max_lat is not -1:
+        if self.box.min_lat is not -1 and self.box.max_lat is not -1 and \
+           self.box.min_lon is not -1 and self.box.max_lat is not -1:
             name = '{0}_{1}'.format(Box.get_lat_str(self.box), Box.get_lon_str(self.box))
 
             masks[name] = mesh.get_region_mask(self.box.min_lat,
@@ -179,13 +179,11 @@ class RegionSum(Diagnostic):
         else:
             self._sum_2d_var(data, mesh, masks)
 
-
     def _sum_2d_var(self, data, mesh, masks):
         areacello = mesh.get_areacello(cell_point=self.grid_point)
         varsum = regsum.compute_regsum_2D(data.data, masks, areacello)
         self._save_result_2D('sum', varsum, data)
 
-
     def _sum_3d_variable(self, data, mesh, masks):
         areacello = mesh.get_areacello(cell_point=self.grid_point)
         tmask = iris.load_cube('mesh_hgr.nc', '{0}mask'.format(self.grid_point))
@@ -206,7 +204,6 @@ class RegionSum(Diagnostic):
                                                     volcello, tmask)
             self._save_result_3D('sum', varsum3d, data)
 
-
     def _try_load_cube(self, number):
         try:
             cube = iris.load_cube('mesh_hgr.nc', 'e{0}{1}'.format(number, self.grid_point))
@@ -239,7 +236,6 @@ class RegionSum(Diagnostic):
         data = iris.load_cube(self.variable_file.local_file)
         return self._rename_depth(data)
 
-
     def _rename_depth(self, data):
         for coord_name in ('model_level_number', 'Vertical T levels', 'lev'):
             if data.coords(coord_name):
@@ -249,7 +245,6 @@ class RegionSum(Diagnostic):
                 break
         return data
 
-
     def _fix_file_metadata(self):
         handler = Utils.open_cdf(self.variable_file.local_file)
         var = handler.variables[self.variable]
@@ -314,7 +309,6 @@ class RegionSum(Diagnostic):
         self.declared[final_name] = self.declare_chunk(ModelingRealms.ocean, final_name, self.startdate, self.member,
                                                        self.chunk, box=box_save, region=self.basins, grid=self.grid)
 
-
     def _save_result_2D(self, var, result, data):
         final_name = '{1}{0}'.format(var, self.variable)
         temp = TempFile.get()
@@ -334,7 +328,6 @@ class RegionSum(Diagnostic):
         handler_temp.close()
         self.declared[final_name].set_local_file(temp, diagnostic=self)
 
-
     def _save_result_3D(self, var, result, data):
         final_name = '{1}3d{0}'.format(var, self.variable)
         temp = TempFile.get()
diff --git a/earthdiagnostics/ocean/siasiesiv.py b/earthdiagnostics/ocean/siasiesiv.py
index 5cf4a905..f623f483 100644
--- a/earthdiagnostics/ocean/siasiesiv.py
+++ b/earthdiagnostics/ocean/siasiesiv.py
@@ -150,7 +150,8 @@ class Siasiesiv(Diagnostic):
         mesh = Nemo('mesh_hgr.nc', 'mask_regions.nc')
         areacello = mesh.get_areacello(cell_point='T')
         gphit = mesh.get_grid_latitude(cell_point='T')
-        self.results['siextentn'], self.results['siextents'], self.results['siarean'], self.results['siareas'] = siasie.compute(gphit, areacello, sic_slices, self.masks)
+        self.results['siextentn'], self.results['siextents'], self.results['siarean'], self.results['siareas'] = \
+            siasie.compute(gphit, areacello, sic_slices, self.masks)
 
         self.save()
 
@@ -165,7 +166,6 @@ class Siasiesiv(Diagnostic):
             ' '.join(set(coordinates) | add_coordinates)
         handler.close()
 
-
     def save(self):
         for var in self.results.keys():
             res = self.results[var]
@@ -175,15 +175,15 @@ class Siasiesiv(Diagnostic):
             Utils.copy_variable(handler_source, handler_temp, 'time', True, True)
             handler_temp.createDimension('region', len(self.masks))
             handler_temp.createDimension('region_length', 50)
-            var_region = handler_temp.createVariable('region', 'S1',
-                                                 ('region', 'region_length'))
-            var_res = handler_temp.createVariable('{0}'.format(var), float,
-                                                  ('time', 'region',))
+            var_region = handler_temp.createVariable(
+                'region', 'S1', ('region', 'region_length')
+            )
+            var_res = handler_temp.createVariable(
+                '{0}'.format(var), float, ('time', 'region',)
+            )
             var_res.units = 'm^2'
             for i, basin in enumerate(self.masks):
                 var_region[i, ...] = netCDF4.stringtoarr(str(basin), 50)
                 var_res[..., i] = res[i, ...]
             handler_temp.close()
             self.generated[var].set_local_file(temp, diagnostic=self)
-
-
diff --git a/test/unit/data_convention/test_primavera.py b/test/unit/data_convention/test_primavera.py
index f2e6739b..52a06faa 100644
--- a/test/unit/data_convention/test_primavera.py
+++ b/test/unit/data_convention/test_primavera.py
@@ -32,6 +32,7 @@ class TestPrimaveraConvention(TestCase):
 
         self.config.cmor.activity = 'activity'
         self.config.cmor.append_startdate = False
+        self.config.cmor.append_startdate_year_only = False
         self.config.cmor.initialization_number = 1
         self.config.cmor.version = 'version'
         self.config.cmor.default_ocean_grid = 'ocean_grid'
@@ -148,6 +149,19 @@ class TestPrimaveraConvention(TestCase):
         self.assertEqual(file_path,
                          'var_Omon_model_experiment_name_s19900101-r2i1p1f1_ocean_grid_199001-199001.nc')
 
+    def test_get_filename_with_startdate_year_only(self):
+        """Test get_filename with startdate"""
+        self.config.cmor.append_startdate = True
+        self.config.cmor.append_startdate_year_only = True
+        cmor_var = Mock()
+        omon = Mock()
+        omon.name = 'Omon'
+        cmor_var.get_table.return_value = omon
+        file_path = self.convention.get_file_name('19900101', 1, ModelingRealms.ocean, 'var', cmor_var,
+                                                  Frequencies.monthly, 1, None, None, None)
+        self.assertEqual(file_path,
+                         'var_Omon_model_experiment_name_s1990-r2i1p1f1_ocean_grid_199001-199001.nc')
+
     def test_get_filename_no_cmor_var(self):
         """Test get_filename not passing cmor_var"""
         cmor_var = Mock()
diff --git a/test/unit/ocean/test_region_mean.py b/test/unit/ocean/test_region_mean.py
index 966ed536..35f8a09c 100644
--- a/test/unit/ocean/test_region_mean.py
+++ b/test/unit/ocean/test_region_mean.py
@@ -38,23 +38,23 @@ class TestRegionMean(TestCase):
         jobs = RegionMean.generate_jobs(self.diags, ['diagnostic', 'ocean', 'var'])
         self.assertEqual(len(jobs), 2)
         self.assertEqual(jobs[0], RegionMean(self.data_manager, '20010101', 0, 0, ModelingRealms.ocean, 'var',
-                                             box, True, 'weights', False, Basins().Global))
+                                             box, True, False, Basins().Global, 't'))
         self.assertEqual(jobs[1], RegionMean(self.data_manager, '20010101', 0, 1, ModelingRealms.ocean, 'var',
-                                             box, True, 'weights', False, Basins().Global))
+                                             box, True, False, Basins().Global, 't'))
 
         jobs = RegionMean.generate_jobs(self.diags, ['diagnostic', 'ocean', 'var', 'U'])
         self.assertEqual(len(jobs), 2)
         self.assertEqual(jobs[0], RegionMean(self.data_manager, '20010101', 0, 0, ModelingRealms.ocean, 'var',
-                                             box, True, 'weights', False, Basins().Global))
+                                             box, True, False, Basins().Global, 'u'))
         self.assertEqual(jobs[1], RegionMean(self.data_manager, '20010101', 0, 1, ModelingRealms.ocean, 'var',
-                                             box, True, 'weights', False, Basins().Global))
+                                             box, True, False, Basins().Global, 'u'))
 
         jobs = RegionMean.generate_jobs(self.diags, ['diagnostic', 'ocean', 'var', 'U', 'global'])
         self.assertEqual(len(jobs), 2)
         self.assertEqual(jobs[0], RegionMean(self.data_manager, '20010101', 0, 0, ModelingRealms.ocean, 'var',
-                                             box, True, 'weights', False, Basins().Global))
+                                             box, True, False, Basins().Global, 'u'))
         self.assertEqual(jobs[1], RegionMean(self.data_manager, '20010101', 0, 1, ModelingRealms.ocean, 'var',
-                                             box, True, 'weights', False, Basins().Global))
+                                             box, True, False, Basins().Global, 'u'))
 
         box = Box()
         box.min_depth = 1.0
@@ -63,27 +63,34 @@ class TestRegionMean(TestCase):
         jobs = RegionMean.generate_jobs(self.diags, ['diagnostic', 'ocean', 'var', 'U', 'global', '1', '10'])
         self.assertEqual(len(jobs), 2)
         self.assertEqual(jobs[0], RegionMean(self.data_manager, '20010101', 0, 0, ModelingRealms.ocean, 'var',
-                                             box, True, 'weights', False, Basins().Global))
+                                             box, True, False, Basins().Global, 'u'))
         self.assertEqual(jobs[1], RegionMean(self.data_manager, '20010101', 0, 1, ModelingRealms.ocean, 'var',
-                                             box, True, 'weights', False, Basins().Global))
+                                             box, True, False, Basins().Global, 'u'))
 
-        jobs = RegionMean.generate_jobs(self.diags, ['diagnostic', 'ocean', 'var', 'U', 'global', '1', '10', 'false'])
+        jobs = RegionMean.generate_jobs(
+            self.diags,
+            ['diagnostic', 'ocean', 'var', 'U', 'global', '1', '10', '', '', '', '', 'false']
+        )
         self.assertEqual(len(jobs), 2)
         self.assertEqual(jobs[0], RegionMean(self.data_manager, '20010101', 0, 0, ModelingRealms.ocean, 'var',
-                                             box, False, Basins().Global, False, ''))
+                                             box, False, False, Basins().Global, 'u'))
         self.assertEqual(jobs[1], RegionMean(self.data_manager, '20010101', 0, 1, ModelingRealms.ocean, 'var',
-                                             box, False, Basins().Global, False, ''))
+                                             box, False, False, Basins().Global, 'u'))
 
-        jobs = RegionMean.generate_jobs(self.diags, ['diagnostic', 'ocean', 'var', 'U', 'global', '1', '10', 'false',
-                                                     'True'])
+        jobs = RegionMean.generate_jobs(
+            self.diags,
+            ['diagnostic', 'ocean', 'var', 'U', 'global', '1', '10', '', '', '', '', 'false', 'True']
+        )
         self.assertEqual(len(jobs), 2)
         self.assertEqual(jobs[0], RegionMean(self.data_manager, '20010101', 0, 0, ModelingRealms.ocean, 'var',
-                                             box, False, Basins().Global, True, ''))
+                                             box, False, True,  Basins().Global, 'u'))
         self.assertEqual(jobs[1], RegionMean(self.data_manager, '20010101', 0, 1, ModelingRealms.ocean, 'var',
-                                             box, False, Basins().Global, True, ''))
+                                             box, False, True, Basins().Global, 'u'))
 
-        jobs = RegionMean.generate_jobs(self.diags, ['diagnostic', 'ocean', 'var', 'U', 'global', '1', '10', 'false',
-                                                     'True', 'grid'])
+        jobs = RegionMean.generate_jobs(
+            self.diags,
+            ['diagnostic', 'ocean', 'var', 'U', 'global', '1', '10', '', '', '', '', 'false', 'True', 'grid']
+        )
         self.assertEqual(len(jobs), 2)
         self.assertEqual(jobs[0], RegionMean(self.data_manager, '20010101', 0, 0, ModelingRealms.ocean, 'var',
                                              box, False, Basins().Global, True, 'grid'))
@@ -94,8 +101,11 @@ class TestRegionMean(TestCase):
             RegionMean.generate_jobs(self.diags, ['diagnostic'])
 
         with self.assertRaises(DiagnosticOptionError):
-            RegionMean.generate_jobs(self.diags, ['diagnostic', 'ocean', 'var', 'U', 'global', '1', '10', 'false',
-                                                  'True', 'grid', 'extra'])
+            RegionMean.generate_jobs(
+                self.diags,
+                ['diagnostic', 'ocean', 'var', 'U', 'global', '1', '10', '', '', '', '', 'false',
+                 'True', 'grid', 'extra']
+            )
 
     def test_str(self):
         box = Box()
-- 
GitLab


From 78192263c681828d18dae565e3b1b38bad411c59 Mon Sep 17 00:00:00 2001
From: sloosvel <saskia.loosveldt@bsc.es>
Date: Fri, 26 Jul 2019 14:54:16 +0200
Subject: [PATCH 29/30] Convert units to % to work with diagonals

---
 earthdiagnostics/ocean/siasiesiv.py | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/earthdiagnostics/ocean/siasiesiv.py b/earthdiagnostics/ocean/siasiesiv.py
index cc65228e..94c65d3c 100644
--- a/earthdiagnostics/ocean/siasiesiv.py
+++ b/earthdiagnostics/ocean/siasiesiv.py
@@ -143,8 +143,10 @@ class Siasiesiv(Diagnostic):
             self.sic.local_file, self.sic_varname
         )
         sic = iris.load_cube(self.sic.local_file)
-        if sic.units.origin == '%' and sic.data.max() < 2:
-            sic.units = '1.0'
+        #if sic.units.origin == '%' and sic.data.max() < 2:
+        #    sic.units = '1.0'
+        if sic.units.origin is not '%': 
+            sic.convert_units('%')
 
         sic_slices = []
         for sic_data in sic.slices_over('time'):
-- 
GitLab


From 598e79d2c1c3cfa502f4204c9cf22e462f5d596f Mon Sep 17 00:00:00 2001
From: sloosvel <saskia.loosveldt@bsc.es>
Date: Fri, 6 Sep 2019 16:44:40 +0200
Subject: [PATCH 30/30] Adapt to last version of diagonals

---
 earthdiagnostics/ocean/heatcontentlayer.py | 142 ++++++++++++---------
 1 file changed, 79 insertions(+), 63 deletions(-)

diff --git a/earthdiagnostics/ocean/heatcontentlayer.py b/earthdiagnostics/ocean/heatcontentlayer.py
index 023909a0..bf0cc580 100644
--- a/earthdiagnostics/ocean/heatcontentlayer.py
+++ b/earthdiagnostics/ocean/heatcontentlayer.py
@@ -4,10 +4,19 @@ import numpy as np
 
 from earthdiagnostics.box import Box
 from earthdiagnostics.constants import Basins
-from earthdiagnostics.diagnostic import Diagnostic, DiagnosticIntOption, DiagnosticBasinOption
+from earthdiagnostics.diagnostic import Diagnostic, DiagnosticIntOption, DiagnosticBasinListOption
 from earthdiagnostics.modelingrealm import ModelingRealms
 from earthdiagnostics.utils import Utils, TempFile
+from bscearth.utils.log import Log
 
+import diagonals
+from diagonals import ohc
+from diagonals.mesh_helpers.nemo import Nemo
+
+import iris
+import iris.coords
+
+import netCDF4
 
 class HeatContentLayer(Diagnostic):
     """
@@ -36,15 +45,17 @@ class HeatContentLayer(Diagnostic):
     alias = 'ohclayer'
     "Diagnostic alias for the configuration file"
 
-    def __init__(self, data_manager, startdate, member, chunk, box, weight, min_level, max_level, data_convention):
+    def __init__(self, data_manager, startdate, member, chunk, box, areas,
+                 weight, layers, basins, data_convention):
         Diagnostic.__init__(self, data_manager)
         self.startdate = startdate
         self.member = member
         self.chunk = chunk
         self.box = box
+        self.areas = areas
         self.weight = weight
-        self.min_level = min_level
-        self.max_level = max_level
+        self.layers = layers
+        self.basins = basins
         self.required_vars = ['so', 'mlotst']
         self.generated_vars = ['scvertsum']
         self.data_convention = data_convention
@@ -65,62 +76,47 @@ class HeatContentLayer(Diagnostic):
         """
         options_available = (DiagnosticIntOption('min_depth'),
                              DiagnosticIntOption('max_depth'),
-                             DiagnosticBasinOption('basin', Basins().Global))
+                             DiagnosticBasinListOption('basins',
+                             Basins().Global))
         options = cls.process_options(options, options_available)
 
         box = Box(True)
         box.min_depth = options['min_depth']
         box.max_depth = options['max_depth']
+        basins = options['basins']
+        if not basins:
+            Log.error('Basins not recognized')
+            return ()
+
+        mesh = Nemo('mesh_hgr.nc', 'mask_regions.nc')
+
+        layers = ((box.min_depth, box.max_depth),)
+
+        masks = {}
+        basins.sort()
+        for basin in basins:
+            masks[basin] = Utils.get_mask(basin)
+
+        areacello = mesh.get_areacello()
+        areas = ohc.get_basin_area(areacello, masks)
+
         job_list = list()
 
-        max_level, min_level, weight = cls._compute_weights(box)
+        e3t = mesh.get_k_length()
+        mask = mesh.get_landsea_mask()
+        depth = mesh.get_depth(cell_point='W')
+        weight = ohc.get_weights(layers, mask, e3t, depth)
+
+        del mask, depth, e3t
 
         for startdate, member, chunk in diags.config.experiment.get_chunk_list():
             job_list.append(HeatContentLayer(
                 diags.data_manager, startdate, member, chunk, box,
-                weight, min_level, max_level,
+                areas, weight, layers, basins,
                 diags.config.data_convention
             ))
         return job_list
 
-    @classmethod
-    def _compute_weights(cls, box):
-        depth, e3t, mask = cls._get_mesh_info()
-
-        def calculate_weight(e3t_point, depth_point, mask_point):
-            """Calculate the weight for each cell"""
-            if not mask_point:
-                return 0
-            top = depth_point
-            bottom = top + e3t_point
-            if bottom < box.min_depth or top > box.max_depth:
-                return 0
-            else:
-                if top < box.min_depth:
-                    top = box.min_depth
-                if bottom > box.max_depth:
-                    bottom = box.max_depth
-
-                return (bottom - top) * 1020 * 4000
-
-        calc = np.vectorize(calculate_weight, otypes='f')
-        weight = calc(e3t, depth, mask)
-        max_level, min_level = cls._get_used_levels(weight)
-        weight = weight[:, min_level:max_level, :]
-        return max_level, min_level, weight
-
-    @classmethod
-    def _get_used_levels(cls, weight):
-        # Now we will reduce to the levels with any weight != 0 to avoid loading too much data on memory
-        levels = weight.shape[1]
-        min_level = 0
-        while min_level < levels and not weight[:, min_level, :].any():
-            min_level += 1
-        max_level = min_level
-        while max_level < (levels - 1) and weight[:, max_level + 1, :].any():
-            max_level += 1
-        return max_level, min_level
-
     @classmethod
     def _get_mesh_info(cls):
         handler = Utils.open_cdf('mesh_zgr.nc')
@@ -150,42 +146,62 @@ class HeatContentLayer(Diagnostic):
 
     def request_data(self):
         """Request data required by the diagnostic"""
-        self.thetao = self.request_chunk(ModelingRealms.ocean, 'thetao', self.startdate, self.member, self.chunk)
+        self.thetao = self.request_chunk(ModelingRealms.ocean, 'thetao',
+                                         self.startdate, self.member,
+                                         self.chunk)
 
     def declare_data_generated(self):
         """Declare data to be generated by the diagnostic"""
-        self.heatc = self.declare_chunk(ModelingRealms.ocean, 'heatc', self.startdate, self.member, self.chunk,
-                                        box=self.box)
+        self.heatc = self.declare_chunk(ModelingRealms.ocean, 'heatc',
+                                        self.startdate, self.member,
+                                        self.chunk, box=self.box)
+        self.heatcsum = self.declare_chunk(ModelingRealms.ocean, 'heatcsum',
+                                           self.startdate, self.member,
+                                           self.chunk, box=self.box)
 
     def compute(self):
         """Run the diagnostic"""
-        nco = Utils.nco()
+        
         thetao_file = TempFile.get()
         results = TempFile.get()
-
+        results1D = TempFile.get()
         Utils.copy_file(self.thetao.local_file, thetao_file)
 
         handler = Utils.open_cdf(thetao_file)
         Utils.convert_units(handler.variables['thetao'], 'K')
-        heatc_sl = np.sum(handler.variables['thetao'][:, self.min_level:self.max_level, :] * self.weight, 1)
+        heatc_sl, heatc_sl1D = ohc.compute(self.layers, self.weight,
+                                           handler.variables['thetao'][:],
+                                           self.areas)
         handler.sync()
         handler.renameVariable('thetao', 'heatc_sl')
-        handler.close()
 
-        nco.ncks(
-            input=thetao_file,
-            output=results,
-            options=(
-                '-O -v {0.lon_name},{0.lat_name},time'.format(self.data_convention),
-            )
-        )
-        Utils.rename_variables(results, {'x': 'i', 'y': 'j'}, False)
-        handler_results = Utils.open_cdf(results)
-        var = handler_results.createVariable('heatc', float, ('time', 'j', 'i'), fill_value=1.e20)
+        results = TempFile.get()
+        handler_results = Utils.open_cdf(results, 'w')
+        Utils.copy_variable(handler, handler_results, 'time', True, True)
+        Utils.copy_variable(handler, handler_results, 'i', True, True)
+        Utils.copy_variable(handler, handler_results, 'j', True, True)
+#        Utils.rename_variables(results, {'x': 'i', 'y': 'j'}, False)
+        var = handler_results.createVariable('heatc', float,
+                                             ('time', 'j', 'i'),
+                                             fill_value=1.e20)
         var.units = 'J m-2'
         handler_results.sync()
-        handler_results.variables['heatc'][:] = heatc_sl
+        handler_results.variables['heatc'][:] = heatc_sl[0] # temporary fix, needs to loop over layers
         handler_results.close()
 
+        results1D = TempFile.get()
+        handler_results1D = Utils.open_cdf(results1D, 'w')
+        Utils.copy_variable(handler, handler_results1D, 'time', True, True)
+        handler_results1D.createDimension('region', len(self.basins))
+        handler_results1D.createDimension('region_length', 50)
+        var_region = handler_results1D.createVariable('region', 'S1', ('region', 'region_length'))
+        var_ohc1D = handler_results1D.createVariable('heatcsum', float, ('time', 'region',),)
+        handler_results1D.sync()
+        for i, basin in enumerate(self.basins):
+            var_region[i, ...] = netCDF4.stringtoarr(basin.name, 50)
+            var_ohc1D[..., i] = heatc_sl1D[i]
+        handler_results1D.close()
+
         Utils.setminmax(results, 'heatc')
         self.heatc.set_local_file(results)
+        self.heatcsum.set_local_file(results1D)
-- 
GitLab