diff --git a/.github/workflows/conda_test.yml b/.github/workflows/conda_test.yml index ec389ee..046ad9c 100644 --- a/.github/workflows/conda_test.yml +++ b/.github/workflows/conda_test.yml @@ -28,21 +28,21 @@ jobs: cd .. pytest --pyargs gridded - lint: - name: Flake8 linting - runs-on: "ubuntu-latest" - steps: - - uses: actions/checkout@v2 - - uses: conda-incubator/setup-miniconda@v2 - with: - activate-environment: lint - python-version: 3.9 - auto-activate-base: false - - shell: bash -l {0} - run: | - conda config --add channels conda-forge - - name: Lint - shell: bash -l {0} - run: | - conda install flake8 - python -m flake8 --statistics gridded/ + # lint: + # name: Flake8 linting + # runs-on: "ubuntu-latest" + # steps: + # - uses: actions/checkout@v2 + # - uses: conda-incubator/setup-miniconda@v2 + # with: + # activate-environment: lint + # python-version: 3.9 + # auto-activate-base: false + # - shell: bash -l {0} + # run: | + # conda config --add channels conda-forge + # - name: Lint + # shell: bash -l {0} + # run: | + # conda install flake8 + # python -m flake8 --statistics gridded/ diff --git a/conda_requirements.txt b/conda_requirements.txt index 0eddcc9..5b50588 100644 --- a/conda_requirements.txt +++ b/conda_requirements.txt @@ -3,8 +3,8 @@ # requires python, but if you leave it out, it will use python already in # the conda environment, rather than upgrading or downgrading -# python>=3.9.* -numpy>=1.15.* +python>=3.9 +numpy>=1.15 scipy -netcdf4>=1.4.* -cell_tree2d>=0.3.* +netcdf4>=1.4 +cell_tree2d>=0.3 diff --git a/conda_requirements_test.txt b/conda_requirements_test.txt index f5045de..926ab40 100644 --- a/conda_requirements_test.txt +++ b/conda_requirements_test.txt @@ -8,7 +8,7 @@ # conda install --file conda_requirements.txt --file conda_requirements_test.txt -# this is the minimal that will support travis tests +# this is the minimal that will support CI tests pytest pytest-cov progressbar diff --git a/docs/source/index.rst b/docs/source/index.rst index dcdf737..ff7f6c2 100644 --- a/docs/source/index.rst +++ b/docs/source/index.rst @@ -10,7 +10,17 @@ Welcome to ``gridded``'s documentation! :align: center -``gridded`` is a single API for accessing / working with gridded model results on multiple grid types +``gridded`` is a single Python API for accessing / working with gridded model results on multiple grid types. + +The source code for ``gridded`` can be found on gitHub at: + +https://github.com/NOAA-ORR-ERD/gridded + +It is also available as a conda package on conda-forge: + +``conda install -c conda-forge gridded`` + +It is not currently available on PyPi, due to complex compiled dependences -- help accepted. Contents diff --git a/gridded/__init__.py b/gridded/__init__.py index 42a73d3..06a626f 100644 --- a/gridded/__init__.py +++ b/gridded/__init__.py @@ -1,6 +1,6 @@ #!/usr/bin/env python -__version__ = "0.5.9" +__version__ = "0.6.3" VALID_SGRID_LOCATIONS = (None, 'center', 'edge1', 'edge2', 'node') VALID_UGRID_LOCATIONS = (None, 'node', 'face', 'edge', 'boundary') @@ -11,10 +11,13 @@ from .variable import Variable, VectorVariable from .time import Time -from .depth import DepthBase +from .depth import DepthBase, L_Depth, S_Depth DepthBase._default_component_types['variable'] = Variable +S_Depth._default_component_types['variable'] = Variable +S_Depth._default_component_types['bathymetry'] = Variable +S_Depth._default_component_types['zeta'] = Variable - +from .time import OutOfTimeRangeError __all__ = ["Variable", "VectorVariable", @@ -23,5 +26,7 @@ "Grid_S", "Grid_U", "Dataset", - "Time"] + "Time", + "OutOfTimeRangeError", + ] diff --git a/gridded/depth.py b/gridded/depth.py index fc8a1ad..1e2dc01 100644 --- a/gridded/depth.py +++ b/gridded/depth.py @@ -1,5 +1,6 @@ import warnings +import os import numpy as np from gridded.time import Time @@ -8,27 +9,44 @@ search_dataset_for_any_long_name, search_dataset_for_variables_by_longname, search_dataset_for_variables_by_varname, - merge_var_search_dicts + merge_var_search_dicts, + search_netcdf_vars, + can_create_class, + parse_filename_dataset_args ) class DepthBase(object): + _def_count = 0 _default_component_types = {'time': Time, 'grid': Grid, 'variable': None} + #'variable' is None here to avoid import issues. It is set in the __init__.py def __init__(self, + name=None, surface_index=None, bottom_index=None, + default_surface_boundary_condition='extrapolate', + default_bottom_boundary_conditon='mask', **kwargs): ''' :param surface_index: array index of 'highest' level (closest to sea level) :param bottom_index: array index of 'lowest' level (closest to seafloor) - :param positive_down: orientation of points coordinates - (for interpolation functions) ''' + self.name=name self.surface_index = surface_index self.bottom_index = bottom_index + self.default_surface_boundary_condition = default_surface_boundary_condition + self.default_bottom_boundary_condition = default_bottom_boundary_conditon + + @classmethod + def _can_create_from_netCDF(cls, + filename=None, + dataset=None, + grid_file=None, + data_file=None,): + return True @classmethod def from_netCDF(cls, @@ -116,37 +134,41 @@ def _gen_varname(cls, class L_Depth(DepthBase): - _def_count=0 - default_terms = [('depth_levels', 'depth')] - default_names = {'depth': ['depth','depth_levels']} - cf_names = {'depth': 'depth'} + + default_terms = [('depth_levels')] + default_names = {'depth_levels': ['depth', 'depth_levels']} + cf_names = {'depth_levels': 'depth'} def __init__(self, - name=None, terms=None, - surface_index=None, - bottom_index=None, - default_surface_boundary_condition='extrapolate', - default_bottom_boundary_conditon='mask', **kwargs): super(L_Depth, self).__init__(**kwargs) - self.name=name - self.terms={} + self.terms = {} if terms is None: raise ValueError('Must provide terms for level depth coordinate') else: - self.terms=terms + self.terms = terms for k, v in terms.items(): setattr(self, k, v) - self.surface_index = surface_index - self.bottom_index = bottom_index - self.default_surface_boundary_condition = default_surface_boundary_condition - self.default_bottom_boundary_condition = default_bottom_boundary_conditon @classmethod + def _can_create_from_netCDF(cls, + filename=None, + dataset=None, + grid_file=None, + data_file=None,): + ds, dg = parse_filename_dataset_args(filename=filename, + dataset=dataset, + grid_file=grid_file, + data_file=data_file) + + return can_create_class(cls, ds, dg) + @classmethod def from_netCDF(cls, filename=None, dataset=None, + grid_file=None, + data_file=None, name=None, topology=None, terms=None, @@ -154,21 +176,17 @@ def from_netCDF(cls, bottom_index=None, **kwargs ): - df = dataset if filename is None else get_dataset(filename, dataset) - - if df is None: - raise ValueError('No filename or dataset provided') + df, dg = parse_filename_dataset_args(filename=filename, + dataset=dataset, + grid_file=grid_file, + data_file=data_file) + nc_vars = search_netcdf_vars(cls, df, dg) if name is None: name = cls.__name__ + str(cls._def_count) if terms is None: - terms={} - for tn, tln in cls.default_terms: - vname=tn - #if tn not in dataset.variables.keys(): # 2023 - if tn not in df.variables.keys(): - vname = cls._gen_varname(filename, dataset, [tn], [tln]) - #terms[tn] = dataset[vname][:] # 2023 - terms[tn] = df[vname][:] + terms = {} + for term, tvar in nc_vars.items(): + terms[term] = tvar[:] if surface_index is None: surface_index = np.argmin(terms['depth_levels']) if bottom_index is None: @@ -182,7 +200,6 @@ def from_netCDF(cls, surface_index=surface_index, bottom_index=bottom_index, **kwargs) - def interpolation_alphas(self, points, @@ -270,31 +287,16 @@ def _find_required_depth_attrs(cls, filename, dataset=None, topology=None): return init_args, gt - class S_Depth(DepthBase): ''' Represents the ocean s-coordinates as implemented by ROMS, It may or may not be useful for other systems. ''' - _def_count = 0 - default_names = {'Cs_r': ['Cs_r'], - 'Cs_w': ['Cs_w'], - 's_rho': ['s_rho'], - 's_w': ['s_w'], - 'hc': ['hc'], - 'bathymetry': ['h'], - 'zeta': ['zeta'] - } - - # NOTE: I don't think these are CFnames :-( - cf_names = {'Cs_r': ['S-coordinate stretching curves at RHO-points'], - 'Cs_w': ['S-coordinate stretching curves at W-points'], - 's_rho': ['S-coordinate at RHO-points'], - 's_w': ['S-coordinate at W-points'], - 'hc': ['S-coordinate parameter, critical depth'], - 'bathymetry': ['bathymetry at RHO-points', 'Final bathymetry at RHO-points'], - 'zeta': ['free-surface'] - } + _default_component_types = {'time': Time, + 'grid': Grid, + 'variable': None, + 'bathymetry': None, + 'zeta': None,} def __init__(self, name=None, @@ -304,7 +306,6 @@ def __init__(self, zeta=None, terms=None, vtransform=2, - positive_down=True, surface_boundary_condition='extrapolate', bottom_boundary_condition='mask', **kwargs): @@ -331,11 +332,6 @@ def __init__(self, :param vtransform: S-coordinate transform type. 1 = Old, 2 = New :type vtransform: int (default 2) - :param positive_down: Flag for interpreting depth values as positive or negative - This applies to any coordinates passed into the various functions, NOT the values - of the S-Coordinates that this object represents - :type positive_down: boolean - :param surface_boundary_condition: Determines how to handle points above the surface :type surface_boundary_condition: string ('extrapolate' or 'mask') @@ -347,11 +343,10 @@ def __init__(self, self.name = name self.time = time self.grid = grid - self.bathymetry = bathymetry + self.bathymetry = bathymetry # Nodal bathymetry only. self.zeta = zeta self.terms = {} self.vtransform = vtransform - self.positive_down = positive_down if terms is None: raise ValueError('Must provide terms for sigma coordinate') else: @@ -359,7 +354,7 @@ def __init__(self, for k, v in terms.items(): setattr(self, k, v) if self.surface_index is None: - self.surface_index = self.num_w_levels - 1 + self.surface_index = self.num_levels - 1 if self.bottom_index is None: self.bottom_index = 0 @@ -374,9 +369,7 @@ def from_netCDF(cls, varnames=None, grid_topology=None, name=None, - # units=None, time=None, - # time_origin=None, grid=None, dataset=None, data_file=None, @@ -385,9 +378,7 @@ def from_netCDF(cls, bathymetry=None, zeta=None, terms=None, - # fill_value=0, vtransform=2, - positive_down=True, **kwargs ): ''' @@ -419,75 +410,61 @@ def from_netCDF(cls, Grid = cls._default_component_types['grid'] Time = cls._default_component_types['time'] Variable = cls._default_component_types['variable'] - if filename is not None: - data_file = filename - grid_file = filename - - ds = None - dg = None - if dataset is None: - if grid_file == data_file: - ds = dg = get_dataset(grid_file) - else: - ds = get_dataset(data_file) - dg = get_dataset(grid_file) - else: - if grid_file is not None: - dg = get_dataset(grid_file) - else: - dg = dataset - ds = dataset + Bathymetry = cls._default_component_types['bathymetry'] + Zeta = cls._default_component_types['zeta'] + + ds, dg = parse_filename_dataset_args(filename=filename, + dataset=dataset, + grid_file=grid_file, + data_file=data_file) + if data_file is None: + data_file = os.path.split(ds.filepath())[-1] if grid is None: - grid = Grid.from_netCDF(grid_file, - dataset=dg, + grid = Grid.from_netCDF(dataset=dg, grid_topology=grid_topology) if name is None: name = cls.__name__ + str(cls._def_count) cls._def_count += 1 - varnames = cls.default_names.copy() # Do a comprehensive search for netCDF4 Variables all at once - vn_search = search_dataset_for_variables_by_varname(ds, varnames) - ds_search = search_dataset_for_variables_by_longname(ds, cls.cf_names) - varnames = merge_var_search_dicts(ds_search, vn_search) - if ds != dg: - dg_search = search_dataset_for_variables_by_longname(dg, cls.cf_names) - varnames = merge_var_search_dicts(varnames, dg_search) + nc_vars = search_netcdf_vars(cls, ds, dg) if bathymetry is None: - bathy_var = varnames.get('bathymetry', None) + bathy_var = nc_vars.get('bathymetry', None) if bathy_var is None: err = 'Unable to locate bathymetry in data file' if dg: raise ValueError(err + ' or grid file') raise ValueError(err) - bathymetry = Variable.from_netCDF(dataset=bathy_var._grp, + bathymetry = Bathymetry.from_netCDF(dataset=bathy_var._grp, varname=bathy_var.name, grid=grid, + time=time, name='bathymetry', ) if zeta is None: - zeta_var = varnames.get('zeta', None) + zeta_var = nc_vars.get('zeta', None) if zeta_var is None: warn = 'Unable to locate zeta in data file' if dg: warnings.warn(warn + ' or grid file.') warn += ' Generating constant (0) zeta.' warnings.warn(err) - zeta = Variable.constant(0) + zeta = Zeta.constant(0) else: - zeta = Variable.from_netCDF(dataset=zeta_var._grp, + zeta = Zeta.from_netCDF(dataset=zeta_var._grp, varname=zeta_var.name, grid=grid, + time=time, name='zeta') if time is None: time = zeta.time if terms is None: terms = {} - for term, tvar in varnames.items(): + for term, tvar in nc_vars.items(): if term in ['bathymetry', 'zeta']: # skip these because they're done separately... continue @@ -503,52 +480,32 @@ def from_netCDF(cls, zeta=zeta, terms=terms, vtransform=vtransform, - positive_down=positive_down, **kwargs) - @property - def num_w_levels(self): - return len(self.s_w) + def num_levels(self): + raise NotImplementedError('num_levels not implemented for S_Depth, required in subclasses') @property - def num_r_levels(self): - return len(self.s_rho) + def num_layers(self): + raise NotImplementedError('num_layers not implemented for S_Depth, required in subclasses') - def __len__(self): - return self.num_w_levels + @classmethod + def _can_create_from_netCDF(cls, + filename=None, + dataset=None, + grid_file=None, + data_file=None,): + ds, dg = parse_filename_dataset_args(filename=filename, + dataset=dataset, + grid_file=grid_file, + data_file=data_file) + + return can_create_class(cls, ds, dg) - def compute_coordinates(self, rho_or_w): - """ - Uses zeta, bathymetry, and the s-coordinate terms to produce the - time-dependent coordinate variable. + def __len__(self): + return self.num_levels - :return: numpy array of s-coordinates - """ - if rho_or_w == 'rho': - s_c = self.s_rho - Cs = self.Cs_r - elif rho_or_w == 'w': - s_c = self.s_w - Cs = self.Cs_w - else: - raise ValueError('invalid rho_or_w argument (must be "rho" or "w")') - if len(self.zeta.data.shape) == 2: #need this check to workaround a 'constant' zeta object - zeta = np.zeros((len(self.time.data), ) + self.bathymetry.data.shape) - else: - zeta = self.zeta.data[:] - h = self.bathymetry.data[:] - hc = self.hc - hCs = h * Cs[:,np.newaxis, np.newaxis] - s_coord = None - if self.vtransform == 1: - S = (hc * s_c)[:, np.newaxis, np.newaxis] + hCs - (hc * Cs)[:, np.newaxis, np.newaxis] - s_coord = -(S + zeta[:, np.newaxis, :, :] * (1 + S / h)) - elif self.vtransform == 2: - S = ((hc * s_c)[:, np.newaxis, np.newaxis] + hCs) / (hc + h) - s_coord = -(zeta[:, np.newaxis, :, :] + (zeta[:, np.newaxis, :, :] + h) * S) - return s_coord - - def get_transect(self, points, time, rho_or_w='w', _hash=None, **kwargs): + def get_transect(self, points, time, data_shape, _hash=None, **kwargs): ''' :param points: array of points to interpolate to :type points: numpy array of shape (n, 3) @@ -556,65 +513,28 @@ def get_transect(self, points, time, rho_or_w='w', _hash=None, **kwargs): :param time: time to interpolate to :type time: datetime.datetime - :param rho_or_w: 'rho' or 'w' to interpolate to rho or w points - :type rho_or_w: string + :param data_shape: Shape of the variable to be interpolated. The first dimension is expected to be depth + :type rho_or_w: tuple of int - :return: numpy array of shape (n, num_w_levels) of interpolated values + :return: numpy array of shape (n, data_shape[0]) of n depth level transects ''' - - s_c = self.s_rho if rho_or_w == 'rho' else self.s_w - C_s = self.Cs_r if rho_or_w == 'rho' else self.Cs_w - h = self.bathymetry.at(points, time, unmask=False, _hash=_hash, **kwargs) - zeta = self.zeta.at(points, time, unmask=False, _hash=_hash, **kwargs) - hc = self.hc - hCs = h * C_s[np.newaxis, :] - if self.vtransform == 1: - S = (hc* s_c) + hCs - (hc * C_s)[np.newaxis, :] - s_coord = -(S + zeta * (1 + S / h)) - elif self.vtransform == 2: - S = ((hc * s_c) + hCs) / (hc + h) - s_coord = -(zeta + (zeta + h) * S) - return s_coord + raise NotImplementedError('get_transect not implemented for S_Depth, required in subclasses') - def _L_Depth_given_bathymetry_t1(self, bathy, zeta, lvl, rho_or_w): - """ - computes the depth (positive up) of a level given the bathymetry - and zeta. Produces an output array same shape as bathy and zeta - Output is always 'positive down' - This implements transform 1 (https://www.myroms.org/wiki/Vertical_S-coordinate) - """ - if rho_or_w == 'rho': - s = self.s_rho[lvl] - Cs = self.Cs_r[lvl] - elif rho_or_w == 'w': - s = self.s_w[lvl] - Cs = self.Cs_w[lvl] - else: - raise ValueError('invalid rho_or_w argument (must be "rho" or "w")') + def get_surface_depth(self, points, time, data_shape, _hash=None, **kwargs): + ''' + :param points: array of points to interpolate to + :type points: numpy array of shape (n, 3) - hc = self.hc - S = hc * s + (bathy - hc) * Cs - return -(S + zeta * (1 + S / bathy)) + :param time: time to interpolate to + :type time: datetime.datetime - def _L_Depth_given_bathymetry_t2(self, bathy, zeta, lvl, rho_or_w): - """ - computes the depth (positive up) of a level given the bathymetry - and zeta. Produces an output array same shape as bathy and zeta - Output is always 'positive down' - This implements transform 2 (https://www.myroms.org/wiki/Vertical_S-coordinate) - """ - if rho_or_w == 'rho': - s = self.s_rho[lvl] - Cs = self.Cs_r[lvl] - elif rho_or_w == 'w': - s = self.s_w[lvl] - Cs = self.Cs_w[lvl] - else: - raise ValueError('invalid rho_or_w argument (must be "rho" or "w")') + :param data_shape: Shape of the variable to be interpolated. The first dimension is expected to be depth + :type rho_or_w: tuple of int + + :return: numpy array of shape (n, 1) of n surface level values + ''' + raise NotImplementedError('get_surface_depth not implemented for S_Depth, required in subclasses') - hc = self.hc - S = (hc * s + bathy * Cs) / (hc + bathy) - return -(zeta + (zeta + bathy) * S) def interpolation_alphas(self, points, @@ -638,39 +558,20 @@ def interpolation_alphas(self, 'extrapolate' will set the index to the surface or bottom index, and the alpha to 0 or 1 depending on the orientation of the layers ''' - depths = points[:, 2] - surface_boundary_condition = (self.default_surface_boundary_condition - if surface_boundary_condition is None else - surface_boundary_condition) - bottom_boundary_condition = (self.default_bottom_boundary_condition - if bottom_boundary_condition is None else - bottom_boundary_condition) - - rho_or_w = surface_index = None # parameter necessary for ldgb - if data_shape[0] == self.num_w_levels: - # data assumed to be on w points - rho_or_w = 'w' + depths = points[:,2] + surface_boundary_condition = self.default_surface_boundary_condition if surface_boundary_condition is None else surface_boundary_condition + bottom_boundary_condition = self.default_bottom_boundary_condition if bottom_boundary_condition is None else bottom_boundary_condition + + surface_index = None + if data_shape[0] == self.num_levels: surface_index = self.surface_index - elif data_shape[0] == self.num_r_levels: - # data assumed to be on rho points - rho_or_w = 'rho' + elif data_shape[0] == self.num_layers: surface_index = self.surface_index-1 else: raise ValueError('Cannot get depth interpolation alphas for data shape specified; does not fit r or w depth axis') - # ldbg = 'level depth given bathymetry'. This is the depth of the layer - # at each x/y point. - if self.vtransform == 2: - ldgb = self._L_Depth_given_bathymetry_t2 - elif self.vtransform == 1: - ldgb = self._L_Depth_given_bathymetry_t1 - else: - raise ValueError('invalid vtransform attribute on depth object') - zetas = self.zeta.at(points, time, _hash=_hash, unmask=False, extrapolate=extrapolate).reshape(-1) - bathy = self.bathymetry.at(points, time, _hash=_hash, unmask=False, - extrapolate=extrapolate).reshape(-1) - - surface_depth = ldgb(bathy, zetas, surface_index, rho_or_w) + transects = self.get_transect(points, time, data_shape=data_shape, _hash=_hash, extrapolate=extrapolate) + surface_depth = transects[:, surface_index] indices = np.ma.MaskedArray(data=-np.ones((len(points)), dtype=np.int64), mask=np.zeros((len(points)), dtype=bool)) alphas = np.ma.MaskedArray(data=np.empty((len(points)), dtype=np.float64) * np.nan, mask=np.zeros((len(points)), dtype=bool)) @@ -686,8 +587,7 @@ def interpolation_alphas(self, alphas = np.ma.masked_invalid(alphas, 0) return indices, alphas - # use np.digitize to bin the depths into the layers - transects = self.get_transect(points, time, rho_or_w=rho_or_w, _hash=_hash, extrapolate=extrapolate) + #use np.digitize to bin the depths into the layers vf = np.vectorize(np.digitize, signature='(),(n)->()', excluded=['right']) indices = vf(depths, transects, right=True) - 1 @@ -716,69 +616,153 @@ def interpolation_alphas(self, raise ValueError('Some alphas are still unmasked and NaN. Please file a bug report') return indices, alphas - def get_section(self, time, coord='w', x_coord=None, y_coord=None, vtransform=2): + +class ROMS_Depth(S_Depth): + ''' + Sigma coordinate depth object for ROMS style output + ''' + _def_count = 0 + default_names = {'Cs_r': ['Cs_r'], + 'Cs_w': ['Cs_w'], + 's_rho': ['s_rho'], + 's_w': ['s_w'], + 'hc': ['hc'], + 'bathymetry': ['h'], + 'zeta': ['zeta'] + } + + cf_names = {'Cs_r': ['S-coordinate stretching curves at RHO-points'], + 'Cs_w': ['S-coordinate stretching curves at W-points'], + 's_rho': ['S-coordinate at RHO-points'], + 's_w':['S-coordinate at W-points'], + 'hc': ['S-coordinate parameter, critical depth'], + 'bathymetry': ['bathymetry at RHO-points', 'Final bathymetry at RHO-points'], + 'zeta': ['free-surface'] + } + + @property + def num_levels(self): + return len(self.s_w) + + @property + def num_layers(self): + return len(self.s_rho) + + def get_transect(self, points, time, data_shape=None, _hash=None, **kwargs): ''' - Returns a section of the z-level space in time. All s-levels are - returned in the data. By providing a x_coord or y_coord you can - get cross sections in the direction specified, or both for the level - depths at a single point. + :param points: array of points to interpolate to + :type points: numpy array of shape (n, 3) + + :param time: time to interpolate to + :type time: datetime.datetime + + :param data_shape: shape of the variable to be interpolated. This param is used to determine + whether to index on the sigma layers or levels + :type data_shape: tuple of int + + :return: numpy array of shape (n, num_w_levels) of n s-coordinate transects ''' - if coord not in ['w', 'rho']: - raise ValueError('Can only specify "w" or "rho" for coord kwarg') - ldgb = None - z_shp = None - if vtransform == 2: - ldgb = self._L_Depth_given_bathymetry_t2 - else: - ldgb = self._L_Depth_given_bathymetry_t1 - - if coord == 'w': - z_shp = (self.num_w_levels, self.zeta.data.shape[-2], self.zeta.data.shape[-1]) - if coord == 'rho': - z_shp = (self.num_r_levels, self.zeta.data.shape[-2], self.zeta.data.shape[-1]) - time_idx = self.time.index_of(time) - time_alpha = self.time.interp_alpha(time) - z0 = self.zeta.data[time_idx] - if time_idx == len(self.time.data) - 1 or time_idx == 0: - zeta = self.zeta.data[time_idx] - else: - z1 = self.zeta.data[time_idx + 1] - zeta = (z1 - z0)*(1-time_alpha) + z0 - bathy = self.bathymetry.data[:] - z_data=np.empty(z_shp) - for lvl in range(0,len(z_data)): - z_data[lvl] = ldgb(bathy, zeta, lvl, coord) - return z_data + s_c = self.s_rho if data_shape[0] == self.num_layers else self.s_w + C_s = self.Cs_r if data_shape[0] == self.num_layers else self.Cs_w + h = self.bathymetry.at(points, time, unmask=False, _hash=_hash, **kwargs) + zeta = self.zeta.at(points, time, unmask=False, _hash=_hash, **kwargs) + hc = self.hc + hCs = h * C_s[np.newaxis, :] + if self.vtransform == 1: + S = (hc* s_c) + hCs - (hc * C_s)[np.newaxis, :] + s_coord = -(S + zeta * (1 + S / h)) + elif self.vtransform == 2: + S = ((hc * s_c) + hCs) / (hc + h) + s_coord = -(zeta + (zeta + h) * S) + return s_coord + + +class FVCOM_Depth(S_Depth): + _def_count = 0 + default_names = { + 'siglay': ['siglay'], # mid layer depth coordinate on nodes + 'siglay_center': ['siglev_center'], # mid layer depth coordinate on centers + 'siglev': ['siglev'], # layer depth coordinate on nodes + 'siglev_center': ['siglev_center'], # layer depth coordinate on centers + 'bathymetry': ['h'], # bathymetry on nodes + 'h_center': ['h_center'], # bathymetry on centers + 'zeta': ['zeta'], # free surface + } + + cf_names = { + 'siglay': ['ocean_sigma/general_coordinate'], + 'siglay_center': ['ocean_sigma/general_coordinate'], + 'siglev': ['ocean_sigma/general_coordinate'], + 'siglev_center': ['ocean_sigma/general_coordinate'], + 'bathymetry': ['sea_floor_depth_below_geoid'], + 'h_center': ['sea_floor_depth_below_geoid'], + 'zeta': ['sea_surface_height_above_geoid'], + } + + @property + def num_levels(self): + return len(self.siglev) + + @property + def num_layers(self): + return len(self.siglay) + + def get_transect(self, points, time, data_shape=None, _hash=None, **kwargs): + ''' + :param points: array of points to interpolate to + :type points: numpy array of shape (n, 3) + + :param time: time to interpolate to + :type time: datetime.datetime + + :param rho_or_w: 'rho' or 'w' to interpolate to rho or w points + :type rho_or_w: string + + :return: numpy array of shape (n, num_w_levels) of n s-coordinate transects + ''' + + #because FVCOM sigma is defined for every node separately. + sigma = self.grid.interpolate_var_to_points(points[:, 0:2], self.siglev[:].T, location='node') + + bathy = self.bathymetry.at(points, time, unmask=False, _hash=_hash, **kwargs) + zeta = self.zeta.at(points, time, unmask=False, _hash=_hash, **kwargs) + + transects = -(zeta + (zeta + bathy) * sigma) + return transects + + # elif self.vtransform == 2: + # S = ((hc * s_c) + hCs) / (hc + h) + # s_coord = -(zeta + (zeta + h) * S) + # if no stretching or crit depth (hc, Cs_r, Cs_w) then S = s_c class Depth(object): ''' Factory class that generates depth objects. Also handles common loading and parsing operations ''' - ld_names = ['level', 'levels', 'L_Depth', 'depth_levels' 'depth_level'] - sd_names = ['sigma'] - surf_names = ['surface', 'surface only', 'surf', 'none'] + ld_types = [L_Depth] + sd_types = [ROMS_Depth, FVCOM_Depth] + surf_types = [DepthBase] def __init__(self): raise NotImplementedError("Depth is not meant to be instantiated. " "Please use the 'from_netCDF' or 'surface_only' function") @staticmethod - def surface_only(*args, **kwargs): + def surface_only(surface_index=-1, + **kwargs): ''' If instantiated directly, this will always return a DepthBase It is assumed index -1 is the surface index ''' - return DepthBase(*args,**kwargs) + return DepthBase(surface_index=surface_index, + **kwargs) - @staticmethod - def from_netCDF(filename=None, + @classmethod + def from_netCDF(cls, + filename=None, dataset=None, data_file=None, grid_file=None, - depth_type=None, - varname=None, - topology=None, - _default_types=(('level', L_Depth), ('sigma', S_Depth), ('surface', DepthBase)), **kwargs): ''' :param filename: File containing a depth @@ -791,86 +775,21 @@ def from_netCDF(filename=None, :type depth_type: string :returns: Instance of L_Depth or S_Depth ''' - if filename is not None: - data_file = filename - grid_file = filename - - ds = None - dg = None - if dataset is None: - if grid_file == data_file: - ds = dg = get_dataset(grid_file) - else: - ds = get_dataset(data_file) - dg = get_dataset(grid_file) - else: - if grid_file is not None: - dg = get_dataset(grid_file) - else: - dg = dataset - ds = dataset - - df = dataset if filename is None else get_dataset(filename, dataset) - if df is None: - raise ValueError('No filename or dataset provided') - - cls = depth_type - if (depth_type is None or isinstance(depth_type, str) or - not issubclass(depth_type, DepthBase)): - cls = Depth._get_depth_type(df, depth_type, topology, _default_types) - - return cls.from_netCDF(filename=filename, - dataset=dataset, - grid_file=grid_file, - data_file=data_file, - topology=topology, - **kwargs) - - # @staticmethod - # def depth_type_of_var(dataset, varname): - # ''' - # Given a varname or netCDF variable and dataset, try to determine the - # depth type from dimensions, attributes, or other metadata information, - # and return a grid_type and topology - # ''' - # var = dataset[varname] - - @staticmethod - def _get_depth_type(dataset, depth_type=None, topology=None, _default_types=None): - if _default_types is None: - _default_types = dict() + ds, dg = parse_filename_dataset_args(filename=filename, + dataset=dataset, + data_file=data_file, + grid_file=grid_file) + + typs = cls.sd_types + cls.ld_types + available_to_create = [typ._can_create_from_netCDF(grid_file=dg, data_file=ds) for typ in typs] + if not any(available_to_create): + warnings.warn('''Unable to automatically determine depth system so + reverting to surface-only mode. Manually check the + (depth_object).surface_index attribute and set it + to the appropriate array index for your model data''', RuntimeWarning) + return cls.surf_types[0].from_netCDF(data_file=ds, grid_file=dg, **kwargs) else: - _default_types = dict(_default_types) - - S_Depth = _default_types.get('sigma', None) - L_Depth = _default_types.get('level', None) - Surface_Depth = _default_types.get('surface', None) - if depth_type is not None: - if depth_type.lower() in Depth.ld_names: - return L_Depth - elif depth_type.lower() in Depth.sd_names: - return S_Depth - elif depth_type.lower() in Depth.surf_names: - return Surface_Depth - else: - raise ValueError('Specified depth_type not recognized/supported') - if topology is not None: - if ('faces' in topology.keys() - or topology.get('depth_type', 'notype').lower() in ld_names): - return L_Depth - elif topology.get('depth_type', 'notype').lower() in Depth.sd_names: - return S_Depth - else: - return Surface_Depth - - else: - t1 = search_dataset_for_any_long_name(dataset, S_Depth.cf_names) - if t1: - return S_Depth - else: - try: - L_Depth.from_netCDF(dataset=dataset) - return L_Depth - except: - warnings.warn('''Unable to automatically determine depth system so reverting to surface-only mode. Please verify the index of the surface is correct.''', RuntimeWarning) - return Surface_Depth + typ = typs[np.argmax(available_to_create)] + if sum(available_to_create) > 1: + warnings.warn('''Multiple depth systems detected. Using the first one found: {0}'''.format(typ.__repr__), RuntimeWarning) + return typ.from_netCDF(data_file=ds, grid_file=dg, **kwargs) diff --git a/gridded/grids.py b/gridded/grids.py index 427fee7..efc1a75 100644 --- a/gridded/grids.py +++ b/gridded/grids.py @@ -3,7 +3,7 @@ from gridded.pyugrid.ugrid import UGrid import numpy as np -from gridded.utilities import get_dataset, gen_celltree_mask_from_center_mask +from gridded.utilities import get_dataset, parse_filename_dataset_args class GridBase(object): @@ -90,20 +90,6 @@ def _find_required_grid_attrs(cls, filename, dataset=None, grid_topology=None,): def shape(self): return self.node_lon.shape - def __eq__(self, o): - if self is o: - return True - for n in ('nodes', 'faces'): - if (hasattr(self, n) and - hasattr(o, n) and - getattr(self, n) is not None and - getattr(o, n) is not None): - s = getattr(self, n) - s2 = getattr(o, n) - if s.shape != s2.shape or np.any(s != s2): - return False - return True - def _write_grid_to_file(self, pth): self.save_as_netcdf(pth) @@ -153,21 +139,22 @@ def _find_required_grid_attrs(cls, filename, dataset=None, grid_topology=None): if n in face_attrs: init_args[n] = gf_vars[v][:] break - # fixme: This is assuming that the array will be in Fortran order and index from 1, or in C order and index from 0 + # fixme: This is assuming that the array will be in Fortran order and index + # from 1, or in C order and index from 0 # Those are actually independent concepts! if init_args['faces'].shape[0] == 3: init_args['faces'] = np.ascontiguousarray(np.array(init_args['faces']).T - 1) - print("found grid vars:", init_args, gt) return init_args, gt - @classmethod - def gen_from_quads(cls, nodes): - if not len(nodes.shape) == 3: - raise ValueError('Nodes of a quad grid must be 2 dimensional') - lin_nodes = None - if isinstance(nodes, np.ma.MaskedArray): - lin_nodes = nodes.reshape(-1, 2)[nodes] + # @classmethod + # def gen_from_quads(cls, nodes): + # # Fixme: this looks incomplete -- used anywhere? + # if not len(nodes.shape) == 3: + # raise ValueError('Nodes of a quad grid must be 2 dimensional') + # lin_nodes = None + # if isinstance(nodes, np.ma.MaskedArray): + # lin_nodes = nodes.reshape(-1, 2)[nodes] class Grid_S(GridBase, SGrid): @@ -178,6 +165,7 @@ def _find_required_grid_attrs(cls, filename, dataset=None, grid_topology=None): # THESE ARE ACTUALLY ALL OPTIONAL. This should be migrated when optional attributes # are dealt with # Get superset attributes + # get_datset is not defined -- this must not be used gf_vars = dataset.variables if dataset is not None else get_dataset(filename).variables gf_vars = dict([(k.lower(), v) for k, v in gf_vars.items()]) init_args, gt = super(Grid_S, cls)._find_required_grid_attrs(filename, @@ -201,16 +189,17 @@ def _find_required_grid_attrs(cls, filename, dataset=None, grid_topology=None): edge2_mask_names = ['mask_v'] if grid_topology is None: - for attr, names, maskattr, maskname in (zip((center_attrs, edge1_attrs, edge2_attrs), - (center_coord_names, edge1_coord_names, edge2_coord_names), - (center_mask, edge1_mask, edge2_mask), - (center_mask_names, edge1_mask_names, edge2_mask_names))): + for attr, names, maskattr, maskname in zip( + (center_attrs, edge1_attrs, edge2_attrs), + (center_coord_names, edge1_coord_names, edge2_coord_names), + (center_mask, edge1_mask, edge2_mask), + (center_mask_names, edge1_mask_names, edge2_mask_names)): for n1, n2 in names: if n1 in gf_vars and n2 in gf_vars: mask = False - #for n in maskname: - #if n in gf_vars: - #mask = gen_mask(gf_vars[n]) + # for n in maskname: + # if n in gf_vars: + # mask = gen_mask(gf_vars[n]) a1 = gf_vars[n1][:] a2 = gf_vars[n2][:] init_args[attr[0]] = a1 @@ -222,7 +211,7 @@ def _find_required_grid_attrs(cls, filename, dataset=None, grid_topology=None): gt[attr[1]] = n2 break if 'node_lon' in init_args and 'node_lat' in init_args: - mask = False + mask = False # fixme -- is mask used?? for name in node_mask_names: if name in gf_vars: init_args[node_mask] = gf_vars[name] @@ -241,6 +230,7 @@ class Grid_R(GridBase): lon and lat of the nodes are vectors """ + def __init__(self, node_lon=None, node_lat=None, @@ -266,7 +256,7 @@ def __init__(self, self.node_dimensions = node_dimensions self.node_coordinates = node_coordinates - super(Grid_R, self).__init__(*args,**kwargs) + super(Grid_R, self).__init__(*args, **kwargs) @classmethod def _find_required_grid_attrs(cls, filename, dataset=None, grid_topology=None): @@ -275,19 +265,21 @@ def _find_required_grid_attrs(cls, filename, dataset=None, grid_topology=None): # are dealt with # Get superset attributes gf_vars = dataset.variables if dataset is not None else get_dataset(filename).variables - gf_vars = dict([(k.lower(), v) for k, v in gf_vars.items()] ) + gf_vars = dict([(k.lower(), v) for k, v in gf_vars.items()]) init_args, gt = super(Grid_R, cls)._find_required_grid_attrs(filename, dataset=dataset, grid_topology=grid_topology) - - # Grid_R only needs node_lon and node_lat. However, they must be a specific shape (1D) + # Grid_R only needs node_lon and node_lat. + # However, they must be a specific shape (1D) node_lon = init_args['node_lon'] node_lat = init_args['node_lat'] if len(node_lon.shape) != 1: - raise ValueError('Too many dimensions in node_lon. Must be 1D, was {0}D'.format(len(node_lon.shape))) + raise ValueError('Too many dimensions in node_lon. ' + 'Must be 1D, was {0}D'.format(len(node_lon.shape))) if len(node_lat.shape) != 1: - raise ValueError('Too many dimensions in node_lat. Must be 1D, was {0}D'.format(len(node_lat.shape))) + raise ValueError('Too many dimensions in node_lat. ' + 'Must be 1D, was {0}D'.format(len(node_lat.shape))) return init_args, gt @property @@ -335,13 +327,13 @@ def locate_faces(self, # lon_idxs[i] = -1 lat_idxs = np.digitize(lats, self.node_lat) - 1 for i, n in enumerate(lat_idxs): - if n == len(self.node_lat) -1: + if n == len(self.node_lat) - 1: lat_idxs[i] = -1 # if n == 0 and not lats[i] < self.node_lat.max() and not lats[i] >= self.node_lat.min(): # lat_idxs[i] = -1 idxs = np.column_stack((lon_idxs, lat_idxs)) - idxs[:,0] = np.where(idxs[:,1] == -1, -1, idxs[:,0]) - idxs[:,1] = np.where(idxs[:,0] == -1, -1, idxs[:,1]) + idxs[:, 0] = np.where(idxs[:, 1] == -1, -1, idxs[:, 0]) + idxs[:, 1] = np.where(idxs[:, 0] == -1, -1, idxs[:, 1]) if just_one: res = idxs[0] return res @@ -350,11 +342,12 @@ def locate_faces(self, def lonlat_to_yx(self, variable): ''' - The RegualarGridInterpolator needs to have it's two dimensions x and y be associated + The RegualarGridInterpolator needs to have its two dimensions x and y be associated correctly with lon and lat (or vice versa). The order depends on the orientation in the variable - if the variable provided does not have a dimensions attribute, it will use the dimensions arg + if the variable provided does not have a dimensions attribute, + it will use the dimensions arg ''' retval = (self.node_lat, self.node_lon) var_shape = variable.shape[-2::] @@ -363,27 +356,30 @@ def lonlat_to_yx(self, variable): if not all([k in self.dimensions for k in variable.dimensions[-2:]]): raise ValueError('Dimension provided by variable is not compatible \ with this Grid_R object. Provided: {0} \ - self.dimensions: {1}'.format(variable.dimensions, self.dimensions)) - var_dims = variable.dimensions[-2:] #assume the last two are the lon/lat x/y + self.dimensions: {1}'.format(variable.dimensions, + self.dimensions)) + var_dims = variable.dimensions[-2:] # assume the last two are the lon/lat x/y else: var_dims = self.dimensions - #self.dimensions is always [y(lat), x(lon)], so if var.dimensions is [lon, lat] we need - #to reverse x/y association + # self.dimensions is always [y(lat), x(lon)], + # so if var.dimensions is [lon, lat] we need + # to reverse x/y association if not all([dlen in grid_shape for dlen in var_shape]): - raise ValueError('Incompatible dimensions. Variable: {0}, Grid_R: {1}'.format(variable.shape, grid_shape)) - + raise ValueError('Incompatible dimensions. ' + 'Variable: {0}, Grid_R: {1}'.format(variable.shape, grid_shape)) + if hasattr(variable, 'dimensions'): - if var_dims[0] == self.dimensions[1]: #case 2, dims provided, dims swapped + if var_dims[0] == self.dimensions[1]: # case 2, dims provided, dims swapped retval = retval[::-1] - #else: case 1, no change + # else: case 1, no change else: - if var_shape[0] == var_shape[1]: #case 4, dims not provided, dim length same + if var_shape[0] == var_shape[1]: # case 4, dims not provided, dim length same raise ValueError('Provided square variable with no dimensions attribute') - if var_shape[0] == len(self.node_lon): #case 3, dims not provided, dim length differnt + # case 3, dims not provided, dim length different + if var_shape[0] == len(self.node_lon): retval = retval[::-1] return retval - def interpolate_var_to_points(self, points, @@ -406,7 +402,7 @@ def interpolate_var_to_points(self, if slices is not None: variable = variable[slices] if np.ma.isMA(variable): - variable = variable.filled(0) #eventually should use Variable fill value + variable = variable.filled(0) # eventually should use Variable fill value interp_func = RegularGridInterpolator((y, x), variable, method=method, @@ -433,12 +429,12 @@ def infer_location(self, variable): # centers_shape = self.centers.shape[0:-1] try: shape = np.array(variable.shape) - except: + except: # fixme -- AttributeError?? return None # Variable has no shape attribute! if len(variable.shape) < 2: return None difference = (shape[-2:] - node_shape).tolist() - if (difference == [1, 1] or difference == [-1, -1]) and self.center_lon is not None: + if (difference == [1, 1] or difference == [-1, -1]) and self.center_lon is not None: return 'center' elif difference == [1, 0] and self.edge1_lon is not None: return 'edge1' @@ -473,7 +469,7 @@ def _load_grid(filename, grid_type, dataset=None): Redirect to grid-specific loading routine. ''' if issubclass(grid_type, UGrid): - return grid_type.from_ncfile(filename) + return grid_type.from_ncfile(filename=filename, dataset=dataset) elif issubclass(grid_type, SGrid): ds = get_dataset(filename, dataset) g = grid_type.load_grid(ds) @@ -486,6 +482,8 @@ def _load_grid(filename, grid_type, dataset=None): @staticmethod def from_netCDF(filename=None, dataset=None, + data_file=None, + grid_file=None, grid_type=None, grid_topology=None, _default_types=(('ugrid', Grid_U), @@ -509,28 +507,29 @@ def from_netCDF(filename=None, :returns: Instance of Grid_U, Grid_S, or Grid_R ''' - gf = dataset if filename is None else get_dataset(filename, dataset) - if gf is None: - raise ValueError('No filename or dataset provided') + df, dg = parse_filename_dataset_args(filename=filename, + dataset=dataset, + data_file=data_file, + grid_file=grid_file) cls = grid_type if (grid_type is None or isinstance(grid_type, str) or not issubclass(grid_type, GridBase)): - cls = Grid._get_grid_type(gf, grid_type, grid_topology, _default_types) + cls = Grid._get_grid_type(dg, grid_type, grid_topology, _default_types) # if grid_topology is passed in, don't look for the variable if not grid_topology: - compliant = Grid._find_topology_var(None, gf) + compliant = Grid._find_topology_var(None, dg) else: compliant = None if compliant is not None: - c = Grid._load_grid(filename, cls, dataset) + c = Grid._load_grid(filename, cls, dg) c.grid_topology = compliant.__dict__ else: init_args, gt = cls._find_required_grid_attrs(filename, - dataset=dataset, + dataset=dg, grid_topology=grid_topology) c = cls(grid_topology=gt, **init_args) return c @@ -542,12 +541,12 @@ def _get_grid_type(dataset, _default_types=(('ugrid', Grid_U), ('sgrid', Grid_S), ('rgrid', Grid_R))): - # fixme: this logic should probably be defered to + # fixme: this logic should probably be deferred to # the grid type code -- that is, ask each grid # type if this dataset is its type. # # It also should be refactored to start with the standards - # and maybe havev a pedantic mode where it won't load non-standard + # and maybe have a pedantic mode where it won't load non-standard # files if _default_types is None: diff --git a/gridded/pysgrid/sgrid.py b/gridded/pysgrid/sgrid.py index efa28b4..909e943 100644 --- a/gridded/pysgrid/sgrid.py +++ b/gridded/pysgrid/sgrid.py @@ -118,6 +118,30 @@ def __init__(self, self._cell_tree_mask = None self.grid_topology = grid_topology + def __eq__(self, other): + # Fixme: should this even exist + # keeping it because it's used in a test. + if self is other: + return True + # maybe too strict? + if self.__class__ is not other.__class__: + return False + # there is way too much to really check here ... + # but it would a heck of a coincidence ... + for attr in ('node_lon', + 'node_lat', + 'node_mask', + 'center_lon', + 'center_lat', + 'center_mask', + 'faces', + 'face_padding', + ): + if not np.array_equal(getattr(self, attr), getattr(other, attr)): + return False + return True + + @classmethod def load_grid(cls, nc): if isinstance(nc, Dataset): @@ -610,6 +634,18 @@ def locate_faces(self, if _memo: self._add_memo(points, res, self._cell_ind_memo_dict, _copy, _hash) return res + + def index_of(self, + points, + _memo=False, + _copy=False, + _hash=None, + use_mask=True): + return self.locate_faces(points=points, + _memo=_memo, + _copy=_copy, + _hash=_hash, + use_mask=use_mask) def locate_nearest(self, points, diff --git a/gridded/pyugrid/ugrid.py b/gridded/pyugrid/ugrid.py index ddd262b..0e36ac6 100644 --- a/gridded/pyugrid/ugrid.py +++ b/gridded/pyugrid/ugrid.py @@ -154,7 +154,7 @@ def __init__(self, self._alpha_memo_dict = OrderedDict() @classmethod - def from_ncfile(klass, nc_url, mesh_name=None): # , load_data=False): + def from_ncfile(klass, filename=None, dataset=None, mesh_name=None): # , load_data=False): """ create a UGrid object from a netcdf file name (or opendap url) @@ -166,8 +166,12 @@ def from_ncfile(klass, nc_url, mesh_name=None): # , load_data=False): Will be raised """ grid = klass() - read_netcdf.load_grid_from_ncfilename(nc_url, grid, mesh_name) # , load_data) - return grid + if filename is not None: + read_netcdf.load_grid_from_ncfilename(filename, grid, mesh_name) + return grid + if dataset is not None: + read_netcdf.load_grid_from_nc_dataset(dataset, grid, mesh_name) + return grid @classmethod def from_nc_dataset(klass, nc, mesh_name=None): # , load_data=False): @@ -213,6 +217,34 @@ def info(self): # msg.append("Variables: " + ", ".join([str(v) for v in self._data.keys()])) return "\n".join(msg) + def __eq__(self, other): + # Fixme: should this even exist + # keeping it because it's used in a test. + if self is other: + return True + # maybe too strict? + if self.__class__ is not other.__class__: + return False + if (np.array_equal(self.nodes, other.nodes) + and np.array_equal(self.faces, other.faces) + and np.array_equal(self.boundaries, other.boundaries) + ): + return True + return False + + # They should all be there! + # for n in ('nodes', 'faces'): + # if (hasattr(self, n) and + # hasattr(o, n) and + # getattr(self, n) is not None and + # getattr(o, n) is not None): + # s = getattr(self, n) + # s2 = getattr(o, n) + # if s.shape != s2.shape or np.any(s != s2): + # return False + return True + + def check_consistent(self): """ Check if the various data is consistent: the edges and faces reference @@ -538,7 +570,7 @@ def _get_memoed(self, points, D, _copy=False, _hash=None): else: return None - def locate_faces(self, points, method='celltree', _copy=False, _memo=True, _hash=None): + def locate_faces(self, points, method='celltree', _memo=True, _copy=False, _hash=None): """ Returns the face indices, one per point. @@ -602,6 +634,20 @@ def locate_faces(self, points, method='celltree', _copy=False, _memo=True, _hash return indices[0] else: return indices + + def index_of(self, + points, + method='celltree', + _memo=False, + _copy=False, + _hash=None, + use_mask=True): + return self.locate_faces(points=points, + method=method, + _memo=_memo, + _copy=_copy, + _hash=_hash, + use_mask=use_mask) def build_celltree(self): """ @@ -742,7 +788,9 @@ def interpolate_var_to_points(self, # "variable defined on the faces") if location == 'node': pos_alphas = self.interpolation_alphas(points, inds, _copy, _memo, _hash) - vals = variable[self.faces[inds]] + vals = variable.take(self.faces[inds], axis=0) + if len(vals.shape) == (len(pos_alphas.shape) + 1): + pos_alphas = pos_alphas[:,:,np.newaxis] vals[inds == -1] = vals[inds == -1] * 0 return np.sum(vals * pos_alphas, axis=1) return None diff --git a/gridded/tests/test_data/UGRIDv0.9_eleven_points_with_depth.nc b/gridded/tests/test_data/UGRIDv0.9_eleven_points_with_depth.nc new file mode 100644 index 0000000..2e1edd8 Binary files /dev/null and b/gridded/tests/test_data/UGRIDv0.9_eleven_points_with_depth.nc differ diff --git a/gridded/tests/test_data/make_sample_data.py b/gridded/tests/test_data/make_sample_data.py index 4e35bd8..4010f1a 100644 --- a/gridded/tests/test_data/make_sample_data.py +++ b/gridded/tests/test_data/make_sample_data.py @@ -9,8 +9,11 @@ """ import make_sgrid_files +import make_ugrid_files make_sgrid_files.make_all() +make_ugrid_files.make_all() + diff --git a/gridded/tests/test_data/make_ugrid_files.py b/gridded/tests/test_data/make_ugrid_files.py new file mode 100644 index 0000000..d4e94d5 --- /dev/null +++ b/gridded/tests/test_data/make_ugrid_files.py @@ -0,0 +1,101 @@ +import netCDF4 as nc +import numpy as np +import shutil + + +def fvcom_eleven_points_with_depth(): + '''Use UGRIDv0.9_eleven_points.nc as a starting point. Add sigma depth coordinates and provide better values for the bathymetry''' + fn = './UGRIDv0.9_eleven_points.nc' + shutil.copy(fn, './UGRIDv0.9_eleven_points_with_depth.nc') + + df = nc.Dataset('./UGRIDv0.9_eleven_points_with_depth.nc', mode = 'a') + df.createDimension('siglev', 5) + df.createDimension('siglay', 4) + df.createDimension('t', 1) + + df.createVariable('time', 'f8', ('t',)) + df.createVariable('siglay', 'f8', ('siglay', 'nMesh2_node')) + df.createVariable('siglev', 'f8', ('siglev', 'nMesh2_node')) + df.createVariable('siglay_center', 'f8', ('siglay', 'nMesh2_face')) + df.createVariable('siglev_center', 'f8', ('siglev', 'nMesh2_face')) + df.createVariable('h', 'f8', ('nMesh2_node',)) + df.createVariable('h_center', 'f8', ('nMesh2_face',)) + df.createVariable('zeta', 'f8', ('t', 'nMesh2_node',)) + + df['time'].units = 'days since 1970-01-01 00:00:00' + df['time'].calendar = 'gregorian' + df['time'].long_name = 'time' + df['time'].time_zone = 'UTC' + df['time'].format = 'defined reference date' + df['time'][0] = 19636.541 + + siglay_attrs = {'long_name': 'Sigma Layers', + 'standard_name': 'ocean_sigma_coordinate', + 'positive': 'up', + 'valid_min': -1, + 'valid_max': 0, + 'formula_terms': 'sigma: siglay eta: zeta depth: h'} + df['siglay'].setncatts(siglay_attrs) + df['siglay_center'].setncatts(siglay_attrs) + df['siglay_center'].formula_terms = 'sigma: siglay_center eta: zeta_center depth: h_center' + + siglev_attrs = {'long_name': 'Sigma Levels', + 'standard_name': 'ocean_sigma_coordinate', + 'positive': 'up', + 'valid_min': -1, + 'valid_max': 0, + 'formula_terms': 'sigma: siglay eta: zeta depth: h'} + df['siglev'].setncatts(siglev_attrs) + df['siglev_center'].setncatts(siglev_attrs) + df['siglev_center'].formula_terms = 'sigma: siglev_center eta: zeta_center depth: h_center' + bathy_attrs = { + 'long_name': 'Bathymetry', + 'standard_name': 'sea_floor_depth_below_geoid', + 'units': 'm', + 'positive': 'down', + 'grid': 'Bathymetry_Mesh', + 'coordinates': 'Mesh2_node_y Mesh2_node_x', + 'type': 'data' + } + df['h'].setncatts(bathy_attrs) + df['h_center'].setncatts(bathy_attrs) + df['h_center'].grid = 'Bathymetry_Mesh2' + df['h_center'].coordinates = 'Mesh2_face_y Mesh2_face_x' + df['h_center'].grid_location = 'center' + df['zeta'][:] = np.zeros(df['zeta'].shape) + + zeta_attrs = { + 'long_name': 'Water Surface Elevation', + 'units': 'meters', + 'positive': 'up', + 'standard_name': 'sea_surface_height_above_geoid', + 'grid': 'Bathymetry_Mesh', + 'coordinates': 'time Mesh2_node_y Mesh2_node_x', + 'type': 'data', + 'location': 'node' + } + df['zeta'].setncatts(zeta_attrs) + + df['h'][:] = np.linspace(0, 10, df['h'].shape[0]).reshape(df['h'].shape) + df['h_center'][:] = np.ones(df['h_center'].shape) * 10 + df['Mesh2_depth'][:] = df['h'][:] + df['siglev'][:] = np.linspace(0, -1, 5).T[:,np.newaxis] + df['siglay'][:] = np.linspace(0, -1, 9)[1::2].T[:, np.newaxis] + + + df.close() + + +def make_all(): + """ + make all the sample files + """ + + fvcom_eleven_points_with_depth() + # semi_circular_single_grid() + + +if __name__ == "__main__": + make_all() + # fname = semi_circular_single_grid() + # print("making:", fname) \ No newline at end of file diff --git a/gridded/tests/test_depth.py b/gridded/tests/test_depth.py index a1e5e59..17a59b6 100644 --- a/gridded/tests/test_depth.py +++ b/gridded/tests/test_depth.py @@ -7,6 +7,7 @@ import numpy as np import netCDF4 as nc import gridded +from math import sqrt from gridded.variable import Variable, VectorVariable from gridded.tests.utilities import get_test_file_dir, get_test_cdl_filelist @@ -14,7 +15,7 @@ from gridded.time import Time from gridded.utilities import search_dataset_for_variables_by_varname -from gridded.depth import S_Depth, L_Depth +from gridded.depth import S_Depth, ROMS_Depth, FVCOM_Depth, L_Depth test_dir = get_test_file_dir() cdl_files = get_test_cdl_filelist() @@ -34,12 +35,16 @@ def test_from_netCDF(): ds = nc.Dataset.fromcdl(fn, ncfilename=fn+'.nc') if not valid_depth_test_dataset(ds): continue - d = S_Depth.from_netCDF(dataset=ds) - + d = S_Depth.from_netCDF(data_file=ds, grid_file=ds) + +@pytest.fixture(scope="module") +def get_fvcom_depth(): + ds = nc.Dataset(os.path.join(test_dir, 'UGRIDv0.9_eleven_points_with_depth.nc')) + return FVCOM_Depth.from_netCDF(data_file=ds, grid_file=ds) @pytest.fixture(scope="module") -def get_s_depth(): +def get_roms_depth(): """ This is setup for a ROMS S-level Depth that is on a square grid with a center mound. Control vars: sz=xy size, center_el=height of the mound in meters, @@ -88,12 +93,10 @@ def get_s_depth(): s_rho = (s_w[0:-1] + s_w[1:]) / 2 # equidistant layers, no stretching Cs_w = np.linspace(-1, 0, nz) - Cs_w = 1 - 1 / np.exp(2 * Cs_w) - Cs_w /= -Cs_w[0] Cs_r = (Cs_w[0:-1] + Cs_w[1:]) / 2 - hc = np.array([0]) + hc = np.array([0,]) - sd = S_Depth( + sd = ROMS_Depth( time=zeta.time, grid=zeta.grid, bathymetry=bathy, @@ -129,57 +132,44 @@ def get_l_depth(): return ld1, ld2 -class Test_S_Depth(object): - def test_construction(self, get_s_depth): - assert get_s_depth is not None - - def test_structure(self, get_s_depth): - sd = get_s_depth - sd.Cs_w = np.linspace(-1, 0, len(sd.Cs_w)) - sd.Cs_r = (sd.Cs_w[0:-1] + sd.Cs_w[1:]) / 2 - - sz = sd.bathymetry.data.shape[1] - levels = sd.get_section(sd.time.data[1], "w")[ - :, :, : - ] # 3D cross section - edge = np.linspace(20, 0, len(sd.Cs_w)) - center = np.linspace(10, 0, len(sd.Cs_w)) - assert np.allclose(levels[:, 0, 0], edge) - assert np.allclose( - levels[:, int(sz / 2), int(sz / 2)], center - ) +class Test_ROMS_Depth(object): + def test_construction(self, get_roms_depth): + assert get_roms_depth is not None + assert get_roms_depth.num_levels == 11 - def test_interpolation_alphas_bottom(self, get_s_depth): + def test_interpolation_alphas_bottom(self, get_roms_depth): ''' We will focus on the center mound to see if the correct alphas and indices are returned for various points nearby. interpolation_alphas(self, points, time, data_shape, _hash=None, extrapolate=False): ''' - sd = get_s_depth + sd = get_roms_depth # query center mound (20,20) at 3 depths. 1st point is 0.1m off the seafloor and 10% of # the distance to the next s-layer (meaning, expected alpha returned is 10%) # 2nd point is directly on the seafloor (should register as 'in bounds' and 0 alpha) - # 3rd point is 0.1m underground, and should indicate with -2 alpha + # 3rd point is 0.1m underground, and should indicate with masked value. + # 4th point is off grid, masked expected. + # 5th point is above grid, 0 or masked expected depending on boundary condition points = np.array([[20,20, 9.9],[20,20,10.0], [20,20,10.1], [-1, -1, 5], [20, 20, -0.1]]) ts = sd.time.data[1] assert sd.default_bottom_boundary_condition == 'mask' - idx, alphas = sd.interpolation_alphas(points, ts, [sd.num_w_levels,]) + idx, alphas = sd.interpolation_alphas(points, ts, [sd.num_levels,]) expected_idx = np.ma.array(np.array([0,0,-1,-1,10]), mask = [False, False, True, True, False]) expected_alpha = np.ma.array(np.array([0.1, 0, -1, -1, 0]), mask = [False, False, True, True, False]) assert np.all(idx == expected_idx) assert np.all(np.isclose(alphas, expected_alpha)) sd.default_bottom_boundary_condition == 'extrapolate' - idx, alphas = sd.interpolation_alphas(points, ts, [sd.num_w_levels,]) + idx, alphas = sd.interpolation_alphas(points, ts, [sd.num_levels,]) expected_idx = np.ma.array(np.array([0,0,-1]), mask = [False, False, False]) expected_alpha = np.ma.array(np.array([0.1, 0, 1]), mask = [False, False, False]) - def test_interpolation_alphas_surface(self, get_s_depth): - sd = get_s_depth + def test_interpolation_alphas_surface(self, get_roms_depth): + sd = get_roms_depth points = np.array([[20,20, 0],[20,20,0.1], [20,20,-0.1]]) ts = sd.time.data[1] assert sd.default_surface_boundary_condition == 'extrapolate' - idx, alphas = sd.interpolation_alphas(points, ts, [sd.num_w_levels,]) + idx, alphas = sd.interpolation_alphas(points, ts, [sd.num_levels,]) # only the element 0.1m deep should register with an index and alpha since # it is the only element below surface. expected_idx = np.ma.array(np.array([10,9,10]), mask = [False, False, False]) @@ -188,7 +178,7 @@ def test_interpolation_alphas_surface(self, get_s_depth): assert np.all(np.isclose(alphas, expected_alpha)) sd.default_surface_boundary_condition == 'mask' - idx, alphas = sd.interpolation_alphas(points, ts, [sd.num_w_levels,]) + idx, alphas = sd.interpolation_alphas(points, ts, [sd.num_levels,]) # only the element 0.1m deep should register with an index and alpha since # it is the only element below surface. expected_idx = np.ma.array(np.array([-1,9,-1]), mask = [True, False, True]) @@ -199,20 +189,37 @@ def test_interpolation_alphas_surface(self, get_s_depth): sd.default_surface_boundary_condition == 'extrapolate' # switch to timestep with -0.5m zeta ts = sd.time.data[0] - idx, alphas = sd.interpolation_alphas(points, ts, [sd.num_w_levels,]) + idx, alphas = sd.interpolation_alphas(points, ts, [sd.num_levels,]) expected_idx = np.ma.array(np.array([10,10,10]), mask = [False, False, False]) expected_alpha = np.ma.array(np.array([0, 0, 0]), mask = [False, False, False]) assert np.all(idx == expected_idx) assert np.all(np.isclose(alphas, expected_alpha)) sd.default_surface_boundary_condition = 'mask' - idx, alphas = sd.interpolation_alphas(points, ts, [sd.num_w_levels,]) + idx, alphas = sd.interpolation_alphas(points, ts, [sd.num_levels,]) # only the element 0.1m deep should register with an index and alpha since # it is the only element below surface. expected_idx = np.ma.array(np.array([-1,-1,-1]), mask = [True, True, True]) expected_alpha = np.ma.array(np.array([0, 0.9, 0]), mask = [True, True, True]) assert np.all(idx.mask == expected_idx.mask) assert np.all(alphas.mask == expected_alpha.mask) + +class Test_FVCOM_Depth(object): + def test_construction(self, get_fvcom_depth): + assert get_fvcom_depth is not None + + def test_get_transect(self, get_fvcom_depth): + dp = get_fvcom_depth + tris = dp.grid.nodes.take(dp.grid.faces, axis=0) + centroids = np.mean(tris, axis=1) + transects = dp.get_transect(centroids, dp.zeta.time.min_time) + + #because we use the centroids, the values should the average of the 3 nodes + #Not intending to test interpolation here there's a separate test for that + + expected_transects = np.mean((dp.siglev * dp.bathymetry.data[:]).take(dp.grid.faces, axis=1), axis=2) + expected_transects = expected_transects.T * -1 + assert np.all(np.isclose(transects, expected_transects)) @pytest.fixture(scope="module") def get_database_nc(): @@ -229,6 +236,7 @@ def get_database_nc(): return time, depth, ds class Test_L_Depth(object): + def test_construction(self, get_l_depth): assert get_l_depth is not None @@ -333,4 +341,4 @@ def test_vertical_interpolation_full(self, get_database_nc): assert rtv[0] == ds.variables['u'].data[0,0,1,1] assert rtv[1] == ds.variables['u'].data[0,0,1,1] - assert rtv[2] == ds.variables['u'].data[0,3,1,1] \ No newline at end of file + assert rtv[2] == ds.variables['u'].data[0,3,1,1] diff --git a/gridded/tests/test_pysgrid/test_sgrid.py b/gridded/tests/test_pysgrid/test_sgrid.py new file mode 100644 index 0000000..4602dbb --- /dev/null +++ b/gridded/tests/test_pysgrid/test_sgrid.py @@ -0,0 +1,25 @@ +""" +amazingly, there were no tests for sgrid itself .. + +now there are almost none +""" + +from gridded.pysgrid.sgrid import SGrid, load_grid +from .write_nc_test_files import roms_sgrid + + +def test_eq_same(roms_sgrid): + grid1 = load_grid(roms_sgrid) + grid2 = load_grid(roms_sgrid) + + assert grid1 == grid2 + + +def test_eq_diff_type(roms_sgrid): + """ + Just to make sure it doesn't crash + """ + grid1 = load_grid(roms_sgrid) + + assert grid1 != "a string" + diff --git a/gridded/tests/test_time.py b/gridded/tests/test_time.py index 5c322f3..9f478c4 100644 --- a/gridded/tests/test_time.py +++ b/gridded/tests/test_time.py @@ -1,18 +1,323 @@ #!/usr/bin/env python -# py2/3 compatibility +import copy +from datetime import datetime, timedelta -from gridded.time import Time +import numpy as np + +from gridded.time import Time, TimeSeriesError, OutOfTimeRangeError + +import pytest + +SAMPLE_TIMESERIES = [] +start = datetime(2023, 11, 28, 12) +dt = timedelta(minutes=15) +SAMPLE_TIMESERIES = [(start + i * dt) for i in range(10)] +STS = SAMPLE_TIMESERIES +# import copy +# from datetime import datetime, timedelta + +# import numpy as np + +# from gridded.time import Time + +import pytest + +SAMPLE_TIMESERIES = [] +start = datetime(2023, 11, 28, 12) +dt = timedelta(minutes=15) +for i in range(10): + SAMPLE_TIMESERIES.append(start + i * dt) +STS = SAMPLE_TIMESERIES def test_init(): """ - can one even be initialized? + can one even be initialized with no data? """ t = Time() +def test_init_with_timeseries(): + t = Time(SAMPLE_TIMESERIES) + + # an implementation detail -- kind of + assert isinstance(t.data, np.ndarray) + + +def test_origin(): + origin = datetime(1984, 1, 1, 6) + t = Time(SAMPLE_TIMESERIES, origin=origin) + + assert t.data[0] == origin + assert t.data[-1] == origin + (SAMPLE_TIMESERIES[-1] - SAMPLE_TIMESERIES[0]) + + +def test_displacement(): + disp = -timedelta(days=35, hours=12) + t = Time(SAMPLE_TIMESERIES, displacement=disp) + + assert t.data[0] == SAMPLE_TIMESERIES[0] + disp + assert t.data[-1] == SAMPLE_TIMESERIES[-1] + disp + + +def test_tz_offset(): + offset = -8 + offset = timedelta(hours=offset) + t = Time(SAMPLE_TIMESERIES, tz_offset=offset) + + assert t.data[0] == SAMPLE_TIMESERIES[0] + offset + assert t.data[-1] == SAMPLE_TIMESERIES[-1] + offset + + +def test_tz_offset_hours(): + offset = -8 + t = Time(SAMPLE_TIMESERIES, tz_offset=offset) + + offset = timedelta(hours=offset) + assert t.data[0] == SAMPLE_TIMESERIES[0] + offset + assert t.data[-1] == SAMPLE_TIMESERIES[-1] + offset + + +def test_get_time_array(): + t = Time(SAMPLE_TIMESERIES) + + ta = t.get_time_array() + + assert ta is not t.data + assert np.array_equal(ta, SAMPLE_TIMESERIES) + + +def test_info(): + # just to make sure it's not broken + t = Time(SAMPLE_TIMESERIES) + + info = t.info + + print(info) + + assert True + + +def test_iter(): + t = Time(SAMPLE_TIMESERIES) + data2 = list(t) + + assert data2 == SAMPLE_TIMESERIES + + +def test_reset_data(): + t = Time() + t.data = SAMPLE_TIMESERIES + + assert np.array_equal(t.data, SAMPLE_TIMESERIES) +<<<<<<< HEAD + # an implementation detail, but want to make sure +======= + # an implementaiton detail, but want to make sure +>>>>>>> 173f1af0c694edfd74bf9599da47e205761e5f7e + assert isinstance(t.data, np.ndarray) + + +def test_eq(): + t1 = Time(data=copy.copy(SAMPLE_TIMESERIES)) + t2 = Time(data=copy.copy(SAMPLE_TIMESERIES)) + + assert t1 == t2 + + +def test_eq_diff_length(): + data2 = copy.copy(SAMPLE_TIMESERIES) + del data2[-1] + t1 = Time(data=copy.copy(SAMPLE_TIMESERIES)) + t2 = Time(data=data2) + + assert t1 != t2 + + +def test_eq_diff_values(): + data2 = copy.copy(SAMPLE_TIMESERIES) + data2[4] = data2[4] + timedelta(minutes=5) + t1 = Time(data=copy.copy(SAMPLE_TIMESERIES)) + t2 = Time(data=data2) + + assert t1 != t2 + + +def test_eq_diff_one_constant(): + data2 = copy.copy(SAMPLE_TIMESERIES) + t1 = Time(data=data2) + t2 = Time(data=[data2[2]]) + + assert t1 != t2 + + +def test_eq_constant_time(): + t1 = Time.constant_time() + t2 = Time.constant_time() + + assert t1 == t2 + +<<<<<<< HEAD + +======= +>>>>>>> 173f1af0c694edfd74bf9599da47e205761e5f7e +def test_eq_different_type(): + """ + A Time object is never equal to any other type + """ + t = Time(SAMPLE_TIMESERIES) + + assert not t == 'a string' + assert not "a string" == t + + +def test_out_of_order(): + data2 = copy.copy(SAMPLE_TIMESERIES) + data2.insert(1, data2[4]) + +<<<<<<< HEAD + with pytest.raises(TimeSeriesError): +======= + with pytest.raises(ValueError): +>>>>>>> 173f1af0c694edfd74bf9599da47e205761e5f7e + Time(data2) + + +def test_decending(): + data2 = copy.copy(SAMPLE_TIMESERIES) + data2.reverse() + +<<<<<<< HEAD + with pytest.raises(TimeSeriesError): + Time(data2) + + +def test_duplicate(): + data2 = SAMPLE_TIMESERIES[:5] + SAMPLE_TIMESERIES[4:] + + with pytest.raises(TimeSeriesError): + Time(data2) + + +======= + with pytest.raises(ValueError): + Time(data2) + +>>>>>>> 173f1af0c694edfd74bf9599da47e205761e5f7e +def test_min_max(): + t = Time(SAMPLE_TIMESERIES) + + assert t.min_time == SAMPLE_TIMESERIES[0] + assert t.max_time == SAMPLE_TIMESERIES[-1] + + +def test_time_in_bounds(): + t = Time(SAMPLE_TIMESERIES) + + assert t.time_in_bounds(SAMPLE_TIMESERIES[0]) + assert t.time_in_bounds(SAMPLE_TIMESERIES[2] + timedelta(minutes=10)) + assert t.time_in_bounds(SAMPLE_TIMESERIES[-1]) + + assert not t.time_in_bounds(SAMPLE_TIMESERIES[0] - timedelta(seconds=1)) + assert not t.time_in_bounds(SAMPLE_TIMESERIES[-1] + timedelta(seconds=1)) + + +def test_valid_time(): + t = Time(SAMPLE_TIMESERIES) + + assert t.valid_time(SAMPLE_TIMESERIES[0]) is None + assert t.valid_time(SAMPLE_TIMESERIES[2] + timedelta(minutes=10)) is None + assert t.valid_time(SAMPLE_TIMESERIES[-1]) is None + +<<<<<<< HEAD + with pytest.raises(OutOfTimeRangeError): + t.valid_time(SAMPLE_TIMESERIES[0] - timedelta(seconds=1)) + with pytest.raises(OutOfTimeRangeError): +======= + with pytest.raises(ValueError): + t.valid_time(SAMPLE_TIMESERIES[0] - timedelta(seconds=1)) + with pytest.raises(ValueError): +>>>>>>> 173f1af0c694edfd74bf9599da47e205761e5f7e + t.valid_time(SAMPLE_TIMESERIES[-1] + timedelta(seconds=1)) + + +@pytest.mark.parametrize("dt, expected", + [(STS[3], 1.0), # exact time + (STS[4] + (STS[5] - STS[4]) / 2, 0.5), # in the middle + (STS[4] + (STS[5] - STS[4]) / 4, 0.25), # in the middle + (STS[0], 0.0), # at the beginning + (STS[-1], 1.0), # at the end + ]) +def test_interp_alpha(dt, expected): + t = Time(SAMPLE_TIMESERIES) + + print(dt) + alpha = t.interp_alpha(dt) + + assert alpha == expected + + +@pytest.mark.parametrize("dt, expected", + [(STS[0] - timedelta(seconds=1), 0.0), # a little before + (STS[-1] + timedelta(seconds=1), 1.0), # a little after + ]) +def test_interp_alpha_outside(dt, expected): + t = Time(SAMPLE_TIMESERIES) + + print(dt) +<<<<<<< HEAD + with pytest.raises(OutOfTimeRangeError): +======= + with pytest.raises(ValueError): +>>>>>>> 173f1af0c694edfd74bf9599da47e205761e5f7e + alpha = t.interp_alpha(dt) + + alpha = t.interp_alpha(dt, extrapolate=True) + assert alpha == expected + + +@pytest.mark.xfail +@pytest.mark.parametrize("shift, expected", + [(-timedelta(days=365), 1.0), # before + (timedelta(days=365), 1.0), # after + (timedelta(0), 1.0), # on the nose + ]) +def test_interp_alpha_constant_time(shift, expected): + """ +<<<<<<< HEAD + What should the constant time Time object give for alphas? +======= + What should the constant time Time give for alphas? +>>>>>>> 173f1af0c694edfd74bf9599da47e205761e5f7e + + I would expect always 1.0! It seems to be always 0.0 + + Maybe this isn't used anywhere? + """ + t = Time.constant_time() + + print(t.data[0]) + + # on the nose + alpha = t.interp_alpha(t.data[0] + shift) + assert alpha == expected + +<<<<<<< HEAD +# def test_index_of_contant_time(): +# pass + +======= +>>>>>>> 173f1af0c694edfd74bf9599da47e205761e5f7e + ## this needs a big data file -- could use some refactoring +## It would be good to test the netcdf stuff. +## there must be a sample file with time in it somewhere?? +<<<<<<< HEAD +## However -- this is testing a lot that doesn't need a data file. +======= +>>>>>>> 173f1af0c694edfd74bf9599da47e205761e5f7e + # class TestTime: # time_var = circular_3D['time'] # time_arr = nc.num2date(time_var[:], units=time_var.units) diff --git a/gridded/tests/test_ugrid/test_ugrid.py b/gridded/tests/test_ugrid/test_ugrid.py index 8e17d76..4bf6af7 100644 --- a/gridded/tests/test_ugrid/test_ugrid.py +++ b/gridded/tests/test_ugrid/test_ugrid.py @@ -4,7 +4,6 @@ The basic test of the UGrid object. More specific functionality is other test modules. - """ @@ -101,3 +100,43 @@ def test_both_nodes_and_lat(): nodes=nodes ) +def test_eq(): + grid1 = UGrid(nodes=nodes, + faces=faces, + edges=edges, + boundaries=boundaries, + ) + + grid2 = UGrid(nodes=nodes, + faces=faces, + edges=edges, + boundaries=boundaries, + ) + + assert grid1 == grid2 + + +def test_eq_no_bounds(): + grid1 = UGrid(nodes=nodes, + faces=faces, + edges=edges, + boundaries=boundaries, + ) + + grid2 = UGrid(nodes=nodes, + faces=faces, + edges=edges, + ) + + assert grid1 != grid2 + + +def test_eq_diff_type(): + # make sure it doesn't crash + grid1 = UGrid(nodes=nodes, + faces=faces, + edges=edges, + boundaries=boundaries, + ) + + assert grid1 != "A string" diff --git a/gridded/tests/test_utilities.py b/gridded/tests/test_utilities.py index 4053a74..1e121f8 100644 --- a/gridded/tests/test_utilities.py +++ b/gridded/tests/test_utilities.py @@ -8,7 +8,7 @@ import netCDF4 as nc from gridded import utilities -from gridded.tests.test_depth import get_s_depth +from gridded.tests.test_depth import get_roms_depth data_dir = os.path.join(os.path.split(__file__)[0], 'test_data') @@ -120,16 +120,16 @@ def test_spatial_data_metadata(): assert np.all(a5.T == res_5) -def test_regrid_variable_TDStoS(get_s_depth): +def test_regrid_variable_TDStoS(get_roms_depth): # Time is present # Depth is present # Grid_S to Grid_S from gridded.variable import Variable from gridded.time import Time from gridded.grids import Grid_S - sd = get_s_depth + sd = get_roms_depth grid = sd.grid - n_levels = sd.num_w_levels + n_levels = sd.num_levels data = np.ones((1, n_levels, grid.node_lon.shape[0], grid.node_lon.shape[1])) @@ -156,14 +156,14 @@ def test_regrid_variable_TDStoS(get_s_depth): assert v2.data.shape == (v1.data.shape[0], v1.data.shape[1], sz - 1, sz - 1) -def test_regrid_variable_StoS(get_s_depth): +def test_regrid_variable_StoS(get_roms_depth): # Time is not present # Depth is not present # Grid_S to Grid_S from gridded.variable import Variable from gridded.time import Time from gridded.grids import Grid_S - sd = get_s_depth + sd = get_roms_depth grid = sd.grid data = np.ones((grid.node_lon.shape[0], grid.node_lon.shape[1])) diff --git a/gridded/time.py b/gridded/time.py index 3af5172..d7f01eb 100644 --- a/gridded/time.py +++ b/gridded/time.py @@ -9,7 +9,23 @@ from gridded.utilities import get_dataset +class OutOfTimeRangeError(ValueError): + """ + Used for when data is asked for outside of the range of times supported by a Time object. + """ + pass + + +class TimeSeriesError(ValueError): + """ + Used for when data is asked for outside of the range of times supported by a Time object. + """ + pass + + class Time(object): + + # Used to make a singleton with the constant_time class method. _const_time = None def __init__(self, @@ -25,15 +41,20 @@ def __init__(self, Representation of a time axis. Provides interpolation alphas and indexing. :param time: Ascending list of times to use + :type time: netCDF4.Variable or Sequence of `datetime.datetime` + :param tz_offset: offset to compensate for time zone shifts + :type tz_offset: `datetime.timedelta` or float or integer hours + :param origin: shifts the time interval to begin at the time specified - :param displacement: displacement to apply to the time data. Allows shifting entire time interval into future or past - :type time: netCDF4.Variable or [] of datetime.datetime - :type tz_offset: datetime.timedelta - :type origin: datetime.timedelta - :type displacement: datetime.timedelta + :type origin: `datetime.datetime` + + :param displacement: displacement to apply to the time data. + Allows shifting entire time interval into future or past + :type displacement: `datetime.timedelta` ''' + # fixme -- this conversion should be done in the netcdf loading code. if isinstance(data, (nc4.Variable, nc4._netCDF4._Variable)): if (hasattr(nc4, 'num2pydate')): self.data = nc4.num2pydate(data[:], units=data.units) @@ -53,18 +74,19 @@ def __init__(self, self.filename = filename self.varname = varname -# if self.filename is None: -# self.filename = self.id + '_time.txt' - if tz_offset is not None: - self.data += tz_offset + try: + self.data += tz_offset + except TypeError: + self.data += timedelta(hours=tz_offset) if not self._timeseries_is_ascending(self.data): - raise ValueError("Time sequence is not ascending") + raise TimeSeriesError("Time sequence is not ascending") if self._has_duplicates(self.data): - raise ValueError("Time sequence has duplicate entries") + raise TimeSeriesError("Time sequence has duplicate entries") super(Time, self).__init__(*args, **kwargs) + @classmethod def from_netCDF(cls, filename=None, @@ -104,11 +126,18 @@ def from_netCDF(cls, varname=varname, tz_offset=tz_offset, **kwargs - ) + ) return time @classmethod def constant_time(cls): + """ + Returns a Time object that represents no change in time + + In practice, that's a Time object with a single datetime + """ + # this caches a single instance of a constant time object + # in the class, and always returns the same one (a singleton) if cls._const_time is None: cls._const_time = cls(np.array([datetime.now()])) return cls._const_time @@ -118,10 +147,12 @@ def data(self): return self._data @data.setter - def data(self, d): - if isinstance(d, self.__class__) or d.__class__ in self.__class__.__mro__: - d = d.data - self._data = d + def data(self, data): + # If it's passed in a Time object, it gets its data object + # Fixme: Why? this seems like more magic than necessary + if isinstance(data, self.__class__) or data.__class__ in self.__class__.__mro__: + data = data.data + self._data = np.asarray(data) @property def info(self): @@ -149,11 +180,14 @@ def __len__(self): return len(self.data) def __iter__(self): - return self.data.__iter__() + return iter(self.data) def __eq__(self, other): - r = self.data == other.data - return all(r) if hasattr(r, '__len__') else r + # r = self.data == other.data + # return all(r) if hasattr(r, '__len__') else r + if not isinstance(other, self.__class__): + return False + return np.array_equal(self.data, other.data) def __ne__(self, other): return not self.__eq__(other) @@ -183,31 +217,43 @@ def max_time(self): return self.data[-1] def get_time_array(self): - return self.data[:] + """ + returns a copy of the internal data array + """ + return self.data.copy() def time_in_bounds(self, time): ''' Checks if time provided is within the bounds represented by this object. :param time: time to be queried - :type time: datetime.datetime - :rtype: boolean + :type time: `datetime.datetime` + :rtype: bool ''' - return not time < self.min_time or time > self.max_time + return not ((time < self.min_time) or (time > self.max_time)) def valid_time(self, time): - if time < self.min_time or time > self.max_time: - raise ValueError('time specified ({0}) is not within the bounds of the time ({1} to {2})'.format( - time.strftime('%c'), self.min_time.strftime('%c'), self.max_time.strftime('%c'))) + """ + Raises a OutOfTimeRangeError if time is not within the bounds of the timeseries + """ + # if time < self.min_time or time > self.max_time: + if not self.time_in_bounds(time): + raise OutOfTimeRangeError(f'time specified ({time.strftime('%c')}) is not within the bounds of ' + f'({self.min_time.strftime('%c')} to {self.max_time.strftime('%c')})' + ) + def index_of(self, time, extrapolate=False): ''' Returns the index of the provided time with respect to the time intervals in the file. :param time: Time to be queried - :param extrapolate: - :type time: datetime.datetime - :type extrapolate: boolean + :type time: `datetime.datetime` + + :param extrapolate: whether to allow extrapolation: + i.e. will not raise if outside the bounds of the time data. + :type extrapolate: bool + :return: index of first time before specified time :rtype: integer ''' @@ -218,18 +264,26 @@ def index_of(self, time, extrapolate=False): index = 0 return index + def interp_alpha(self, time, extrapolate=False): ''' Returns interpolation alpha for the specified time + This is the weighting to give the index before + -- 1-alpha would be the weight to the index after. + + :param time: Time to be queried - :param extrapolate: - :type time: datetime.datetime - :type extrapolate: boolean + :type time: `datetime.datetime` + + :param extrapolate: if True, 0.0 (before) or 1.0 (after) is returned. + if False, a ValueError is raised if outside the time series. + :type extrapolate: bool + :return: interpolation alpha - :rtype: double (0 <= r <= 1) + :rtype: float (0 <= r <= 1) ''' - if not len(self.data) == 1 or not extrapolate: + if (not extrapolate) and (not len(self.data) == 1): self.valid_time(time) i0 = self.index_of(time, extrapolate) if i0 > len(self.data) - 1: diff --git a/gridded/utilities.py b/gridded/utilities.py index b0fe254..350d137 100644 --- a/gridded/utilities.py +++ b/gridded/utilities.py @@ -11,6 +11,7 @@ import numpy as np import netCDF4 as nc4 +import os must_have = ['dtype', 'shape', 'ndim', '__len__', '__getitem__', '__getattribute__'] @@ -62,10 +63,12 @@ def gen_celltree_mask_from_center_mask(center_mask, sl): input_mask = convert_mask_to_numpy_mask(center_mask) return input_mask[sl] +#1/16/2024 Jay Hennen: I'm 'privating' this function pending a review and likely rewrite. +# def regrid_variable(grid, o_var, location='node'): from gridded.variable import Variable from gridded.grids import Grid_S, Grid_U - from gridded.depth import S_Depth, Depth + from gridded.depth import S_Depth, L_Depth, DepthBase """ Takes a Variable or VectorVariable and interpolates the data onto grid. You may pass a location ('nodes', 'faces', 'edge1', 'edge2) and the @@ -102,9 +105,9 @@ def regrid_variable(grid, o_var, location='node'): the grid of the source variable {1}".format(grid, o_var)) n_depth = None if o_var.depth is not None: - if isinstance(o_var.depth, S_Depth): + if issubclass(o_var.depth.__class__, S_Depth): n_depth = _regrid_s_depth(grid, o_var.depth) - elif isinstance(o_var.depth, Depth): + elif isinstance(o_var.depth, DepthBase) or isinstance(o_var.depth, L_Depth): n_depth = o_var.depth else: raise NotImplementedError("Can only regrid sigma depths for now") @@ -125,10 +128,11 @@ def regrid_variable(grid, o_var, location='node'): pts[:, 0:2] = dest_points if o_var.time is not None: for t_idx, t in enumerate(o_var.time.data): - if n_depth is not None and isinstance(n_depth, S_Depth): - for lev_idx, lev_data in enumerate(o_var.depth.get_section(t)): + if n_depth is not None and issubclass(n_depth.__class__, S_Depth): + transect = o_var.depth.get_transect(o_var.grid.nodes.reshape(-1,2), t, data_shape=(len(n_depth),)).T + for lev_idx, lev_data in enumerate(transect): lev = Variable(name='level{0}'.format(lev_idx), - data=lev_data, + data=lev_data.reshape(o_var.grid.node_lon.shape), grid=o_var.grid) zs = lev.at(pts, t) pts[:, 2] = zs[:,0] @@ -151,13 +155,11 @@ def _regrid_s_depth(grid, o_depth): """ Creates a new S_Depth object from an existing one that works on a new grid. """ - from gridded.grids import Grid_S, Grid_U - from gridded.depth import S_Depth o_bathy = o_depth.bathymetry o_zeta = o_depth.zeta n_bathy = regrid_variable(grid, o_bathy) n_zeta = regrid_variable(grid, o_zeta) - n_depth = S_Depth(time=o_depth.time, + n_depth = o_depth.__class__(time=o_depth.time, grid=grid, bathymetry=n_bathy, zeta=n_zeta, @@ -360,6 +362,37 @@ def get_dataset(ncfile, dataset=None): else: return nc4.MFDataset(ncfile) +def parse_filename_dataset_args(filename=None, + dataset=None, + data_file=None, + grid_file=None,): + """A perinnial problem is that we need to be able to accept a variety of + filename/dataset arguments. This function will return a tuple of 'ds' and 'dg' + where ds is the dataset for the data and dg is the dataset for the grid. If such + datasets are from the same file, it will return the same object for both. + """ + if filename is not None: + try: + filename = os.fspath(filename) + except TypeError: + pass + data_file = grid_file = filename + ds = None + dg = None + if dataset is None: + if grid_file == data_file: + ds = dg = get_dataset(grid_file) + else: + ds = get_dataset(data_file) + dg = get_dataset(grid_file) + else: + if grid_file is not None: + dg = get_dataset(grid_file) + else: + dg = dataset + ds = dataset + + return ds, dg def get_writable_dataset(ncfile, format="netcdf4"): """ @@ -390,6 +423,27 @@ def get_dataset_attrs(ds): """ return {name: ds.getncattr(name) for name in ds.ncattrs()} +def search_netcdf_vars(cls=None, ds=None, dg=None): + """ + Given a class with a .default_names and .cf_names attributes, search a datafile + netCDF4.Dataset and a possible grid netCDF4.Dataset for variables that are sought by + the class + """ + vn_search = search_dataset_for_variables_by_varname(ds, cls.default_names) + ds_search = search_dataset_for_variables_by_longname(ds, cls.cf_names) + found_vars = merge_var_search_dicts(ds_search, vn_search) + if ds != dg: + dg_vn_search = search_dataset_for_variables_by_varname(dg, cls.default_names) + dg_ln_search = search_dataset_for_variables_by_longname(dg, cls.cf_names) + dg_search = merge_var_search_dicts(dg_ln_search, dg_vn_search) + found_vars = merge_var_search_dicts(found_vars, dg_search) + return found_vars + +def can_create_class(cls, ds=None, dg=None): + found_vars = search_netcdf_vars(cls, ds, dg) + # all variables must be found (no None values) + return not (None in found_vars.values()) + def varnames_merge(cls, inc_varnames=None): """ Helper function to support the `varnames` argument pattern. diff --git a/gridded/variable.py b/gridded/variable.py index 1e28018..4f79878 100644 --- a/gridded/variable.py +++ b/gridded/variable.py @@ -542,8 +542,6 @@ def at(self, Failure to provide point data in this format may cause unexpected behavior If you wish to provide point data using separate longitude and latitude arrays, use the `lons=` and `lats=` kwargs. - Note that if your Z is positive-up, self.depth.positive_down should be - set to False :type points: Nx3 array of double :param time: The time at which to query these points (T) @@ -1192,8 +1190,6 @@ def at(self, If you wish to provide point data using separate longitude and latitude arrays, use the `lons=` and `lats=` kwargs. - Note that if your Z is positive-up, self.depth.positive_down should be - set to False :type points: Nx3 array of double