Source code for opendrift.models.basemodel.environment

import logging
from typing import OrderedDict, Dict, List
import traceback
import numpy as np
import pyproj

from opendrift.timer import Timeable
from opendrift.config import CONFIG_LEVEL_BASIC, CONFIG_LEVEL_ADVANCED
from opendrift.readers.basereader import BaseReader, standard_names
from opendrift.readers.reader_lazy import Reader
from opendrift.readers import reader_from_url, reader_global_landmask
from opendrift.errors import NotCoveredError
from opendrift.models import physics_methods as pm

from opendrift.config import Configurable

logger = logging.getLogger(__name__)


[docs] class Environment(Timeable, Configurable): readers: OrderedDict priority_list: OrderedDict required_variables: Dict discarded_readers: Dict proj_latlon = pyproj.Proj('+proj=latlong') __finalized__ = False def __init__(self, required_variables, _config): super().__init__() self.readers = OrderedDict() self.priority_list = OrderedDict() self.discarded_readers = {} self.required_variables = required_variables self._config = _config # reference to simulation config # Add constant and fallback environment variables to config c = {} for v in self.required_variables: minval = maxval = units = None description_constant = 'Use constant value for %s' % v description_fallback = 'Fallback value for %s if not available from any reader' % v if v in standard_names: if 'valid_min' in standard_names[v]: minval = standard_names[v]['valid_min'] if 'valid_max' in standard_names[v]: maxval = standard_names[v]['valid_max'] if 'long_name' in standard_names[v]: description_constant = description_fallback = standard_names[ v]['long_name'] if 'units' in standard_names[v]: units = standard_names[v]['units'] c['environment:constant:%s' % v] = { 'type': 'float', 'min': minval, 'max': maxval, 'units': units, 'default': None, 'level': CONFIG_LEVEL_BASIC, 'description': description_constant } c['environment:fallback:%s' % v] = { 'type': 'float', 'min': minval, 'max': maxval, 'units': units, 'default': self.required_variables[v]['fallback'] if 'fallback' in self.required_variables[v] else None, 'level': CONFIG_LEVEL_BASIC, 'description': description_fallback } self._add_config(c) self._add_config({ 'general:use_auto_landmask': { 'type': 'bool', 'default': True, 'description': 'A built-in GSHHG global landmask is used if True, ' 'otherwise landmask is taken from reader or fallback value.', 'level': CONFIG_LEVEL_ADVANCED }, 'drift:current_uncertainty': { 'type': 'float', 'default': 0, 'min': 0, 'max': 5, 'units': 'm/s', 'description': 'Add gaussian perturbation with this standard deviation to current components at each time step', 'level': CONFIG_LEVEL_ADVANCED }, 'drift:current_uncertainty_uniform': { 'type': 'float', 'default': 0, 'min': 0, 'max': 5, 'units': 'm/s', 'description': 'Add gaussian perturbation with this standard deviation to current components at each time step', 'level': CONFIG_LEVEL_ADVANCED }, 'drift:max_speed': { 'type': 'float', 'default': 1, 'min': 0, 'max': np.inf, 'units': 'seconds', 'description': 'Typical maximum speed of elements, used to estimate reader buffer size', 'level': CONFIG_LEVEL_ADVANCED }, 'readers:max_number_of_fails': { 'type': 'int', 'default': 1, 'min': 0, 'max': 1e6, 'units': 'number', 'description': 'Readers are discarded if they fail (e.g. corrupted data, og hanging servers) move than this number of times', 'level': CONFIG_LEVEL_ADVANCED }, }) # Find variables which require profiles self.required_profiles = [ var for var in self.required_variables if 'profiles' in self.required_variables[var] and self.required_variables[var]['profiles'] is True ] # Find variables which are desired, but not required self.desired_variables = [ var for var in self.required_variables if 'important' in self.required_variables[var] and self.required_variables[var]['important'] is False ]
[docs] def finalize(self, simulation_extent=None, start=None, end=None): """ Prepare environment for simulation. Args: simulation_extent: The expected extent of the simulation. start: Expected start time of simulation. end: Expected end time of simulation. """ if self.__finalized__ is False: # TODO: discard irrelevant readers self.__generate_constant_readers__() self.__add_auto_landmask__() self.__assert_no_missing_variables__() if simulation_extent is not None: self.prepare_readers(simulation_extent, start, end) self.__finalized__ = True
[docs] def prepare_readers(self, extent, start_time, end_time): if extent is not None: self.simulation_extent = extent for reader in self.readers.values(): logger.debug(f'\tPreparing {reader.name} for extent {extent}') reader.prepare(extent=extent, start_time=start_time, end_time=end_time, max_speed=self.get_config('drift:max_speed'))
[docs] def __generate_constant_readers__(self): # Make constant readers if config environment:constant:<var> is set c = self.get_configspec('environment:constant:') mr = {} for var in list(c): if c[var]['value'] is not None: mr[var.split(':')[-1]] = c[var]['value'] if len(mr) > 0: from opendrift.readers import reader_constant rc = reader_constant.Reader(mr) self.add_reader(rc, first=True)
[docs] def __add_auto_landmask__(self): ############################################################## # If no landmask has been added, we determine it dynamically ############################################################## # TODO: some more error checking here # If landmask is requested, it shall not be obtained from other readers if self.get_config('general:use_auto_landmask', False) is True: if 'land_binary_mask' in self.priority_list: if 'global_landmask' in self.priority_list['land_binary_mask']: self.priority_list['land_binary_mask'] = [ 'global_landmask' ] else: if self.get_config('environment:constant:land_binary_mask') is None: del self.priority_list['land_binary_mask'] if self.get_config('general:use_auto_landmask', False) is True and \ ('land_binary_mask' in self.required_variables and \ 'land_binary_mask' not in self.priority_list): if self.get_config('environment:constant:land_binary_mask') is None: logger.info( 'Adding a dynamical landmask with max. priority based on ' 'assumed maximum speed of %s m/s. ' 'Adding a customised landmask may be faster...' % self.get_config('drift:max_speed')) self.timer_start('preparing main loop:making dynamical landmask') reader_landmask = reader_global_landmask.Reader() self.add_reader(reader_landmask) self.timer_end('preparing main loop:making dynamical landmask') else: logger.warning('Using constant reader for land_binary_mask, ' 'although config general:use_auto_landmask is True')
[docs] def __assert_no_missing_variables__(self): missing_variables = self.missing_variables() if len(missing_variables) > 0: has_fallback = [ var for var in missing_variables if (self.get_config(f'environment:fallback:{var}') is not None or self.get_config(f'environment:constant:{var}') is not None) ] has_no_fallback = [ var for var in missing_variables if (self.get_config(f'environment:fallback:{var}') is None and self.get_config(f'environment:constant:{var}') is None) ] if len(has_fallback) > 0: # == missing_variables: logger.info('Fallback values will be used for the following ' 'variables which have no readers: ') for var in has_fallback: logger.info('\t%s: %f' % (var, self.get_config(f'environment:fallback:{var}'))) if len(has_no_fallback) > 0 and len( self._lazy_readers()) == 0: # == missing_variables: logger.warning( 'No readers added for the following variables: ' + str(has_no_fallback)) raise ValueError('Readers must be added for the ' 'following required variables: ' + str(has_no_fallback))
[docs] def add_readers_from_file(self, filename, timeout=10, lazy=True): fp = open(filename, 'r') sources = fp.readlines() sources = [line.strip() for line in sources if line[0] != '#'] self.add_readers_from_list(sources, timeout, lazy=lazy)
[docs] def add_readers_from_list(self, urls, timeout=10, lazy=True, variables=None): '''Make readers from a list of URLs or paths to netCDF datasets''' if isinstance(urls, str): urls = [urls] if lazy is True: from opendrift.readers.reader_lazy import Reader readers = [Reader(u) for u in urls] self.add_reader(readers, variables=variables) return readers = [reader_from_url(u, timeout) for u in urls] self.add_reader([r for r in readers if r is not None], variables=variables)
[docs] def add_reader(self, readers, variables=None, first=False): """Add one or more readers providing variables used by this model. Method may be called subsequently to add more readers for other variables. Args: readers: one or more (list) Reader objects. variables (optional): list of strings of standard_name of variables to be provided by this/these reader(s). first: Set to True if this reader should be set as first option """ # Convert any strings to lists, for looping if isinstance(variables, str): variables = [variables] if isinstance(readers, BaseReader): readers = [readers] if isinstance(readers, Reader): readers = [readers] for reader in readers: # Check if input class is of correct type if not isinstance(reader, BaseReader) and \ not hasattr(reader, '_lazyname'): raise TypeError('Please provide Reader object') # Check that reader class contains the requested variables if variables is not None: missingVariables = set(variables) - set(reader.variables) if missingVariables: raise ValueError( 'Reader %s does not provide variables: %s' % (reader.name, list(missingVariables))) # Finally add new reader to list if reader.name in self.readers: # Reader names must be unique, adding integer for n in range(99999): tmp_name = reader.name + '_%d' % n if tmp_name not in self.readers: reader.name = tmp_name break self.readers[reader.name] = reader logger.debug('Added reader ' + reader.name) # Add this reader for each of the given variables if reader.is_lazy is False: for variable in variables if variables else reader.variables: if variable in list(self.priority_list): if reader.name not in self.priority_list[variable]: if first is True: self.priority_list[variable].insert( 0, reader.name) else: self.priority_list[variable].append( reader.name) else: self.priority_list[variable] = [reader.name] # Remove/hide variables not needed by the current trajectory model for variable in list(self.priority_list): if variable not in self.required_variables: del self.priority_list[variable]
[docs] def list_environment_variables(self): """Return list of all variables provided by the added readers.""" variables = [] for reader in self.readers: variables.extend(self.readers[reader].variables) return variables
[docs] def get_reader_groups(self, variables=None): """Find which groups of variables are provided by the same readers. This function loops through 'priority_list' (see above) and groups all variables returned by the same readers in the same order. This allows asking readers for several variables simultaneously, improving performance. Used by method 'get_environment'. Returns: variable_groups: list of lists of (environment) variables. reader_groups: list of list of reader names, corresponding to each of the variable_groups. """ if variables is None: variables = list(self.required_variables) reader_groups = [] # Find all unique reader groups for variable, readers in self.priority_list.items(): if (variable in variables) and (readers not in reader_groups): reader_groups.append(readers) # Find all variables returned by the same reader group variable_groups = [None] * len(reader_groups) for variable, readers in self.priority_list.items(): if variable not in variables: continue for i, readerGroup in enumerate(reader_groups): if readers == readerGroup: if variable_groups[i]: variable_groups[i].append(variable) else: variable_groups[i] = [variable] missing_variables = list( set(variables) - set(self.priority_list.keys())) return variable_groups, reader_groups, missing_variables
[docs] def _lazy_readers(self): return [r for r in self.readers if self.readers[r].is_lazy is True]
[docs] def _unlazy_readers(self): return [r for r in self.readers if self.readers[r].is_lazy is False]
[docs] def _initialise_next_lazy_reader(self): '''Returns reader if successful and None if no more readers''' lazy_readers = self._lazy_readers() if len(lazy_readers) == 0: return None lazyname = lazy_readers[0] reader = self.readers[lazyname] try: reader.initialise() except Exception as e: logger.debug(e) logger.warning('Reader could not be initialised, and is' ' discarded: ' + lazyname) self.discard_reader(reader, reason='could not be initialized') return self._initialise_next_lazy_reader() # Call self reader.set_buffer_size(max_speed=self.get_config('drift:max_speed')) # Update reader lazy name with actual name self.readers[reader.name] = \ self.readers.pop(lazyname) for var in reader.variables: if var in list(self.priority_list): self.priority_list[var].append(reader.name) else: self.priority_list[var] = [reader.name] # Remove variables not needed for variable in list(self.priority_list): if variable not in self.required_variables: del self.priority_list[variable] return reader
[docs] def discard_reader_if_not_relevant(self, reader, time): if reader.is_lazy: return False if reader.start_time is not None and reader.always_valid is False: if hasattr(self, 'expected_end_time' ) and reader.start_time > self.expected_end_time: self.discard_reader(reader, 'starts after simulation end') return True if hasattr(self, 'start_time') and reader.end_time < self.start_time: self.discard_reader(reader, 'ends before simuation start') return True if time is not None and reader.end_time < time: self.discard_reader(reader, 'ends before simuation is finished') return True if len(set(self.required_variables) & set(reader.variables)) == 0: self.discard_reader( reader, reason='does not contain any relevant variables') return True if not hasattr(reader, 'checked_for_overlap'): if not reader.global_coverage(): if not hasattr(self, 'simulation_extent'): logger.warning( 'Simulation has no simulation_extent, cannot check reader coverage' ) return False # TODO # need a better coverage/overlap check below corners = reader.xy2lonlat( [reader.xmin, reader.xmin, reader.xmax, reader.xmax], [reader.ymax, reader.ymin, reader.ymax, reader.ymin]) rlonmin = np.min(corners[0]) rlonmax = np.max(corners[0]) rlatmin = np.min(corners[1]) rlatmax = np.max(corners[1]) if hasattr( reader, 'proj4' ) and 'stere' in reader.proj4 and 'lat_0=90' in reader.proj4: rlatmax = 90 if hasattr( reader, 'proj4' ) and 'stere' in reader.proj4 and 'lat_0=-90' in reader.proj4: rlatmin = -90 if rlatmin > self.simulation_extent[3]: self.discard_reader(reader, reason='too far north') return True if rlatmax < self.simulation_extent[1]: self.discard_reader(reader, reason='too far south') return True # Disabling below checks, as +/-360 deg is not considered #if rlonmax < self.simulation_extent[0]: # self.discard_reader(reader, reason='too far west') # return True #if rlonmin > self.simulation_extent[2]: # self.discard_reader(reader, reason='too far east') # return True reader.checked_for_overlap = True return False # Reader is not discarded
[docs] def discard_reader(self, reader, reason): readername = reader.name logger.debug('Discarding reader (%s): %s' % (reason, readername)) del self.readers[readername] self.discarded_readers[readername] = reason # Remove from priority list for var in self.priority_list.copy(): self.priority_list[var] = [ r for r in self.priority_list[var] if r != readername ] if len(self.priority_list[var]) == 0: del self.priority_list[var]
[docs] def missing_variables(self): """Return list of all variables for which no reader has been added.""" return [ var for var in self.required_variables if var not in self.priority_list ]
[docs] def get_environment(self, variables, time, lon, lat, z, profiles=None, profiles_depth=None): '''Retrieve environmental variables at requested positions. Args: variables: list of variable names time: time to get environment for lon: array of longitudes lat: array of latitudes z: depth to get value for profiles: list of variables for which profiles are needed profiles_depth: depth of profiles in meters, as a positive number Updates: Buffer (raw data blocks) for each reader stored for performance: [readers].var_block_before (last before requested time) [readers].var_block_after (first after requested time) - lists of one ReaderBlock per variable group: time, x, y, [vars] Returns: environment: recarray with variables as named attributes, interpolated to requested positions/time. ''' assert self.__finalized__ is True, 'The environment has not been finalized.' self.timer_start('main loop:readers') # Initialise ndarray to hold environment variables dtype = [(var, np.float32) for var in variables] env = np.ma.array(np.zeros(len(lon)) * np.nan, dtype=dtype) num_elements_active = len(lon) # Discard any existing readers which are not relevant for readername, reader in self.readers.copy().items(): self.discard_reader_if_not_relevant(reader, time) if profiles_depth is None: profiles_depth = np.abs(z).max() if 'drift:truncate_ocean_model_below_m' in self._config: truncate_depth = self.get_config( 'drift:truncate_ocean_model_below_m') if truncate_depth is not None: logger.debug('Truncating ocean models below %s m' % truncate_depth) z = z.copy() z[z < -truncate_depth] = -truncate_depth profiles_depth = np.minimum(profiles_depth, truncate_depth) # Initialise more lazy readers if necessary missing_variables = ['missingvar'] while (len(missing_variables) > 0 and len(self._lazy_readers()) > 0): variable_groups, reader_groups, missing_variables = \ self.get_reader_groups(variables) if hasattr(self, 'desired_variables'): missing_variables = list( set(missing_variables) - set(self.desired_variables)) if len(missing_variables) > 0: logger.debug('Variables not covered by any reader: ' + str(missing_variables)) reader = 'NotNone' while reader is not None: reader = self._initialise_next_lazy_reader() if reader is not None: if self.discard_reader_if_not_relevant(reader, time): reader = None if reader is not None: if (reader.covers_time(time) and len( reader.covers_positions(lon, lat)[0]) > 0): missing_variables = list( set(missing_variables) - set(reader.variables)) if len(missing_variables) == 0: break # We cover now all variables # For each variable/reader group: variable_groups, reader_groups, missing_variables = \ self.get_reader_groups(variables) for variable in variables: # Fill with fallback value if no reader co = self.get_config('environment:fallback:%s' % variable) if co is not None: env[variable] = np.ma.ones(env[variable].shape) * co for variable_group, reader_group in zip(variable_groups, reader_groups): logger.debug('----------------------------------------') logger.debug('Variable group %s' % (str(variable_group))) logger.debug('----------------------------------------') missing_indices = np.array(range(len(lon))) # For each reader: for reader_name in reader_group: logger.debug('Calling reader ' + reader_name) logger.debug('----------------------------------------') self.timer_start('main loop:readers:' + reader_name.replace(':', '<colon>')) reader = self.readers[reader_name] if not reader.covers_time(time): logger.debug('\tOutside time coverage of reader.') if reader_name == reader_group[-1]: if self._initialise_next_lazy_reader() is not None: logger.debug( 'Missing variables: calling get_environment recursively' ) return self.get_environment( variables, time, lon, lat, z, profiles, profiles_depth) continue # Fetch given variables at given positions from current reader try: logger.debug('Data needed for %i elements' % len(missing_indices)) # Check if vertical profiles are requested from reader if profiles is not None: profiles_from_reader = list( set(variable_group) & set(profiles)) if profiles_from_reader == []: profiles_from_reader = None else: profiles_from_reader = None env_tmp, env_profiles_tmp = \ reader.get_variables_interpolated( variable_group, profiles = profiles_from_reader, profiles_depth = profiles_depth, time = time, lon=lon[missing_indices], lat=lat[missing_indices], z=z[missing_indices], rotate_to_proj=self.proj_latlon) except NotCoveredError as e: logger.info(e) self.timer_end('main loop:readers:' + reader_name.replace(':', '<colon>')) if reader_name == reader_group[-1]: if self._initialise_next_lazy_reader() is not None: logger.debug( 'Missing variables: calling get_environment recursively' ) return self.get_environment( variables, time, lon, lat, z, profiles, profiles_depth) continue except Exception as e: # Unknown error # TODO: # This could e.g. be due to corrupted files or # hangig thredds-servers. A reader could be discarded # after e.g. 3 such failed attempts logger.info('========================') logger.exception(e) logger.debug(traceback.format_exc()) logger.info('========================') reader.number_of_fails = reader.number_of_fails + 1 max_fails = self.get_config('readers:max_number_of_fails') if reader.number_of_fails > max_fails: logger.warning( f'Reader {reader.name} is discarded after failing ' f'more times than allowed ({max_fails})') self.discard_reader( reader, reason=f'failed more than {max_fails} times') self.timer_end('main loop:readers:' + reader_name.replace(':', '<colon>')) if reader_name == reader_group[-1]: if self._initialise_next_lazy_reader() is not None: logger.debug( 'Missing variables: calling get_environment recursively' ) return self.get_environment( variables, time, lon, lat, z, profiles, profiles_depth) continue # Copy retrieved variables to env array, and mask nan-values for var in variable_group: if var not in self.required_variables: logger.debug('Not returning env-variable: ' + var) continue if var not in env.dtype.names: continue # Skipping variables that are only used to derive needed variables env[var][missing_indices] = np.ma.masked_invalid( env_tmp[var][0:len(missing_indices)]).astype('float32') if profiles_from_reader is not None and var in profiles_from_reader: if 'env_profiles' not in locals(): env_profiles = env_profiles_tmp # TODO: fix to be checked if var in env_profiles and var in env_profiles_tmp: # If one profile has fewer vertical layers than # the other, we use only the overlapping part if len(env_profiles['z']) != len( env_profiles_tmp['z']): logger.debug('Warning: different number of ' ' vertical layers: %s and %s' % (len(env_profiles['z']), len(env_profiles_tmp['z']))) z_ind = np.arange( np.minimum( len(env_profiles['z']) - 1, len(env_profiles_tmp['z']) - 1)) # len(missing_indices) since 2 points might have been added and not removed env_profiles_tmp[var] = np.ma.atleast_2d( env_profiles_tmp[var]) env_profiles[var][np.ix_(z_ind, missing_indices)] = \ np.ma.masked_invalid(env_profiles_tmp[var][z_ind,0:len(missing_indices)]).astype('float32') # For profiles with different numbers of layers, we extrapolate if env_profiles[var].shape[0] > 1: missingbottom = np.isnan( env_profiles[var][-1, :]) env_profiles[var][ -1, missingbottom] = env_profiles[var][ -2, missingbottom] # Detect elements with missing data, for present reader group if not hasattr(env_tmp[variable_group[0]], 'mask'): env_tmp[variable_group[0]] = np.ma.masked_invalid( env_tmp[variable_group[0]]) try: del combined_mask except: pass for var in variable_group: tmp_var = np.ma.masked_invalid(env_tmp[var]) # Changed 13 Oct 2016, but uncertain of effect # TODO: to be checked #tmp_var = env_tmp[var] if 'combined_mask' not in locals(): combined_mask = np.ma.getmask(tmp_var) else: combined_mask = \ np.ma.mask_or(combined_mask, np.ma.getmask(tmp_var), shrink=False) try: if len(missing_indices) != len(combined_mask): # TODO: mask mismatch due to 2 added points raise ValueError('Mismatch of masks') missing_indices = missing_indices[combined_mask] except Exception as ex: # Not sure what is happening here logger.info('Problems setting mask on missing_indices!') logger.exception(ex) if (type(missing_indices) == np.int64) or (type(missing_indices) == np.int32): missing_indices = [] self.timer_end('main loop:readers:' + reader_name.replace(':', '<colon>')) if len(missing_indices) == 0: logger.debug('Obtained data for all elements.') break else: logger.debug('Data missing for %i elements.' % (len(missing_indices))) if len(self._lazy_readers()) > 0: if self._initialise_next_lazy_reader() is not None: logger.warning( 'Missing variables: calling get_environment recursively' ) return self.get_environment( variables, time, lon, lat, z, profiles) logger.debug('---------------------------------------') logger.debug('Finished processing all variable groups') self.timer_start('main loop:readers:postprocessing') #for var in self.fallback_values: # if (var not in variables) and (profiles is None # or var not in profiles): # continue for var in variables: if self.get_config(f'environment:fallback:{var}') is None: continue mask = env[var].mask fallback = self.get_config(f'environment:fallback:{var}') if any(mask == True): logger.debug( ' Using fallback value %s for %s for %s elements' % (fallback, var, np.sum(mask == True))) env[var][mask] = fallback # Profiles if profiles is not None and var in profiles: if 'env_profiles' not in locals(): logger.debug('Creating empty dictionary for profiles not ' 'profided by any reader: ' + str(self.required_profiles)) env_profiles = {'z': [0, -profiles_depth]} if var not in env_profiles: logger.debug( ' Using fallback value %s for %s for all profiles' % (fallback, var)) env_profiles[var] = fallback*\ np.ma.ones((len(env_profiles['z']), num_elements_active)) else: mask = env_profiles[var].mask num_masked_values_per_element = np.sum(mask == True) num_missing_profiles = np.sum(num_masked_values_per_element == len(env_profiles['z'])) env_profiles[var][mask] = fallback logger.debug( ' Using fallback value %s for %s for %s profiles' % ( fallback, var, num_missing_profiles, )) num_missing_individual = np.sum( num_masked_values_per_element > 0) - num_missing_profiles if num_missing_individual > 0: logger.debug( ' ...plus %s individual points in other profiles' % num_missing_individual) ####################################################### # Some extra checks of units and realistic magnitude ####################################################### if 'sea_water_temperature' in variables: t_kelvin = np.where(env['sea_water_temperature'] > 100)[0] if len(t_kelvin) > 0: logger.warning( 'Converting temperatures from Kelvin to Celcius') env['sea_water_temperature'][ t_kelvin] = env['sea_water_temperature'][t_kelvin] - 273.15 if 'env_profiles' in locals( ) and 'sea_water_temperature' in env_profiles.keys(): env_profiles['sea_water_temperature'][:,t_kelvin] = \ env_profiles['sea_water_temperature'][:,t_kelvin] - 273.15 ############################################################ # Parameterisation of unavailable variables # TODO: use instead "environment mapping" mechanism for this ############################################################# if 'drift:use_tabularised_stokes_drift' in self._config and self.get_config( 'drift:use_tabularised_stokes_drift') is True: if 'x_wind' in variables: if 'sea_surface_wave_stokes_drift_x_velocity' not in variables or ( env['sea_surface_wave_stokes_drift_x_velocity'].max() == 0 and env['sea_surface_wave_stokes_drift_y_velocity'].max() == 0): logger.debug('Calculating parameterised stokes drift') env['sea_surface_wave_stokes_drift_x_velocity'], \ env['sea_surface_wave_stokes_drift_y_velocity'] = \ pm.wave_stokes_drift_parameterised((env['x_wind'], env['y_wind']), self.get_config('drift:tabularised_stokes_drift_fetch')) if (env['sea_surface_wave_significant_height'].max() == 0): logger.debug( 'Calculating parameterised significant wave height') env['sea_surface_wave_significant_height'] = \ pm.wave_significant_height_parameterised((env['x_wind'], env['y_wind']), self.get_config('drift:tabularised_stokes_drift_fetch')) ############################# # Add uncertainty/diffusion ############################# # Current if 'x_sea_water_velocity' in variables and \ 'y_sea_water_velocity' in variables: std = self.get_config('drift:current_uncertainty') if std > 0: logger.debug('Adding uncertainty for current: %s m/s' % std) env['x_sea_water_velocity'] += np.random.normal( 0, std, num_elements_active) env['y_sea_water_velocity'] += np.random.normal( 0, std, num_elements_active) std = self.get_config('drift:current_uncertainty_uniform') if std > 0: logger.debug('Adding uncertainty for current: %s m/s' % std) env['x_sea_water_velocity'] += np.random.uniform( -std, std, num_elements_active) env['y_sea_water_velocity'] += np.random.uniform( -std, std, num_elements_active) # Wind if 'x_wind' in variables and 'y_wind' in variables: std = self.get_config('drift:wind_uncertainty') if std > 0: logger.debug('Adding uncertainty for wind: %s m/s' % std) env['x_wind'] += np.random.normal(0, std, num_elements_active) env['y_wind'] += np.random.normal(0, std, num_elements_active) ##################### # Diagnostic output ##################### if len(env) > 0: logger.debug('------------ SUMMARY -------------') for var in variables: logger.debug(' %s: %g (min) %g (max)' % (var, env[var].min(), env[var].max())) logger.debug('---------------------------------') # Prepare array indiciating which elements contain any invalid values missing = np.ma.masked_invalid(env[variables[0]]).mask for var in variables[1:]: missing = np.ma.mask_or(missing, np.ma.masked_invalid(env[var]).mask, shrink=False) # Convert dictionary to recarray and return if 'env_profiles' not in locals(): env_profiles = None # Convert masked arrays to regular arrays for increased performance env = np.array(env) if env_profiles is not None: for var in env_profiles: env_profiles[var] = np.array(env_profiles[var]) self.timer_end('main loop:readers:postprocessing') self.timer_end('main loop:readers') return env.view(np.recarray), env_profiles, missing
[docs] def get_variables_along_trajectory(self, variables, lons, lats, times, z=0): self.finalize() data = {'time': times, 'lon': lons, 'lat': lats} for var in variables: data[var] = np.zeros(len(times)) for i, time in enumerate(times): self.time = time d = self.get_environment(lon=np.atleast_1d(lons[i]), lat=np.atleast_1d(lats[i]), z=np.atleast_1d(z), time=time, variables=variables, profiles=None) for var in variables: data[var][i] = d[0][var][0] return data