Source code for opendrift.readers.interpolation.structured

import numpy as np
from scipy.ndimage import map_coordinates
import scipy.ndimage as ndimage
from scipy.interpolate import interp1d, LinearNDInterpolator

from .interpolators import Nearest2DInterpolator, fill_NaN_towards_seafloor, horizontal_interpolation_methods, vertical_interpolation_methods
from opendrift.readers.basereader import variables

import logging
logger = logging.getLogger(__name__)

[docs] class ReaderBlock(): """Class to store and interpolate the output from a reader with data on a regular (structured) grid.""" def __init__(self, data_dict, interpolation_horizontal='linearNDFast', interpolation_vertical='linear'): # Make pointers to data values, for convenience self.x = data_dict['x'] self.y = data_dict['y'] self.time = data_dict['time'] self.data_dict = data_dict del self.data_dict['x'] del self.data_dict['y'] del self.data_dict['time'] try: self.z = data_dict['z'] del self.data_dict['z'] except: self.z = None # Mask any extremely large values, e.g. if missing netCDF _Fill_value filled_variables = set() for var in self.data_dict: self.data_dict[var] = variables.Variables.__check_variable_array__(var, self.data_dict[var]) if isinstance(self.data_dict[var], np.ma.core.MaskedArray): self.data_dict[var] = np.ma.masked_outside( np.ma.masked_invalid(self.data_dict[var]), -1E+9, 1E+9) # Convert masked arrays to numpy arrays self.data_dict[var] = np.ma.filled(self.data_dict[var], fill_value=np.nan) # Fill missing data towards seafloor if 3D if isinstance(self.data_dict[var], (list,)): logger.warning('Ensemble data currently not extrapolated towards seafloor') elif self.data_dict[var].ndim == 3: filled = fill_NaN_towards_seafloor(self.data_dict[var]) if filled is True: filled_variables.add(var) if len(filled_variables) > 0: logger.debug('Filled NaN-values toward seafloor for :' + str(list(filled_variables))) # Set 1D (vertical) and 2D (horizontal) interpolators try: self.Interpolator2DClass = \ horizontal_interpolation_methods[interpolation_horizontal] except Exception: raise NotImplementedError( 'Valid interpolation methods are: ' + str(horizontal_interpolation_methods.keys())) try: self.Interpolator1DClass = \ vertical_interpolation_methods[interpolation_vertical] except Exception: raise NotImplementedError( 'Valid interpolation methods are: ' + str(vertical_interpolation_methods.keys())) if 'land_binary_mask' in self.data_dict.keys() and \ interpolation_horizontal != 'nearest': logger.debug('Nearest interpolation will be used ' 'for landmask, and %s for other variables' % interpolation_horizontal)
[docs] def _initialize_interpolator(self, x, y, z=None): logger.debug('Initialising interpolator.') self.interpolator2d = self.Interpolator2DClass(self.x, self.y, x, y) if self.z is not None and len(np.atleast_1d(self.z)) > 1: self.interpolator1d = self.Interpolator1DClass(self.z, z)
[docs] def interpolate(self, x, y, z=None, variables=None, profiles=[], profiles_depth=None): self._initialize_interpolator(x, y, z) env_dict = {} if profiles is not []: profiles_dict = {'z': self.z} for varname, data in self.data_dict.items(): nearest = False if varname == 'land_binary_mask': nearest = True self.interpolator2d_nearest = Nearest2DInterpolator(self.x, self.y, x, y) if type(data) is list: num_ensembles = len(data) logger.debug('Interpolating %i ensembles for %s' % (num_ensembles, varname)) if data[0].ndim == 2: horizontal = np.zeros(x.shape)*np.nan else: horizontal = np.zeros((len(self.z), len(x)))*np.nan ensemble_number = np.remainder(range(len(x)), num_ensembles) for en in range(num_ensembles): elnum = ensemble_number == en int_full = self._interpolate_horizontal_layers(data[en], nearest=nearest) if int_full.ndim == 1: horizontal[elnum] = int_full[elnum] else: horizontal[:, elnum] = int_full[:, elnum] else: horizontal = self._interpolate_horizontal_layers(data, nearest=nearest) if profiles is not None and varname in profiles: profiles_dict[varname] = horizontal if horizontal.ndim > 1: env_dict[varname] = self.interpolator1d(horizontal) else: env_dict[varname] = horizontal if 'z' in profiles_dict: profiles_dict['z'] = np.atleast_1d(profiles_dict['z']) return env_dict, profiles_dict
[docs] def _interpolate_horizontal_layers(self, data, nearest=False): '''Interpolate all layers of 3d (or 2d) array.''' if nearest is True: interpolator2d = self.interpolator2d_nearest else: interpolator2d = self.interpolator2d if data.ndim == 2: return interpolator2d(data) if data.ndim == 3: num_layers = data.shape[0] # Allocate output array result = np.ma.empty((num_layers, len(interpolator2d.x))) for layer in range(num_layers): result[layer, :] = self.interpolator2d(data[layer, :, :]) return result
[docs] def covers_positions(self, x, y, z=None): '''Check if given positions are covered by this reader block.''' indices = np.where((x >= self.x.min()) & (x <= self.x.max()) & (y >= self.y.min()) & (y <= self.y.max()))[0] if len(indices) == len(x): return True else: return False