Source code for opendrift.models.basemodel

# This file is part of OpenDrift.
#
# OpenDrift is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, version 2
#
# OpenDrift is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with OpenDrift.  If not, see <https://www.gnu.org/licenses/>.
#
# Copyright 2015, Knut-Frode Dagestad, MET Norway
# Copyright 2020, Gaute Hope, MET Norway

import sys
import os
import glob
import types
import traceback
import inspect
import logging; logging.captureWarnings(True); logger = logging.getLogger(__name__)
import warnings
from datetime import datetime, timedelta
from collections import OrderedDict
from abc import ABCMeta, abstractmethod, abstractproperty
import geojson
import netCDF4
import nc_time_axis
import xarray as xr
from netCDF4 import Dataset, date2num

import numpy as np
import scipy
import pyproj
try:
    import matplotlib
    matplotlib.rcParams['legend.numpoints'] = 1
    matplotlib.rcParams['legend.scatterpoints'] = 1
    if ('DISPLAY' not in os.environ and
            'PYCHARM_HOSTED' not in os.environ and
            os.name != 'nt'):
        logger.info('No display found. Using non-interactive Agg backend')
        matplotlib.use('agg')
    import matplotlib.pyplot as plt
    from matplotlib import animation
    from matplotlib.patches import Polygon
    from matplotlib.path import Path
    import cartopy
    import cartopy.crs as ccrs
    import cartopy.feature as cfeature
except ImportError:
    print('matplotlib and/or cartopy is not available, can not make plots')

import opendrift
from opendrift.timer import Timeable
from opendrift.readers.basereader import BaseReader, vector_pairs_xy, standard_names
from opendrift.readers import reader_from_url, reader_global_landmask
from opendrift.models.physics_methods import PhysicsMethods

[docs]class OpenDriftSimulation(PhysicsMethods, Timeable): """Generic trajectory model class, to be extended (subclassed). This as an Abstract Base Class, meaning that only subclasses can be initiated and used. Any specific subclass ('model') must contain its own (or shared) specific type of particles (ElementType), whose properties are updated at each time_step using method update() on basis of model physics/chemistry/biology and 'required_variables' (environment) which are provided by one or more Reader objects. Attributes: ElementType: the type (class) of particles to be used by this model elements: object of the class ElementType, storing the specific particle properties (ndarrays and scalars) of all active particles as named attributes. Elements are added by seeding-functions (presently only one implemented: seed_elements). elements_deactivated: ElementType object containing particles which have been deactivated (and removed from 'elements') elements_scheduled: ElementType object containing particles which have been scheduled, but not yet activated required_variables: list of strings of CF standard_names which is needed by this model (update function) to update properties of particles ('elements') at each time_step. This core class has no required_elements, this is implemented by subclasses/modules. environment: recarray storing environment variables (wind, waves, current etc) as named attributes. Attribute names follow standard_name from CF-convention, allowing any OpenDriftSimulation module/subclass using environment data from any readers which can provide the requested variables. Used in method 'update' to update properties of elements every time_step. time_step: timedelta object, time interval at which element properties are updated (including advection). time_step_output: timedelta object, time interval at which element properties are stored in memory and eventually written to file readers: Dictionary where values are Reader objects, and names are unique reference keywords used to access a given reader (typically filename or URL) priority_list: OrderedDict where names are variable names, and values are lists of names (kewywords) of the reader, in the order of priority (user defined) of which the readers shall be called to retrieve environmental data. """ __metaclass__ = ABCMeta status_categories = ['active'] # Particles are active by default # Default plotting colors of trajectory endpoints status_colors_default = {'initial': 'green', 'active': 'blue', 'missing_data': 'gray'} CONFIG_LEVEL_ESSENTIAL=1 CONFIG_LEVEL_BASIC=2 CONFIG_LEVEL_ADVANCED=3 max_speed = 1 # Assumed max average speed of any element required_profiles_z_range = None # [min_depth, max_depth] plot_comparison_colors = ['k', 'r', 'g', 'b', 'm', 'c', 'y'] proj_latlon = pyproj.Proj('+proj=latlong')
[docs] @classmethod def SRS(cls): return cls.proj_latlon
def __init__(self, seed=0, iomodule='netcdf', loglevel=logging.DEBUG, logtime='%H:%M:%S', logfile=None): """Initialise OpenDriftSimulation Args: seed: integer or None. A given integer will yield identical random numbers drawn each simulation. Random numbers are e.g. used to distribute particles spatially when seeding, and may be used by modules (subclasses) for e.g. diffusion. Specifying a fixed value (default: 0) is useful for sensitivity tests. With seed = None, different random numbers will be drawn for subsequent runs, even with identical configuration/input. iomodule: name of module used to export data default: netcdf, see :py:mod:`opendrift.io` for more alternatives. `iomodule` is module/filename without preceeding `io_` loglevel: set to 0 (default) to retrieve all debug information. Provide a higher value (e.g. 20) to receive less output. Use the string 'custom' to configure logging from outside. logtime: if True, a time stamp is given for each logging line. logtime can also be given as a python time specifier (e.g. '%H:%M:%S') """ self.show_continuous_performance = False self.origin_marker = None # Dictionary to store named seeding locations self.minvals = {} # Dicionaries to store minimum and maximum values of variables self.maxvals = {} # List to store GeoJSON dicts of seeding commands self.seed_geojson = [] # Dict to store readers self.readers = OrderedDict() # Dictionary, key=name, value=reader object self.priority_list = OrderedDict() # Make copies of dictionaries so that they are private to each instance self.status_categories = ['active'] # Particles are active by default self.status_colors_default = self.status_colors_default.copy() if hasattr(self, 'status_colors'): # Append model specific colors to (and override) default colors self.status_colors_default.update(self.status_colors) self.status_colors = self.status_colors_default else: self.status_colors = self.status_colors_default # Using a fixed seed will generate the same random numbers # each run, useful for sensitivity tests # Use seed = None to get different random numbers each time np.random.seed(seed) self.steps_calculation = 0 # Increase for each simulation step self.steps_output = 0 self.elements_deactivated = self.ElementType() # Empty array self.elements = self.ElementType() # Empty array if loglevel != 'custom': format = '%(levelname)-7s %(name)s: %(message)s' datefmt = None if logtime is not False: format = '%(asctime)s ' + format if logtime is not True: datefmt = logtime formatter = logging.Formatter(format, datefmt=datefmt) if loglevel < 10: # 0 is NOTSET, giving no output loglevel=10 od_loggers = [ logging.getLogger('opendrift'), logging.getLogger('opendrift_landmask_data') ] if logfile is not None: handler = logging.FileHandler(logfile, mode='w') handler.setFormatter(formatter) for l in od_loggers: l.setLevel(loglevel) l.handlers = [] l.addHandler(handler) else: import coloredlogs fields = coloredlogs.DEFAULT_FIELD_STYLES fields['levelname']['color'] = 'magenta' # coloredlogs does not create duplicate handlers for l in od_loggers: coloredlogs.install(level = loglevel, logger = l, fmt = format, datefmt = datefmt, field_styles = fields) # Prepare outfile try: io_module = __import__('opendrift.export.io_' + iomodule, fromlist=['init', 'write_buffer', 'close', 'import_file']) except ImportError: logger.info('Could not import iomodule ' + iomodule) self.io_init = types.MethodType(io_module.init, self) self.io_write_buffer = types.MethodType(io_module.write_buffer, self) self.io_close = types.MethodType(io_module.close, self) self.io_import_file = types.MethodType(io_module.import_file, self) self.io_import_file_xarray = types.MethodType(io_module.import_file_xarray, self) # Set configuration options self._add_config({ # type, default, min, max, enum, important, value, units, description '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': self.CONFIG_LEVEL_ADVANCED}, 'general:coastline_action': {'type': 'enum', 'enum': ['none', 'stranding', 'previous'], 'default': 'stranding', 'level': self.CONFIG_LEVEL_BASIC, 'description': 'None means that objects may also move over land. ' 'stranding means that objects are deactivated if they hit land. ' 'previous means that objects will move back to the previous location ' 'if they hit land'}, 'general:time_step_minutes': {'type': 'float', 'min': .01, 'max': 1440, 'default': 60, 'units': 'minutes', 'level': self.CONFIG_LEVEL_BASIC, 'description': 'Calculation time step used for the simulation. The output time step may ' 'be equal or larger than this.'}, 'general:time_step_output_minutes': {'type': 'float', 'min': 1, 'max': 1440, 'default': None, 'units': 'minutes', 'level': self.CONFIG_LEVEL_BASIC, 'description': 'Output time step, i.e. the interval at which output is saved. This must be larger than ' 'the calculation time step, and be an integer multiple of this.'}, 'seed:ocean_only': {'type': 'bool', 'default': True, 'description': 'If True, elements seeded on land will be moved to the closest ' 'position in ocean', 'level': self.CONFIG_LEVEL_ADVANCED}, 'seed:number': {'type': 'int', 'default': 1, 'min': 1, 'max': 100000000, 'units': 1, 'description': 'The number of elements for the simulation.', 'level': self.CONFIG_LEVEL_BASIC}, 'drift:max_age_seconds': {'type': 'float', 'default': None, 'min': 0, 'max': np.inf, 'units': 'seconds', 'description': 'Elements will be deactivated when this age is reached', 'level': self.CONFIG_LEVEL_ADVANCED}, 'drift:advection_scheme': {'type': 'enum', 'enum': ['euler', 'runge-kutta', 'runge-kutta4'], 'default': 'euler', 'level': self.CONFIG_LEVEL_ADVANCED, 'description': 'Numerical advection scheme for ocean current advection'}, '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': self.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': self.CONFIG_LEVEL_ADVANCED}, 'drift:horizontal_diffusivity': {'type': 'float', 'default': 0, 'min': 0, 'max': 100000, 'units': 'm2/s', 'description': 'Add horizontal diffusivity (random walk)', 'level': self.CONFIG_LEVEL_ADVANCED}, 'drift:wind_uncertainty': {'type': 'float', 'default': 0, 'min': 0, 'max': 5, 'units': 'm/s', 'description': 'Add gaussian perturbation with this standard deviation to wind components at each time step.', 'level': self.CONFIG_LEVEL_ADVANCED}, 'drift:relative_wind': {'type': 'bool', 'default': False, 'description': 'If True, wind drift is calculated for absolute wind (wind vector minus ocean surface current vector).', 'level': self.CONFIG_LEVEL_ADVANCED}, 'drift:deactivate_north_of': {'type': 'float', 'default': None, 'min': -90, 'max': 90, 'units': 'degrees', 'description': 'Elements are deactivated if the move further north than this limit', 'level': self.CONFIG_LEVEL_ADVANCED}, 'drift:deactivate_south_of': {'type': 'float', 'default': None, 'min': -90, 'max': 90, 'units': 'degrees', 'description': 'Elements are deactivated if the move further south than this limit', 'level': self.CONFIG_LEVEL_ADVANCED}, 'drift:deactivate_east_of': {'type': 'float', 'default': None, 'min': -360, 'max': 360, 'units': 'degrees', 'description': 'Elements are deactivated if the move further east than this limit', 'level': self.CONFIG_LEVEL_ADVANCED}, 'drift:deactivate_west_of': {'type': 'float', 'default': None, 'min': -360, 'max': 360, 'units': 'degrees', 'description': 'Elements are deactivated if the move further west than this limit', 'level': self.CONFIG_LEVEL_ADVANCED}, }) # Add default element properties to config c = {} for p in self.ElementType.variables: v = self.ElementType.variables[p] if 'seed' in v and v['seed'] is False: continue # Properties which may not be provided by user minval = v['min'] if 'min' in v else None maxval = v['max'] if 'max' in v else None units = v['units'] if 'units' in v else None c['seed:%s' % p] = { 'type': v['type'] if 'type' in v else 'float', 'min': v['min'] if 'min' in v else None, 'max': v['max'] if 'max' in v else None, 'units': v['units'] if 'units' in v else None, 'default': v['default'] if 'default' in v else None, 'description': v['description'] if 'description' in v else 'Seeding value of %s' % p, 'level': v['level'] if 'level' in v else self.CONFIG_LEVEL_ADVANCED} self._add_config(c) # 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': OpenDriftSimulation.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': OpenDriftSimulation.CONFIG_LEVEL_BASIC, 'description': description_fallback} self._add_config(c) # 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] self.timer_start('total time') self.timer_start('configuration') self.add_metadata('opendrift_version', opendrift.__version__) logger.info('OpenDriftSimulation initialised (version %s)' % opendrift.version.version_or_git()) # Check if dependencies are outdated import importlib if importlib.util.find_spec("cmocean") is None: logger.warning('#'*82) logger.warning('Dependencies are outdated, please update with: conda env update -f environment.yml') logger.warning('#'*82)
[docs] def list_config(self, prefix=''): """List all possible configuration settings with values""" str = '\n=============================================\n' for key in self._config: if key.startswith(prefix): str += '%s [%s]\n' % (key, self.get_config(key)) str += '=============================================\n' logger.info(str)
[docs] def list_configspec(self, prefix=''): """Readable formatting of config specification""" for c, i in self._config.items(): if c.startswith(prefix): val = i['value'] if 'value' in i else None if i['type'] == 'bool': rang='' elif i['type'] in ['float', 'int']: rang = 'min: %s, max: %s [%s]' % (i['min'], i['max'], i['units']) elif i['type'] == 'enum': rang = i['enum'] print('%-35s [%s] %-5s %s %s...' % (c, val, i['type'], rang, i['description'][0:20]))
[docs] def get_configspec(self, prefix='', level=[1,2,3]): if not isinstance(level, list): level = [level] configspec = {k:v for (k,v) in self._config.items() if k.startswith(prefix) and self._config[k]['level'] in level} return configspec
[docs] def _add_config(self, config, overwrite=True): """Add configuration settings config is a dictionary where keys are configuration keywords, and values are dictionaries with the following contents: type (string): 'float', 'int', 'bool' or 'enum' min, max (float/int/None): (only when type is 'float' or 'int') The minimum and maximum allowed values for this setting. May also be None if there are no upper/lowe limits. units (string): (only when type is 'float' or 'int') The units of this config setting. enum (list): (only when type is 'enum') A list of possible values for this setting. default (number/bool/string/None): The default value for this setting. value (number/bool/string/None): The actual value for this setting. This is updated with self.set_config(key, value) and retrieved with self.get_config(key) description (string): A description of this config setting, for users/documentation/GUIs. level (int): A parameter to determine the level of exposure in GUIs 1 self.CONFIG_LEVEL_ESSENTIAL: important setting which user has to consider 2 self.CONFIG_LEVEL_BASIC: setting which many users may consider 3 self.CONFIG_LEVEL_ADVANCED: setting relevant only to advanced users """ caller = inspect.stack()[1] caller = os.path.splitext(os.path.basename(caller.filename))[0] logger.debug('Adding %i config items from %s' % (len(config), caller)) if not hasattr(self, '_config'): self._config = {} remove = [] for c, i in config.items(): # Check that provided config is conistent if c in self._config: if overwrite is False: logger.debug(' Config item %s is already specified, not overwriting' % c) remove.append(c) else: logger.debug(' Overwriting config item %s' % c) for p in ['type', 'description', 'level']: if p not in i: raise ValueError('"%s" must be specified for config item %s' % (p, c)) if i['level'] != self.CONFIG_LEVEL_ESSENTIAL and 'default' not in i:#or i['default'] is None: raise ValueError('A default value must be provided for config item %s' % c) if i['type'] == 'enum': if 'enum' not in i or not isinstance(i['enum'], list): raise ValueError('"enum" of type list must be provided for config item %s' % (c)) elif i['type'] in ['float', 'int']: for p in ['min', 'max', 'units']: if p not in i: raise ValueError('"%s" not provided for config item %s' % (p, c)) elif i['type'] == 'bool': pass # no check for bool else: raise ValueError('Config type "%s" (%s) is not defined. Valid options are: ' 'float, int, enum, bool' % (i['type'], c)) if 'default' in i: i['value'] = i['default'] for r in remove: del config[r] self._config.update(config)
[docs] def set_config(self, key, value): if not key in self._config: self.list_config() raise ValueError('No config setting named %s' % key) i = self._config[key] if i['type'] == 'bool': if value not in [True, False]: raise ValueError('Config value %s must be True or False' % key) elif i['type'] in ['float', 'int'] and value is not None: if (i['min'] is not None and value < i['min']) or (i['max'] is not None and value > i['max']): raise ValueError('Config value %s must be between %s and %s' % (key, i['min'], i['max'])) if i['type'] == 'float' and value is not None: value = float(value) elif i['type'] == 'int' and value is not None: value = int(value) elif i['type'] == 'enum': if value not in i['enum']: if len(i['enum']) > 5: import difflib matches = difflib.get_close_matches(value, i['enum'], n=20, cutoff=.3) containing = [e for e in i['enum'] if value in e] matches = list(set(matches) | set(containing)) if len(matches) > 0: matches.sort() suggestion = '\nDid you mean any of these?\n%s' % str(matches) else: suggestion = '' raise ValueError('Wrong configuration, possible values are:\n\t%s\n%s' % (i['enum'], suggestion)) self._config[key]['value'] = value
[docs] def _set_config_default(self, key, value): """Update both default and actual value of a config setting""" self.set_config(key, value) self._config[key]['default'] = self.get_config(key)
[docs] def get_config(self, key): if not key in self._config: raise ValueError('No config setting named %s' % key) return(self._config[key]['value'])
[docs] def add_metadata(self, key, value): """Add item to metadata dictionary, for export as netCDF global attributes""" if not hasattr(self, 'metadata_dict'): from collections import OrderedDict self.metadata_dict = OrderedDict() self.metadata_dict[key] = value
[docs] def prepare_run(self): pass # to be overloaded when needed
[docs] def store_present_positions(self, IDs=None, lons=None, lats=None): """Store present element positions, in case they shall be moved back""" if self.get_config('general:coastline_action') == 'previous' or ('general:seafloor_action' in self._config and self.get_config('general:seafloor_action') == 'previous'): if not hasattr(self, 'previous_lon'): self.previous_lon = np.ma.masked_all(self.num_elements_total()) self.previous_lat = np.ma.masked_all(self.num_elements_total()) if IDs is None: IDs = self.elements.ID lons = self.elements.lon lats = self.elements.lat self.newly_seeded_IDs = None else: # to check if seeded on land if len (IDs) > 0: self.newly_seeded_IDs = np.copy(IDs) else: self.newly_seeded_IDs = None self.previous_lon[IDs-1] = np.copy(lons) self.previous_lat[IDs-1] = np.copy(lats)
[docs] def store_previous_variables(self): """Store some environment variables, for access at next time step""" if not hasattr(self, 'store_previous'): return if not hasattr(self, 'variables_previous'): # Create ndarray to store previous variables dtype = [(var, np.float32) for var in self.store_previous] self.variables_previous = np.array( np.full(self.num_elements_total(), np.nan), dtype=dtype) # Copying variables_previous to environment_previous self.environment_previous = self.variables_previous[self.elements.ID-1] # Use new values for new elements which have no previous value for var in self.store_previous: undefined = np.isnan(self.environment_previous[var]) self.environment_previous[var][undefined] = getattr(self.environment, var)[undefined] self.environment_previous = self.environment_previous.view(np.recarray) for var in self.store_previous: self.variables_previous[var][self.elements.ID-1] = getattr(self.environment, var)
[docs] def interact_with_coastline(self, final=False): """Coastline interaction according to configuration setting""" if self.num_elements_active() == 0: return i = self.get_config('general:coastline_action') if not hasattr(self, 'environment') or not hasattr(self.environment, 'land_binary_mask'): return if i == 'none': # Do nothing return if final is True: # Get land_binary_mask for final location en, en_prof, missing = \ self.get_environment(['land_binary_mask'], self.time, self.elements.lon, self.elements.lat, self.elements.z, None) self.environment.land_binary_mask = en.land_binary_mask if i == 'stranding': # Deactivate elements on land self.deactivate_elements( self.environment.land_binary_mask == 1, reason='stranded') elif i == 'previous': # Go back to previous position (in water) if self.newly_seeded_IDs is not None: self.deactivate_elements( (self.environment.land_binary_mask == 1) & (self.elements.age_seconds == self.time_step.total_seconds()), reason='seeded_on_land') on_land = np.where(self.environment.land_binary_mask == 1)[0] if len(on_land) == 0: logger.debug('No elements hit coastline.') else: logger.debug('%s elements hit coastline, ' 'moving back to water' % len(on_land)) on_land_ID = self.elements.ID[on_land] self.elements.lon[on_land] = \ np.copy(self.previous_lon[on_land_ID - 1]) self.elements.lat[on_land] = \ np.copy(self.previous_lat[on_land_ID - 1]) self.environment.land_binary_mask[on_land] = 0
[docs] def interact_with_seafloor(self): """Seafloor interaction according to configuration setting""" if self.num_elements_active() == 0: return if 'sea_floor_depth_below_sea_level' not in self.priority_list: return sea_floor_depth = self.sea_floor_depth() below = np.where(self.elements.z < -sea_floor_depth)[0] if len(below) == 0: logger.debug('No elements hit seafloor.') return i = self.get_config('general:seafloor_action') if i == 'lift_to_seafloor': 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') elif i == 'previous': # Go back to previous position (in water) logger.warning('%s elements hit seafloor, ' 'moving back ' % len(below)) below_ID = self.elements.ID[below] self.elements.lon[below] = \ np.copy(self.previous_lon[below_ID - 1]) self.elements.lat[below] = \ np.copy(self.previous_lat[below_ID - 1])
[docs] @abstractmethod def update(self): """Any trajectory model implementation must define an update method. This method must/can use environment data (self.environment) to update properties (including position) of its particles (self.elements) """
@abstractproperty def ElementType(self): """Any trajectory model implementation must define an ElementType.""" @abstractproperty def required_variables(self): """Any trajectory model implementation must list needed variables."""
[docs] def test_data_folder(self): import opendrift return os.path.abspath( os.path.join(os.path.dirname(opendrift.__file__), '..', 'tests', 'test_data')) + os.path.sep
[docs] def performance(self): '''Report the time spent on various tasks''' outStr = '--------------------\n' outStr += 'Reader performance:\n' for r in self.readers: reader = self.readers[r] if reader.is_lazy: continue outStr += '--------------------\n' outStr += r + '\n' outStr += reader.performance() outStr += '--------------------\n' outStr += 'Performance:\n' for category, time in self.timing.items(): timestr = str(time)[0:str(time).find('.') + 2] for i, c in enumerate(timestr): if c in '123456789.': timestr = timestr[i:] # Strip leading 0 and : if c == '.': timestr = '0' + timestr break parts = category.split(':') indent = ' '*(len(parts) - 1) category = parts[-1] category = category.replace('<colon>', ':') outStr += '%s%7s %s\n' % (indent, timestr, category) outStr += '--------------------\n' return outStr
[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] 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 # Horizontal buffer of reader must be large enough to cover # the distance possibly covered by elements within a time step if not reader.is_lazy: reader.set_buffer_size(max_speed=self.max_speed) 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 add_readers_from_list(self, urls, timeout=10, lazy=True): '''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) return readers = [reader_from_url(u, timeout) for u in urls] self.add_reader([r for r in readers if r is not None])
[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 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 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_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) return self._initialise_next_lazy_reader() # Call self reader.set_buffer_size(max_speed=self.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 earliest_time(self): return min(self.start_time, self.expected_end_time)
[docs] def latest_time(self): return min(self.start_time, self.expected_end_time)
[docs] def discard_reader_if_not_relevant(self, reader): if hasattr(self, 'expected_endtime') and reader.start_time is not None and ( (reader.start_time > self.latest_time() or reader.end_time < self.earliest_time())): logger.debug('Reader does not cover simulation period') self.discard_reader(reader) return True if len(set(self.required_variables) & set(reader.variables)) == 0: logger.debug('Reader does not contain any relevant variables') self.discard_reader(reader) return True return False
[docs] def discard_reader(self, reader): readername = reader.name logger.debug('Discarding reader: ' + readername) del self.readers[readername] if not hasattr(self, 'discarded_readers'): self.discarded_readers = [readername] else: self.discarded_readers.append(readername) # Remove from priority list for var in self.priority_list: 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 discard_irrelevant_readers(self): for readername in self.readers: reader = self.readers[readername] if reader.is_lazy: continue if self.discard_reader_if_not_relevant(reader): logger.debug('DISCARDED: ' + readername)
[docs] def get_environment(self, variables, time, lon, lat, z, profiles): '''Retrieve environmental variables at requested positions. 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. ''' 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) if not hasattr(self, 'fallback_values'): self.set_fallback_values(refresh=False) # Discard any existing readers which are not relevant self.discard_irrelevant_readers() 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 if self.required_profiles_z_range is not None: self.required_profiles_z_range = np.array( self.required_profiles_z_range) self.required_profiles_z_range[self.required_profiles_z_range<-truncate_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): reader = None if reader is not None: if (reader.covers_time(self.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 i, variable_group in enumerate(variable_groups): logger.debug('----------------------------------------') logger.debug('Variable group %s' % (str(variable_group))) logger.debug('----------------------------------------') reader_group = reader_groups[i] 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 reader.is_lazy: logger.warning('Reader is lazy, should not happen') import sys; sys.exit('Should not happen') 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) 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_from_reader, self.required_profiles_z_range, time, lon[missing_indices], lat[missing_indices], z[missing_indices], self.proj_latlon) except Exception as e: logger.info('========================') logger.info('Exception:') logger.info(e) logger.debug(traceback.format_exc()) logger.info('========================') 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) 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 hasattr(env_tmp[variable_group[0]], 'mask'): 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) else: missing_indices = [] # temporary workaround 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 mask = env[var].mask if any(mask==True): logger.debug(' Using fallback value %s for %s for %s elements' % (self.fallback_values[var], var, np.sum(mask==True))) env[var][mask] = self.fallback_values[var] # 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 = {} env_profiles['z'] = \ np.array(self.required_profiles_z_range)[::-1] if var not in env_profiles: logger.debug(' Using fallback value %s for %s for all profiles' % (self.fallback_values[var], var)) env_profiles[var] = self.fallback_values[var]*\ np.ma.ones((len(env_profiles['z']), self.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] = self.fallback_values[var] logger.debug(' Using fallback value %s for %s for %s profiles' % (self.fallback_values[var], 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 ####################################################### if 'drift:use_tabularised_stokes_drift' in self._config and self.get_config('drift:use_tabularised_stokes_drift') is True: if 'x_wind' not in variables: logger.debug('No wind available to calculate Stokes drift') else: 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'] = \ self.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'] = \ self.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, self.num_elements_active()) env['y_sea_water_velocity'] += np.random.normal( 0, std, self.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, self.num_elements_active()) env['y_sea_water_velocity'] += np.random.uniform( -std, std, self.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, self.num_elements_active()) env['y_wind'] += np.random.normal( 0, std, self.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('---------------------------------') logger.debug('\t\t%s active elements' % self.num_elements_active()) if self.num_elements_active() > 0: lonmin = self.elements.lon.min() lonmax = self.elements.lon.max() latmin = self.elements.lat.min() latmax = self.elements.lat.max() zmin = self.elements.z.min() zmax = self.elements.z.max() if latmin == latmax: logger.debug('\t\tlatitude = %s' % (latmin)) else: logger.debug('\t\t%s <- latitude -> %s' % (latmin, latmax)) if lonmin == lonmax: logger.debug('\t\tlongitude = %s' % (lonmin)) else: logger.debug('\t\t%s <- longitude -> %s' % (lonmin, lonmax)) if zmin == zmax: logger.debug('\t\tz = %s' % (zmin)) else: logger.debug('\t\t%s <- z -> %s' % (zmin, zmax)) 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 num_elements_active(self): """The number of active elements.""" if hasattr(self, 'elements'): return len(self.elements) else: return 0
[docs] def num_elements_deactivated(self): """The number of deactivated elements.""" if hasattr(self, 'elements_deactivated'): return len(self.elements_deactivated) else: return 0
[docs] def num_elements_scheduled(self): if hasattr(self, 'elements_scheduled'): return len(self.elements_scheduled) else: return 0
[docs] def num_elements_total(self): """The total number of scheduled, active and deactivated elements.""" return self.num_elements_activated() + self.num_elements_scheduled()
[docs] def num_elements_activated(self): """The total number of active and deactivated elements.""" return self.num_elements_active() + self.num_elements_deactivated()
[docs] def schedule_elements(self, elements, time): """Schedule elements to be seeded during runtime. Also assigns a unique ID to each particle, monotonically increasing.""" # prepare time if isinstance(time, np.ndarray): time = list(time) if not isinstance(time, list): time = [time] if len(time) == 1 and len(elements) > 1: time = time*len(elements) if not hasattr(self, 'elements_scheduled'): self.elements_scheduled = elements self.elements_scheduled_time = np.array(time) # We start simulation at time of release of first element: self.start_time = time[0] self.elements_scheduled.ID = np.arange(1, len(elements) + 1) else: elements.ID = np.arange(self.num_elements_scheduled() + 1, self.num_elements_scheduled() + 1 + len(elements)) # Increase ID successively self.elements_scheduled.extend(elements) self.elements_scheduled_time = np.append( self.elements_scheduled_time, np.array(time)) min_time = np.min(time) if hasattr(self, 'start_time'): if min_time < self.start_time: self.start_time = min_time logger.debug('Setting simulation start time to %s' % str(min_time)) else: self.start_time = min_time logger.debug('Setting simulation start time to %s' % str(min_time))
[docs] def release_elements(self): """Activate elements which are scheduled within following timestep.""" logger.debug('to be seeded: %s, already seeded %s' % ( len(self.elements_scheduled), self.num_elements_activated())) if len(self.elements_scheduled) == 0: return if self.time_step.days >= 0: indices = (self.elements_scheduled_time >= self.time) & \ (self.elements_scheduled_time < self.time + self.time_step) else: indices = (self.elements_scheduled_time <= self.time) & \ (self.elements_scheduled_time > self.time + self.time_step) self.store_present_positions( self.elements_scheduled.ID[indices], self.elements_scheduled.lon[indices], self.elements_scheduled.lat[indices]) self.elements_scheduled.move_elements(self.elements, indices) self.elements_scheduled_time = self.elements_scheduled_time[~indices] logger.debug('Released %i new elements.' % np.sum(indices))
[docs] def closest_ocean_points(self, lon, lat): """Return the closest ocean points for given lon, lat""" deltalon = 0.01 # grid deltalat = 0.01 numbuffer = 10 lonmin = lon.min() - deltalon*numbuffer lonmax = lon.max() + deltalon*numbuffer latmin = lat.min() - deltalat*numbuffer latmax = lat.max() + deltalat*numbuffer if not 'land_binary_mask' in self.priority_list: logger.info('No land reader added, ' 'making a temporary landmask reader') from opendrift.models.oceandrift import OceanDrift reader_landmask = reader_global_landmask.Reader( extent = [ np.maximum(-360, self.elements_scheduled.lon.min() - deltalon), np.maximum(-89, self.elements_scheduled.lat.min() - deltalat), np.minimum(720, self.elements_scheduled.lon.max() + deltalon), np.minimum(89, self.elements_scheduled.lat.max() + deltalat) ]) o = OceanDrift(loglevel='custom') o.add_reader(reader_landmask) land_reader = reader_landmask else: logger.info('Using existing reader for land_binary_mask') land_reader_name = self.priority_list['land_binary_mask'][0] land_reader = self.readers[land_reader_name] o = self land = o.get_environment(['land_binary_mask'], lon=lon, lat=lat, z=0*lon, time=land_reader.start_time, profiles=None)[0]['land_binary_mask'] if land.max() == 0: logger.info('All points are in ocean') return lon, lat logger.info('Moving %i out of %i points from land to water' % (np.sum(land!=0), len(lon))) landlons = lon[land!=0] landlats = lat[land!=0] longrid = np.arange(lonmin, lonmax, deltalon) latgrid = np.arange(latmin, latmax, deltalat) longrid, latgrid = np.meshgrid(longrid, latgrid) longrid = longrid.ravel() latgrid = latgrid.ravel() # Remove grid-points not covered by this reader latgrid_covered = land_reader.covers_positions(longrid, latgrid)[0] longrid = longrid[latgrid_covered] latgrid = latgrid[latgrid_covered] landgrid = o.get_environment( ['land_binary_mask'], lon=longrid, lat=latgrid, z=0*longrid, time=land_reader.start_time, profiles=None)[0]['land_binary_mask'] if landgrid.min() == 1 or np.isnan(landgrid.min()): logger.warning('No ocean pixels nearby, cannot move elements.') return lon, lat oceangridlons = longrid[landgrid==0] oceangridlats = latgrid[landgrid==0] from scipy import spatial tree = scipy.spatial.cKDTree( np.dstack([oceangridlons, oceangridlats])[0]) landpoints = np.dstack([landlons, landlats]) dist, indices = tree.query(landpoints) indices = indices.ravel() lon[land!=0] = oceangridlons[indices] lat[land!=0] = oceangridlats[indices] return lon, lat
[docs] def seed_elements(self, lon, lat, time, radius=0, number=None, radius_type='gaussian', **kwargs): """Seed elements with given position(s), time and properties. Arguments: lon: scalar or array central longitude(s). lat: scalar or array central latitude(s). radius: scalar or array radius in meters around each lon-lat pair, within which particles will be randomly seeded. number: integer, total number of particles to be seeded If number is None, the number of elements is the length of lon/lat or time if these are arrays. Otherwise the number of elements are obtained from the config-default. time: datenum or list The time at which particles are seeded/released. If time is a list with two elements, elements are seeded continously from start/first to end/last time. If time is a list with more than two elements, the number of elements is equal to len(time) and are seeded as a time series. radius_type: string If 'gaussian' (default), the radius is the standard deviation in x-y-directions. If 'uniform', elements are spread evenly and always inside a circle with the given radius. kwargs: keyword arguments containing properties/attributes and values corresponding to the actual particle type (ElementType). These are forwarded to the ElementType class. All properties for which there are no default value must be specified. """ if 'cone' in kwargs: raise ValueError('Keyword *cone* for seed_elements is deprecated, use seed_cone() instead.') if self.origin_marker is None: self.origin_marker = {} if 'origin_marker' in kwargs: origin_marker = kwargs['origin_marker'] else: origin_marker = len(self.origin_marker) if 'origin_marker_name' in kwargs: origin_marker_name = kwargs['origin_marker_name'] del kwargs['origin_marker_name'] else: origin_marker_name = 'Seed %d' % len(self.origin_marker) if not 'origin_marker' in kwargs: kwargs['origin_marker'] = origin_marker if '_' in origin_marker_name: raise ValueError('Underscore (_) not allowed in origin_marker_name') self.origin_marker[str(origin_marker)] = origin_marker_name.replace(' ', '_') lon = np.atleast_1d(lon).ravel() lat = np.atleast_1d(lat).ravel() radius = np.atleast_1d(radius).ravel() time = np.atleast_1d(time) if lat.max() > 90 or lat.min() < -90: raise ValueError('Latitude must be between -90 and 90 degrees') if len(lon) != len(lat): raise ValueError('Lon and lat must have same lengths') if len(lon) > 1: if number is not None and number != len(lon): raise ValueError('Lon and lat have length %s, but number is %s' % (len(lon), number)) number = len(lon) else: if number is None: if len(time) > 2: number = len(time) # Interpreting as time series else: number = self.get_config('seed:number') lon = lon*np.ones(number) lat = lat*np.ones(number) if len(time) != number and len(time) > 1: if len(time) == 2: # start -> end td = (time[1]-time[0])/(number-1) # timestep between points time = [time[0] + i*td for i in range(number)] else: raise ValueError('Time array has length %s, must be 1, 2 or %s' % (len(time), number)) # Add radius / perturbation if radius.max() > 0: geod = pyproj.Geod(ellps='WGS84') ones = np.ones(np.sum(number)) if radius_type == 'gaussian': x = np.random.randn(np.sum(number))*radius y = np.random.randn(np.sum(number))*radius az = np.degrees(np.arctan2(x, y)) dist = np.sqrt(x*x+y*y) elif radius_type == 'uniform': az = np.random.randn(np.sum(number))*360 dist = np.sqrt(np.random.uniform(0, 1, np.sum(number)))*radius lon, lat, az = geod.fwd(lon, lat, az, dist, radians=False) # If z is 'seafloor' if not 'z' in kwargs or kwargs['z'] is None: if 'seed:seafloor' in self._config: if self.get_config('seed:seafloor') is True: kwargs['z'] = 'seafloor' logger.debug('Seafloor is selected, neglecting z') if 'z' in kwargs and isinstance(kwargs['z'], str) \ and kwargs['z'][0:8] == 'seafloor': # We need to fetch seafloor depth from reader seafloor_constant = self.get_config('environment:constant:sea_floor_depth_below_sea_level') seafloor_fallback = self.get_config('environment:fallback:sea_floor_depth_below_sea_level') if seafloor_constant is not None: env = {'sea_floor_depth_below_sea_level': np.array(seafloor_constant)} elif ('sea_floor_depth_below_sea_level' in self.priority_list ) or len(self._lazy_readers()): if not hasattr(self, 'time'): self.time = time[0] env, env_profiles, missing = \ self.get_environment(['sea_floor_depth_below_sea_level'], time=time[0], lon=lon, lat=lat, z=0*lon, profiles=None) elif seafloor_fallback is not None: env = {'sea_floor_depth_below_sea_level': np.array(seafloor_fallback)} else: raise ValueError('A reader providing the variable ' 'sea_floor_depth_below_sea_level must be ' 'added before seeding elements at seafloor.') # Add M meters if given as 'seafloor+M' if len(kwargs['z']) > 8 and kwargs['z'][8] == '+': meters_above_seafloor = float(kwargs['z'][9::]) logger.info('Seeding elements %f meters above seafloor' % meters_above_seafloor) else: meters_above_seafloor = 0 kwargs['z'] = \ -env['sea_floor_depth_below_sea_level'].astype('float32') + meters_above_seafloor # Creating and scheduling elements elements = self.ElementType(lon=lon, lat=lat, **kwargs) time_array = np.array(time) self.schedule_elements(elements, time)
[docs] def seed_cone(self, lon, lat, time, radius=0, number=None, **kwargs): """Seed elements along a transect/cone between two points/times Arguments: lon: scalar or list with 2 elements [lon0, lon1] lat: scalar or list with 2 elements [lat0, lat] time: datetime or list with 2 elements [t0, t1] radius: scalar or list with 2 elements [r0, r1] Unit: meters number (int): The number of elements. If this is None, the number of elements is taken from configuration. Elements are seeded along a transect from (lon0, lat0) with uncertainty radius r0 at time t0, towards (lon1, lat1) with uncertainty radius r1 at time t1. If r0 != r1, the unceetainty radius is linearly changed along the transect, thus outlining a "cone". """ if number is None: number = self.get_config('seed:number') if number == 1: raise ValueError('For a cone, the number of elements must be at least 2 or more, given is 1') lon = np.atleast_1d(lon).ravel() lat = np.atleast_1d(lat).ravel() radius = np.atleast_1d(radius).ravel() if len(lon) != len(lat): raise ValueError('Lon and lat must have same length (1 or 2)') elif len(lon) > 2: raise ValueError('Lon and lat must have length 1 or 2, given length is %s' % (len(lon))) elif len(lon) == 1: lon = lon*np.ones(number) lat = lat*np.ones(number) elif len(lon) == 2: # Segment from lon0,lat1 to lon1,lat2 geod = pyproj.Geod(ellps='WGS84') lonin = lon latin = lat # Note that npts places points in-between start and end, and does not include these conelonlats = geod.npts(lon[0], lat[0], lon[1], lat[1], number, radians=False) lon, lat = zip(*conelonlats) if len(radius) > 2: raise ValueError('Seed radius must have length 1 or 2') elif len(radius) == 2: # Linear increase from r0 to r1 radius = np.linspace(radius[0], radius[1], number) if isinstance(time, list) and len(time)==1: time = time[0] if hasattr(time, '__len__'): timespan = [time[0], time[-1]] else: timespan = [time, time] radius = radius.astype(np.float32) lonin = lonin if 'lonin' in locals() else [lon.min(), lon.max()] latin = latin if 'latin' in locals() else [lat.min(), lat.max()] self.seed_cone_arguments = {'lon': lonin, 'lat': latin, 'radius': [float(radius[0]), float(radius[-1])], 'time': timespan, 'number': number} # Make GeoJson seeding dict to be saved in netCDF metadata geo = geojson.LineString([(float(lonin[0]), float(latin[0])), (float(lonin[1]), float(latin[1]))]) seed_defaults = self.get_configspec('seed') default_seed = {k.split(':')[-1]:seed_defaults[k]['value'] for k in seed_defaults} if 'seafloor' in default_seed and default_seed['seafloor'] is True: default_seed['z'] = 'seafloor' default_seed = {**default_seed, **kwargs} # Overwrite with explicitly provided values properties = {**default_seed, 'time': [str(timespan[0]), str(timespan[1])], 'radius': [float(radius[0]), float(radius[-1])], 'number': number} f = geojson.Feature(geometry=geo, properties=properties) self.seed_geojson.append(f) # Forwarding calculated cone points/radii to seed_elements self.seed_elements(lon=lon, lat=lat, time=time, radius=radius, number=number, **kwargs)
[docs] def seed_from_geojson(self, gjson): """Under development""" try: gj = geojson.loads(gjson) except: raise ValueError('Could not load GeoJSON string: %s' % gjson) if not gj.is_valid: raise ValueError('GeoJSON string is not valid: %s' % gj.errors()) # Assuming temporally that g is a Feature, and not a FeatureCollection properties = gj['properties'] if 'time' not in properties: raise ValueError('Property "time" is not available') kwargs = {} for prop in properties: if prop == 'time': t = properties['time'] if isinstance(t, list): time = [datetime.fromisoformat(t[0].replace("Z", "+00:00")), datetime.fromisoformat(t[1].replace("Z", "+00:00"))] else: time = datetime.fromisoformat(t.replace("Z", "+00:00")) else: kwargs[prop] = properties[prop] geometry = gj['geometry'] if geometry['type'] == 'Polygon': coords = list(geojson.utils.coords(gj)) lon, lat = zip(*[(c[0], c[1]) for c in coords]) self.seed_within_polygon(lons=lon, lats=lat, time=time, **kwargs) elif geometry['type'] == 'LineString': coords = list(geojson.utils.coords(gj)) lon, lat = zip(*[(c[0], c[1]) for c in coords]) self.seed_cone(lon=lon, lat=lat, time=time, **kwargs) elif geometry['type'] == 'Point': coords = list(geojson.utils.coords(gj)) lon, lat = zip(*[(c[0], c[1]) for c in coords]) self.seed_elements(lon=lon, lat=lat, time=time, **kwargs) else: raise ValueError('Not yet implemented')
[docs] def seed_repeated_segment(self, lons, lats, start_time, end_time, time_interval=None, number_per_segment=None, total_number=None, **kwargs): """Seed elements repeatedly in time along a segment. The segment goes from lon[0],lat[0] to lon[1],lat[1]. The number of elements should be proved as either: 1) number_per_segment, in which case total number of elements is number_per_segment * len(times), or 2) total_number, in which case the number of elements per segment is: total_number / len(times). Any extra elements are duplicated along at the first segment. """ numtimes = int((end_time-start_time).total_seconds()/ time_interval.total_seconds() + 1) times = [start_time+i*time_interval for i in range(numtimes)] geod = pyproj.Geod(ellps='WGS84') if number_per_segment is None: number_per_segment = int(np.floor(total_number/numtimes)) s_lonlats= geod.npts(lons[0], lats[0], lons[1], lats[1], number_per_segment, radians=False) slon, slat = list(zip(*s_lonlats)) slon = np.atleast_1d(slon) slat = np.atleast_1d(slat) lon, time = np.meshgrid(slon, times) lat, time = np.meshgrid(slat, times) lon = lon.ravel() lat = lat.ravel() time = time.ravel() if total_number is not None: additional_elements = total_number - len(lon.ravel()) print('Repeating the %d last points, to obtain %d elements' % (additional_elements, total_number)) lon = np.concatenate((lon, lon[-additional_elements::])) lat = np.concatenate((lat, lat[-additional_elements::])) time = np.concatenate((time, time[-additional_elements::])) self.seed_elements(lon=lon, lat=lat, time=time, **kwargs)
[docs] def seed_within_polygon(self, lons, lats, number=None, **kwargs): """Seed a number of elements within given polygon. Arguments: lon: array of longitudes lat: array of latitudes number: int, number of elements to be seeded kwargs: keyword arguments containing properties/attributes and values corresponding to the actual particle type (ElementType). These are forwarded to method seed_elements(). All properties for which there are no default value must be specified. """ if number == 0: return if number is None: number = self.get_config('seed:number') lons = np.asarray(lons) lats = np.asarray(lats) if len(lons) < 3: logger.info('At least three points needed to make a polygon') return if len(lons) != len(lats): raise ValueError('lon and lat arrays must have same length.') poly = Polygon(list(zip(lons, lats)), closed=True) # Place N points within the polygons proj = pyproj.Proj('+proj=aea +lat_1=%f +lat_2=%f +lat_0=%f ' '+lon_0=%f +R=6370997.0 +units=m +ellps=WGS84' % (lats.min(), lats.max(), (lats.min()+lats.max())/2, (lons.min()+lons.max())/2)) lonlat = poly.get_xy() lon = lonlat[:, 0] lat = lonlat[:, 1] x, y = proj(lon, lat) area = 0.0 for i in range(-1, len(x)-1): area += x[i] * (y[i+1] - y[i-1]) area = abs(area) / 2 # Make points, evenly distributed deltax = np.sqrt(area/number) lonpoints = np.array([]) latpoints = np.array([]) lonlat = poly.get_xy() lon = lonlat[:, 0] lat = lonlat[:, 1] x, y = proj(lon, lat) xvec = np.linspace(x.min() + deltax/2, x.max() - deltax/2, int((x.max()-x.min())/deltax)) yvec = np.linspace(y.min() + deltax/2, y.max() - deltax/2, int((y.max()-y.min())/deltax)) x, y = np.meshgrid(xvec, yvec) lon, lat = proj(x, y, inverse=True) lon = lon.ravel() lat = lat.ravel() points = np.c_[lon, lat] ind = Path(poly.xy).contains_points(points) if not any(ind): # No elements are inside, we seed on border lonpoints = np.append(lonpoints, lons[0:number]) latpoints = np.append(latpoints, lats[0:number]) else: lonpoints = np.append(lonpoints, lon[ind]) latpoints = np.append(latpoints, lat[ind]) if len(ind) == 0: logger.info('Small or irregular polygon, using center point.') lonpoints = np.atleast_1d(np.mean(lons)) latpoints = np.atleast_1d(np.mean(lats)) # Truncate if too many # NB: should also repeat some points, if too few lonpoints = lonpoints[0:number] latpoints = latpoints[0:number] if len(lonpoints) < number: # If number of positions is smaller than requested, # we duplicate the first ones missing = number - len(lonpoints) lonpoints = np.append(lonpoints, lonpoints[0:missing]) latpoints = np.append(latpoints, latpoints[0:missing]) # Finally seed at calculated positions self.seed_elements(lonpoints, latpoints, number=number, **kwargs)
[docs] def seed_from_wkt(self, wkt, number=None, **kwargs): """Seeds elements within (multi)polygons from WKT""" try: from osgeo import ogr, osr except Exception as e: logger.warning(e) raise ValueError('OGR library is needed to parse WKT') if number is None: number = self.get_config('seed:number') geom = ogr.CreateGeometryFromWkt(wkt) total_area = 0 for i in range(0, geom.GetGeometryCount()): g = geom.GetGeometryRef(i) total_area += g.GetArea() logger.info('Total area of all polygons: %s m2' % total_area) num_seeded = 0 for i in range(0, geom.GetGeometryCount()): g = geom.GetGeometryRef(i) num_elements = int(number*g.GetArea()/total_area) if i == geom.GetGeometryCount()-1: # For the last feature we seed the remaining number, # avoiding difference due to rounding: num_elements = number - num_seeded logger.info('\tSeeding %s elements within polygon number %s' % (num_elements, str(i))) try: g.Transform(coordTrans) except: pass b = g.GetBoundary() if b is not None: points = b.GetPoints() lons = [p[0] for p in points] lats = [p[1] for p in points] else: # Alternative if OGR is not built with GEOS support r = g.GetGeometryRef(0) lons = [r.GetX(j) for j in range(r.GetPointCount())] lats = [r.GetY(j) for j in range(r.GetPointCount())] self.seed_within_polygon(lons=lons, lats=lats, number=num_elements, **kwargs) num_seeded += num_elements
[docs] def seed_from_shapefile(self, shapefile, number, layername=None, featurenum=None, **kwargs): """Seeds elements within contours read from a shapefile""" try: from osgeo import ogr, osr except Exception as e: logger.warning(e) raise ValueError('OGR library is needed to read shapefiles.') if 'timeformat' in kwargs: # Recondstructing time from filename, where 'timeformat' # is forwarded to datetime.strptime() kwargs['time'] = datetime.strptime(os.path.basename(shapefile), kwargs['timeformat']) del kwargs['timeformat'] num_seeded_before = self.num_elements_scheduled() targetSRS = osr.SpatialReference() targetSRS.ImportFromEPSG(4326) try: s = ogr.Open(shapefile) except: s = shapefile for layer in s: if layername is not None and layer.GetName() != layername: logger.info('Skipping layer: ' + layer.GetName()) continue else: logger.info('Seeding for layer: %s (%s features)' % (layer.GetDescription(), layer.GetFeatureCount())) coordTrans = osr.CoordinateTransformation(layer.GetSpatialRef(), targetSRS) if featurenum is None: featurenum = range(1, layer.GetFeatureCount() + 1) else: featurenum = np.atleast_1d(featurenum) if max(featurenum) > layer.GetFeatureCount(): raise ValueError('Only %s features in layer.' % layer.GetFeatureCount()) # Loop first through all features to determine total area layer.ResetReading() area_srs = osr.SpatialReference() area_srs.ImportFromEPSG(3857) areaTransform = osr.CoordinateTransformation(layer.GetSpatialRef(), area_srs) areas = np.zeros(len(featurenum)) for i, f in enumerate(featurenum): feature = layer.GetFeature(f - 1) # Note 1-indexing, not 0 if feature is not None: gom = feature.GetGeometryRef().Clone() gom.Transform(areaTransform) areas[i] = gom.GetArea() total_area = np.sum(areas) layer.ResetReading() # Rewind to first layer logger.info('Total area of all polygons: %s m2' % total_area) # Find number of points per polygon numbers = np.round(number*areas/total_area).astype(int) numbers[numbers.argmax()] += int(number-sum(numbers)) for i, f in enumerate(featurenum): feature = layer.GetFeature(f - 1) if feature is None: continue num_elements = numbers[i] geom = feature.GetGeometryRef() logger.info('\tSeeding %s elements within polygon number %s' % (num_elements, featurenum[i])) try: geom.Transform(coordTrans) except Exception as e: logger.warning('Could not transform coordinates:') logger.warning(e) pass #b = geom.GetBoundary() #if b is not None: # points = b.GetPoints() # lons = [p[0] for p in points] # lats = [p[1] for p in points] #else: # Alternative if OGR is not built with GEOS support r = geom.GetGeometryRef(0) lons = [r.GetY(j) for j in range(r.GetPointCount())] lats = [r.GetX(j) for j in range(r.GetPointCount())] self.seed_within_polygon(lons=lons, lats=lats, number=num_elements, **kwargs)
[docs] def seed_letters(self, text, lon, lat, time, number, scale=1.2): """Seed elements within text polygons""" from matplotlib.font_manager import FontProperties fp = FontProperties(family='Bitstream Vera Sans', weight='bold') pol = matplotlib.textpath.TextPath((lon, lat), text, size=1, prop=fp) pol._vertices *= np.array([scale, scale]) patch = matplotlib.patches.PathPatch( pol, facecolor='none', edgecolor='black', transform=ccrs.PlateCarree()) po = patch.get_path().to_polygons() for p in po: self.seed_within_polygon(lons=p[:,0], lats=p[:,1], number=number, time=time)
[docs] def seed_from_ladim(self, ladimfile, roms): """Seed elements from ladim \\*.rls text file: [time, x, y, z, name]""" data = np.loadtxt(ladimfile, dtype={'names': ('time', 'x', 'y', 'z'), 'formats': ('S20', 'f4', 'f4', 'f4')}, usecols=(0,1,2,3)) time = [datetime.strptime(t, "%Y-%m-%dT%H") for t in data['time']] time = np.array(time) lon, lat = roms.xy2lonlat(data['x'], data['y']) z = -data['z'] logger.info('Seeding %i elements from %s:' % (len(lon), ladimfile)) logger.info(' Lons: %f to %f' % (lon.min(), lon.max())) logger.info(' Lats: %f to %f' % (lat.min(), lat.max())) logger.info(' Depths: %f to %f' % (z.min(), z.max())) logger.info(' Time: %s to %s' % (time.min(), time.max())) elements = self.ElementType(lon=lon, lat=lat, z=-z) self.schedule_elements(elements, time)
[docs] def horizontal_diffusion(self): """Move elements with random walk according to given horizontal diffuivity.""" D = self.get_config('drift:horizontal_diffusivity') if D == 0: logger.debug('Horizontal diffusivity is 0, no random walk.') return dt = np.abs(self.time_step.total_seconds()) x_vel = self.elements.moving * np.sqrt(2*D/dt)*np.random.normal(scale=1, size=self.num_elements_active()) y_vel = self.elements.moving * np.sqrt(2*D/dt)*np.random.normal(scale=1, size=self.num_elements_active()) speed = np.sqrt(x_vel*x_vel+y_vel*y_vel) logger.debug('Moving elements according to horizontal diffusivity of %s, with speeds between %s and %s m/s' % (D, speed.min(), speed.max())) self.update_positions(x_vel, y_vel)
[docs] def deactivate_elements(self, indices, reason='deactivated'): """Schedule deactivated particles for deletion (at end of step)""" if any(indices) is False: return if reason not in self.status_categories: self.status_categories.append(reason) logger.debug('Added status %s' % (reason)) reason_number = self.status_categories.index(reason) #if not hasattr(self.elements.status, "__len__"): if len(np.atleast_1d(self.elements.status)) == 1: status = self.elements.status.item() self.elements.status = np.zeros(self.num_elements_active()) self.elements.status.fill(status) # Deactivate elements, if they have not already been deactivated self.elements.status[indices & (self.elements.status ==0)] = \ reason_number logger.debug('%s elements scheduled for deactivation (%s)' % (np.sum(indices), reason)) logger.debug('\t(z: %f to %f)' % (self.elements.z[indices].min(), self.elements.z[indices].max()))
[docs] def remove_deactivated_elements(self): """Moving deactivated elements from self.elements to self.elements_deactivated.""" # All particles scheduled for deletion indices = (self.elements.status != 0) #try: # len(indices) #except: if len(indices) == 0 or np.sum(indices) == 0: logger.debug('No elements to deactivate') return # No elements scheduled for deactivation # Basic, but some more housekeeping will be required later self.elements.move_elements(self.elements_deactivated, indices) logger.debug('Removed %i elements.' % (np.sum(indices))) if hasattr(self, 'environment'): self.environment = self.environment[~indices] logger.debug('Removed %i values from environment.' % (np.sum(indices))) if hasattr(self, 'environment_profiles') and \ self.environment_profiles is not None: for varname, profiles in self.environment_profiles.items(): logger.debug('remove items from profile for '+varname) if varname != 'z': self.environment_profiles[varname] = \ profiles[:, ~indices] logger.debug('Removed %i values from environment_profiles.' % (np.sum(indices)))
#if self.num_elements_active() == 0: # raise ValueError('No more active elements.') # End simulation
[docs] def set_fallback_values(self, refresh=False): if hasattr(self, 'fallback_values') and refresh is False: raise ValueError('Manually editing fallback_values dict is deprecated, please use set_config()') else: c = self.get_configspec('environment:fallback:') self.fallback_values = {} for var in list(c): if c[var]['value'] is not None: self.fallback_values[var.split(':')[-1]] = c[var]['value']
[docs] def run(self, time_step=None, steps=None, time_step_output=None, duration=None, end_time=None, outfile=None, export_variables=None, export_buffer_length=100, stop_on_error=False): """Start a trajectory simulation, after initial configuration. Performs the main loop: - Obtain environment data for positions of all particles. - Call method 'update' to update (incl advect) particle properties. until one of the following conditions are met: - Maximum number of steps are reached - A needed variable can not be obtained by any reader (outside spatial/temporal domain) and has no fallback (default) value. - All particles have been deactivated (e.g. by stranding) - Occurance of any error, whose trace will be output to terminal. Before starting a model run, readers must be added for all required variables, unless fallback values have been specified. Some particles/elements must have been scheduled for seeding, and the run will start at the time when the first element has been scheduled.. Arguments: time_step: interval between particles updates, in seconds or as timedelta. Default: 3600 seconds (1 hour) time_step_output: Time step at which element properties are stored and eventually written to file. Timedelta object or seconds. Default: same as time_step, meaning that all steps are stored The length of the simulation is specified by defining one (and only one) of the following parameters: - steps: integer, maximum number of steps. End of simulation will be self.start_time + steps*self.time_step - duration: timedelta defining the length of the simulation - end_time: datetime object defining the end of the simulation export_variables: list of variables and parameter names to be saved to file. Default is None (all variables are saved) """ # Exporting software and hardware specification, for possible debugging logger.debug(opendrift.versions()) self.timer_end('configuration') self.timer_start('preparing main loop') if self.num_elements_scheduled() == 0: raise ValueError('Please seed elements before starting a run.') self.elements = self.ElementType() # Export seed_geojson as FeatureCollection string self.add_metadata('seed_geojson', geojson.FeatureCollection(self.seed_geojson)) # Collect fallback values from config into dict self.set_fallback_values(refresh=True) if outfile is None and export_buffer_length is not None: logger.debug('No output file is specified, ' 'neglecting export_buffer_length') export_buffer_length = None # Make constant readers if config environment:constant:<var> is 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) missing_variables = self.missing_variables() missing_variables = [m for m in missing_variables if m != 'land_binary_mask'] if len(missing_variables) > 0: has_fallback = [var for var in missing_variables if var in self.fallback_values] has_no_fallback = [var for var in missing_variables if var not in self.fallback_values] #if has_fallback == missing_variables: 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.fallback_values[var])) #else: 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)) # Some cleanup needed if starting from imported state if self.steps_calculation >= 1: self.steps_calculation = 0 if hasattr(self, 'history'): # Delete history matrix before new run delattr(self, 'history') # Renumbering elements from 0 to num_elements, necessary fix when # importing from file, where elements may have been deactivated # TODO: should start from 1? self.elements.ID = np.arange(0, self.num_elements_active()) ######################## # Simulation time step ######################## if time_step is None: time_step = timedelta(minutes=self.get_config('general:time_step_minutes')) if type(time_step) is not timedelta: # Time step may be given in seconds, as alternative to timedelta time_step = timedelta(seconds=time_step) self.time_step = time_step if time_step_output is None: time_step_output = self.get_config('general:time_step_output_minutes') if time_step_output is None: self.time_step_output = self.time_step else: self.time_step_output = timedelta(minutes=time_step_output) else: if type(time_step_output) is timedelta: self.time_step_output = time_step_output else: self.time_step_output = timedelta(seconds=time_step_output) if self.time_step_output.days >= 0 and self.time_step.days < 0: self.time_step_output = -self.time_step_output time_step_ratio = self.time_step_output.total_seconds() / \ self.time_step.total_seconds() if time_step_ratio < 1: raise ValueError('Output time step must be equal or larger ' 'than calculation time step.') if not time_step_ratio.is_integer(): raise ValueError('Ratio of calculation and output time steps ' 'must be an integer - given ratio is %s' % time_step_ratio) ######################## # Simulation duration ######################## if time_step.days < 0: logger.info('Backwards simulation, starting from last seeded element') self.start_time = self.elements_scheduled_time.max() if (duration is not None and end_time is not None) or \ (duration is not None and steps is not None) or \ (steps is not None and end_time is not None): raise ValueError('Only one of "steps", "duration" and "end_time" ' 'may be provided simultaneously') if duration is None and end_time is None: if steps is not None: duration = steps*self.time_step else: for reader in self.readers.values(): if reader.end_time is not None: if end_time is None: end_time = reader.end_time else: end_time = min(end_time, reader.end_time) logger.info('Duration, steps or end time not specified, ' 'running until end of first reader: %s' % (end_time)) if duration is None: duration = end_time - self.start_time if time_step.days < 0 and duration.days >= 0: # Duration shall also be negative for backwards run duration = -duration if np.sign(duration.total_seconds()) * np.sign(time_step.total_seconds()) < 0: raise ValueError("Time step must be negative if duration is negative.") self.expected_steps_output = duration.total_seconds() / \ self.time_step_output.total_seconds() + 1 # Includes start and end self.expected_steps_calculation = duration.total_seconds() / \ self.time_step.total_seconds() self.expected_steps_output = int(self.expected_steps_output) self.expected_steps_calculation = int(self.expected_steps_calculation) self.expected_end_time = self.start_time + self.expected_steps_calculation*self.time_step ############################################################## # Prepare readers for the requested simulation domain/time ############################################################## max_distance = \ self.max_speed*self.expected_steps_calculation * \ np.abs(self.time_step.total_seconds()) deltalat = max_distance/111000. deltalon = deltalat/np.cos( np.radians(np.mean(self.elements_scheduled.lat))) # TODO: extent should ideally be a general polygon, not only lon/lat-min/max # TODO: Should also take into account eventual lifetime of elements simulation_extent = [ np.maximum(-360, self.elements_scheduled.lon.min() - deltalon), np.maximum(-89, self.elements_scheduled.lat.min() - deltalat), np.minimum(720, self.elements_scheduled.lon.max() + deltalon), np.minimum(89, self.elements_scheduled.lat.max() + deltalat)] logger.debug('Preparing readers for simulation coverage (%s) and time (%s to %s)' % (simulation_extent, self.start_time, self.expected_end_time)) for reader in self.readers.values(): logger.debug('\tPreparing %s' % reader.name) reader.prepare( extent=simulation_extent, start_time=self.start_time, end_time = self.expected_end_time) ############################################################## # 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') is True: if 'land_binary_mask' in self.priority_list: if 'basemap_landmask' in self.priority_list['land_binary_mask']: self.priority_list['land_binary_mask'] = ['basemap_landmask'] elif 'global_landmask' in self.priority_list['land_binary_mask']: self.priority_list['land_binary_mask'] = ['global_landmask'] else: del self.priority_list['land_binary_mask'] if self.get_config('general:use_auto_landmask') is True and \ ('land_binary_mask' in self.required_variables and \ 'land_binary_mask' not in self.priority_list \ and 'land_binary_mask' not in self.fallback_values): 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.max_speed) self.timer_start('preparing main loop:making dynamical landmask') reader_landmask = reader_global_landmask.Reader(extent = simulation_extent) self.add_reader(reader_landmask) self.timer_end('preparing main loop:making dynamical landmask') # Move point seeded on land to ocean if self.get_config('seed:ocean_only') is True and \ ('land_binary_mask' in self.required_variables): #('land_binary_mask' not in self.fallback_values) and \ self.timer_start('preparing main loop:moving elements to ocean') self.elements_scheduled.lon, self.elements_scheduled.lat = \ self.closest_ocean_points(self.elements_scheduled.lon, self.elements_scheduled.lat) self.timer_end('preparing main loop:moving elements to ocean') #################################################################### # Preparing history array for storage in memory and eventually file #################################################################### if export_buffer_length is None: self.export_buffer_length = self.expected_steps_output else: self.export_buffer_length = export_buffer_length if self.time_step.days < 0: # For backwards simulation, we start at last seeded element logger.info('Backwards simulation, starting at ' 'time of last seeded element') self.time = self.elements_scheduled_time.max() # Flipping ID array, so that lowest IDs are released first self.elements_scheduled.ID = \ np.flipud(self.elements_scheduled.ID) else: # Forward simulation, start time has been set when seeding self.time = self.start_time # Add the output variables which are always required if export_variables is not None: export_variables = list(set(export_variables + ['lon', 'lat', 'ID', 'status'])) self.export_variables = export_variables # Initialise array to hold history (element properties and environment) # for export to file. history_dtype_fields = [(name, self.ElementType.variables[name]['dtype']) for name in self.ElementType.variables] # Add environment variables self.history_metadata = self.ElementType.variables.copy() for env_var in self.required_variables: history_dtype_fields.append((env_var, np.dtype('float32'))) self.history_metadata[env_var] = {} # Remove variables from output array, if only subset is requested if self.export_variables is not None: history_dtype_fields = [f for f in history_dtype_fields if f[0] in self.export_variables] for m in list(self.history_metadata): if m not in self.export_variables: del self.history_metadata[m] history_dtype = np.dtype(history_dtype_fields) self.history = np.ma.array(np.zeros((len(self.elements_scheduled), self.export_buffer_length)), dtype=history_dtype) self.history.mask = True self.steps_exported = 0 if outfile is not None: self.io_init(outfile) else: self.outfile = None ############################# # Check validity domain ############################# validity_domain = [ self.get_config('drift:deactivate_west_of'), self.get_config('drift:deactivate_east_of'), self.get_config('drift:deactivate_south_of'), self.get_config('drift:deactivate_north_of')] if validity_domain == [None, None, None, None]: self.validity_domain = None else: self.validity_domain = validity_domain ############################# # Model specific preparation ############################# self.prepare_run() ########################## # Main loop ########################## self.add_metadata('simulation_time', datetime.now()) self.timer_end('preparing main loop') self.timer_start('main loop') for i in range(self.expected_steps_calculation): try: # Release elements self.release_elements() if self.num_elements_active() == 0 and self.num_elements_scheduled() > 0: self.steps_calculation += 1 logger.info('No active but %s scheduled elements, skipping timestep %s (%s)' % (self.num_elements_scheduled(), self.steps_calculation, self.time)) self.state_to_buffer() # Append status to history array if self.time is not None: self.time = self.time + self.time_step continue self.increase_age_and_retire() self.interact_with_seafloor() if self.show_continuous_performance is True: logger.info(self.performance()) # Display time to terminal logger.debug('==================================='*2) logger.info('%s - step %i of %i - %i active elements ' '(%i deactivated)' % (self.time, self.steps_calculation + 1, self.expected_steps_calculation, self.num_elements_active(), self.num_elements_deactivated())) logger.debug('%s elements scheduled.' % self.num_elements_scheduled()) logger.debug('==================================='*2) self.environment, self.environment_profiles, missing = \ self.get_environment(list(self.required_variables), self.time, self.elements.lon, self.elements.lat, self.elements.z, self.required_profiles) self.store_previous_variables() self.calculate_missing_environment_variables() if any(missing): self.report_missing_variables() self.interact_with_coastline() self.interact_with_seafloor() self.deactivate_elements(missing, reason='missing_data') self.state_to_buffer() # Append status to history array self.remove_deactivated_elements() # Propagate one timestep forwards self.steps_calculation += 1 if self.num_elements_active() == 0 and self.num_elements_scheduled() == 0: raise ValueError('No more active or scheduled elements, quitting.') # Store location, in case elements shall be moved back self.store_present_positions() ##################################################### if self.num_elements_active() > 0: logger.debug('Calling %s.update()' % type(self).__name__) self.timer_start('main loop:updating elements') self.update() self.timer_end('main loop:updating elements') else: logger.info('No active elements, skipping update() method') ##################################################### self.horizontal_diffusion() if self.num_elements_active() == 0 and self.num_elements_scheduled() == 0: raise ValueError('No active or scheduled elements, quitting simulation') logger.debug('%s active elements (%s deactivated)' % (self.num_elements_active(), self.num_elements_deactivated())) # Updating time if self.time is not None: self.time = self.time + self.time_step except Exception as e: message = ('The simulation stopped before requested ' 'end time was reached.') logger.warning(message) self.store_message(message) logger.info('========================') logger.info('End of simulation:') logger.info(e) logger.info(traceback.format_exc()) logger.info(self.get_messages()) if not hasattr(self, 'environment'): sys.exit('Simulation aborted. ' + self.get_messages()) logger.info('========================') if stop_on_error is True: sys.exit('Stopping on error. ' + self.get_messages()) if self.steps_calculation <= 1: raise ValueError('Simulation stopped within ' 'first timestep. ' + self.get_messages()) break self.timer_end('main loop') self.timer_start('cleaning up') logger.debug('Cleaning up') self.interact_with_coastline(final=True) self.state_to_buffer() # Append final status to buffer ############################# # Add some metadata ############################# for var in self.required_variables: keyword = 'reader_' + var if var not in self.priority_list: self.add_metadata(keyword, self.fallback_values[var]) else: readers = self.priority_list[var] if readers[0].startswith('constant_reader') and var in self.readers[readers[0]]._parameter_value_map: self.add_metadata(keyword, self.readers[readers[ 0]]._parameter_value_map[var][0]) else: self.add_metadata(keyword, self.priority_list[var]) if outfile is not None: logger.debug('Writing and closing output file: %s' % outfile) # Write buffer to outfile, and close if self.steps_output >= self.steps_exported: # Write last lines, if needed self.io_write_buffer() self.io_close() # Remove any elements scheduled for deactivation during last step #self.remove_deactivated_elements() if export_buffer_length is None: # Remove columns for unseeded elements in history array if self.num_elements_scheduled() > 0: logger.info('Removing %i unseeded elements from history array' % self.num_elements_scheduled()) mask = np.ones(self.history.shape[0], dtype=bool) mask[self.elements_scheduled.ID-1] = False self.history = self.history[mask, :] # Remove rows for unreached timsteps in history array self.history = self.history[:, range(self.steps_output)] else: # If output has been flushed to file during run, we # need to reimport from file to get all data in memory del self.environment if hasattr(self, 'environment_profiles'): del self.environment_profiles self.io_import_file(outfile) self.timer_end('cleaning up') self.timer_end('total time')
[docs] def increase_age_and_retire(self): """Increase age of elements, and retire if older than config setting.""" # Increase age of elements self.elements.age_seconds += self.time_step.total_seconds() # Deactivate elements that exceed a certain age if self.get_config('drift:max_age_seconds') is not None: self.deactivate_elements(self.elements.age_seconds >= self.get_config('drift:max_age_seconds'), reason='retired') # Deacticate any elements outside validity domain set by user if self.validity_domain is not None: W, E, S, N = self.validity_domain if W is not None: self.deactivate_elements(self.elements.lon < W, reason='outside') if E is not None: self.deactivate_elements(self.elements.lon > E, reason='outside') if S is not None: self.deactivate_elements(self.elements.lat < S, reason='outside') if N is not None: self.deactivate_elements(self.elements.lat > N, reason='outside')
[docs] def state_to_buffer(self): """Append present state (elements and environment) to recarray.""" steps_calculation_float = \ (self.steps_calculation * self.time_step.total_seconds() / self.time_step_output.total_seconds()) + 1 if self.time_step <= timedelta(seconds=1): self.steps_output = int(np.round(steps_calculation_float)) else: self.steps_output = int(np.floor(steps_calculation_float)) ID_ind = self.elements.ID - 1 time_ind = self.steps_output - 1 - self.steps_exported if steps_calculation_float.is_integer() or self.time_step < timedelta(seconds=1): element_ind = range(len(ID_ind)) # We write all elements else: deactivated = np.where(self.elements.status != 0)[0] if len(deactivated) == 0: return # No deactivated elements this sub-timestep # We write history for deactivated elements only: logger.debug('Writing history for %s deactivated elements' % len(deactivated)) ID_ind = ID_ind[deactivated] element_ind = deactivated time_ind = np.minimum(time_ind + 1, self.history.shape[1] - 1) # TODO: storing of variables and environment below should be collected in a single loop # Store present state in history recarray for i, var in enumerate(self.elements.variables): if self.export_variables is not None and \ var not in self.export_variables: continue # Temporarily assuming elements numbered # from 0 to num_elements_active() # Does not hold when importing ID from a saved file, where # some elements have been deactivated self.history[var][ID_ind, time_ind] = \ getattr(self.elements, var)[element_ind] if len(ID_ind) > 0: newmin = np.min(self.history[var][ID_ind, time_ind]) newmax = np.max(self.history[var][ID_ind, time_ind]) if var not in self.minvals: self.minvals[var] = newmin self.maxvals[var] = newmax else: self.minvals[var] = np.minimum(self.minvals[var], newmin) self.maxvals[var] = np.maximum(self.maxvals[var], newmax) # Copy environment data to history array for i, var in enumerate(self.environment.dtype.names): if self.export_variables is not None and \ var not in self.export_variables: continue self.history[var][ID_ind, time_ind] = \ getattr(self.environment, var)[element_ind] if len(ID_ind) > 0: newmin = np.min(self.history[var][ID_ind, time_ind]) newmax = np.max(self.history[var][ID_ind, time_ind]) if var not in self.minvals: self.minvals[var] = newmin self.maxvals[var] = newmax else: self.minvals[var] = np.minimum(self.minvals[var], newmin) self.maxvals[var] = np.maximum(self.maxvals[var], newmax) # Call writer if buffer is full if (self.outfile is not None) and \ ((self.steps_output - self.steps_exported) == self.export_buffer_length): self.io_write_buffer()
[docs] def report_missing_variables(self): """Issue warning if some environment variables missing.""" missing_variables = [] for var in self.required_variables: if np.isnan(getattr(self.environment, var).min()): missing_variables.append(var) if len(missing_variables) > 0: logger.warning('Missing variables: ' + str(missing_variables)) self.store_message('Missing variables: ' + str(missing_variables))
[docs] def index_of_activation_and_deactivation(self): """Return the indices when elements were seeded and deactivated.""" firstlast = np.ma.notmasked_edges(self.history['lon'], axis=1) index_of_activation = firstlast[0][1] index_of_deactivation = firstlast[1][1] if len(index_of_deactivation) < self.history['lon'].shape[0]: missingind = np.setdiff1d( np.arange(0, self.history['lon'].shape[0]), firstlast[0][0]) logger.warning('%s elements were never seeded, removing from history array (this is probably caused by importing an old file)' % len(missingind)) self.history = self.history[firstlast[0][0], :] return index_of_activation, index_of_deactivation
[docs] def set_up_map(self, corners=None, buffer=.1, delta_lat=None, lscale=None, fast=False, hide_landmask=False, **kwargs): """ Generate Figure instance on which trajectories are plotted. :param hide_landmask: do not plot landmask (default False) :type hide_landmask: bool provide corners=[lonmin, lonmax, latmin, latmax] for specific map selection """ # Initialise map if hasattr(self, 'ds'): # If dataset is lazily imported lons = self.ds.lon lats = self.ds.lat logger.debug('Finding min longitude...') self.lonmin = np.nanmin(self.ds.lon) logger.debug('Finding max longitude...') self.lonmax = np.nanmax(self.ds.lon) logger.debug('Finding min latitude...') self.latmin = np.nanmin(self.ds.lat) logger.debug('Finding max latitude...') self.latmax = np.nanmax(self.ds.lat) if not hasattr(self, 'af'): if os.path.exists(self.analysis_file): self.af = Dataset(self.analysis_file, 'a') else: self.af = Dataset(self.analysis_file, 'w') self.af.lonmin = self.lonmin self.af.lonmax = self.lonmax self.af.latmin = self.latmin self.af.latmax = self.latmax self.af.close() else: lons, lats = self.get_lonlats() # TODO: to be removed if corners is not None: # User provided map corners lonmin = corners[0] lonmax = corners[1] latmin = corners[2] latmax = corners[3] elif hasattr(self, 'lonmin'): # if dataset is lazily imported lonmin = self.lonmin lonmax = self.lonmax latmin = self.latmin latmax = self.latmax else: lons, lats = self.get_lonlats() lonmin = np.nanmin(lons) - buffer*2 lonmax = np.nanmax(lons) + buffer*2 latmin = np.nanmin(lats) - buffer latmax = np.nanmax(lats) + buffer if fast is True: logger.warning('Plotting fast. This will make your plots less accurate.') import matplotlib.style as mplstyle mplstyle.use(['fast']) # use a spherical earth axis = 57.29577951308232 # something to do with pi globe = ccrs.Globe(ellipse = None, semimajor_axis = axis, semiminor_axis = axis) crs = ccrs.Mercator(globe = globe) if lscale is None: lscale = 'c' else: crs = ccrs.Mercator() if lscale is None: lscale = 'auto' meanlat = (latmin + latmax)/2 aspect_ratio = float(latmax-latmin) / (float(lonmax-lonmin)) aspect_ratio = aspect_ratio / np.cos(np.radians(meanlat)) if aspect_ratio > 1: fig = plt.figure(figsize=(11./aspect_ratio, 11.)) else: fig = plt.figure(figsize=(11., 11.*aspect_ratio)) ax = fig.add_subplot(111, projection=crs) # need '111' for Python 2 ax.set_extent([lonmin, lonmax, latmin, latmax], crs=ccrs.PlateCarree()) if 'ocean_color' in kwargs: ax.patch.set_facecolor(kwargs['ocean_color']) ocean_color = kwargs['ocean_color'] else: ocean_color = 'white' if 'land_color' in kwargs: land_color = kwargs['land_color'] else: if fast is True: land_color = 'gray' else: land_color = cfeature.COLORS['land'] if 'text' in kwargs: if not isinstance(kwargs['text'], list): text = list(kwargs['text']) else: text = kwargs['text'] for te in text: plt.text(transform=ccrs.Geodetic(), **te) if 'box' in kwargs: if not isinstance(kwargs['box'], list): box = list(kwargs['box']) else: box = kwargs['box'] for bx in box: lonmn = bx['lon'][0] lonmx = bx['lon'][1] latmn = bx['lat'][0] latmx = bx['lat'][1] del bx['lon'] del bx['lat'] if 'text' in bx: plt.text(x=lonmn, y=latmx, s=bx['text'], transform=ccrs.Geodetic()) del bx['text'] patch = matplotlib.patches.Rectangle(xy=[lonmn, latmn], width=lonmx-lonmn, height=latmx-latmn, transform=ccrs.Geodetic(), **bx) ax.add_patch(patch) if not hide_landmask: if 'land_binary_mask' in self.priority_list and self.priority_list['land_binary_mask'][0] == 'shape': logger.debug('Using custom shapes for plotting land..') ax.add_geometries(self.readers['shape'].polys, ccrs.PlateCarree(), facecolor=land_color, edgecolor='black') else: reader_global_landmask.plot_land(ax, lonmin, latmin, lonmax, latmax, fast, ocean_color, land_color, lscale) gl = ax.gridlines(ccrs.PlateCarree(), draw_labels=True) if cartopy.__version__ < '0.18.0': gl.xlabels_top = False # Cartopy < 0.18 else: gl.top_labels = None # Cartopy >= 0.18 fig.canvas.draw() fig.set_tight_layout(True) if not hasattr(self, 'ds'): try: firstlast = np.ma.notmasked_edges(lons, axis=1) index_of_first = firstlast[0][1] index_of_last = firstlast[1][1] except: index_of_last = 0 else: index_of_first = None index_of_last = None try: # Activate figure zooming mng = plt.get_current_fig_manager() mng.toolbar.zoom() except: pass try: # Maximise figure window size mng.resize(*mng.window.maxsize()) except: pass return fig, ax, crs, lons.T, lats.T, index_of_first, index_of_last
[docs] def get_lonlats(self): if hasattr(self, 'history'): lons = self.history['lon'] lats = self.history['lat'] else: if self.steps_output > 0: lons = np.ma.array(np.reshape(self.elements.lon, (1, -1))).T lats = np.ma.array(np.reshape(self.elements.lat, (1, -1))).T else: lons = np.ma.array( np.reshape(self.elements_scheduled.lon, (1, -1))).T lats = np.ma.array( np.reshape(self.elements_scheduled.lat, (1, -1))).T return lons, lats
[docs] def animation(self, buffer=.2, corners=None, filename=None, compare=None, background=None, bgalpha=.5, vmin=None, vmax=None, drifter=None, skip=5, scale=10, color=False, clabel=None, colorbar=True, cmap=None, density=False, show_elements=True, show_trajectories=False, trajectory_alpha=.1, hide_landmask=False, density_pixelsize_m=1000, unitfactor=1, lcs=None, surface_only=False, markersize=20, origin_marker=None, legend=None, legend_loc='best', fps=8, lscale=None, fast=False, **kwargs): """Animate last run.""" if self.num_elements_total() == 0 and not hasattr(self, 'ds'): raise ValueError('Please run simulation before animating') markersizebymass = False if isinstance(markersize, str): if markersize == 'mass': markersizebymass = True markersize = 20 start_time = datetime.now() if cmap is None: cmap = 'jet' if isinstance(cmap, str): cmap = matplotlib.cm.get_cmap(cmap) if color is False and background is None and lcs is None and density is False: colorbar = False markercolor = self.plot_comparison_colors[0] if isinstance(density, str): # Density field is weighted by this variable # TODO: not yet implemented! density_weight = density density = True else: if density is True: density_weight = None elif density is not False: density_weight=density density=True if density is True: # Get density arrays if hasattr(self, 'ds'): # opened with Xarray if origin_marker is None: origin_marker = 0 per_origin_marker = False else: per_origin_marker = True H, H_om, lon_array, lat_array = self.get_density_xarray(pixelsize_m=density_pixelsize_m, weights=density_weight) if per_origin_marker is True: H = H_om[:,:,:,origin_marker] else: if origin_marker is not None: raise ValueError('Separation by origin_marker is only active when imported from file with ' 'open_xarray: https://opendrift.github.io/gallery/example_huge_output.html') H, H_submerged, H_stranded, lon_array, lat_array = \ self.get_density_array(pixelsize_m=density_pixelsize_m, weight=density_weight) H = H + H_submerged + H_stranded # x, y are longitude, latitude -> i.e. in a PlateCarree CRS gcrs = ccrs.PlateCarree() def plot_timestep(i): """Sub function needed for matplotlib animation.""" ax.set_title('%s\n%s UTC' % (self._figure_title(), times[i])) if background is not None: if isinstance(background, xr.DataArray): scalar = background[i,:,:].values elif isinstance(background, str): map_x, map_y, scalar, u_component, v_component, qmap_x, qmap_y = \ self.get_map_background(ax, background, time=times[i]) # https://stackoverflow.com/questions/18797175/animation-with-pcolormesh-routine-in-matplotlib-how-do-i-initialize-the-data bg.set_array(scalar[:-1,:-1].ravel()) if type(background) is list: bg_quiv.set_UVC(u_component[::skip, ::skip], v_component[::skip, ::skip]) if lcs is not None: ax.pcolormesh( lcs['lon'], lcs['lat'], lcs['ALCS'][i,:,:], alpha=bgalpha, vmin=vmin, vmax=vmax, cmap=cmap, transform = gcrs) if density is True: # Update density plot pm.set_array(H[i,:,:].ravel()) # Move points if show_elements is True: points.set_offsets(np.c_[x[i, range(x.shape[1])], y[i, range(x.shape[1])]]) points_deactivated.set_offsets(np.c_[ x_deactive[index_of_last_deactivated < i], y_deactive[index_of_last_deactivated < i]]) if markersizebymass: points.set_sizes(100*(self.history['mass'][:,i] / (self.history['mass'][:,i] + self.history['mass_degraded'][:,i] + self.history['mass_volatilized'][:,i]))) if color is not False: # Update colors points.set_array(colorarray[:, i]) if isinstance(color, str) or hasattr(color, '__len__'): points_deactivated.set_array( colorarray_deactivated[ index_of_last_deactivated < i]) if drifter is not None: drifter_pos.set_offsets(np.c_[drifter['x'][i], drifter['y'][i]]) if show_elements is True: if compare is not None: for cd in compare_list: cd['points_other'].set_offsets(np.c_[ cd['x_other'][range(cd['x_other'].shape[0]), i], cd['y_other'][range(cd['x_other'].shape[0]), i]]) cd['points_other_deactivated'].set_offsets(np.c_[ cd['x_other_deactive'][ cd['index_of_last_deactivated_other'] < i], cd['y_other_deactive'][ cd['index_of_last_deactivated_other'] < i]]) return points, cd['points_other'] else: return points # Find map coordinates and plot points with empty data fig, ax, crs, x, y, index_of_first, index_of_last = \ self.set_up_map(buffer=buffer, corners=corners, lscale=lscale, fast=fast, hide_landmask=hide_landmask, **kwargs) if surface_only is True: z = self.get_property('z')[0] x[z<0] = np.nan y[z<0] = np.nan if show_trajectories is True: ax.plot(x, y, color='gray', alpha=trajectory_alpha, transform = gcrs) if color is not False and show_elements is True: if isinstance(color, str): colorarray = self.get_property(color)[0].T colorarray = colorarray*unitfactor colorarray_deactivated = \ self.get_property(color)[0][ index_of_last[self.elements_deactivated.ID-1], self.elements_deactivated.ID-1].T elif hasattr(color, '__len__'): # E.g. array/list of ensemble numbers colorarray_deactivated = color[self.elements_deactivated.ID-1] colorarray = np.tile(color, (self.steps_output, 1)).T else: colorarray = color if vmin is None: vmin = colorarray.min() vmax = colorarray.max() if background is not None: if isinstance(background, xr.DataArray): map_x = background.coords['lon_bin'] map_y = background.coords['lat_bin'] scalar = background[0,:,:] map_y, map_x = np.meshgrid(map_y, map_x) else: map_x, map_y, scalar, u_component, v_component, qmap_x, qmap_y = \ self.get_map_background(ax, background, time=self.start_time) bg = ax.pcolormesh(map_x, map_y, scalar[:-1,:-1], alpha=bgalpha, antialiased=True, linewidth=0.0, rasterized=True, vmin=vmin, vmax=vmax, cmap=cmap, transform = gcrs) if type(background) is list: bg_quiv = ax.quiver(qmap_x[::skip, ::skip], qmap_y[::skip, ::skip], u_component[::skip, ::skip], v_component[::skip, ::skip], scale=scale, transform = gcrs) if lcs is not None: if vmin is None: vmin = lcs['ALCS'].min() vmax = lcs['ALCS'].max() lcsh = ax.pcolormesh(lcs['lon'], lcs['lat'], lcs['ALCS'][0,:,:], vmin=vmin, vmax=vmax, cmap=cmap, transform = gcrs) times = self.get_time_array()[0] if show_elements is True: index_of_last_deactivated = \ index_of_last[self.elements_deactivated.ID-1] if legend is None: legend = [''] if color is False: c = markercolor else: c = [] if markersizebymass: points = ax.scatter([], [], c=c, zorder=10, edgecolor=[], cmap=cmap, alpha=.4, vmin=vmin, vmax=vmax, label=legend[0], transform = gcrs) else: points = ax.scatter([], [], c=c, zorder=10, edgecolor=[], cmap=cmap, s=markersize, vmin=vmin, vmax=vmax, label=legend[0], transform = gcrs) if (compare is None) and (legend != ['']): markers=[] for legend_index in np.arange(len(legend)): markers.append(matplotlib.lines.Line2D([0], [0], marker='o', color='w', markerfacecolor=cmap(legend_index/(len(legend)-1)), markersize=10, label=legend[legend_index])) ax.legend(markers, legend, loc=legend_loc) # Plot deactivated elements, with transparency if markersizebymass: points_deactivated = ax.scatter([], [], c=c, zorder=9, vmin=vmin, vmax=vmax, s=markersize, cmap=cmap, edgecolor=[], alpha=0, transform = gcrs) else: points_deactivated = ax.scatter([], [], c=c, zorder=9, vmin=vmin, vmax=vmax, s=markersize, cmap=cmap, edgecolor=[], alpha=.3, transform = gcrs) x_deactive, y_deactive = (self.elements_deactivated.lon, self.elements_deactivated.lat) if compare is not None: compare_list = self._get_comparison_xy_for_plots(compare) for cn, cd in enumerate(compare_list): if legend != ['']: legstr = legend[cn+1] else: legstr = None cd['points_other'] = \ ax.scatter([], [], c=self.plot_comparison_colors[cn+1], s=markersize, label=legstr, zorder=10, transform = gcrs) # Plot deactivated elements, with transparency cd['points_other_deactivated'] = \ ax.scatter([], [], alpha=.3, zorder=9, c=self.plot_comparison_colors[cn+1], s=markersize, transform = gcrs) if legend != ['', '']: plt.legend(markerscale=2, loc=legend_loc) if density is True: cmap.set_under('w') H = np.ma.masked_where(H==0, H) lat_array, lon_array = np.meshgrid(lat_array, lon_array) if vmax is None: vmax=H.max() pm = ax.pcolormesh(lon_array, lat_array, H[0,:,:], vmin=0.1, vmax=vmax, cmap=cmap, transform=gcrs) if drifter is not None: # Interpolate drifter time series onto simulation times sts = np.array([t.total_seconds() for t in np.array(times) - times[0]]) dts = np.array([t.total_seconds() for t in np.array(drifter['time']) - times[0]]) drifter['x'] = np.interp(sts, dts, drifter['lon']) drifter['y'] = np.interp(sts, dts, drifter['lat']) drifter['x'][sts<dts[0]] = np.nan drifter['x'][sts>dts[-1]] = np.nan drifter['y'][sts<dts[0]] = np.nan drifter['y'][sts>dts[-1]] = np.nan dlabel = drifter['label'] if 'label' in drifter else 'Drifter' dcolor = drifter['color'] if 'color' in drifter else 'r' dlinewidth = drifter['linewidth'] if 'linewidth' in drifter else 2 dmarkersize = drifter['markersize'] if 'markersize' in drifter else 20 drifter_pos = ax.scatter([], [], c=dcolor, zorder=15, s=dmarkersize, label=dlabel, transform=gcrs) ax.plot(drifter['x'], drifter['y'], color=dcolor, linewidth=dlinewidth, zorder=14, transform=gcrs) plt.legend() fig.canvas.draw() fig.set_tight_layout(True) if colorbar is True: if color is not False: if isinstance(color, str) or clabel is not None: if clabel is None: clabel = color item = points elif density is not False: item = pm if clabel is None: clabel = 'density' elif lcs is not None: item = lcsh if clabel is None: clabel = 'LCS' elif background is not None: item = bg if clabel is None: if isinstance(background, xr.DataArray): clabel = background.name else: clabel = background cb = fig.colorbar(item, orientation='horizontal', pad=.05, aspect=30, shrink=.8, drawedges=False) cb.set_label(clabel) anim = animation.FuncAnimation( plt.gcf(), plot_timestep, blit=False, frames=x.shape[0], interval=50) if filename is not None or 'sphinx_gallery' in sys.modules: self._save_animation(anim, filename, fps) logger.debug('Time to make animation: %s' % (datetime.now()-start_time)) else: try: plt.show() except AttributeError: pass
[docs] def animation_profile(self, filename=None, compare=None, legend=['', ''], markersize=5, fps=20, color=None, cmap=None, vmin=None, vmax=None, legend_loc=None): """Animate vertical profile of the last run.""" def plot_timestep(i): """Sub function needed for matplotlib animation.""" #plt.gcf().gca().set_title(str(i)) ax.set_title('%s UTC' % times[i]) if PlotColors: points.set_offsets(np.array( [x[range(x.shape[0]), i].T,z[range(x.shape[0]), i].T] ).T) points.set_array(colorarray[:, i]) else: points.set_data(x[range(x.shape[0]), i], z[range(x.shape[0]), i]) points_deactivated.set_data( x_deactive[index_of_last_deactivated < i], z_deactive[index_of_last_deactivated < i]) if compare is not None: points_other.set_data(x_other[range(x_other.shape[0]), i], z_other[range(x_other.shape[0]), i]) points_other_deactivated.set_data( x_other_deactive[index_of_last_deactivated_other < i], z_other_deactive[index_of_last_deactivated_other < i]) return points, points_other else: return points PlotColors=(compare is None) and (legend != ['','']) if PlotColors: if cmap is None: cmap = 'jet' if isinstance(cmap, str): cmap = matplotlib.cm.get_cmap(cmap) if color is not False: if isinstance(color, str): colorarray = self.get_property(color)[0].T # Set up plot index_of_first, index_of_last = \ self.index_of_activation_and_deactivation() z = self.get_property('z')[0].T x = self.get_property('lon')[0].T #seafloor_depth = \ # -self.get_property('sea_floor_depth_below_sea_level')[0].T fig = plt.figure(figsize=(10, 6.)) # Suitable aspect ratio ax = fig.gca() plt.xlabel('Longitude [degrees]') plt.ylabel('Depth [m]') times = self.get_time_array()[0] index_of_last_deactivated = \ index_of_last[self.elements_deactivated.ID-1] if PlotColors: points = ax.scatter([], [], c=[], zorder=10, edgecolor=[], cmap=cmap, s=markersize, vmin=vmin, vmax=vmax) markers=[] for legend_index in np.arange(len(legend)): markers.append(matplotlib.lines.Line2D([0], [0], marker='o', linewidth=0, markeredgewidth=0, markerfacecolor=cmap(legend_index/(len(legend)-1)), markersize=10, label=legend[legend_index])) leg=ax.legend(markers, legend, loc=legend_loc) leg.set_zorder(20) else: points = plt.plot([], [], '.k', label=legend[0], markersize=markersize)[0] # Plot deactivated elements, with transparency points_deactivated = plt.plot([], [], '.k', alpha=.3)[0] x_deactive = self.elements_deactivated.lon z_deactive = self.elements_deactivated.z if compare is not None: if type(compare) is str: # Other is given as filename other = self.__class__(loglevel=0) other.io_import_file(compare) else: # Other is given as an OpenDrift object other = compare z_other = other.get_property('z')[0].T x_other = self.get_property('lon')[0].T points_other = plt.plot(x_other[0, 0], z_other[0, 0], '.r', label=legend[1], markersize=markersize)[0] # Plot deactivated elements, with transparency points_other_deactivated = plt.plot([], [], '.r', alpha=.3)[0] x_other_deactive = other.elements_deactivated.lon z_other_deactive = other.elements_deactivated.z firstlast = np.ma.notmasked_edges(x_other, axis=1) index_of_last_other = firstlast[1][1] index_of_last_deactivated_other = \ index_of_last_other[other.elements_deactivated.ID-1] xmax = np.maximum(x.max(), x_other.max()) xmin = np.minimum(x.min(), x_other.min()) zmax = np.maximum(z.max(), z_other.max()) zmin = np.minimum(z.min(), z_other.min()) else: xmin = x.min() xmax = x.max() zmin = z.min() zmax = z.max() # Set figure limits sky = (zmax-zmin)*.1 # Sky height is 10% of water depth plt.xlim([xmin, xmax]) plt.ylim([zmin, sky]) ax.add_patch(plt.Rectangle((xmin, 0), xmax-xmin, sky, color='lightsteelblue')) ax.add_patch(plt.Rectangle((xmin, zmin), xmax-xmin, -zmin, color='cornflowerblue')) if legend != ['', ''] and PlotColors is False: plt.legend(loc=4) anim = animation.FuncAnimation(plt.gcf(), plot_timestep, blit=False, frames=x.shape[1], interval=150) if filename is not None or 'sphinx_gallery' in sys.modules: self._save_animation(anim, filename, fps) else: try: plt.show() except AttributeError: pass
[docs] def _get_comparison_xy_for_plots(self, compare): if not type(compare) is list: compare = [compare] compare_list = [{}]*len(compare) for cn, comp in enumerate(compare): compare_list[cn] = {} cd = compare_list[cn] # pointer to dict with data if type(comp) is str: # Other is given as filename other = self.__class__(loglevel=0) other.io_import_file(comp) else: # Other is given as an OpenDrift object other = comp # Find map coordinates of comparison simulations cd['x_other'], cd['y_other'] = \ (other.history['lon'].copy(), other.history['lat'].copy()) cd['x_other_deactive'], cd['y_other_deactive'] = \ (other.elements_deactivated.lon.copy(), other.elements_deactivated.lat.copy()) cd['firstlast'] = np.ma.notmasked_edges( cd['x_other'], axis=1) cd['index_of_last_other'] = cd['firstlast'][1][1] cd['index_of_last_deactivated_other'] = \ cd['index_of_last_other'][other.elements_deactivated.ID-1] return compare_list
[docs] def plot(self, background=None, buffer=.2, corners=None, linecolor=None, filename=None, show=True, vmin=None, vmax=None, compare=None, cmap='jet', lvmin=None, lvmax=None, skip=2, scale=10, show_scalar=True, contourlines=False, trajectory_dict=None, colorbar=True, linewidth=1, lcs=None, show_particles=True, show_initial=True, density_pixelsize_m=1000, bgalpha=1, clabel=None, surface_color=None, submerged_color=None, markersize=20, title='auto', legend=True, legend_loc='best', lscale=None, fast=False, hide_landmask=False, **kwargs): """Basic built-in plotting function intended for developing/debugging. Plots trajectories of all particles. Positions marked with colored stars: - green: all start positions - red: deactivated particles - blue: particles still active at end of simulation Requires availability of Cartopy. Arguments: background: string, name of variable (standard_name) which will be plotted as background of trajectories, provided that it can be read with one of the available readers. buffer: float; spatial buffer of plot in degrees of longitude/latitude around particle collection. background: name of variable to be plotted as background field. Use two element list for vector fields, e.g. ['x_wind', 'y_wind'] vmin, vmax: minimum and maximum values for colors of background. linecolor: name of variable to be used for coloring trajectories, or matplotlib color string. lvmin, lvmax: minimum and maximum values for colors of trajectories. lscale (string): resolution of land feature ('c', 'l', 'i', 'h', 'f', 'auto'). default is 'auto'. fast (bool): use some optimizations to speed up plotting at the cost of accuracy :param hide_landmask: do not plot landmask (default False). :type hide_landmask: bool """ mappable = None if self.num_elements_total() == 0: raise ValueError('Please run simulation before plotting') start_time = datetime.now() # x, y are longitude, latitude -> i.e. in a PlateCarree CRS gcrs = ccrs.PlateCarree() fig, ax, crs, x, y, index_of_first, index_of_last = \ self.set_up_map(buffer=buffer, corners=corners, lscale=lscale, fast=fast, hide_landmask=hide_landmask, **kwargs) markercolor = self.plot_comparison_colors[0] # The more elements, the more transparent we make the lines min_alpha = 0.1 max_elements = 5000.0 alpha = min_alpha**(2*(self.num_elements_total()-1)/(max_elements-1)) alpha = np.max((min_alpha, alpha)) if legend is False: legend = None if hasattr(self, 'history') and linewidth != 0: # Plot trajectories from matplotlib.colors import is_color_like if linecolor is None or is_color_like(linecolor) is True: if is_color_like(linecolor): linecolor = linecolor else: linecolor = 'gray' if compare is not None and legend is not None: if legend is True: if hasattr(compare, 'len'): numleg = len(compare) else: numleg = 2 legend = ['Simulation %d' % (i+1) for i in range(numleg)] ax.plot(x[:,0], y[:,0], color=linecolor, alpha=alpha, label=legend[0], linewidth=linewidth, transform = gcrs) ax.plot(x, y, color=linecolor, alpha=alpha, label='_nolegend_', linewidth=linewidth, transform = gcrs) else: ax.plot(x, y, color=linecolor, alpha=alpha, linewidth=linewidth, transform = gcrs) else: #colorbar = True # Color lines according to given parameter try: if isinstance(linecolor, str): param = self.history[linecolor] elif hasattr(linecolor, '__len__'): param = np.tile(linecolor, (self.steps_output, 1)).T else: param = linecolor except: raise ValueError( 'Available parameters to be used for linecolors: ' + str(self.history.dtype.fields)) from matplotlib.collections import LineCollection for i in range(x.shape[1]): vind = np.arange(index_of_first[i], index_of_last[i] + 1) points = np.array( [x[vind, i].T, y[vind, i].T]).T.reshape(-1, 1, 2) segments = np.concatenate([points[:-1], points[1:]], axis=1) if lvmin is None: lvmin = param.min() lvmax = param.max() lc = LineCollection(segments, #cmap=plt.get_cmap('Spectral'), cmap=cmap, norm=plt.Normalize(lvmin, lvmax), transform = gcrs) #lc.set_linewidth(3) lc.set_array(param.T[vind, i]) mappable = ax.add_collection(lc) #axcb = fig.colorbar(lc, ax = ax, orientation = 'horizontal') #try: # Add unit to colorbar if available # colorbarstring = linecolor + ' [%s]' % \ # (self.history_metadata[linecolor]['units']) #except: # colorbarstring = linecolor ##axcb.set_label(colorbarstring) #axcb.set_label(colorbarstring, size=14) #axcb.ax.tick_params(labelsize=14) if compare is None: label_initial = 'initial (%i)' % x.shape[1] label_active = 'active (%i)' % (x.shape[1] - self.num_elements_deactivated()) color_initial = self.status_colors['initial'] color_active = self.status_colors['active'] else: label_initial = None label_active = None color_initial = 'gray' color_active = 'gray' if show_particles is True: if show_initial is True: ax.scatter(x[index_of_first, range(x.shape[1])], y[index_of_first, range(x.shape[1])], s=markersize, zorder=10, edgecolor=markercolor, linewidths=.2, c=color_initial, label=label_initial, transform = gcrs) if surface_color is not None: color_active = surface_color label_active = 'surface' ax.scatter(x[index_of_last, range(x.shape[1])], y[index_of_last, range(x.shape[1])], s=markersize, zorder=3, edgecolor=markercolor, linewidths=.2, c=color_active, label=label_active, transform = gcrs) #if submerged_color is not None: # map.scatter(x[range(x.shape[0]), index_of_last], # y[range(x.shape[0]), index_of_last], s=markersize, # zorder=3, edgecolor=markercolor, linewidths=.2, # c=submerged_color, label='submerged') x_deactivated, y_deactivated = (self.elements_deactivated.lon, self.elements_deactivated.lat) # Plot deactivated elements, labeled by deactivation reason for statusnum, status in enumerate(self.status_categories): if status == 'active': continue # plotted above if status not in self.status_colors: # If no color specified, pick an unused one for color in ['red', 'blue', 'green', 'black', 'gray', 'cyan', 'DarkSeaGreen', 'brown']: if color not in self.status_colors.values(): self.status_colors[status] = color break indices = np.where(self.elements_deactivated.status == statusnum) if len(indices[0]) > 0: if (status == 'seeded_on_land' or status == 'seeded_at_nodata_position'): zorder = 11 else: zorder = 3 if compare is not None: legstr = None else: legstr = '%s (%i)' % (status, len(indices[0])) if compare is None: color_status = self.status_colors[status] else: color_status = 'gray' ax.scatter(x_deactivated[indices], y_deactivated[indices], s=markersize, zorder=zorder, edgecolor=markercolor, linewidths=.1, c=color_status, label=legstr, transform = gcrs) if compare is not None: cd = self._get_comparison_xy_for_plots(compare) for i, c in enumerate(cd): if legend != None: legstr = legend[i+1] else: legstr = None ax.plot(c['x_other'].T[:,0], c['y_other'].T[:,0], self.plot_comparison_colors[i+1] + '-', label=legstr, linewidth=linewidth, transform = gcrs) ax.plot(c['x_other'].T, c['y_other'].T, self.plot_comparison_colors[i+1] + '-', label='_nolegend_', linewidth=linewidth, transform = gcrs) ax.scatter(c['x_other'][range(c['x_other'].shape[0]), c['index_of_last_other']], c['y_other'][range(c['y_other'].shape[0]), c['index_of_last_other']], s=markersize, zorder=3, edgecolor=markercolor, linewidths=.2, c=self.plot_comparison_colors[i+1], transform = gcrs) try: handles, labels = ax.get_legend_handles_labels() if legend is not None and len(handles)>0: plt.legend(loc=legend_loc, markerscale=2) except Exception as e: logger.warning('Cannot plot legend, due to bug in matplotlib:') logger.warning(traceback.format_exc()) if background is not None: if hasattr(self, 'time'): time = self.time - self.time_step_output else: time = None if background == 'residence': scalar,lon_res,lat_res = self.get_residence_time( pixelsize_m=density_pixelsize_m) scalar[scalar==0] = np.nan lon_res, lat_res = np.meshgrid(lon_res[0:-1], lat_res[0:-1]) lon_res = lon_res.T lat_res = lat_res.T map_x, map_y = (lon_res, lat_res) else: map_x, map_y, scalar, u_component, v_component, qmap_x, qmap_y = \ self.get_map_background(ax, background, time=time) #self.time_step_output) if show_scalar is True: if contourlines is False: scalar = np.ma.masked_invalid(scalar) mappable = ax.pcolormesh(map_x, map_y, scalar, alpha=bgalpha, vmin=vmin, vmax=vmax, cmap=cmap, transform = gcrs) else: if contourlines is True: CS = ax.contour(map_x, map_y, scalar, colors='gray', transform = gcrs) else: # contourlines is an array of values CS = ax.contour(map_x, map_y, scalar, contourlines, colors='gray', transform = gcrs) plt.clabel(CS, fmt='%g') if mappable is not None and colorbar is True: cb = fig.colorbar(mappable, orientation='horizontal', pad=.05, aspect=30, shrink=.8, drawedges=False) # TODO: need better control of colorbar content if clabel is not None: cb.set_label(clabel) elif isinstance(linecolor, str) and linecolor != 'gray': cb.set_label(str(linecolor)) if background is not None: cb.set_label(str(background)) if type(background) is list: delta_x = (map_x[1,2] - map_x[1,1])/2. delta_y = (map_y[2,1] - map_y[1,1])/2. ax.quiver(qmap_x[::skip, ::skip] + delta_x, qmap_y[::skip, ::skip] + delta_y, u_component[::skip, ::skip], v_component[::skip, ::skip], scale=scale, transform = gcrs) if lcs is not None: map_x_lcs, map_y_lcs = (lcs['lon'], lcs['lat']) ax.pcolormesh( map_x_lcs, map_y_lcs, lcs['ALCS'][0,:,:], alpha=1, vmin=vmin, vmax=vmax, cmap=cmap, transform = gcrs) if title is not None: if title == 'auto': if hasattr(self, 'time'): plt.title('%s\n%s to %s UTC (%i steps)' % ( self._figure_title(), self.start_time.strftime('%Y-%m-%d %H:%M'), self.time.strftime('%Y-%m-%d %H:%M'), self.steps_output)) else: plt.title('%s\n%i elements seeded at %s UTC' % ( self._figure_title(), self.num_elements_scheduled(), self.elements_scheduled_time[0].strftime( '%Y-%m-%d %H:%M'))) else: plt.title(title) if trajectory_dict is not None: self._plot_trajectory_dict(ax, trajectory_dict) #plt.gca().tick_params(labelsize=14) #fig.canvas.draw() #fig.set_tight_layout(True) if filename is not None: plt.savefig(filename) logger.info('Time to make plot: ' + str(datetime.now() - start_time)) else: if show is True: plt.show() return ax, plt
[docs] def _substance_name(self): return None
[docs] def _figure_title(self): if self._substance_name() is None: return 'OpenDrift - ' + type(self).__name__ else: return 'OpenDrift - ' + type(self).__name__ + ' (%s)' % self._substance_name()
[docs] def _plot_trajectory_dict(self, ax, trajectory_dict): '''Plot provided trajectory along with simulated''' time = trajectory_dict['time'] time = np.array(time) i = np.where((time>=self.start_time) & (time<=self.time))[0] x, y = (np.atleast_1d(trajectory_dict['lon'])[i], np.atleast_1d(trajectory_dict['lat'])[i]) ls = trajectory_dict['linestyle'] gcrs = ccrs.PlateCarree() ax.plot(x, y, ls, linewidth=2, transform = gcrs) ax.plot(x[0], y[0], 'ok', transform = gcrs) ax.plot(x[-1], y[-1], 'xk', transform = gcrs)
[docs] def get_map_background(self, ax, background, time=None): # Get background field for plotting on map or animation # TODO: this method should be made more robust if type(background) is list: variable = background[0] # A vector is requested else: variable = background # A scalar is requested for readerName in self.readers: reader = self.readers[readerName] if variable in reader.variables: if time is None or reader.start_time is None or ( time>= reader.start_time and time <= reader.end_time) or ( reader.always_valid is True): break if time is None: if hasattr(self, 'elements_scheduled_time'): # Using time of first seeded element time = self.elements_scheduled_time[0] # Get reader coordinates covering given map area axisproj = pyproj.Proj(ax.projection.proj4_params) xmin, xmax, ymin, ymax = ax.get_extent(ccrs.PlateCarree()) cornerlons = np.array([xmin, xmin, xmax, xmax]) cornerlats = np.array([ymin, ymax, ymin, ymax]) reader_x, reader_y = reader.lonlat2xy(cornerlons, cornerlats) if sum(~np.isfinite(reader_x+reader_y)) > 0: # Axis corner points are not within reader domain reader_x = np.array([reader.xmin, reader.xmax]) reader_y = np.array([reader.ymin, reader.ymax]) else: reader_x = np.linspace(reader_x.min(), reader_x.max(), 10) reader_y = np.linspace(reader_y.min(), reader_y.max(), 10) data = reader.get_variables( background, time, reader_x, reader_y, None) reader_x, reader_y = np.meshgrid(data['x'], data['y']) if type(background) is list: # Ensemble reader, using first member u_component = data[background[0]] v_component = data[background[1]] if isinstance(u_component, list): u_component = u_component[0] v_component = v_component[0] with np.errstate(invalid='ignore'): scalar = np.sqrt(u_component**2 + v_component**2) # NB: rotation not completed! u_component, v_component = reader.rotate_vectors( reader_x, reader_y, u_component, v_component, reader.proj, ccrs.PlateCarree( globe=ccrs.Globe(datum='WGS84', ellipse='WGS84') ).proj4_init) else: scalar = data[background] if isinstance(scalar, list): # Ensemble reader, using first member scalar = scalar[0] u_component = v_component = None # Shift one pixel for correct plotting reader_x = reader_x - reader.delta_x reader_y = reader_y - reader.delta_y if reader.projected is False: reader_y[reader_y<0] = 0 reader_x[reader_x<0] = 0 rlons, rlats = reader.xy2lonlat(reader_x, reader_y) qrlons, qrlats = reader.xy2lonlat(reader_x+reader.delta_x/2, reader_y+reader.delta_y/2) if rlons.max() > 360: rlons = rlons - 360 qrlons = qrlons - 360 map_x, map_y = (rlons, rlats) qmap_x, qmap_y = (qrlons, qrlats) scalar = np.ma.masked_invalid(scalar) return map_x, map_y, scalar, u_component, v_component, qmap_x, qmap_y
[docs] def get_density_timeseries(self, lon, lat): """Get timeseries at a point from pre-calculated density field""" if not hasattr(self, 'analysis_file'): raise ValueError('Density has to be calculated with get_density_array or animation') self.af = xr.open_dataset(self.analysis_file) lon_bin = self.af.lon_bin.values lat_bin = self.af.lat_bin.values lon = np.atleast_1d(lon) lat = np.atleast_1d(lat) if lon.min()<lon_bin.min() or lon.max()>lon_bin.max(): raise ValueError('Longitude %s is outside range %s to %s' % (lon, lon_bin.min(), lon_bin.max())) if lat.min()<lat_bin.min() or lat.max()>lat_bin.max(): raise ValueError('Latitude %s is outside range %s to %s' % (lat, lat_bin.min(), lat_bin.max())) if len(lon) == 1: t_total = self.af.density.sel(lon_bin=lon, lat_bin=lat, method='nearest') t_origin_marker = self.af.density_origin_marker.sel(lon_bin=lon, lat_bin=lat, method='nearest') elif len(lon) == 2: t_total = self.af.density.sel(lon_bin=slice(lon[0], lon[1]), lat_bin=slice(lat[0], lat[1])) t_total = t_total.sum(('lon_bin', 'lat_bin')) t_origin_marker = self.af.density_origin_marker.sel(lon_bin=slice(lon[0], lon[1]), lat_bin=slice(lat[0], lat[1])) t_origin_marker = t_origin_marker.sum(('lon_bin', 'lat_bin')) else: raise ValueError('Lon and lat must have length 1 (point) or 2 (min, max)') return t_total, t_origin_marker
[docs] def get_density_xarray(self, pixelsize_m, weights=None, per_origin_marker=True, dominating_origin_marker=False): from xhistogram.xarray import histogram if os.path.exists(self.analysis_file): self.ads = xr.open_dataset(self.analysis_file) # Import from file logger.info('Importing from saved analysis file') lon_bin = self.ads.lon_bin.values lat_bin = self.ads.lat_bin.values density = self.ads.density.values if 'density_origin_marker' in self.ads: density_origin_marker = self.ads.density_origin_marker.values else: density_origin_marker = None self.ads.close() return density, density_origin_marker, lon_bin, lat_bin else: # Calculate self.ads = xr.Dataset() # Dataset to store all data, to be saved to self.analysis_file start = datetime.now() logger.info('Calculating density array, this may take some time...') deltalat = pixelsize_m/111000.0 # m to degrees if hasattr(self, 'lonmin'): lonmin = self.lonmin lonmax = self.lonmax latmin = self.latmin latmax = self.latmax else: logger.debug('Finding min and max of lon and lat...') lonmin = np.nanmin(self.ds.lon) lonmax = np.nanmax(self.ds.lon) latmin = np.nanmin(self.ds.lat) latmax = np.nanmax(self.ds.lat) self.ads.attrs['lonmin'] = lonmin self.ads.attrs['lonmax'] = lonmax self.ads.attrs['latmin'] = latmin self.ads.attrs['latmax'] = latmax self.ads.coords['time'] = self.ds.coords['time'] deltalon = deltalat/np.cos(np.radians((latmin+latmax)/2)) latbin = np.arange(latmin-deltalat, latmax+deltalat, deltalat) lonbin = np.arange(lonmin-deltalon, lonmax+deltalon, deltalon) if weights is not None: if isinstance(weights, str): # element property # We do not add 2D weights to analysis file, to save space self.ads.attrs['weights'] = weights weights = self.get_property(weights)[0] elif not hasattr(weights, '__len__'): # scalar weights = xr.DataArray(weights*xr.ones_like(self.ds.lon), dims=['trajectory'], coords={'trajectory': self.ds.coords['trajectory']}) self.ads['weights'] = weights else: # provided array weights = xr.DataArray(weights, dims=['trajectory'], coords={'trajectory': self.ds.coords['trajectory']}) self.ads['weights'] = weights if per_origin_marker is True: max_om = self.ds.origin_marker.max().compute().values origin_marker = np.arange(max_om+1) self.ads.coords['origin_marker'] = origin_marker density_om = xr.DataArray(np.zeros((len(self.ds.coords['time']), len(lonbin)-1, len(latbin)-1, len(origin_marker))), name='density_origin_marker', dims=('time', 'lon_bin', 'lat_bin', 'origin_marker')) for om in range(max_om+1): logger.info('\tcalculating for origin_marker %s...' % om) h = histogram(self.ds.lon.where(self.ds.origin_marker==om), self.ds.lat.where(self.ds.origin_marker==om), bins=[lonbin, latbin], weights=weights, dim=['trajectory']) lon_bin = h.coords['lon_bin'] lat_bin = h.coords['lat_bin'] density_om[:, :, :, om] = h.copy() self.ads['density_origin_marker'] = density_om # Computing total density historgram from om-components density = density_om.sum(dim='origin_marker') self.ads['density'] = density if dominating_origin_marker is True: logger.info('Calculating dominating origin_marker') self.ads['dominating_origin_marker'] = density_om.argmax( dim='origin_marker', skipna=True) self.ads['dominating_origin_marker'] = \ self.ads['dominating_origin_marker'].where(density>0) else: # Adding total density historgram density = histogram(self.ds.lon, self.ds.lat, bins=[lonbin, latbin], weights=weights, dim=['trajectory']) lon_bin = density.coords['lon_bin'] lat_bin = density.coords['lat_bin'] self.ads['density'] = density density_om = None self.ads.coords['lon_bin'] = lon_bin self.ads.coords['lat_bin'] = lat_bin self.ads.to_netcdf(self.analysis_file) logger.info('Time to calculate density array: %s' % (datetime.now() - start)) return density, density_om, lonbin, latbin
[docs] def get_density_array(self, pixelsize_m, weight=None): lon = self.get_property('lon')[0] lat = self.get_property('lat')[0] times = self.get_time_array()[0] deltalat = pixelsize_m/111000.0 # m to degrees deltalon = deltalat/np.cos(np.radians((np.nanmin(lat) + np.nanmax(lat))/2)) lat_array = np.arange(np.nanmin(lat)-deltalat, np.nanmax(lat)+deltalat, deltalat) lon_array = np.arange(np.nanmin(lon)-deltalat, np.nanmax(lon)+deltalon, deltalon) bins=(lon_array, lat_array) z = self.get_property('z')[0] if weight is not None: weight_array = self.get_property(weight)[0] status = self.get_property('status')[0] lon_submerged = lon.copy() lat_submerged = lat.copy() lon_stranded = lon.copy() lat_stranded = lat.copy() lon_submerged[z>=0] = 1000 lat_submerged[z>=0] = 1000 lon[z<0] = 1000 lat[z<0] = 1000 H = np.zeros((len(times), len(lon_array) - 1, len(lat_array) - 1))#.astype(int) H_submerged = H.copy() H_stranded = H.copy() try: strandnum = self.status_categories.index('stranded') lon_stranded[status!=strandnum] = 1000 lat_stranded[status!=strandnum] = 1000 contains_stranded = True except ValueError: contains_stranded = False for i in range(len(times)): if weight is not None: weights = weight_array[i,:] else: weights = None H[i,:,:], dummy, dummy = \ np.histogram2d(lon[i,:], lat[i,:], weights=weights, bins=bins) H_submerged[i,:,:], dummy, dummy = \ np.histogram2d(lon_submerged[i,:], lat_submerged[i,:], weights=weights, bins=bins) if contains_stranded is True: H_stranded[i,:,:], dummy, dummy = \ np.histogram2d(lon_stranded[i,:], lat_stranded[i,:], weights=weights, bins=bins) return H, H_submerged, H_stranded, lon_array, lat_array
[docs] def get_density_array_proj(self, pixelsize_m, density_proj=None, llcrnrlon=None,llcrnrlat=None,urcrnrlon=None,urcrnrlat=None, weight=None): # # TODO: should be merged with get_density_array # KFD Jan 2021 # lon = self.get_property('lon')[0] lat = self.get_property('lat')[0] times = self.get_time_array()[0] #deltalat = pixelsize_m/111000.0 # m to degrees #deltalon = deltalat/np.cos(np.radians((np.nanmin(lat) + # np.nanmax(lat))/2)) #lat_array = np.arange(np.nanmin(lat)-deltalat, # np.nanmax(lat)+deltalat, deltalat) #lon_array = np.arange(np.nanmin(lon)-deltalat, # np.nanmax(lon)+deltalon, deltalon) #bins=(lon_array, lat_array) if density_proj is None: # add default projection with equal-area property density_proj = pyproj.Proj('+proj=moll +ellps=WGS84 +lon_0=0.0') density_proj = pyproj.Proj('+proj=longlat +a=6371229 +no_defs') # create a grid in the specified projection x,y = density_proj(lon, lat) if llcrnrlon is not None: llcrnrx,llcrnry = density_proj(llcrnrlon,llcrnrlat) urcrnrx,urcrnry = density_proj(urcrnrlon,urcrnrlat) else: llcrnrx,llcrnry = x.min()-pixelsize_m, y.min()-pixelsize_m urcrnrx,urcrnry = x.max()+pixelsize_m, y.max()+pixelsize_m x_array = np.arange(llcrnrx,urcrnrx, pixelsize_m) y_array = np.arange(llcrnry,urcrnry, pixelsize_m) bins=(x_array, y_array) outsidex, outsidey = max(x_array)*1.5, max(y_array)*1.5 z = self.get_property('z')[0] if weight is not None: weight_array = self.get_property(weight)[0] status = self.get_property('status')[0] #lon_submerged = lon.copy() #lat_submerged = lat.copy() #lon_stranded = lon.copy() #lat_stranded = lat.copy() #lon_submerged[z>=0] = 1000 #lat_submerged[z>=0] = 1000 #lon[z<0] = 1000 #lat[z<0] = 1000 #H = np.zeros((len(times), len(lon_array) - 1, # len(lat_array) - 1))#.astype(int) x_submerged = x.copy() y_submerged = y.copy() x_stranded = x.copy() y_stranded = y.copy() x_submerged[z>=0] = outsidex y_submerged[z>=0] = outsidey x[z<0] = outsidex y[z<0] = outsidey H = np.zeros((len(times), len(x_array) - 1, len(y_array) - 1))#.astype(int) H_submerged = H.copy() H_stranded = H.copy() try: strandnum = self.status_categories.index('stranded') #lon_stranded[status!=strandnum] = 1000 #lat_stranded[status!=strandnum] = 1000 x_stranded[status!=strandnum] = outsidex y_stranded[status!=strandnum] = outsidey contains_stranded = True except ValueError: contains_stranded = False for i in range(len(times)): if weight is not None: weights = weight_array[i,:] else: weights = None H[i,:,:], dummy, dummy = \ np.histogram2d(x[i,:], y[i,:], weights=weights, bins=bins) H_submerged[i,:,:], dummy, dummy = \ np.histogram2d(x_submerged[i,:], y_submerged[i,:], weights=weights, bins=bins) if contains_stranded is True: H_stranded[i,:,:], dummy, dummy = \ np.histogram2d(x_stranded[i,:], y_stranded[i,:], weights=weights, bins=bins) if density_proj is not None: Y,X = np.meshgrid(y_array, x_array) lon_array, lat_array = density_proj(X,Y,inverse=True) return H, H_submerged, H_stranded, lon_array, lat_array
[docs] def get_residence_time(self, pixelsize_m): H,H_sub, H_str,lon_array,lat_array = \ self.get_density_array(pixelsize_m) residence = np.sum(H, axis=0) return residence, lon_array, lat_array
[docs] def write_netcdf_density_map(self, filename, pixelsize_m='auto'): '''Write netCDF file with map of particles densities''' if pixelsize_m == 'auto': lon, lat = self.get_lonlats() latspan = lat.max()-lat.min() pixelsize_m=30 if latspan > .05: pixelsize_m = 50 if latspan > .1: pixelsize_m = 300 if latspan > .3: pixelsize_m = 500 if latspan > .7: pixelsize_m = 1000 if latspan > 2: pixelsize_m = 2000 if latspan > 5: pixelsize_m = 4000 H, H_submerged, H_stranded, lon_array, lat_array = \ self.get_density_array(pixelsize_m) lon_array = (lon_array[0:-1] + lon_array[1::])/2 lat_array = (lat_array[0:-1] + lat_array[1::])/2 nc = Dataset(filename, 'w') nc.createDimension('lon', len(lon_array)) nc.createDimension('lat', len(lat_array)) nc.createDimension('time', H.shape[0]) times = self.get_time_array()[0] timestr = 'seconds since 1970-01-01 00:00:00' nc.createVariable('time', 'f8', ('time',)) nc.variables['time'][:] = date2num(times, timestr) nc.variables['time'].units = timestr nc.variables['time'].standard_name = 'time' # Projection nc.createVariable('projection_lonlat', 'i8') nc.variables['projection_lonlat'].grid_mapping_name = \ 'latitude_longitude' nc.variables['projection_lonlat'].earth_radius = 6371229. nc.variables['projection_lonlat'].proj4 = \ '+proj=longlat +a=6371229 +no_defs' # Coordinates nc.createVariable('lon', 'f8', ('lon',)) nc.createVariable('lat', 'f8', ('lat',)) nc.variables['lon'][:] = lon_array nc.variables['lon'].long_name = 'longitude' nc.variables['lon'].short_name = 'longitude' nc.variables['lon'].units = 'degrees_east' nc.variables['lat'][:] = lat_array nc.variables['lat'].long_name = 'latitude' nc.variables['lat'].short_name = 'latitude' nc.variables['lat'].units = 'degrees_north' # Density nc.createVariable('density_surface', 'u1', ('time','lat', 'lon')) H = np.swapaxes(H, 1, 2).astype('uint8') H = np.ma.masked_where(H==0, H) nc.variables['density_surface'][:] = H nc.variables['density_surface'].long_name = 'Detection probability' nc.variables['density_surface'].grid_mapping = 'projection_lonlat' nc.variables['density_surface'].units = '1' # Density submerged nc.createVariable('density_submerged', 'u1', ('time','lat', 'lon')) H_sub = np.swapaxes(H_submerged, 1, 2).astype('uint8') H_sub = np.ma.masked_where(H_sub==0, H_sub) nc.variables['density_submerged'][:] = H_sub nc.variables['density_submerged'].long_name = 'Detection probability submerged' nc.variables['density_submerged'].grid_mapping = 'projection_lonlat' nc.variables['density_submerged'].units = '1' # Density stranded nc.createVariable('density_stranded', 'u1', ('time','lat', 'lon')) H_stranded = np.swapaxes(H_stranded, 1, 2).astype('uint8') H_stranded = np.ma.masked_where(H_stranded==0, H_stranded) nc.variables['density_stranded'][:] = H_stranded nc.variables['density_stranded'].long_name = 'Detection probability stranded' nc.variables['density_stranded'].grid_mapping = 'projection_lonlat' nc.variables['density_stranded'].units = '1' nc.close()
[docs] def write_netcdf_density_map_proj(self, filename, pixelsize_m='auto', density_proj=None, llcrnrlon=None, llcrnrlat=None, urcrnrlon=None, urcrnrlat=None): '''Write netCDF file with map of particles densities for a given projection or area''' # # TODO: should be merged with write_netcdf_density_map_proj # KFD Jan 2021 # if pixelsize_m == 'auto': lon, lat = self.get_lonlats() latspan = lat.max()-lat.min() pixelsize_m=30 if latspan > .05: pixelsize_m = 50 if latspan > .1: pixelsize_m = 300 if latspan > .3: pixelsize_m = 500 if latspan > .7: pixelsize_m = 1000 if latspan > 2: pixelsize_m = 2000 if latspan > 5: pixelsize_m = 4000 if density_proj is None: # add default projection with equal-area property density_proj = pyproj.Proj('+proj=moll +ellps=WGS84 +lon_0=0.0') H, H_submerged, H_stranded, lon_array, lat_array = \ self.get_density_array_proj(pixelsize_m=pixelsize_m, density_proj=density_proj, llcrnrlon=llcrnrlon, llcrnrlat=llcrnrlat, urcrnrlon=urcrnrlon, urcrnrlat=urcrnrlat) # calculate center coordinates print(lon_array.shape, lat_array.shape) lon_array = (lon_array[:-1,:-1] + lon_array[1:,1:])/2. lat_array = (lat_array[:-1,:-1] + lat_array[1:,1:])/2. from netCDF4 import Dataset, date2num nc = Dataset(filename, 'w') nc.createDimension('x', lon_array.shape[0]) nc.createDimension('y', lon_array.shape[1]) nc.createDimension('time', H.shape[0]) times = self.get_time_array()[0] timestr = 'seconds since 1970-01-01 00:00:00' nc.createVariable('time', 'f8', ('time',)) nc.variables['time'][:] = date2num(times, timestr) nc.variables['time'].units = timestr nc.variables['time'].standard_name = 'time' # Projection nc.createVariable('projection', 'i8') nc.variables['projection'].proj4 = density_proj.definition_string() # Coordinates nc.createVariable('lon', 'f8', ('y','x')) nc.createVariable('lat', 'f8', ('y','x')) nc.variables['lon'][:] = lon_array.T nc.variables['lon'].long_name = 'longitude' nc.variables['lon'].short_name = 'longitude' nc.variables['lon'].units = 'degrees_east' nc.variables['lat'][:] = lat_array.T nc.variables['lat'].long_name = 'latitude' nc.variables['lat'].short_name = 'latitude' nc.variables['lat'].units = 'degrees_north' # Density nc.createVariable('density_surface', 'u1', ('time','y', 'x')) H = np.swapaxes(H, 1, 2).astype('uint8') H = np.ma.masked_where(H==0, H) nc.variables['density_surface'][:] = H nc.variables['density_surface'].long_name = 'Detection probability' nc.variables['density_surface'].grid_mapping = 'projection' nc.variables['density_surface'].units = '1' # Density submerged nc.createVariable('density_submerged', 'u1', ('time','y', 'x')) H_sub = np.swapaxes(H_submerged, 1, 2).astype('uint8') H_sub = np.ma.masked_where(H_sub==0, H_sub) nc.variables['density_submerged'][:] = H_sub nc.variables['density_submerged'].long_name = 'Detection probability submerged' nc.variables['density_submerged'].grid_mapping = 'projection' nc.variables['density_submerged'].units = '1' # Density stranded nc.createVariable('density_stranded', 'u1', ('time','y', 'x')) H_stranded = np.swapaxes(H_stranded, 1, 2).astype('uint8') H_stranded = np.ma.masked_where(H_stranded==0, H_stranded) nc.variables['density_stranded'][:] = H_stranded nc.variables['density_stranded'].long_name = 'Detection probability stranded' nc.variables['density_stranded'].grid_mapping = 'projection' nc.variables['density_stranded'].units = '1' nc.close()
[docs] def write_geotiff(self, filename, pixelsize_km=.2): '''Write one GeoTiff image per timestep. filename should contain date identifiers, e.g. 'img_%Y%m%d_%H%M.tif' https://docs.python.org/2/library/datetime.html#strftime-and-strptime-behavior ''' try: from osgeo import gdal, osr except: raise ValueError('GDAL is needed to write geotiff images.') import matplotlib.pyplot as plt driver = gdal.GetDriverByName('GTiff') srs = osr.SpatialReference() srs.ImportFromEPSG(4326) colortable = gdal.ColorTable() colortable.SetColorEntry(0,(255,255,255, 0)) colortable.SetColorEntry(1,(0,0,0, 255)) colortable.SetColorEntry(2,(255,0,0, 255)) colortable.SetColorEntry(3,(0,255,0, 255)) colortable.SetColorEntry(4,(0,0,255, 255)) lon = self.get_property('lon')[0] lat = self.get_property('lat')[0] status = self.get_property('status')[0] times = self.get_time_array()[0] deltalat = pixelsize_km/111.0 # km to degrees deltalon = deltalat/np.cos(np.radians((lat.min() + lat.max())/2)) lat_array = np.arange(lat.min()-deltalat, lat.max()+deltalat, deltalat) lon_array = np.arange(lon.min()-deltalat, lon.max()+deltalon, deltalon) ilon = (np.round((lon-lon.min())/deltalon)).astype(int) ilat = (np.round((lat-lat.min())/deltalat)).astype(int) # Setting masked values to zero, for use as indices ilon[ilon.mask] = 0 ilat[ilat.mask] = 0 status[ilon.mask] = 0 image = np.zeros((len(times), len(lon_array), len(lat_array))).astype(int) geotransform = [lon_array.min(), deltalon, 0, lat_array.max(), 0, -deltalat] for i, t in enumerate(times): image[i, ilon[i,:], ilat[i,:]] = status[i, :] + 1 filename_i = t.strftime(filename) ds = driver.Create(filename_i, len(lon_array), len(lat_array), 1, gdal.GDT_Byte, ) ds.SetProjection(srs.ExportToWkt()) ds.SetGeoTransform(geotransform) outband=ds.GetRasterBand(1) outband.SetNoDataValue(0) outband.WriteArray(np.fliplr(image[i, :, :]).transpose()) outband.SetColorTable(colortable) ds = None
[docs] def get_time_array(self): """Return a list of output times of last run.""" # Making sure start_time is datetime, and not cftime object self.start_time = datetime(self.start_time.year, self.start_time.month, self.start_time.day, self.start_time.hour, self.start_time.minute, self.start_time.second) td = self.time_step_output time_array = [self.start_time + td*i for i in range(self.steps_output)] time_array_relative = [td*i for i in range(self.steps_output)] return time_array, time_array_relative
[docs] def plot_environment(self, filename=None, ax=None, show=True): """Plot mean wind and current velocities of element of last run.""" x_wind = self.get_property('x_wind')[0] y_wind = self.get_property('y_wind')[0] wind = np.sqrt(x_wind**2 + y_wind**2) x_sea_water_velocity = self.get_property('x_sea_water_velocity')[0] y_sea_water_velocity = self.get_property('y_sea_water_velocity')[0] current = np.sqrt(x_sea_water_velocity**2 + y_sea_water_velocity**2) wind = np.ma.mean(wind, axis=1) current = np.ma.mean(current, axis=1) time, time_relative = self.get_time_array() time = np.array([t.total_seconds()/3600. for t in time_relative]) if ax is None: fig = plt.figure() ax = fig.add_subplot(111) ax.plot(time, wind, 'b', label='Wind speed') ax.set_ylabel('Wind speed [m/s]', color='b') ax.set_xlim([0, time[-1]]) ax.set_ylim([0, wind.max()*1.1]) ax2 = ax.twinx() ax2.plot(time, current, 'r', label='Current speed') ax2.set_ylabel('Current speed [m/s]', color='r') ax2.set_xlim([0, time[-1]]) ax2.set_ylim([0, current.max()*1.1]) for tl in ax.get_yticklabels(): tl.set_color('b') for tl in ax2.get_yticklabels(): tl.set_color('r') ax.set_xlabel('Time [hours]') ax.legend(loc='upper left') ax2.legend(loc='lower right') if filename is None: if show is True: plt.show() else: plt.savefig(filename)
[docs] def plot_property(self, prop, filename=None, mean=False): """Basic function to plot time series of any element properties.""" import matplotlib.pyplot as plt from matplotlib import dates hfmt = dates.DateFormatter('%d %b %Y %H:%M') fig = plt.figure() ax = fig.gca() ax.xaxis.set_major_formatter(hfmt) plt.xticks(rotation='vertical') start_time = self.start_time # In case start_time is unsupported cftime start_time = datetime(start_time.year, start_time.month, start_time.day, start_time.hour, start_time.minute, start_time.second) times = [start_time + n*self.time_step_output for n in range(self.steps_output)] data = self.history[prop].T[0:len(times), :] if mean is True: # Taking average over elements data = np.mean(data, axis=1) plt.plot(times, data) plt.title(prop) plt.xlabel('Time [UTC]') try: plt.ylabel('%s [%s]' % (prop, self.elements.variables[prop]['units'])) except: plt.ylabel(prop) plt.subplots_adjust(bottom=.3) plt.grid() if filename is None: plt.show() else: plt.savefig(filename)
[docs] def get_property(self, propname): """Get property from history, sorted by status.""" index_of_first, index_of_last = \ self.index_of_activation_and_deactivation() prop = self.history[propname].copy() status = self.history['status'].copy() j = np.arange(status.shape[1]) # Fill arrays with last value before deactivation for i in range(status.shape[0]): status[i, j > index_of_last[i]] = status[i, index_of_last[i]] prop[i, j > index_of_last[i]] = prop[i, index_of_last[i]] return prop.T, status.T
[docs] def get_trajectory_lengths(self): """Calculate lengths and speeds along trajectories.""" lons = self.get_property('lon')[0] lats = self.get_property('lat')[0] geod = pyproj.Geod(ellps='WGS84') a1, a2, distances = geod.inv(lons[0:-1,:], lats[0:-1,:], lons[1::,:], lats[1::,:]) distances[np.isnan(distances)] = 0 speeds = distances / self.time_step_output.total_seconds() distances[speeds>100] = 0 # TODO: need better way to mask invalid distances speeds[speeds>100] = 0 # due to masked lons/lats arrays total_length = np.cumsum(distances, 0)[-1,:] return total_length, distances, speeds
[docs] def update_positions(self, x_vel, y_vel): """Move particles according to given velocity components. This method shall account for projection metrics (a distance on a map projection does not necessarily correspond to the same distance over true ground (not yet implemented). Arguments: x_vel and v_vel: floats, velocities in m/s of particle along x- and y-axes of the inherit SRS (proj4). """ geod = pyproj.Geod(ellps='WGS84') azimuth = np.degrees(np.arctan2(x_vel, y_vel)) # Direction of motion velocity = np.sqrt(x_vel**2 + y_vel**2) # Velocity in m/s velocity = velocity * self.elements.moving # Do not move frosen elements # Calculate new positions self.elements.lon, self.elements.lat, back_az = geod.fwd( self.elements.lon, self.elements.lat, azimuth, velocity*self.time_step.total_seconds()) # Check that new positions are valid if (self.elements.lon.min() < -180) or ( self.elements.lon.min() > 360) or ( self.elements.lat.min() < -90) or ( self.elements.lat.max() > 90): logger.info('Invalid new coordinates:') logger.info(self.elements) sys.exit('Quitting')
[docs] def __repr__(self): """String representation providing overview of model status.""" outStr = '===========================\n' if hasattr(self, 'history'): outStr += self.performance() outStr += '===========================\n' outStr += 'Model:\t' + type(self).__name__ + \ ' (OpenDrift version %s)\n' % opendrift.__version__ outStr += '\t%s active %s particles (%s deactivated, %s scheduled)\n'\ % (self.num_elements_active(), self.ElementType.__name__, self.num_elements_deactivated(), self.num_elements_scheduled()) variable_groups, reader_groups, missing = self.get_reader_groups() outStr += '-------------------\n' outStr += 'Environment variables:\n' for i, variableGroup in enumerate(variable_groups): outStr += ' -----\n' readerGroup = reader_groups[i] for variable in sorted(variableGroup): outStr += ' ' + variable + '\n' for i, reader in enumerate(readerGroup): outStr += ' ' + str(i+1) + ') ' + reader + '\n' if len(self.missing_variables()) > 0: outStr += ' -----\n' outStr += 'Readers not added for the following variables:\n' for variable in sorted(self.missing_variables()): outStr += ' ' + variable + '\n' lazy_readers = [r for r in self.readers if self.readers[r].is_lazy is True] if len(lazy_readers) > 0: outStr += '---\nLazy readers:\n' for lr in lazy_readers: outStr += ' ' + lr + '\n' if hasattr(self, 'discarded_readers'): outStr += '\nDiscarded readers:\n' for dr in self.discarded_readers: outStr += ' ' + dr + '\n' if hasattr(self, 'time'): outStr += '\nTime:\n' outStr += '\tStart: %s\n' % (self.start_time) outStr += '\tPresent: %s\n' % (self.time) if hasattr(self, 'time_step'): outStr += '\tCalculation steps: %i * %s - total time: %s\n' % ( self.steps_calculation, self.time_step, self.time-self.start_time) outStr += '\tOutput steps: %i * %s\n' % ( self.steps_output, self.time_step_output) if hasattr(self, 'messages'): outStr += '-------------------\n' outStr += self.get_messages() outStr += '===========================\n' return outStr
[docs] def store_message(self, message): """Store important messages to be displayed to user at end.""" if not hasattr(self, 'messages'): self.messages = [] self.messages.append(message)
[docs] def get_messages(self): """Report any messages stored during simulation.""" if hasattr(self, 'messages'): return str(self.messages).strip('[]') + '\n' else: return ''
[docs] def add_halo_readers(self): """Adding some Thredds and file readers in prioritised order""" self.add_readers_from_file(self.test_data_folder() + '../../opendrift/scripts/data_sources.txt')
[docs] def _save_animation(self, anim, filename, fps): if 'sphinx_gallery' in sys.modules: # This assumes that the calling script is two frames up in the stack. If # _save_animation is called through a more deeply nested method, it will # not give the correct result. caller = inspect.stack()[2] caller = os.path.splitext(os.path.basename(caller.filename))[0] # Calling script is string input (e.g. from ..plot::) if caller == '<string>': caller = 'plot_directive' adir = os.path.realpath('../source/gallery/animations') else: adir = os.path.realpath('../docs/source/gallery/animations') if not hasattr(OpenDriftSimulation, '__anim_no__'): OpenDriftSimulation.__anim_no__ = { } if caller not in OpenDriftSimulation.__anim_no__: OpenDriftSimulation.__anim_no__[caller] = 0 os.makedirs(adir, exist_ok=True) filename = '%s_%d.gif' % (caller, OpenDriftSimulation.__anim_no__[caller]) OpenDriftSimulation.__anim_no__[caller] += 1 filename = os.path.join(adir, filename) logger.info('Saving animation to ' + filename + '...') try: if filename[-4:] == '.gif': # GIF logger.info('Making animated gif...') anim.save(filename, fps=fps, writer='imagemagick') else: # MP4 try: FFwriter=animation.FFMpegWriter(fps=fps, codec='libx264', bitrate=1800, extra_args=['-profile:v', 'baseline', '-vf', 'crop=trunc(iw/2)*2:trunc(ih/2)*2', # cropping 1 pixel if not even '-pix_fmt', 'yuv420p', '-an']) anim.save(filename, writer=FFwriter) except Exception as e: logger.info(e) anim.save(filename, fps=fps) logger.warning('Animation might not be HTML5 compatible.') except Exception as e: logger.info('Could not save animation:') logger.info(e) logger.debug(traceback.format_exc()) if 'sphinx_gallery' in sys.modules: plt.close()
[docs] def calculate_ftle(self, reader=None, delta=None, domain=None, time=None, time_step=None, duration=None, z=0, RLCS=True, ALCS=True): if reader is None: logger.info('No reader provided, using first available:') reader = list(self.readers.items())[0][1] logger.info(reader.name) if isinstance(reader, pyproj.Proj): proj = reader elif isinstance(reader,str): proj = pyproj.Proj(reader) else: proj = reader.proj from opendrift.models.physics_methods import ftle if not isinstance(duration, timedelta): duration = timedelta(seconds=duration) if domain==None: xs = np.arange(reader.xmin, reader.xmax, delta) ys = np.arange(reader.ymin, reader.ymax, delta) else: xmin, xmax, ymin, ymax = domain xs = np.arange(xmin, xmax, delta) ys = np.arange(ymin, ymax, delta) X, Y = np.meshgrid(xs, ys) lons, lats = proj(X, Y, inverse=True) if time is None: time = reader.start_time if not isinstance(time, list): time = [time] # dictionary to hold LCS calculation lcs = {'time': time, 'lon': lons, 'lat':lats} lcs['RLCS'] = np.zeros((len(time), len(ys), len(xs))) lcs['ALCS'] = np.zeros((len(time), len(ys), len(xs))) T = np.abs(duration.total_seconds()) for i, t in enumerate(time): logger.info('Calculating LCS for ' + str(t)) # Forwards if RLCS is True: self.reset() self.seed_elements(lons.ravel(), lats.ravel(), time=t, z=z) self.run(duration=duration, time_step=time_step) f_x1, f_y1 = proj( self.history['lon'].T[-1].reshape(X.shape), self.history['lat'].T[-1].reshape(X.shape)) lcs['RLCS'][i,:,:] = ftle(f_x1-X, f_y1-Y, delta, T) # Backwards if ALCS is True: self.reset() self.seed_elements(lons.ravel(), lats.ravel(), time=t+duration, z=z) self.run(duration=duration, time_step=-time_step) b_x1, b_y1 = proj( self.history['lon'].T[-1].reshape(X.shape), self.history['lat'].T[-1].reshape(X.shape)) lcs['ALCS'][i,:,:] = ftle(b_x1-X, b_y1-Y, delta, T) lcs['RLCS'] = np.ma.masked_invalid(lcs['RLCS']) lcs['ALCS'] = np.ma.masked_invalid(lcs['ALCS']) # Flipping ALCS left-right. Not sure why this is needed lcs['ALCS'] = lcs['ALCS'][:,::-1,::-1] return lcs
[docs] def center_of_gravity(self, onlysurface=False): """ calculate center of mass and variance of all elements returns (lon,lat), variance where (lon,lat) are the coordinates of the center of mass as function of time""" #lon,lat = self.get_property('lon')[0], self.get_property('lat')[0] lon,lat = self.history['lon'], self.history['lat'] x,y = self.proj_latlon(lon, lat) if onlysurface==True: z = self.history['z'] submerged = z < 0 x = np.ma.array(x, mask=submerged) y = np.ma.array(y, mask=submerged) # center of gravity: x_m, y_m = np.ma.mean(x, axis=0), np.ma.mean(y, axis=0) center = self.proj_latlon(x_m, y_m, inverse=True) one = np.ones_like(x) # variance: variance = np.ma.mean((x-x_m*one)**2 + (y-y_m*one)**2, axis=0) return center,variance
[docs] def reset(self): """Preparing OpenDrift object for new run""" if not hasattr(self, 'start_time'): logger.info('Nothing to reset') return for attr in ['start_time', 'history', 'elements']: if hasattr(self, attr): delattr(self, attr) #del self.start_time #del self.history #del self.elements self.elements_deactivated = self.ElementType() # Empty array self.elements = self.ElementType() # Empty array