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/19] 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/19] 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/19] 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/19] 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/19] 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/19] 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/19] 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/19] 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/19] 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/19] 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/19] 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/19] 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/19] 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/19] 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/19] 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/19] 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/19] 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/19] 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/19] 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