Merge pull request #1382 from knutfrode/dev

Cleaned mapping of variables in ROMS native reader. If multiple varia…

opendrift/readers/basereader/__init__.py

@@ -95,7 +95,10 @@ class BaseReader(Variables, Combine, Filter):
                                  'surface_northward_geostrophic_sea_water_velocity_assuming_sea_level_for_geoid'],
         'x_wind': 'eastward_wind', 'y_wind': 'northward_wind',
         'sea_surface_wave_stokes_drift_x_velocity': 'eastward_surface_stokes_drift',
-        'sea_surface_wave_stokes_drift_y_velocity': 'northward_surface_stokes_drift'}
+        'sea_surface_wave_stokes_drift_y_velocity': 'northward_surface_stokes_drift',
+        'sea_ice_x_velocity': 'eastward_sea_ice_velocity',
+        'sea_ice_y_velocity': 'northward_sea_ice_velocity',
+        }
 
     def __init__(self):
         """Common constructor for all readers"""

opendrift/readers/reader_ROMS_native.py

@@ -91,18 +91,16 @@ def __init__(self, filename=None, name=None, gridfile=None, standard_name_mappin
         if filename is None:
             raise ValueError('Need filename as argument to constructor')
 
-        # Map ROMS variable names to CF standard_name
+        # Map ROMS variable names to CF standard_name, for cases where standard_name attribute is missing
         self.ROMS_variable_mapping = {
-            # Removing (temoprarily) land_binary_mask from ROMS-variables,
-            # as this leads to trouble with linearNDFast interpolation
             'mask_rho': 'land_binary_mask',
-            # 'mask_psi': 'land_binary_mask',  # don't want two variables mapping together - raises error now
+            'mask_psi': 'land_binary_mask',
             'h': 'sea_floor_depth_below_sea_level',
             'zeta': 'sea_surface_height',
             'u': 'x_sea_water_velocity',
             'v': 'y_sea_water_velocity',
-            'u_eastward': 'x_sea_water_velocity',
-            'v_northward': 'y_sea_water_velocity',
+            'u_eastward': 'eastward_sea_water_velocity',
+            'v_northward': 'northward_sea_water_velocity',
             'w': 'upward_sea_water_velocity',
             'temp': 'sea_water_temperature',
             'salt': 'sea_water_salinity',
@@ -277,25 +275,61 @@ def drop_non_essential_vars_pop(ds):
         self.precalculate_s2z_coefficients = True
 
         # Find all variables having standard_name
-        self.variables = []
+        self.standard_name_mapping = {}  # Inverse and unique mapping standard_name -> variable name
+        unmapped_variables = []
         for var_name in list(self.Dataset.variables):
             var = self.Dataset.variables[var_name]
             if 'standard_name' in var.attrs and var_name not in self.ROMS_variable_mapping.keys():
-                self.ROMS_variable_mapping[var_name] = var.attrs['standard_name']
-            if var_name in self.ROMS_variable_mapping.keys():
-                self.variables.append(self.ROMS_variable_mapping[var_name])
+                standard_name = var.attrs['standard_name']
+            elif var_name in self.ROMS_variable_mapping:
+                standard_name = self.ROMS_variable_mapping[var_name]
+            else:
+                unmapped_variables.append(var_name)
+                continue  # Variable cannot be mapped to standard_name
+            selected_variable = var_name
+            if standard_name in self.standard_name_mapping:  # We have a duplicate
+                if var_name in standard_name_mapping:  # provided by user
+                    selected_variable = var_name
+                    discarded_variable = self.standard_name_mapping[standard_name]
+                else:
+                    discarded_variable = var_name
+                    selected_variable = self.standard_name_mapping[standard_name]
+                logger.warning(f'Duplicate variables for {standard_name}, selecting {selected_variable}, '
+                            f'and discarding {discarded_variable}')
+            if standard_name in self.variable_aliases:  # Mapping for aliases
+                standard_name = self.variable_aliases[standard_name]
+            self.standard_name_mapping[standard_name] = selected_variable
+
+        if len(unmapped_variables) > 0:
+            logger.info(f'The following variables without standard_name are discarded: {unmapped_variables}')
 
         # A bit hackish solution:
         # If variable names or their standard_name contain "east" or "north", 
         # these should not be rotated from xi-direction to east-direction
         self.do_not_rotate = []
-        for var, stdname in self.ROMS_variable_mapping.items():
+        for stdname, var in self.standard_name_mapping.copy().items():
+            # Renaming east-north-names to x-y, as CRS of ROMS reader is always east-north,
+            # and we want to avoid rotating velocities if not neceassary
+            for xvar, eastnorthvar in self.xy2eastnorth_mapping.items():
+                if stdname in eastnorthvar:
+                    logger.info(f'Mapping {stdname} to {xvar} to avoid unnecessary rotation')
+                    self.standard_name_mapping[xvar] = var
+                    del self.standard_name_mapping[stdname]
+                    self.do_not_rotate.append(xvar)    
+
+        for stdname, var in self.standard_name_mapping.copy().items():
+            # Also avoid rotation of these variables
             if 'east' in var.lower() or 'east' in stdname.lower() or \
                     'north' in var.lower() or 'north' in stdname.lower():
-                self.do_not_rotate.append(stdname)
+                if stdname not in self.do_not_rotate:
+                    self.do_not_rotate.append(stdname)
+
         if len(self.do_not_rotate)>0:
             logger.debug('The following ROMS vectors are considered east-north, and will not be rotated %s' % self.do_not_rotate)
 
+
+        self.variables = list(self.standard_name_mapping)
+
         # Run constructor of parent Reader class
         super(Reader, self).__init__()
 
@@ -528,14 +562,8 @@ def get_mask(mask_name, imask, masks_store):
 
         masks_store = {}  # To store masks for various grids
         for par in requested_variables:
-            varname = [name for name, cf in
-                       self.ROMS_variable_mapping.items() if cf == par]
-            if len(varname) > 1:
-                raise ValueError("Multiple variables exist with standard name and "
-                                 "are present in reader. Either remove the duplicate mapping "
-                                 "or remove the variable from the reader."
-                                 "Variables: " + str(varname))
-            var = self.Dataset.variables[varname[0]]
+            varname = self.standard_name_mapping[par]
+            var = self.Dataset.variables[varname]
 
             if par == 'land_binary_mask':
                 variables[par] = self.land_binary_mask[imask]
@@ -579,7 +607,7 @@ def get_mask(mask_name, imask, masks_store):
             if var.ndim == 4:
                 # Regrid from sigma to z levels
                 if len(np.atleast_1d(indz)) > 1:
-                    logger.debug('sigma to z for ' + varname[0])
+                    logger.debug('sigma to z for ' + varname)
                     if self.precalculate_s2z_coefficients is True:
                         M = self.sea_floor_depth_below_sea_level.shape[0]
                         N = self.sea_floor_depth_below_sea_level.shape[1]

tests/readers/test_roms.py

@@ -50,13 +50,6 @@ def setUp(self):
         )
 
 
-    def test_catch_multiple_variables(self):
-        """Catch if multiple variables are present."""
-        standard_mapping = {"u_eastward": 'x_sea_water_velocity'}
-        reader = reader_ROMS_native.Reader(self.ds, standard_name_mapping=standard_mapping)
-        with pytest.raises(ValueError):
-            reader.get_variables("x_sea_water_velocity", datetime(1970, 1, 1), 0, 0, 0)
-
     def test_get_variables_u_eastward_wetdry_mask_rho(self):
         """get a variable from the dataset and verify correct mask was used
         
@@ -107,8 +100,9 @@ def test_get_variables_u_mask_u(self):
         
         wetdry_mask_u is not available so mask_u should be used.
         """
-        standard_mapping = {'u': 'x_sea_water_velocity'} 
-        reader = reader_ROMS_native.Reader(self.ds, standard_name_mapping=standard_mapping)
+        standard_mapping = {'u': 'x_sea_water_velocity', 'v': 'y_sea_water_velocity'} 
+        reader = reader_ROMS_native.Reader(self.ds.drop_vars(['u_eastward','v_northward']),
+                                           standard_name_mapping=standard_mapping)
         # drop wetdry_mask_rho to test that mask_rho is used
         reader.Dataset = reader.Dataset.drop_vars("wetdry_mask_rho")
         var, key = "u", "x_sea_water_velocity"
@@ -126,7 +120,8 @@ def test_get_variables_u_mask_u(self):
     def test_get_variables_u_depth_zeta_zero(self):
         """get a variable from the dataset and verify in depth."""
         standard_mapping = {'u': 'x_sea_water_velocity'}
-        reader = reader_ROMS_native.Reader(self.ds, standard_name_mapping=standard_mapping)
+        reader = reader_ROMS_native.Reader(self.ds.drop_vars(['u_eastward','v_northward']),
+                                           standard_name_mapping=standard_mapping)
         # drop wetdry_mask_rho to test that mask_rho is used
         reader.Dataset = reader.Dataset.drop_vars("wetdry_mask_rho")
         var, key = "u", "x_sea_water_velocity"
@@ -153,6 +148,24 @@ def test_get_variables_u_depth_zeta_nonzero(self):
 
         output = reader.get_variables(key, datetime(1970, 1, 1), 0, 0, [-5], testing=False)
 
+        u_expected = [1.81428571, 1.91428571]
+        assert key in output
+        assert reader.do_not_rotate == ['x_sea_water_velocity', 'y_sea_water_velocity']
+        assert np.allclose(output[key].data[1,0,:], u_expected)
+
+    def test_get_variables_u_depth_zeta_nonzero_noeastnorth(self):
+        """get a variable from the dataset and verify in depth."""
+        standard_mapping = {'u': 'x_sea_water_velocity'}
+        # Here we remove u_eastward and v_northward, so that u and v are used (and needs rotation)
+        reader = reader_ROMS_native.Reader(self.ds.drop_vars(['u_eastward','v_northward']), standard_name_mapping=standard_mapping)
+        # drop wetdry_mask_rho to test that mask_rho is used
+        reader.Dataset = reader.Dataset.drop_vars("wetdry_mask_rho")
+        var, key = "u", "x_sea_water_velocity"
+        zeta = np.ones((2,2,3))*2
+        reader.Dataset["zeta"] = (["ocean_time", "eta_rho", "xi_rho"], zeta)
+
+        output = reader.get_variables(key, datetime(1970, 1, 1), 0, 0, [-5], testing=False)
+
         u_expected = [0.64285714, 0.74285714]
         assert key in output
         assert reader.do_not_rotate == []
@@ -194,7 +207,7 @@ def test_rotate(self):
         # in opendrift
         standard_name_mapping = {"u": 'eastward_sea_water_velocity', "v": 'northward_sea_water_velocity'}
         reader = reader_ROMS_native.Reader(self.ds, standard_name_mapping=standard_name_mapping)
-        assert reader.do_not_rotate == ['eastward_sea_water_velocity', 'northward_sea_water_velocity']
+        assert reader.do_not_rotate == ['x_sea_water_velocity', 'y_sea_water_velocity']
 
         # currents have unusual clearly east/north variable names and should not
         # be rotated
@@ -206,12 +219,14 @@ def test_rotate(self):
         # analogous for winds
         standard_name_mapping = {"Uwind": 'eastward_wind', "Vwind": 'northward_wind'}
         reader = reader_ROMS_native.Reader(self.ds, standard_name_mapping=standard_name_mapping)
-        assert reader.do_not_rotate == ['eastward_wind', 'northward_wind']
+        assert reader.do_not_rotate == ['x_sea_water_velocity', 'y_sea_water_velocity',
+                                        'x_wind', 'y_wind']
 
         self.ds = self.ds.rename_vars({"Uwind": "Uwind_eastward", "Vwind": "Vwind_northward"})
         standard_name_mapping = {"Uwind_eastward": 'x_wind', "Vwind_northward": 'y_wind'}
         reader = reader_ROMS_native.Reader(self.ds, standard_name_mapping=standard_name_mapping)
-        assert reader.do_not_rotate == ['x_wind', 'y_wind']
+        assert reader.do_not_rotate == ['x_sea_water_velocity', 'y_sea_water_velocity',
+                                        'x_wind', 'y_wind']