diff --git a/opendrift/readers/basereader/finite_volume.py b/opendrift/readers/basereader/finite_volume.py new file mode 100644 index 000000000..e0b0ed0e5 --- /dev/null +++ b/opendrift/readers/basereader/finite_volume.py @@ -0,0 +1,446 @@ +from abc import abstractmethod +import numpy as np +import logging + +logger = logging.getLogger(__name__) + +from .variables import Variables +from opendrift.readers.interpolation.finite_volume import FiniteVolumeReaderBlock +from numba import njit, prange + +class FiniteVolumeReader(Variables): + """ + An unstructured ReaderBlock for Finite Volume ocean models with irregularly gridded data on an Arakawa B-grid + + The initial type of grid that this class supports are `triangular prisms + `_. + Unstructured in xy-coordinates, x and y is constant in z. z might be + non-cartesian (e.g. sigma-levels). + + .. seealso:: + + :py:mod:`opendrift.readers` + + :class:`.structured.StructuredReader` + """ + + # Number of parallel threads to use when possible + boundary = None + + # nodes + x, y = None, None + node_variables = None # list of std-name variables defined at nodes + + # faces + xc, yc = None, None + face_variables = None # list of std-name variables defined at center of faces + + @abstractmethod + def get_variables(self, variables, time=None, x=None, y=None, z=None): + """ + Obtain and return values of the requested variables at all positions + (x, y, z) closest to given time. + """ + + def _get_variables_interpolated_(self, variables, profiles, profiles_depth, time, x, y, z): + """ + The function called by OpenDrift to get an estimate of velocities at the particle at timestep n. + - makes use of thredding to access the next timestep + - treats static and dynamic variables differently + """ + logger.debug(f"Requested variabels: {variables}") + + x, y, z = np.atleast_1d(x), np.atleast_1d(y), np.atleast_1d(z) + + # Isn't this a bit assumptuous? + if len(z) == 1: + z = z[0] * np.ones(x.shape) + + requested_variables, time, x, y, z, _outside = self.check_arguments(variables, time, x, y, z) + + # Figure out where we are in this netCDF files timeline + # ---- + (nearest_time, time_before, time_after, + indx_nearest, _indx_before, _indx_after, + ) = self.nearest_time(time) + logger.debug("Reader time:\n\t\t%s (before)\n\t\t%s (after)" % (time_before, time_after)) + + + # If none of the requested variables are dynamic, we won't need the blocks. + # For time-dependent variables, however, we do. + # ---- + time_dependent = [var for var in requested_variables if var not in self.static_variables] + if not time_dependent and self.block_before is not None: + env_before, env_after = self._interpolate_from_block_to_particles( + x, y, z, time, time_before, time_after, variables + ) + env = env_before + else: + # Get data at timesteps before- and after time (these are added to self) + # --- + logger.debug(f"interpolating {requested_variables} to the particles") + self.timer_start("load_to_blocks") + self._get_var_blocks(variables, time_before, time_after, x, y, z) + self._check_var_blocks(variables, time_before, time_after, x, y, z) + self.timer_end("load_to_blocks") + + # Interpolate data at those timesteps to the particles + self.timer_start("interpolation_in_space") + env_before, env_after = self._interpolate_from_block_to_particles( + x, y, z, time, time_before, time_after, variables + ) + self.timer_end("interpolation_in_space") + + # Interpolation in time + self.timer_start("interpolation_in_time") + if ( + (time_after is not None) + and (time_before != time) + and self.always_valid is False + ): + weight_after = self._find_weight_after(time) + weight_before = 1-weight_after + logger.debug( + ("Interpolating before (%s, weight %.2f) and\n\t\t after (%s, weight %.2f) in time") + % (self.block_before.time, weight_before, self.block_after.time, weight_after) + ) + + # Interpolation procedure + env = {} + for var in variables: + env[var] = env_before[var]*weight_before + env_after[var]*weight_after + else: + logger.debug("No time interpolation needed - right on time.") + env = env_before + self.timer_end("interpolation_in_time") + env_profiles = None + return env, env_profiles + + def _find_weight_after(self, time): + '''Calculates weight for interpolation in time''' + weight_after = ((time - self.block_before.time).total_seconds() + / (self.block_after.time - self.block_before.time).total_seconds()) + + assert weight_after <= 1, f'{weight_after=} is invalid' + assert weight_after >= 0, f'{weight_after=} is invalid' + return weight_after + + # blocks of data for timesteps + # ---- + def _get_var_blocks(self, variables, time_before, time_after, x, y, z): + """ + Make block variables for timestep before and after + - variables -- to be interpolated from the ocean model + - time_before -- model timestep before particle timestep + - time_after -- model timestep after particle timestep + - x, y, z -- Currently does not do anything in particular, since profiles are not supported + and interpolation is done later. + + However get_variables (further down the pipeline) is expected to read x,y,z, so I just kept + it here to keep the call consistent with what OpenDrift expects upstream. + """ + # Timestep (assuming times are regularly spaced) + # ---- + dt = time_after-time_before + + # Swap block when moving over to the next model timestep, see if we can access the next block + # ---- + if self.block_before is not None and self.block_after is not None: + if self.block_before.time != time_before and self.block_after.time == time_before: + self.block_before = self.block_after + next_block = self.get_block_next(time_after) + if next_block is not None: + self.block_after = next_block + logger.debug(f"Fetched env-block for time before {self.block_before.time}") + self.make_block_next(variables, time_after + dt, x, y, z) + + if self.block_after.time != time_after and self.block_before.time == time_after: + self.block_after = self.block_before + next_block = self.get_block_next(time_before) + if next_block is not None: + self.block_before = next_block + logger.debug(f"Fetched env-block for time after {self.block_after.time}") + self.make_block_next(variables, time_before - dt, x, y, z) + + # First access to blocks - when no buffer is available + # ---- + if self.block_before is None or self.block_before.time != time_before: + _ = self.get_block_next(time_before) # make sure to close the next_block thread if active + self.block_before = self._fetch_data(variables, time_before, x, y, z) + logger.debug(f"Fetched env-block for time before {self.block_before.time}") + + if self.block_after is None or self.block_after.time != time_after: + _ = self.get_block_next(time_after) # make sure to close the next_block thread if active + self.block_after = self._fetch_data(variables, time_after, x, y, z) + logger.debug(f"Fetched env-block for time after {self.block_after.time}") + + def _check_var_blocks(self, variables, time_before, time_after, x, y, z): + """ + check that both blocks have the data they need to update the environent + """ + # Check that all requested variables are in the blocks, update both if something is missing + if self.block_before is not None: + missing = [var for var in variables if var not in self.block_before.data_dict.keys()] + updated_variables = [var for var in self.block_before.data_dict.keys() if var != 'time'] + missing + if missing: + logger.info(f'Missing {missing}, reloading block_before with them') + _ = self.get_block_next(time_before) + self.block_before = self._fetch_data(updated_variables, time_before, x, y, z) + + if self.block_after is not None: + missing = [var for var in variables if var not in self.block_after.data_dict.keys()] + updated_variables = [var for var in self.block_before.data_dict.keys() if var != 'time'] + missing + if missing: + logger.info(f'Missing {missing}, reloading block_after with them') + _ = self.get_block_next(time_after) + self.block_after = self._fetch_data(updated_variables, time_after, x, y, z) + + def _fetch_data(self, variables, time, x, y, z): + """ + Get another block of data if not already cached + """ + var_block = FiniteVolumeReaderBlock(self.get_variables(variables, time, x, y, z)) + + # Add some necessary variable info to the object + var_block.node_variables = self.node_variables + var_block.face_variables = self.face_variables + var_block.variable_mapping = self.variable_mapping + var_block.variable_aliases = self.variable_aliases + return var_block + + # Methods needed to use an asynchronous strategy to speed up the reading part + # ---- + def get_block_next(self, time): + ''' + time is the timestep we want to read. + Also used to ensure that the thread is reading data from a netCDF file before accessing more data in the netCDF + ''' + block_next = None + if self.block_next is not None: + block_next = self.block_next.get() + if block_next.time != time: + block_next = None + else: + logger.debug(f"Fetched env-block for next time {block_next.time}") + return block_next + + def make_block_next(self, variables, next_time, x, y, z): + ''' + Start reading the (expected) next block of data + ''' + if self.end_time > next_time: + logger.debug(f"Fetcing env-block for next timestep {next_time}") + _ = self.get_block_next(next_time) # just to make sure that all thredds are closed + self.block_next = self.pool.apply_async(self._fetch_data, args = (variables, next_time, x, y, z)) + else: + self.block_next = None + + # Interpolation + # ---- + def _interpolate_from_block_to_particles(self, x, y, z, time, time_before, time_after, variables): + """ + Interpolate data from block to positions of particles in space + - Warning: Does not interpolate. The routine finds nearest point in space to (x, y, z). + """ + node_variables = [var for var in variables if var in self.node_variables] + cell_variables = [var for var in variables if var in self.face_variables] + + # Find nearest nodes, cells and corresponding sigma layers + # ---- + cells = {} + if cell_variables: + cells["id"] = self._nearest_cell_(x, y) # 0.1 s for 1 millon particles, 2 million cells + if self.use_3d: + cells["sigma_ind"] = self.__nearest_cell_sigma__(cells["id"], z) # 0.4 s (same exp) + (cells["sigma_next"], cells['weight_sigma_ind'], cells['weight_sigma_next'], + ) = self._sigma_for_interpolation(self.ocean_depth_cells, self.siglay_center, cells, z) # 0.2 s (same exp) + + nodes = {} + if node_variables: + nodes["id"] = self._nearest_node_(x, y) + if self.use_3d: + nodes["sigma_ind"] = self.__nearest_node_sigma__(nodes["id"], z) + (nodes["sigma_next"], nodes['weight_sigma_ind'], nodes['weight_sigma_next'], + ) = self._sigma_for_interpolation(self.ocean_depth_nodes, self.siglay, nodes, z) + logger.debug(f"Interpolating before {self.block_before.time} in space {self.interpolation}") + + # Get environment at timestep before + env_before = self.block_before.interpolate(nodes, cells, variables) + + # If relevant, get environment at timestep after + if (time_after is not None) and (time_before != time): + logger.debug(f"Interpolating after {self.block_after.time} in space {self.interpolation}") + env_after = self.block_after.interpolate(nodes, cells, variables) + else: + env_after = env_before + return env_before, env_after + + def __nearest_node_sigma__(self, nodes, z): + """ + Find nearest depth at node. + """ + # For now, we only need siglay + sigmas = self.siglay[:, nodes] + + # Calculating depths from sigmas + depths = self.z_from_sigma(sigmas, self.ocean_depth_nodes[nodes]) + + return self._vector_nearest_(depths, z) + + def __nearest_cell_sigma__(self, cells, z): + """ + Find nearest depth at element. + """ + # add siglay if needed + sigmas = self.siglay_center[:, cells] + + # Calculating depths from sigmas + depths = self.z_from_sigma(sigmas, self.ocean_depth_cells[cells]) + + return self._vector_nearest_(depths, z) + + def _sigma_for_interpolation(self, depth, sigma, grid, z): + ''' + Find the sigma level which closes the "box" bounding nearest sigma and the particle + - nearest --> sigma nearest + - next --> sigma next to nearest that creates a closed box covering the particle + + Using the vector for full depth is not very efficient, and can be rewritten + - This routine has not been properly performance tested, I expect there to be room for rewrites... + ''' + # Get grid depth and depth of the nearest sigma layer + full_depths = self.z_from_sigma(sigma[:, grid['id']], depth[grid['id']]) + sigma_ind_depth = full_depths[grid['sigma_ind'], np.arange(grid['id'].shape[0])] # depth of nearest sigma level + + ## Determine if nearest sigma level is below or above the particle + dz = sigma_ind_depth-z + + # - Negative number if the particle is above nearest sigma, positive if particle is below. + # The other sigma level we wish to extract is therefore either above or below (+1 or -1) + dsig = self.sign(dz.data).astype(int) + sigmas_next = grid['sigma_ind'] + dsig + + # Find levels where the particle is outside of the "sigma range" + particles_over_shallowest_sigma = np.where(sigmas_next<0)[0] + particles_below_deepest_sigma = np.where(sigmas_next==full_depths.shape[0])[0] + + # We need these invalid sigma-layers to have valid indices when we broadcast, set to nearest sigma + sigmas_next[particles_below_deepest_sigma] = grid['sigma_ind'][particles_below_deepest_sigma] + + # We use a linear interpolation scheme in the vertical and set next-sigma + # weights to 0 where we just want to use the nearest sigma + weight_sigma_next = np.abs(sigma_ind_depth - z)/np.abs(sigma_ind_depth - full_depths[sigmas_next, np.arange(grid['id'].shape[0])]) + weight_sigma_next[particles_below_deepest_sigma] = 0 + weight_sigma_next[particles_over_shallowest_sigma] = 0 + + assert weight_sigma_next.min() >= 0, f'invalid vertical interpolation weight: {weight_sigma_next.min()}' + assert weight_sigma_next.max() <= 1, f'invalid vertical interpolation weight: {weight_sigma_next.max()}' + + weight_sigma_ind = 1-weight_sigma_next + + return sigmas_next, weight_sigma_ind, weight_sigma_next + + def _vector_nearest_(self, X, xp): + """ + Find nearest element in vector of vectors `X` for each `xp`. + + Args: + X NxM matrix of levels + xp M vector of positions + + Returns: + i M vector of indices [0..N] of closest element in + X[0..N, i] to xp[i] + """ + xp = np.atleast_2d(xp) + diff = self._abs(X.data - xp) + return self._find_sigma(diff) + + @staticmethod + @njit + def sign(field): + return np.sign(field) + + @staticmethod + @njit + def _abs(difference): + return np.abs(difference) + + @staticmethod + @njit(parallel = True) + def _find_sigma(depth_difference): + """ + Essentially a sped-up version of np.argmin + Benchmarked on 1 mill particles, 34 sigma layers: + np.argmin(depth_difference, axis=0): 330 ms + _find_sigma(depth_difference): 10 ms + + Developed and tested using numba-0.56.4 + """ + out = np.zeros((depth_difference.shape[1],), dtype=np.int32) + for n in prange(depth_difference.shape[1]): + out[n] = np.argmin(depth_difference[:, n]) + return out + + @staticmethod + def z_from_sigma(sigma, depth, elevation=0): + """ + Calculate z-depth from sigma constant. + + Args: + sigma Sigma coefficient(s) + depth Depth below mean sea-surface + elevation Elevation of sea-surface (e.g. tidal) + + Returns: z, depth below sea-surface in meters. + """ + return sigma * (depth + elevation) + + def _nearest_node_(self, x, y): + """ + Return nearest node (id) for x and y + """ + query_points = np.array([x, y], dtype=np.float32).T + _, inds = self.KDTree_node.query(query_points, k=1) + return inds + + def _nearest_cell_(self, x, y): + """ + Return nearest cell or face (id) for x and y + """ + query_points = np.array([x, y], dtype=np.float32).T + _, inds = self.KDTree_cell.query(query_points, k=1) + return inds + + # Rotines expected by OpenDrift + # ---- + def _build_boundary_polygon_(self, x_grid, y_grid): + """ + Build a mesh boundary polygon + """ + from shapely.geometry import Polygon + from shapely.prepared import prep + from scipy.spatial import ConvexHull + + P = np.vstack((x_grid, y_grid)).T + hull = ConvexHull(P) + + boundary = P[hull.vertices, :] + boundary = prep(Polygon(boundary)) + + return boundary + + def covers_positions(self, x, y, z=0): + """ + Check which points are within boundary of mesh. + """ + assert ( + self.boundary is not None + ), "Boundary of mesh has not been prepared by reader" + + # TODO: Check z coordinates + logger.warning("z-coordinates are not bounds-checked") + + from shapely.vectorized import contains + + return contains(self.boundary, x, y) diff --git a/opendrift/readers/reader_FVCOM_netcdf.py b/opendrift/readers/reader_FVCOM_netcdf.py new file mode 100644 index 000000000..d94486f49 --- /dev/null +++ b/opendrift/readers/reader_FVCOM_netcdf.py @@ -0,0 +1,373 @@ +# 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 . +# +# Copyright 2023, Håvard Espenes, University of Oslo / Akvaplan-niva +# Copyright 2020, Gaute Hope, MET Norway +# Copyright 2020, Simon Weppe, MetOcean Solution, MetService New Zealand +# Copyright 2015, Knut-Frode Dagestad, MET Norway + +from datetime import datetime, timedelta +from netCDF4 import Dataset, MFDataset +from functools import cached_property +from pyproj import Proj +from pykdtree.kdtree import KDTree +from multiprocessing.pool import ThreadPool + +import numpy as np +import logging + +logger = logging.getLogger(__name__) + +from opendrift.readers.basereader import BaseReader, FiniteVolumeReader + +# This reader use threads to speed up reading from netCDF files. This, however, means that the +# reader can not access the dataset while OpenDrift is advecting particles (which would not be thread safe) +class Reader(BaseReader, FiniteVolumeReader): + """ + A reader for unstructured (irregularily gridded) `CF compliant `_ netCDF files + from models using the finite volume method. Note, we have just tested it with hydrodynamics from FVCOM 4.x. + + Args: + :param filename: A single netCDF file, or a pattern of files. The netCDF file can also be an URL to an OPeNDAP server. + :type filename: string, required. + + :param name: Name of reader + :type name: string, optional + + :param proj4: PROJ.4 string describing projection of data. + :type proj4: string, optional + + .. seealso:: + + py:mod:`opendrift.readers.basereader.finite_volume`. + """ + + dataset = None + block_before, block_after, block_next = None, None, None + + # We aspire to include more complicated schemes in the future. + interpolation = "nearest" + + def __init__( + self, filename=None, filelist=None, name=None, proj4=None, use_3d=True + ): + """ + The reader is used to read data, either from a single netCDF4 or from a number of files + using the MFDataset class of netCDF4. + + - filename : use a specific file as input + - filelist : use a list of files as input + - name : give a name to the reader + - proj3 : projection of the readers grid coordinates + - use_3d : Switch + - if True (default): use (u,v,w) to advect particles + - if False: use (ua, va) to advect particles + """ + self.use_3d = use_3d + assert ( + filename is not None or filelist is not None + ), "You must either provide a filename or a filelist" + filestr = str(filename) + + name + if name is None: + if filename is not None: + self.name = filestr + elif filelist is not None: + self.name = filelist[0].split('/')[-1].split('_')[0] + else: + self.name = name + + self._open_fvcom_dataset(filestr, filename, filelist) + + self._load_fvcom_grid(proj4) + + self._check_and_prepare_time() + + self._configure_fields_to_read_from_fvcom() + + self._map_fvcom_variables_to_opendrift_internal_names() + + # Run constructor of parent Reader class + super().__init__() + + self.boundary = self._build_boundary_polygon_( + self.x.compressed(), self.y.compressed() + ) + + # Start a threadded pool for reading the next timestep + self.pool = ThreadPool(processes = 1) + + @cached_property + def KDTree_node(self): + return KDTree(np.array([self.x, self.y], dtype=np.float32).T) + + @cached_property + def KDTree_cell(self): + return KDTree(np.array([self.xc, self.yc], dtype=np.float32).T) + + # Get bounds of the model domain + # ---- + @property + def xmin(self): + return np.min(self.x) + + @property + def xmax(self): + return np.max(self.x) + + @property + def ymin(self): + return np.min(self.y) + + @property + def ymax(self): + return np.max(self.y) + + # Routines related to the initialization of the reader + # ---- + def _open_fvcom_dataset(self, filestr, filename, filelist): + """ + Opens fvcom file using netCDF dataset methods + """ + self.timer_start("open dataset") + + if ("*" in filestr) or ("?" in filestr) or ("[" in filestr): + logger.info(f"Opening multiple datasets") + self.dataset = MFDataset(filename) + + elif filelist is not None: + logger.info(f"Opening multiple datasets") + self.dataset = MFDataset(filelist, "r") + + else: + logger.info(f"Opening dataset: {filestr}") + self.dataset = Dataset(filename, "r") + self.timer_end("open dataset") + + def _load_fvcom_grid(self, proj4): + """ + Makes the reader aware of projection and grid positions in this projection + """ + logger.info("Reading grid and coordinate variables..") + self.proj4 = proj4 + if self.dataset.CoordinateSystem == 'GeoReferenced': + project = Proj(self.proj4) + self.x, self.y = project(self.dataset['lon'][:], self.dataset['lat'], inverse = True) + self.xc, self.yc = project(self.dataset['lonc'][:], self.dataset['latc'], inverse = True) + + elif self.dataset.CoordinateSystem == "Cartesian": + self.x, self.y = self.dataset['x'][:], self.dataset['y'][:] + self.xc, self.yc = self.dataset['xc'][:], self.dataset['yc'][:] + + # depth properties + self.siglay, self.siglev = self.dataset['siglay'][:], self.dataset['siglev'][:] + self.ocean_depth_nodes = self.dataset["h"][:] + self.ocean_depth_cells = self.dataset["h_center"][:] + self.siglay_center = self.dataset["siglay_center"][:] + self.siglev_center = self.dataset["siglev_center"][:] + + def _check_and_prepare_time(self): + """ + Checks that the time is in assumed format, defines the datum and updates the reader + with min- and max time for the current file. This is essential information for the blocks later on. + """ + assert self.dataset["time"].time_zone == "UTC" + assert self.dataset["time"].units == "days since 1858-11-17 00:00:00" + assert self.dataset["time"].format == "modified julian day (MJD)" + ref_time = datetime( + 1858, 11, 17, 00, 00, 00 + ) # TODO: include , tzinfo=timezone.utc) + self.times = np.array( + [ref_time + timedelta(days=Itime.item()) + timedelta(milliseconds=Itime2.item()) for (Itime, Itime2) in zip(self.dataset["Itime"][:], self.dataset['Itime2'][:])] + ) + self.start_time = self.times[0] + self.end_time = self.times[-1] + + def _configure_fields_to_read_from_fvcom(self): + """ + FVCOM stores fields with different long_names and variable names than done internally in + opendrift. Furthermore, we want to be able to advect using either 2d or 3d fields. + To do so, we need some dictionaries and lists that correctly references fields. + + Compatible with uk-fvcom branch (not sure how / if this differs from original FVCOM) + """ + # A list of time-independent variables + self.static_variables = [ + 'sea_floor_depth_below_sea_level', + 'sea_surface_height_above_geoid', + ] + + # And then make a list of all fields we might want to read from FVCOM + if self.use_3d: + logger.info("Load 3D (u, v, w) velocity data") + # name_used_in_fvcom : equivalent_CF_name used in opendrift + self.variable_aliases = { + "eastward_sea_water_velocity": "x_sea_water_velocity", + "Northward_sea_water_velocity": "y_sea_water_velocity", + "eastward wind": "x_wind", + "northward wind": "y_wind", + "Upward Water Velocity": "upward_sea_water_velocity", + "sea_floor_depth_below_geoid": "sea_floor_depth_below_sea_level", + "bed stress magnitude from currents": "bed_stress", + } + + self.node_variables = [ + "sea_water_salinity", + "sea_water_temperature", + "sea_floor_depth_below_sea_level", + "sea_surface_height_above_geoid", + ] + + self.face_variables = [ + "x_wind", + "y_wind", + "x_sea_water_velocity", + "y_sea_water_velocity", + "upward_sea_water_velocity", + ] + + assert ( + "u" in self.dataset.variables.keys() + ), f"u not found in {self.dataset.filepath()}" + + assert ( + "v" in self.dataset.variables.keys() + ), f"v not found in {self.dataset.filepath()}" + + else: + logger.info("Load 2D (ua, va) velocity data") + # [name_used_in_fvcom : equivalent_CF_name_for_opendrift] + self.variable_aliases = { + "Vertically Averaged x-velocity": "x_sea_water_velocity", + "Vertically Averaged y-velocity": "y_sea_water_velocity", + "sea_floor_depth_below_geoid": "sea_floor_depth_below_sea_level", + "bed stress magnitude from currents": "q", + } + + self.node_variables = [ + "sea_surface_height_above_geoid", + "sea_floor_depth_below_sea_level", + ] + + self.face_variables = ["x_sea_water_velocity", "y_sea_water_velocity"] + + assert ( + "ua" in self.dataset.variables.keys() + ), f"ua not found in {self.dataset.filepath()}" + + assert ( + "va" in self.dataset.variables.keys() + ), f"va not found in {self.dataset.filepath()}" + + def _map_fvcom_variables_to_opendrift_internal_names(self): + """ + since OpenDrift is a generic tool + """ + self.variable_mapping = {} + for var_name in self.dataset.variables: + # skipping coordinate variables + if var_name in [ + "x", + "y", + "time", + "lon", + "lat", + "lonc", + "latc", + "siglay", + "siglev", + "siglay_center", + "siglev_center", + ]: + continue + + var = self.dataset[var_name] + if "standard_name" in var.ncattrs(): + std_name = var.standard_name + std_name = self.variable_aliases.get(std_name, std_name) + self.variable_mapping[std_name] = str(var_name) + + elif "long_name" in var.ncattrs(): + long_name = var.long_name + long_name = self.variable_aliases.get(long_name, long_name) + self.variable_mapping[long_name] = str(var_name) + self.variables = list(self.variable_mapping.keys()) + + def get_variables(self, requested_variables, time=None, x=None, y=None, z=None): + """ + Routine used to get data from FVCOM results. FVCOM chunks by time, so it makes sense to extract all of it + and dump to a data cacher (ReaderBlock) + + FVCOM Grid: + - Uses triangular prisms + + Variables: + - Face (vector variables): + - Node (scalar variables) + """ + logger.debug("Requested variabels: %s" % (requested_variables)) + requested_variables, time, x, y, z, _outside = self.check_arguments( + requested_variables, time, x, y, z + ) + ( + nearest_time, + _time_before, + _time_after, + indx_nearest, + _indx_before, + _indx_after, + ) = self.nearest_time(time) + + logger.debug("Nearest time: %s" % nearest_time) + node_variables = [ + var for var in requested_variables if var in self.node_variables + ] + face_variables = [ + var for var in requested_variables if var in self.face_variables + ] + assert (len(node_variables) + len(face_variables)) == len( + requested_variables + ), "missing variables requested" + + # Let the routine using this dict know which time was read from FVCOM + # ---- + variables = {"time": nearest_time} + + # Load data from FVCOM + # ---- + if node_variables: + logger.debug("Reading node-variables..") + for var in node_variables: + dvar = self.variable_mapping.get(var) + dvar = self.dataset[dvar] + logger.debug("Reading from netCDF: %s (%s)" % (var, dvar)) + if var == "nodes": + variables[var] = np.empty([]) # Hvorfor prøver den å lese nodes? + elif dvar.ndim == 1: + variables[var] = dvar[:] + else: + variables[var] = dvar[indx_nearest, :] + + if face_variables: + logger.debug("Reading face-variables..") + for var in face_variables: + dvar = self.variable_mapping.get(var) + logger.debug("Reading from netCDF: %s (%s)" % (var, dvar)) + dvar = self.dataset[dvar] + if dvar.ndim == 1: + variables[var] = dvar[:] + else: + variables[var] = dvar[indx_nearest, :] + return variables diff --git a/pyproject.toml b/pyproject.toml index ca69a798b..47149c177 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -30,7 +30,8 @@ dependencies = [ "pykdtree>=1.3", "xhistogram>=0.3", "adios_db>1.2", - "copernicusmarine>=2.0" + "copernicusmarine>=2.0", + "numba>=0.63.1", ] [project.scripts]