From 1feb249c15d12965521ab802ac582d249a4a3899 Mon Sep 17 00:00:00 2001 From: Kristen Thyng Date: Thu, 21 Mar 2024 13:35:22 -0700 Subject: [PATCH 1/9] added sea_surface_height to models that subclass OceanDrift --- opendrift/models/chemicaldrift.py | 1 + opendrift/models/larvalfish.py | 1 + opendrift/models/model_template.py | 1 + opendrift/models/openberg.py | 1 + opendrift/models/openhns.py | 1 + opendrift/models/openoil/openoil.py | 1 + opendrift/models/pelagicegg.py | 1 + opendrift/models/plastdrift.py | 1 + opendrift/models/radionuclides.py | 1 + opendrift/models/sealice.py | 1 + opendrift/models/sedimentdrift.py | 1 + 11 files changed, 11 insertions(+) diff --git a/opendrift/models/chemicaldrift.py b/opendrift/models/chemicaldrift.py index 64d8e154a..9c2951b57 100644 --- a/opendrift/models/chemicaldrift.py +++ b/opendrift/models/chemicaldrift.py @@ -93,6 +93,7 @@ class ChemicalDrift(OceanDrift): required_variables = { 'x_sea_water_velocity': {'fallback': None}, 'y_sea_water_velocity': {'fallback': None}, + 'sea_surface_height': {'fallback': 0}, 'x_wind': {'fallback': 0}, 'y_wind': {'fallback': 0}, 'land_binary_mask': {'fallback': None}, diff --git a/opendrift/models/larvalfish.py b/opendrift/models/larvalfish.py index 61b2f1873..8be3efb01 100644 --- a/opendrift/models/larvalfish.py +++ b/opendrift/models/larvalfish.py @@ -68,6 +68,7 @@ class LarvalFish(OceanDrift): required_variables = { 'x_sea_water_velocity': {'fallback': 0}, 'y_sea_water_velocity': {'fallback': 0}, + 'sea_surface_height': {'fallback': 0}, 'sea_surface_wave_significant_height': {'fallback': 0}, 'x_wind': {'fallback': 0}, 'y_wind': {'fallback': 0}, diff --git a/opendrift/models/model_template.py b/opendrift/models/model_template.py index 96895e160..e197ee835 100644 --- a/opendrift/models/model_template.py +++ b/opendrift/models/model_template.py @@ -83,6 +83,7 @@ class ModelTemplate(OceanDrift): required_variables = { 'x_sea_water_velocity': {'fallback': 0}, 'y_sea_water_velocity': {'fallback': 0}, + 'sea_surface_height': {'fallback': 0}, 'x_wind': {'fallback': 0, 'important': False}, 'y_wind': {'fallback': 0, 'important': False}, 'ocean_vertical_diffusivity': {'fallback': 0.02, 'important': False, 'profiles': True}, diff --git a/opendrift/models/openberg.py b/opendrift/models/openberg.py index 68e29614e..9de178abf 100644 --- a/opendrift/models/openberg.py +++ b/opendrift/models/openberg.py @@ -59,6 +59,7 @@ class OpenBerg(OceanDrift): required_variables = { 'x_sea_water_velocity': {'fallback': None}, 'y_sea_water_velocity': {'fallback': None}, + 'sea_surface_height': {'fallback': 0}, 'sea_floor_depth_below_sea_level': {'fallback': 10000}, 'x_wind': {'fallback': None}, 'y_wind': {'fallback': None}, diff --git a/opendrift/models/openhns.py b/opendrift/models/openhns.py index 0567f9c14..0e8211887 100644 --- a/opendrift/models/openhns.py +++ b/opendrift/models/openhns.py @@ -109,6 +109,7 @@ class OpenHNS(OceanDrift): 'y_sea_water_velocity': { 'fallback': None }, + 'sea_surface_height': {'fallback': 0}, 'x_wind': { 'fallback': None }, diff --git a/opendrift/models/openoil/openoil.py b/opendrift/models/openoil/openoil.py index 9785f9603..19a683031 100644 --- a/opendrift/models/openoil/openoil.py +++ b/opendrift/models/openoil/openoil.py @@ -229,6 +229,7 @@ class OpenOil(OceanDrift): 'y_wind': { 'fallback': None }, + 'sea_surface_height': {'fallback': 0}, 'upward_sea_water_velocity': { 'fallback': 0, 'important': False diff --git a/opendrift/models/pelagicegg.py b/opendrift/models/pelagicegg.py index 6220206b4..452bfc049 100644 --- a/opendrift/models/pelagicegg.py +++ b/opendrift/models/pelagicegg.py @@ -59,6 +59,7 @@ class PelagicEggDrift(OceanDrift): required_variables = { 'x_sea_water_velocity': {'fallback': 0}, 'y_sea_water_velocity': {'fallback': 0}, + 'sea_surface_height': {'fallback': 0}, 'sea_surface_wave_significant_height': {'fallback': 0}, 'sea_ice_area_fraction': {'fallback': 0}, 'x_wind': {'fallback': 0}, diff --git a/opendrift/models/plastdrift.py b/opendrift/models/plastdrift.py index 77520136c..2b2ddf388 100644 --- a/opendrift/models/plastdrift.py +++ b/opendrift/models/plastdrift.py @@ -44,6 +44,7 @@ class PlastDrift(OceanDrift): required_variables = { 'x_sea_water_velocity': {'fallback': 0}, 'y_sea_water_velocity': {'fallback': 0}, + 'sea_surface_height': {'fallback': 0}, 'sea_surface_wave_stokes_drift_x_velocity': {'fallback': 0}, 'sea_surface_wave_stokes_drift_y_velocity': {'fallback': 0}, 'sea_surface_wave_significant_height': {'fallback': 0}, diff --git a/opendrift/models/radionuclides.py b/opendrift/models/radionuclides.py index d60c48fa2..0e0e8c26f 100644 --- a/opendrift/models/radionuclides.py +++ b/opendrift/models/radionuclides.py @@ -75,6 +75,7 @@ class RadionuclideDrift(OceanDrift): required_variables = { 'x_sea_water_velocity': {'fallback': None}, 'y_sea_water_velocity': {'fallback': None}, + 'sea_surface_height': {'fallback': 0}, 'x_wind': {'fallback': 0}, 'y_wind': {'fallback': 0}, 'land_binary_mask': {'fallback': None}, diff --git a/opendrift/models/sealice.py b/opendrift/models/sealice.py index 259db4557..cfa65ee6a 100644 --- a/opendrift/models/sealice.py +++ b/opendrift/models/sealice.py @@ -84,6 +84,7 @@ class SeaLice(OceanDrift): required_variables = { 'x_sea_water_velocity': {'fallback': 0}, 'y_sea_water_velocity': {'fallback': 0}, + 'sea_surface_height': {'fallback': 0}, # 'sea_surface_wave_significant_height': {'fallback': 0}, # 'x_wind': {'fallback': 0}, # 'y_wind': {'fallback': 0}, diff --git a/opendrift/models/sedimentdrift.py b/opendrift/models/sedimentdrift.py index f74880300..553e4d8db 100644 --- a/opendrift/models/sedimentdrift.py +++ b/opendrift/models/sedimentdrift.py @@ -45,6 +45,7 @@ class SedimentDrift(OceanDrift): required_variables = { 'x_sea_water_velocity': {'fallback': 0}, 'y_sea_water_velocity': {'fallback': 0}, + 'sea_surface_height': {'fallback': 0}, 'upward_sea_water_velocity': {'fallback': 0}, 'x_wind': {'fallback': 0}, 'y_wind': {'fallback': 0}, From ad2dae68560f04390e8c788359da0b8c05a6d128 Mon Sep 17 00:00:00 2001 From: Kristen Thyng Date: Thu, 21 Mar 2024 13:37:02 -0700 Subject: [PATCH 2/9] change min sea_floor_depth_below_sea_level make min more negative to accommodate more extreme wet/dry values. --- opendrift/readers/basereader/consts.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opendrift/readers/basereader/consts.py b/opendrift/readers/basereader/consts.py index b54d8d8c2..c073138d0 100644 --- a/opendrift/readers/basereader/consts.py +++ b/opendrift/readers/basereader/consts.py @@ -17,7 +17,7 @@ '(northwards if projection is lonlat/Mercator)'}, 'land_binary_mask': {'valid_min': 0, 'valid_max': 1, 'long_name': '1 is land, 0 is sea'}, - 'sea_floor_depth_below_sea_level': {'valid_min': -10, 'valid_max': 12000, + 'sea_floor_depth_below_sea_level': {'valid_min': -20, 'valid_max': 12000, 'long_name': 'Depth of seafloor'}, 'ocean_vertical_diffusivity': {'valid_min': 0, 'valid_max': 1}} From 0395895835db305111510a936e5e216427276409 Mon Sep 17 00:00:00 2001 From: Kristen Thyng Date: Thu, 21 Mar 2024 13:52:49 -0700 Subject: [PATCH 3/9] main changes here * sea_surface_height added as required_variable to OceanDrift model * general:seafloor_action_dcrit added as a new configuration parameter which is the minimum wet value in a grid cell when a model running wetting and drying * updating vertical_mixing so that bottom checks account for depth of water column and sea_surface_height combined * updated interact_with_seafloor so that it accounts for zeta and new dcrit parameter in bottom checks too * new sea_surface_height function analogous to sea_floor_depth function in physics_methods that is for active drifters --- opendrift/models/basemodel/__init__.py | 22 +++++++++++++++++++--- opendrift/models/oceandrift.py | 14 +++++++++++--- opendrift/models/physics_methods.py | 19 +++++++++++++++++++ 3 files changed, 49 insertions(+), 6 deletions(-) diff --git a/opendrift/models/basemodel/__init__.py b/opendrift/models/basemodel/__init__.py index 53d3ee19d..3ed4e02b1 100644 --- a/opendrift/models/basemodel/__init__.py +++ b/opendrift/models/basemodel/__init__.py @@ -694,8 +694,25 @@ def interact_with_seafloor(self): return if 'sea_floor_depth_below_sea_level' not in self.env.priority_list: return + + if not hasattr(self, 'environment') or not hasattr( + self.environment, 'sea_surface_height'): + logger.warning('Seafloor check not being run because sea_surface_height is missing. ' + 'This will happen the first time the function is run but if it happens ' + 'subsequently there is probably a problem.') + return + + # the shape of these is different than the original arrays + # because it is for active drifters sea_floor_depth = self.sea_floor_depth() - below = np.where(self.elements.z < -sea_floor_depth)[0] + sea_surface_height = self.sea_surface_height() + + # Check if any elements are below sea floor + # But remember that the water column is the sea floor depth + sea surface height + parameter Dcrit + Dcrit = self.get_config('general:seafloor_action_dcrit') + ibelow = self.elements.z < -(sea_floor_depth + sea_surface_height + Dcrit) + below = np.where(ibelow)[0] + if len(below) == 0: logger.debug('No elements hit seafloor.') return @@ -705,8 +722,7 @@ def interact_with_seafloor(self): logger.debug('Lifting %s elements to seafloor.' % len(below)) self.elements.z[below] = -sea_floor_depth[below] elif i == 'deactivate': - self.deactivate_elements(self.elements.z < -sea_floor_depth, - reason='seafloor') + self.deactivate_elements(ibelow, reason='seafloor') self.elements.z[below] = -sea_floor_depth[below] elif i == 'previous': # Go back to previous position (in water) logger.warning('%s elements hit seafloor, ' diff --git a/opendrift/models/oceandrift.py b/opendrift/models/oceandrift.py index b91878a53..bfcc839b5 100644 --- a/opendrift/models/oceandrift.py +++ b/opendrift/models/oceandrift.py @@ -69,6 +69,7 @@ class OceanDrift(OpenDriftSimulation): required_variables = { 'x_sea_water_velocity': {'fallback': 0}, 'y_sea_water_velocity': {'fallback': 0}, + 'sea_surface_height': {'fallback': 0}, 'x_wind': {'fallback': 0}, 'y_wind': {'fallback': 0}, 'upward_sea_water_velocity': {'fallback': 0}, @@ -158,6 +159,10 @@ def __init__(self, *args, **kwargs): 'enum': ['none', 'lift_to_seafloor', 'deactivate', 'previous'], 'description': '"deactivate": elements are deactivated; "lift_to_seafloor": elements are lifted to seafloor level; "previous": elements are moved back to previous position; "none"; seafloor is ignored.', 'level': CONFIG_LEVEL_ADVANCED}, + 'general:seafloor_action_dcrit': {'type': 'float', 'default': 0.0, + 'min': 0.0, 'max': 10000, 'units': 'm', + 'description': 'This parameter represents the amount of water left in a grid cell to keep it wet in a wet/dry simulation for numerical stability. The condition checked for seafloor_action is z < -(sea_floor_depth + sea_surface_height + `general:seafloor_action_dcrit`. It is 0 by default to assume that a wet/dry case is not being run, however, if it is and the correct value is not known 0.1 is a good default to use.', + 'level': CONFIG_LEVEL_ADVANCED}, 'drift:truncate_ocean_model_below_m': {'type': 'float', 'default': None, 'min': 0, 'max': 10000, 'units': 'm', 'description': 'Ocean model data are only read down to at most this depth, and extrapolated below. May be specified to read less data to improve performance.', @@ -489,8 +494,10 @@ def vertical_mixing(self, store_depths=False): dt_mix = self.get_config('vertical_mixing:timestep') - # minimum height/maximum depth for each particle - Zmin = -1.*self.environment.sea_floor_depth_below_sea_level + # minimum height/maximum depth for each particle accouting also for + # the sea surface height and the critical/min wet cell depth + Dcrit = self.get_config('general:seafloor_action_dcrit') + Zmin = -1.*(self.environment.sea_floor_depth_below_sea_level + self.environment.sea_surface_height + Dcrit) # Eventual model specific preparions self.prepare_vertical_mixing() @@ -609,6 +616,7 @@ def vertical_mixing(self, store_depths=False): self.elements.z[reflect] = -self.elements.z[reflect] # Reflect elements going below seafloor + # Zmin accounts for sea_surface_height bottom = np.where(np.logical_and(self.elements.z < Zmin, self.elements.moving == 1)) if len(bottom[0]) > 0: logger.debug('%s elements penetrated seafloor, lifting up' % len(bottom[0])) @@ -629,7 +637,7 @@ def vertical_mixing(self, store_depths=False): # Let particles stick to bottom bottom = np.where(self.elements.z < Zmin) if len(bottom[0]) > 0: - logger.debug('%s elements reached seafloor, set to bottom' % len(bottom[0])) + logger.debug('%s elements reached seafloor, interacting with bottom' % len(bottom[0])) self.interact_with_seafloor() self.bottom_interaction(Zmin) diff --git a/opendrift/models/physics_methods.py b/opendrift/models/physics_methods.py index 21294309a..804debd49 100644 --- a/opendrift/models/physics_methods.py +++ b/opendrift/models/physics_methods.py @@ -1114,6 +1114,25 @@ def sea_floor_depth(self): env['sea_floor_depth_below_sea_level'].astype('float32') return sea_floor_depth + def sea_surface_height(self): + '''Sea surface height (positive/negative) for presently active elements''' + + if hasattr(self, 'environment') and \ + hasattr(self.environment, 'sea_surface_height'): + if len(self.environment.sea_surface_height) == \ + self.num_elements_active(): + sea_surface_height = \ + self.environment.sea_surface_height + if 'sea_surface_height' not in locals(): + env, env_profiles, missing = \ + self.env.get_environment(['sea_surface_height'], + time=self.time, lon=self.elements.lon, + lat=self.elements.lat, + z=0*self.elements.lon, profiles=None) + sea_surface_height = \ + env['sea_surface_height'].astype('float32') + return sea_surface_height + def wind_drag_coefficient(windspeed): '''Large and Pond (1981), J. Phys. Oceanog., 11, 324-336.''' From afd1e928451c43e29f2c09892d3f16a75a8e7180 Mon Sep 17 00:00:00 2001 From: Kristen Thyng Date: Thu, 21 Mar 2024 14:16:54 -0700 Subject: [PATCH 4/9] found another Zmin to update --- opendrift/models/oceandrift.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/opendrift/models/oceandrift.py b/opendrift/models/oceandrift.py index bfcc839b5..79ffe980e 100644 --- a/opendrift/models/oceandrift.py +++ b/opendrift/models/oceandrift.py @@ -416,13 +416,15 @@ def vertical_buoyancy(self): self.elements.z[in_ocean] = np.minimum(0, self.elements.z[in_ocean] + self.elements.terminal_velocity[in_ocean] * self.time_step.total_seconds()) - # check for minimum height/maximum depth for each particle - Zmin = -1.*self.environment.sea_floor_depth_below_sea_level + # check for minimum height/maximum depth for each particle accouting also for + # the sea surface height and the critical/min wet cell depth + Dcrit = self.get_config('general:seafloor_action_dcrit') + Zmin = -1.*(self.environment.sea_floor_depth_below_sea_level + self.environment.sea_surface_height + Dcrit) # Let particles stick to bottom bottom = np.where(self.elements.z < Zmin) if len(bottom[0]) > 0: - logger.debug('%s elements reached seafloor, set to bottom' % len(bottom[0])) + logger.debug('%s elements reached seafloor, interacting with bottom' % len(bottom[0])) self.interact_with_seafloor() self.bottom_interaction(Zmin) From 1b6862b43e3df5aad6dc8441b2e3fe610a9aac4d Mon Sep 17 00:00:00 2001 From: Kristen Thyng Date: Thu, 21 Mar 2024 14:20:56 -0700 Subject: [PATCH 5/9] updated tests * test_run had a test that required sea_surface_height values be added * test_stranding values for some checks changed because zeta was available in the reader and is used now --- tests/models/test_run.py | 6 +++--- tests/models/test_stranding.py | 8 ++++---- 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/tests/models/test_run.py b/tests/models/test_run.py index e670b3a1c..4638ba1f0 100644 --- a/tests/models/test_run.py +++ b/tests/models/test_run.py @@ -390,9 +390,9 @@ def test_vertical_mixing_profiles(self): o.time = time o.time_step = timedelta(hours=2) o.release_elements() - o.environment = np.array(np.ones(N)*100, - dtype=[('sea_floor_depth_below_sea_level', - np.float32)]).view(np.recarray) + o.environment = np.array([(100, 0) for _ in range(N)], + dtype=[('sea_floor_depth_below_sea_level', np.float32), + ('sea_surface_height', np.float32)]).view(np.recarray) o.environment_profiles = { 'z': z, 'ocean_vertical_diffusivity': np.tile(diffusivity, (N, 1)).T} o.env.finalize() diff --git a/tests/models/test_stranding.py b/tests/models/test_stranding.py index 032a29000..f239c81a0 100644 --- a/tests/models/test_stranding.py +++ b/tests/models/test_stranding.py @@ -50,15 +50,15 @@ def test_stranding_3d(self): self.assertEqual(o.elements_deactivated.status.min(), 1) self.assertEqual(o.elements_deactivated.status.max(), 1) self.assertEqual(o.num_elements_scheduled(), 0) - self.assertEqual(o.num_elements_active(), 77) + self.assertEqual(o.num_elements_active(), 82) self.assertEqual(o.num_elements_activated(), 100) - self.assertEqual(o.num_elements_deactivated(), 23) + self.assertEqual(o.num_elements_deactivated(), 18) self.assertEqual(o.num_elements_total(), 100) # Check calculated trajectory lengths and speeds total_length, distances, speeds = o.get_trajectory_lengths() - self.assertAlmostEqual(total_length.max(), 14974.8, 1) - self.assertAlmostEqual(total_length.min(), 1222.3, 1) + self.assertAlmostEqual(total_length.max(), 15719.5, 1) + self.assertAlmostEqual(total_length.min(), 1225.2, 1) self.assertAlmostEqual(speeds.max(), 0.132, 1) self.assertAlmostEqual(distances.max(), 2859.5, 1) From 006cee6f53c966999d9e195c4ae46e4848ae0bdd Mon Sep 17 00:00:00 2001 From: Kristen Thyng Date: Fri, 22 Mar 2024 15:39:59 -0700 Subject: [PATCH 6/9] improved angle handling in ROMS reader also improved handling of rotation of vectors --- opendrift/readers/reader_ROMS_native.py | 93 +++++++++++++++---------- 1 file changed, 58 insertions(+), 35 deletions(-) diff --git a/opendrift/readers/reader_ROMS_native.py b/opendrift/readers/reader_ROMS_native.py index ced5a33cb..6066d42bf 100644 --- a/opendrift/readers/reader_ROMS_native.py +++ b/opendrift/readers/reader_ROMS_native.py @@ -82,11 +82,11 @@ def __init__(self, filename=None, name=None, gridfile=None, standard_name_mappin self._mask_u = None self._mask_v = None self._zeta = None + self._angle = None self.land_binary_mask = None self.sea_floor_depth_below_sea_level = None self.z_rho_tot = None self.s2z_A = None - self.angle_xi_east = None if filename is None: raise ValueError('Need filename as argument to constructor') @@ -153,7 +153,7 @@ def drop_non_essential_vars_pop(ds): dropvars = [v for v in ds.variables if v not in list(self.ROMS_variable_mapping.keys()) + gls_param + ['ocean_time', 'time', 'bulk_time', 's_rho', - 'Cs_r', 'hc', 'angle', 'Vtransform'] + 'Cs_r', 'hc', 'Vtransform'] and v[0:3] not in ['lon', 'lat', 'mas']] logger.debug('Dropping variables: %s' % dropvars) ds = ds.drop_vars(dropvars) @@ -230,7 +230,7 @@ def drop_non_essential_vars_pop(ds): self.lat = self.lat.data if self.lat.ndim == 1: self.lon, self.lat = np.meshgrid(self.lon, self.lat) - self.angle_xi_east = 0 + # self.angle_xi_east = 0 # this was moved to the angle property else: raise ValueError(filename + ' does not contain lon/lat ' 'arrays, please supply a grid-file: "gridfile="') @@ -372,9 +372,23 @@ def zeta(self): self._zeta = np.zeros(self.mask_rho.shape) logger.info("No zeta found, using 0 array for sea surface height") return self._zeta + + @property + def angle(self): + """Grid angle if curvilinear.""" + if self._angle is None: + if 'lat_rho' in self.Dataset.variables and self.lat.ndim == 1: + self._angle = 0 + elif 'angle' in self.Dataset.data_vars: + self._angle = self.Dataset.variables['angle'] + logger.info("Using angle from Dataset.") + + if self._angle is None: + raise ValueError('No angle between xi and east found') + return self._angle def get_variables(self, requested_variables, time=None, - x=None, y=None, z=None): + x=None, y=None, z=None, testing=False): start_time = datetime.now() requested_variables, time, x, y, z, outside = self.check_arguments( requested_variables, time, x, y, z) @@ -483,17 +497,22 @@ def get_variables(self, requested_variables, time=None, # define another set of indices itzxy = (indxTime, indz, indy, indx) - def get_mask(mask_name, imask): - if mask_name in masks_for_loop: - mask = masks_for_loop[mask_name] + def get_mask(mask_name, imask, masks_store): + if mask_name in masks_store: + mask = masks_store[mask_name] else: mask = getattr(self, mask_name)[imask] return mask, mask_name - masks_for_loop = {} # To store maskes for various grids + 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]] if par == 'land_binary_mask': @@ -518,15 +537,15 @@ def get_mask(mask_name, imask): # make sure that var has matching horizontal dimensions with the mask # make sure coord names also match if self.mask_rho.shape[-2:] == var.shape[-2:] and self.mask_rho.dims[-2:] == var.dims[-2:]: - mask, mask_name = get_mask("mask_rho", imask) + mask, mask_name = get_mask("mask_rho", imask, masks_store) elif self.mask_u.shape[-2:] == var.shape[-2:] and self.mask_u.dims[-2:] == var.dims[-2:]: - mask, mask_name = get_mask("mask_u", imask) + mask, mask_name = get_mask("mask_u", imask, masks_store) elif self.mask_v.shape[-2:] == var.shape[-2:] and self.mask_v.dims[-2:] == var.dims[-2:]: - mask, mask_name = get_mask("mask_v", imask) + mask, mask_name = get_mask("mask_v", imask, masks_store) else: raise Exception('No mask found for ' + par) - masks_for_loop[mask_name] = np.asarray(mask) + masks_store[mask_name] = np.asarray(mask) mask = np.asarray(mask) if mask.min() == 0: @@ -659,29 +678,24 @@ def get_mask(mask_name, imask): variables['y'] = variables['y'].astype(np.float32) variables['time'] = nearestTime - if 'x_sea_water_velocity' in variables.keys() or \ - 'sea_ice_x_velocity' in variables.keys() or \ - 'x_wind' in variables.keys() and \ - 'x_sea_water_velocity' not in self.do_not_rotate and \ - 'sea_ice_x_velocity' not in self.do_not_rotate and \ - 'x_wind' not in self.do_not_rotate: - # We must rotate current vectors - if self.angle_xi_east is None: - if 'angle' in self.Dataset.variables: - logger.debug('Reading angle between xi and east...') - self.angle_xi_east = self.Dataset.variables['angle'][:] - if isinstance(self.angle_xi_east, int): - rad = self.angle_xi_east + if 'x_sea_water_velocity' in variables.keys() and \ + 'x_sea_water_velocity' not in self.do_not_rotate: + if isinstance(self.angle, int): + rad = self.angle + else: + rad = np.ma.asarray(self.angle[indy, indx]) + variables['x_sea_water_velocity'], \ + variables['y_sea_water_velocity'] = rotate_vectors_angle( + variables['x_sea_water_velocity'], + variables['y_sea_water_velocity'], rad) + logger.debug('Rotated x_sea_water_velocity and y_sea_water_velocity') + + if 'sea_ice_x_velocity' in variables.keys() and \ + 'sea_ice_x_velocity' not in self.do_not_rotate: + if isinstance(self.angle, int): + rad = self.angle else: - rad = self.angle_xi_east[indy, indx] - rad = np.ma.asarray(rad) - if 'x_sea_water_velocity' in variables.keys() and \ - 'x_sea_water_velocity' not in self.do_not_rotate: - variables['x_sea_water_velocity'], \ - variables['y_sea_water_velocity'] = rotate_vectors_angle( - variables['x_sea_water_velocity'], - variables['y_sea_water_velocity'], rad) - logger.debug('Rotated x_sea_water_velocity and y_sea_water_velocity') + rad = np.ma.asarray(self.angle[indy, indx]) if 'sea_ice_x_velocity' in variables.keys() and \ 'sea_ice_x_velocity' not in self.do_not_rotate: variables['sea_ice_x_velocity'], \ @@ -689,6 +703,12 @@ def get_mask(mask_name, imask): variables['sea_ice_x_velocity'], variables['sea_ice_y_velocity'], rad) logger.debug('Rotated sea_ice_x_velocity and sea_ice_y_velocity') + + if 'x_wind' in variables.keys() and 'x_wind' not in self.do_not_rotate: + if isinstance(self.angle, int): + rad = self.angle + else: + rad = np.ma.asarray(self.angle[indy, indx]) if 'x_wind' in variables.keys() and \ 'x_wind' not in self.do_not_rotate: variables['x_wind'], \ @@ -703,7 +723,10 @@ def get_mask(mask_name, imask): logger.debug('Time for ROMS native reader: ' + str(datetime.now()-start_time)) - return variables + if testing: + return variables, masks_store + else: + return variables def rotate_vectors_angle(u, v, radians): From e3c9ed16780f7b7ff18c122aadfba77c739072a9 Mon Sep 17 00:00:00 2001 From: Kristen Thyng Date: Fri, 22 Mar 2024 15:40:34 -0700 Subject: [PATCH 7/9] added ROMS reader test file plan to add one more set of tests for zeta stuff specifically. --- tests/readers/test_roms.py | 185 +++++++++++++++++++++++++++++++++++++ 1 file changed, 185 insertions(+) create mode 100644 tests/readers/test_roms.py diff --git a/tests/readers/test_roms.py b/tests/readers/test_roms.py new file mode 100644 index 000000000..10a1edf06 --- /dev/null +++ b/tests/readers/test_roms.py @@ -0,0 +1,185 @@ +import unittest +from opendrift.readers import reader_ROMS_native +import xarray as xr +import numpy as np +from datetime import datetime +import pytest + + +class TestROMSReader(unittest.TestCase): + def setUp(self): + + zeta = np.array([[[0.1, 0.2, 0.3], [0.4, 0.5, 0.6]], [[0.7, 0.8, 0.9], [1.0, 1.1, 1.2]]]) + # zeta = np.zeros((2,2,3)) + + u_eastward = np.array([[[[1.0, 1.1, 1.2], [1.3, 1.4, 1.5]], [[1.6, 1.7, 1.8], [1.9, 2.0, 2.1]], [[2.2, 2.3, 2.4], [2.5, 2.6, 2.7]]], [[[2.8, 2.9, 3.0], [3.1, 3.2, 3.3]], [[3.4, 3.5, 3.6], [3.7, 3.8, 3.9]], [[4.0, 4.1, 4.2], [4.3, 4.4, 4.5]]]]) + + u = np.array([[[[0.1, 0.2], [0.3, 0.4]], [[0.5, 0.6], [0.7, 0.8]], [[0.9, 1.0], [1.1, 1.2]]], [[[1.3, 1.4], [1.5, 1.6]], [[1.7, 1.8], [1.9, 2.0]], [[2.1, 2.2], [2.3, 2.4]]]]) + v = np.linspace(0, 1, 2*3*1*3).reshape(2, 3, 1, 3) + Cs_r = np.linspace(-1, 0, 3) + mask_rho = np.ones((2, 3)) + # use odd mask values so I can check when used in tests + wetdry_mask_rho = np.ones((2, 2, 3))*2 + mask_u = np.ones((2, 2)) + mask_v = np.ones((1, 3)) + angle = np.zeros((2, 3)) + + self.ds = xr.Dataset( + { + "lon_rho": (["eta_rho", "xi_rho"], [[-152.0, -151.9, -151.8], [-152.0, -151.9, -151.8]]), + "lat_rho": (["eta_rho", "xi_rho"], [[60.5, 60.5, 60.5], [60.4, 60.4, 60.4]]), + "lon_u": (["eta_u", "xi_u"], [[-151.95, -151.85], [-151.95, -151.85]]), + "lat_u": (["eta_u", "xi_u"], [[60.5, 60.5], [60.4, 60.4]]), + "lon_v": (["eta_v", "xi_v"], [[-152.0, -151.9, -151.8]]), + "lat_v": (["eta_v", "xi_v"], [[60.45, 60.45, 60.45]]), + "mask_rho": (["eta_rho", "xi_rho"], mask_rho), + "wetdry_mask_rho": (["ocean_time", "eta_rho", "xi_rho"], wetdry_mask_rho), + "angle": (["eta_rho", "xi_rho"], angle), + "mask_u": (["eta_u", "xi_u"], mask_u), + "mask_v": (["eta_v", "xi_v"], mask_v), + "ocean_time": ("ocean_time", [0, 1], {"units": "seconds since 1970-01-01"}), + "zeta": (["ocean_time", "eta_rho", "xi_rho"], zeta), + "s_rho": (["s_rho"], Cs_r), + "Cs_r": (["s_rho"], Cs_r), + "hc": 16, + "u_eastward": (("ocean_time", "s_rho", "eta_rho", "xi_rho"), u_eastward), + "v_northward": (("ocean_time", "s_rho", "eta_rho", "xi_rho"), u_eastward), + "u": (("ocean_time", "s_rho", "eta_u", "xi_u"), u), + "v": (("ocean_time", "s_rho", "eta_v", "xi_v"), v), + } + ) + + + 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 + + Since wetdry_mask_rho is available and matches u_eastward dims, it should be used. + """ + standard_mapping = {"u_eastward": 'x_sea_water_velocity', + "v_northward": 'y_sea_water_velocity', + 'u': 'blah', + 'v': 'blah'} # redirect mapping of u/v so these don't overlap + reader = reader_ROMS_native.Reader(self.ds, standard_name_mapping=standard_mapping) + var, key = "u_eastward", "x_sea_water_velocity" + output, masks = reader.get_variables(key, datetime(1970, 1, 1), 0, 0, 0, testing=True) + + # check that variable was handled correctly + assert key in output + assert reader.do_not_rotate == ["x_sea_water_velocity", "y_sea_water_velocity"] + assert np.allclose(output[key].data, reader.Dataset[var][0,-1,0,0:2].values) + # check that wetdry mask was used (can tell by mask value) + assert np.allclose(masks["mask_rho"], reader.Dataset["wetdry_mask_rho"][0,0,0:2].values) + + def test_get_variables_u_eastward_mask_rho(self): + """get a variable from the dataset and verify correct mask was used + + wetdry_mask_rho is not available so mask_rho should be used. + """ + standard_mapping = {"u_eastward": 'x_sea_water_velocity', + "v_northward": 'y_sea_water_velocity', + 'u': 'blah', + 'v': 'blah'} # redirect mapping of u/v so these don't overlap + reader = reader_ROMS_native.Reader(self.ds, 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_eastward", "x_sea_water_velocity" + output, masks = reader.get_variables(key, datetime(1970, 1, 1), 0, 0, 0, testing=True) + + # check that variable was handled correctly + assert key in output + assert reader.do_not_rotate == ["x_sea_water_velocity", "y_sea_water_velocity"] + assert np.allclose(output[key].data, reader.Dataset[var][0,-1,0,0:2].values) + # check that mask was used (can tell by mask value) + # only mask_rho should have been used since u_eastward on rho grid + assert "mask_u" not in masks + assert np.allclose(masks["mask_rho"], reader.Dataset["mask_rho"][0,0:2].values) + + + def test_get_variables_u_mask_u(self): + """get a variable from the dataset and verify correct mask was used + + 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) + # 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" + output, masks = reader.get_variables(key, datetime(1970, 1, 1), 0, 0, 0, testing=True) + + # check that variable was handled correctly + assert key in output + assert reader.do_not_rotate == [] + assert np.allclose(output[key].data, reader.Dataset[var][0,-1,0,0:2].values) + # check that mask was used (can tell by mask value) + # mask_rho should not have been used since u on u-grid + assert "mask_rho" not in masks + + +class TestROMSReaderRotation(unittest.TestCase): + def setUp(self): + + self.ds = xr.Dataset( + data_vars={ + "u": (("ocean_time", "Z", "Y", "X"), np.zeros((2, 3, 2, 3))), + "v": (("ocean_time", "Z", "Y", "X"), np.zeros((2, 3, 2, 3))), + "Uwind": (("ocean_time", "Y", "X"), np.zeros((2, 2, 3))), + "Vwind": (("ocean_time", "Y", "X"), np.zeros((2, 2, 3))), + "Cs_r": (("Z"), np.linspace(-1, 0, 3)), + "hc": 16, + }, + coords={ + "ocean_time": ("ocean_time", [0, 1], {"units": "seconds since 1970-01-01"}), + "s_rho": (("Z"), np.linspace(-1, 0, 3)), + "lon_rho": (("Y", "X"), np.array([[1, 2, 3], [1, 2, 3]])), + "lat_rho": (("Y", "X"), np.array([[1, 1, 1], [2, 2, 2]])), + }, + ) + + + def test_rotate(self): + """test that correct variables are rotated""" + + # currents and winds have normal ROMS names and should be rotated + reader = reader_ROMS_native.Reader(self.ds) + assert reader.do_not_rotate == [] + + # having the "east" or "north" in the standard_name attributes also cause + # a variable to not be rotated + # this is to test the code though this standard mapping would not lead to + # these standard names actually being used outside of the ROMS reader + # 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'] + + # currents have unusual clearly east/north variable names and should not + # be rotated + self.ds = self.ds.rename_vars({"u": "u_eastward", "v": "v_northward"}) + standard_name_mapping = {"u_eastward": 'x_sea_water_velocity', "v_northward": 'y_sea_water_velocity'} + reader = reader_ROMS_native.Reader(self.ds, standard_name_mapping=standard_name_mapping) + assert reader.do_not_rotate == ['x_sea_water_velocity', 'y_sea_water_velocity'] + + # 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'] + + 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'] + + + + + +if __name__ == '__main__': + unittest.main() From 8f4b10a6a4465000baa44c53b9f24b6586ea49b6 Mon Sep 17 00:00:00 2001 From: Kristen Thyng Date: Thu, 28 Mar 2024 08:25:56 -0700 Subject: [PATCH 8/9] transform z_rho to match convention of z --- opendrift/readers/reader_ROMS_native.py | 17 ++++++++++++++++- 1 file changed, 16 insertions(+), 1 deletion(-) diff --git a/opendrift/readers/reader_ROMS_native.py b/opendrift/readers/reader_ROMS_native.py index 6066d42bf..7d29d0b3c 100644 --- a/opendrift/readers/reader_ROMS_native.py +++ b/opendrift/readers/reader_ROMS_native.py @@ -96,7 +96,7 @@ def __init__(self, filename=None, name=None, gridfile=None, standard_name_mappin # 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', + # 'mask_psi': 'land_binary_mask', # don't want two variables mapping together - raises error now 'h': 'sea_floor_depth_below_sea_level', 'zeta': 'sea_surface_height', 'u': 'x_sea_water_velocity', @@ -382,6 +382,9 @@ def angle(self): elif 'angle' in self.Dataset.data_vars: self._angle = self.Dataset.variables['angle'] logger.info("Using angle from Dataset.") + # else: + # self._angle = 0 + # logger.warning("No angle found, using 0 integer for angle.") if self._angle is None: raise ValueError('No angle between xi and east found') @@ -465,11 +468,23 @@ def get_variables(self, requested_variables, time=None, zeta = self.zeta[indxTime] self.z_rho_tot = depth.sdepth(Htot, zeta, self.hc, self.Cs_r, Vtransform=self.Vtransform) + # z_rho is positive relative to mean sea level but z is + # 0 at the surface. + # Transform z_rho to match convention of z. + self.z_rho_tot -= np.asarray(zeta)[np.newaxis] H = self.sea_floor_depth_below_sea_level[indy, indx] zeta = self.zeta[itxy] z_rho = depth.sdepth(H, zeta, self.hc, self.Cs_r, Vtransform=self.Vtransform) + + # z_rho is positive relative to mean sea level but z is + # 0 at the surface. + # Transform z_rho to match convention of z. + z_rho -= np.asarray(zeta)[np.newaxis] + + assert (z_rho <=0).all() + # Element indices must be relative to extracted subset indx_el = np.clip(indx_el - indx.min(), 0, z_rho.shape[2]-1) indy_el = np.clip(indy_el - indy.min(), 0, z_rho.shape[1]-1) From d0a1245e8220558dffda7397b8d397777586e7f8 Mon Sep 17 00:00:00 2001 From: Kristen Thyng Date: Thu, 28 Mar 2024 08:28:47 -0700 Subject: [PATCH 9/9] updated roms tests --- tests/models/test_stranding.py | 4 ++-- tests/readers/test_depth.py | 16 ++++++------- tests/readers/test_netcdf.py | 4 ++++ tests/readers/test_roms.py | 42 +++++++++++++++++++++++++++++++--- 4 files changed, 53 insertions(+), 13 deletions(-) diff --git a/tests/models/test_stranding.py b/tests/models/test_stranding.py index f239c81a0..9b0988267 100644 --- a/tests/models/test_stranding.py +++ b/tests/models/test_stranding.py @@ -57,10 +57,10 @@ def test_stranding_3d(self): # Check calculated trajectory lengths and speeds total_length, distances, speeds = o.get_trajectory_lengths() - self.assertAlmostEqual(total_length.max(), 15719.5, 1) + self.assertAlmostEqual(total_length.max(), 15722.7, 1) self.assertAlmostEqual(total_length.min(), 1225.2, 1) self.assertAlmostEqual(speeds.max(), 0.132, 1) - self.assertAlmostEqual(distances.max(), 2859.5, 1) + self.assertAlmostEqual(distances.max(), 2859.0, 1) def test_stranding_roms(self): o = PelagicEggDrift(loglevel=20) diff --git a/tests/readers/test_depth.py b/tests/readers/test_depth.py index ce6765920..37b82a07f 100644 --- a/tests/readers/test_depth.py +++ b/tests/readers/test_depth.py @@ -25,7 +25,7 @@ def test_sdepth_zeta_zero_Vtransform1(self): expected_output = np.array([-0.98333333, -0.5, -0.01666667]) assert z_rho.shape == (self.NZ, self.NY, self.NX) - np.allclose(z_rho[:,0,0], expected_output) + assert np.allclose(z_rho[:,0,0], expected_output) # zeta 0s array zeta = np.zeros((self.NY, self.NX)) @@ -33,7 +33,7 @@ def test_sdepth_zeta_zero_Vtransform1(self): expected_output = np.array([-0.98333333, -0.5, -0.01666667]) assert z_rho.shape == (self.NZ, self.NY, self.NX) - np.allclose(z_rho[:,0,0], expected_output) + assert np.allclose(z_rho[:,0,0], expected_output) def test_sdepth_zeta_zero_Vtransform2(self): @@ -46,14 +46,14 @@ def test_sdepth_zeta_zero_Vtransform2(self): z_rho = depth.sdepth(self.Htot, zeta, self.hc, self.Cs_r, Vtransform=Vtransform) expected_output = np.array([-0.98484848, -0.5, -0.01515152]) assert z_rho.shape == (self.NZ, self.NY, self.NX) - np.allclose(z_rho[:,0,0], expected_output) + assert np.allclose(z_rho[:,0,0], expected_output) # zeta 0s array zeta = np.zeros((self.NY, self.NX)) z_rho = depth.sdepth(self.Htot, zeta, self.hc, self.Cs_r, Vtransform=Vtransform) expected_output = np.array([-0.98484848, -0.5, -0.01515152]) assert z_rho.shape == (self.NZ, self.NY, self.NX) - np.allclose(z_rho[:,0,0], expected_output) + assert np.allclose(z_rho[:,0,0], expected_output) def test_sdepth_zeta_nonzero_Vtransform1(self): """test nonzero zeta in sdepth, Vtransform 1.""" @@ -65,14 +65,14 @@ def test_sdepth_zeta_nonzero_Vtransform1(self): z_rho = depth.sdepth(self.Htot, zeta, self.hc, self.Cs_r, Vtransform=Vtransform) expected_output = np.array([-0.96666667, 0., 0.96666667]) assert z_rho.shape == (self.NZ, self.NY, self.NX) - np.allclose(z_rho[:,0,0], expected_output) + assert np.allclose(z_rho[:,0,0], expected_output) # zeta 1s array zeta = np.ones((self.NY, self.NX)) z_rho = depth.sdepth(self.Htot, zeta, self.hc, self.Cs_r, Vtransform=Vtransform) expected_output = np.array([-0.96666667, 0., 0.96666667]) assert z_rho.shape == (self.NZ, self.NY, self.NX) - np.allclose(z_rho[:,0,0], expected_output) + assert np.allclose(z_rho[:,0,0], expected_output) def test_sdepth_zeta_nonzero_Vtransform2(self): @@ -85,14 +85,14 @@ def test_sdepth_zeta_nonzero_Vtransform2(self): z_rho = depth.sdepth(self.Htot, zeta, self.hc, self.Cs_r, Vtransform=Vtransform) expected_output = np.array([-0.96969697, 0., 0.96969697]) assert z_rho.shape == (self.NZ, self.NY, self.NX) - np.allclose(z_rho[:,0,0], expected_output) + assert np.allclose(z_rho[:,0,0], expected_output) # zeta 1s array zeta = np.ones((self.NY, self.NX)) z_rho = depth.sdepth(self.Htot, zeta, self.hc, self.Cs_r, Vtransform=Vtransform) expected_output = np.array([-0.96969697, 0., 0.96969697]) assert z_rho.shape == (self.NZ, self.NY, self.NX) - np.allclose(z_rho[:,0,0], expected_output) + assert np.allclose(z_rho[:,0,0], expected_output) diff --git a/tests/readers/test_netcdf.py b/tests/readers/test_netcdf.py index 9a85b5e1a..a25590ca0 100644 --- a/tests/readers/test_netcdf.py +++ b/tests/readers/test_netcdf.py @@ -38,6 +38,10 @@ def test_MFDataset(self): '2Feb2016_Nordic_sigma_3d/Nordic_subset.nc') lon = 14.0 lat = 67.3 + + # nordicMF is missing angle + nordicMF.Dataset["angle"] = nordicMF_all.Dataset.angle + #nordic3d = reader_ROMS_native.Reader(o.test_data_folder() + # '2Feb2016_Nordic_sigma_3d/Nordic-4km_SLEVELS_avg_00_subset2Feb2016.nc') #o.add_reader(nordic3d) # Slightly different results diff --git a/tests/readers/test_roms.py b/tests/readers/test_roms.py index 10a1edf06..5413ee434 100644 --- a/tests/readers/test_roms.py +++ b/tests/readers/test_roms.py @@ -9,8 +9,7 @@ class TestROMSReader(unittest.TestCase): def setUp(self): - zeta = np.array([[[0.1, 0.2, 0.3], [0.4, 0.5, 0.6]], [[0.7, 0.8, 0.9], [1.0, 1.1, 1.2]]]) - # zeta = np.zeros((2,2,3)) + h = np.ones((2, 3))*10 u_eastward = np.array([[[[1.0, 1.1, 1.2], [1.3, 1.4, 1.5]], [[1.6, 1.7, 1.8], [1.9, 2.0, 2.1]], [[2.2, 2.3, 2.4], [2.5, 2.6, 2.7]]], [[[2.8, 2.9, 3.0], [3.1, 3.2, 3.3]], [[3.4, 3.5, 3.6], [3.7, 3.8, 3.9]], [[4.0, 4.1, 4.2], [4.3, 4.4, 4.5]]]]) @@ -33,12 +32,13 @@ def setUp(self): "lon_v": (["eta_v", "xi_v"], [[-152.0, -151.9, -151.8]]), "lat_v": (["eta_v", "xi_v"], [[60.45, 60.45, 60.45]]), "mask_rho": (["eta_rho", "xi_rho"], mask_rho), + "h": (["eta_rho", "xi_rho"], h), "wetdry_mask_rho": (["ocean_time", "eta_rho", "xi_rho"], wetdry_mask_rho), "angle": (["eta_rho", "xi_rho"], angle), "mask_u": (["eta_u", "xi_u"], mask_u), "mask_v": (["eta_v", "xi_v"], mask_v), "ocean_time": ("ocean_time", [0, 1], {"units": "seconds since 1970-01-01"}), - "zeta": (["ocean_time", "eta_rho", "xi_rho"], zeta), + # "zeta": (["ocean_time", "eta_rho", "xi_rho"], zeta), "s_rho": (["s_rho"], Cs_r), "Cs_r": (["s_rho"], Cs_r), "hc": 16, @@ -123,6 +123,42 @@ def test_get_variables_u_mask_u(self): assert "mask_rho" not in masks + 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) + # 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.zeros((2,2,3)) + reader.Dataset["zeta"] = (["ocean_time", "eta_rho", "xi_rho"], zeta) + + output = reader.get_variables(key, datetime(1970, 1, 1), 0, 0, [-5], testing=False) + + assert key in output + assert reader.do_not_rotate == [] + # those expected values at 0.5, 0.6 from the input array + assert np.allclose(output[key].data[1,0,:], self.ds[var][0,1,0,:].values) + + + def test_get_variables_u_depth_zeta_nonzero(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) + # 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 == [] + assert np.allclose(output[key].data[1,0,:], u_expected) + + class TestROMSReaderRotation(unittest.TestCase): def setUp(self):