From 3102764504870519ebac05bb41694c010275c386 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?H=C3=A5vard=20Espenes?= Date: Thu, 18 Dec 2025 11:38:23 +0100 Subject: [PATCH] Adds a FVCOM reader which builds on netCDF4. Reads the next timestep synchronously. In practice, this allows you to lower the model timestep without penalty in terms of simulation time, since you won't need to wait for the next timestep to load as you enter it. (with an appropriately low timestep, that is) --- opendrift/readers/basereader/finite_volume.py | 446 ++++++++++++++++++ opendrift/readers/reader_FVCOM_netcdf.py | 373 +++++++++++++++ pyproject.toml | 3 +- 3 files changed, 821 insertions(+), 1 deletion(-) create mode 100644 opendrift/readers/basereader/finite_volume.py create mode 100644 opendrift/readers/reader_FVCOM_netcdf.py 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]