diff --git a/HISTORY.rst b/HISTORY.rst index e1c1caa0..8df92388 100644 --- a/HISTORY.rst +++ b/HISTORY.rst @@ -4,6 +4,10 @@ History ======= +v0.3.1 (2022-12-14) +------------------- +Fix tarball and wheel for new release (PR 295 by `Filipe Fernandes`_). + v0.3.0 (2022-12-12) ------------------- Migrated master to main (PR 239 by `Mattia Almansi`_). Allow oceandataset to have dict-style access to @@ -27,4 +31,5 @@ Initial release published in the `Journal of Open Source Software`_. .. _`Ali Siddiqui`: https://github.com/asiddi24 .. _`Miguel Jimenez Urias`: https://github.com/Mikejmnez .. _`Wenrui Jiang`: https://github.com/MaceKuailv +.. _`Filipe Fernandes`: https://github.com/ocefpaf .. _`Journal of Open Source Software`: https://joss.theoj.org diff --git a/binder/environment.yml b/binder/environment.yml index 11b30f36..9e6f6f70 100644 --- a/binder/environment.yml +++ b/binder/environment.yml @@ -2,7 +2,7 @@ name: rise-environment channels: - conda-forge dependencies: -- python < 3.10 +- python - numpy - matplotlib - pandas @@ -12,7 +12,7 @@ dependencies: - distributed - bottleneck - netCDF4 -- xarray < 2022.6 +- xarray - cartopy < 0.20 - esmpy - intake-xarray diff --git a/ci/environment.yml b/ci/environment.yml index 8c49b163..52ccdc4d 100644 --- a/ci/environment.yml +++ b/ci/environment.yml @@ -11,7 +11,7 @@ dependencies: - intake-xarray - geopy - xesmf -- esmf < 8.4 +- esmf - xgcm < 0.7 - Ipython - tqdm diff --git a/docs/.DS_Store b/docs/.DS_Store index 2376f41c..8b5b011e 100644 Binary files a/docs/.DS_Store and b/docs/.DS_Store differ diff --git a/docs/api.rst b/docs/api.rst index 98c9b7ba..ab7dfbb1 100644 --- a/docs/api.rst +++ b/docs/api.rst @@ -217,3 +217,25 @@ Functions utils.densjmd95 utils.densmdjwf utils.Coriolis_parameter + +LLC-transformation +================== +.. autosummary:: + :toctree: generated/ + + llc_rearrange + + +Class +--------- +.. autosummary:: + :toctree: generated/ + + llc_rearrange.LLCtransformation + +.. Classmethod +.. --------- +.. .. autosummary:: +.. :toctree: generated/ + +.. llc_rearrange.LLCtransformation.arctic_crown diff --git a/oceanspy/compute.py b/oceanspy/compute.py index 8be29467..82057dd0 100644 --- a/oceanspy/compute.py +++ b/oceanspy/compute.py @@ -89,7 +89,6 @@ geographical_aligned_velocities=["U_zonal", "V_merid"], survey_aligned_velocities=["rot_ang_Vel", "tan_Vel", "ort_Vel"], missing_horizontal_spacing=["dxF", "dxV", "dyF", "dyU"], - missing_cs_sn=["CS", "SN"], ) @@ -2834,36 +2833,6 @@ def salt_budget(od): return _ospy.OceanDataset(ds).dataset -def missing_cs_sn(od): - """ - Compute missing grid orientaioon. - - Parameters - ---------- - od: OceanDataset - oceandataset used to compute - - Returns - ------- - od: OceanDataset - | CS: cosine of the orientation angle - | SN: sine of the orientation angle - """ - xc = _np.deg2rad(_np.array(od._ds.XC)) - yc = _np.deg2rad(_np.array(od._ds.YC)) - cs = _np.zeros_like(xc) - sn = _np.zeros_like(xc) - cs[0], sn[0] = _utils.find_cs_sn(yc[0], xc[0], yc[1], xc[1]) - cs[-1], sn[-1] = _utils.find_cs_sn(yc[-2], xc[-2], yc[-1], xc[-1]) - cs[1:-1], sn[1:-1] = _utils.find_cs_sn(yc[:-2], xc[:-2], yc[2:], xc[2:]) - od._ds["CS"] = od._ds["XC"] - od._ds["CS"].values = cs - - od._ds["SN"] = od._ds["XC"] - od._ds["SN"].values = sn - return od - - def missing_horizontal_spacing(od): """ Compute missing horizontal spacing. diff --git a/oceanspy/llc_rearrange.py b/oceanspy/llc_rearrange.py index 13b6a186..abac3418 100644 --- a/oceanspy/llc_rearrange.py +++ b/oceanspy/llc_rearrange.py @@ -1,9 +1,15 @@ +""" +OceanSpy functionality that transforms a dataset with LLC geometry characterized by +13 faces (or tiles), into one with simple geometry. +""" + import copy as _copy import reprlib import dask import numpy as _np import xarray as _xr +from xgcm import Grid from .utils import _rel_lon, _reset_range, get_maskH @@ -16,7 +22,7 @@ class LLCtransformation: - """A class containing the transformation of LLCgrids""" + """A class containing the transformation types of LLCgrids.""" def __init__( self, @@ -28,15 +34,16 @@ def __init__( faces=None, centered=False, chunks=None, - drop=False, ): self._ds = ds # xarray.DataSet self._varList = varList # variables names to be transformed + self._add_Hbdr = add_Hbdr self._XRange = XRange # lon range of data to retain - self.YRange = YRange # lat range of data to retain. + self._YRange = YRange # lat range of data to retain. self._chunks = chunks # dict. self._faces = faces # faces involved in transformation - self._centered = (centered,) + self._centered = centered + self._chunks = chunks @classmethod def arctic_crown( @@ -49,12 +56,70 @@ def arctic_crown( faces=None, centered=None, chunks=None, - drop=False, ): - """Transforms the dataset in which `face` appears as a dimension into - one without faces, with grids and variables sharing a common grid - orientation. + """This transformation splits the arctic cap (face=6) into four triangular + regions and combines all faces in a quasi lat-lon grid. The triangular + arctic regions form a crown atop faces {7, 10, 2, 5}. The final size of + the transformed dataset depends XRange, YRange or faces. + + Parameters + ---------- + dataset: xarray.Dataset + The multi-dimensional, in memory, array database. e.g., `oceandataset._ds`. + varList: 1D array_like, str, or None + List of variables (strings). + YRange: 1D array_like, scalar, or None + Y axis limits (e.g., latitudes). + If len(YRange)>2, max and min values are used. + XRange: 1D array_like, scalar, or None + X axis limits (e.g., longitudes). Can handle (periodic) discontinuity at + lon=180 deg E. + add_Hbdr: bool, scal + If scalar, add and subtract `add_Hbdr` to the the horizontal range. + of the horizontal ranges. + If True, automatically estimate add_Hbdr. + If False, add_Hbdr is set to zero. + faces: 1D array_like, scalar, or None + List of faces to be transformed. + If None, entire dataset is transformed. + When both [XRange, YRange] and faces are defined, [XRange, YRange] is used. + centered: str or bool. + If 'Atlantic' (default), the transformation creates a dataset in which the + Atlantic Ocean lies at the center of the domain. + If 'Pacific', the transformed data has a layout in which the Pacific Ocean + lies at the center of the domain. + This option is only relevant when transforming the entire dataset. + chunks: bool or dict. + If False (default) - chunking is automatic. + If dict, rechunks the dataset according to the spefications of the + dictionary. See xarray.chunk(). + + Returns + ------- + + ds: xarray.Dataset + face is no longer a dimension of the dataset. + + + Notes + ----- + This functionality is very similar to, takes on similar arguments and is used + internally by subsample.cutout when extracting cutout regions of datasets with + face as a dimension. + + + References + ---------- + https://docs.xarray.dev/en/stable/generated/xarray.Dataset.html + + https://docs.xarray.dev/en/stable/generated/xarray.Dataset.chunk.html + + + See Also + -------- + subsample.cutout """ + print("Warning: This is an experimental feature") if "face" not in ds.dims: raise ValueError("face does not appear as a dimension of the dataset") @@ -82,6 +147,11 @@ def arctic_crown( varList = list(varList) + # store original attributes + attrs = {} + for var in varList: + attrs = {var: ds[var].attrs, **attrs} + # if faces is None: faces = _np.arange(13) @@ -325,8 +395,7 @@ def arctic_crown( if type(DSFacet34) == int: DS = _reorder_ds(DS, dims_c, dims_g).persist() - if drop: - DS = _LLC_check_sizes(DS) + DS = _LLC_check_sizes(DS) # rechunk data. In the ECCO data this is done automatically if chunks: @@ -335,6 +404,14 @@ def arctic_crown( if XRange is not None and YRange is not None: # drop copy var = 'nYg' (line 101) DS = DS.drop_vars(_var_) + + DS = llc_local_to_lat_lon(DS) + + # restore original attrs if lost + for var in varList: + if var in DS.reset_coords().data_vars: + DS[var].attrs = attrs[var] + return DS @@ -392,6 +469,8 @@ def arct_connect(ds, varName, faces=None, masking=False, opt=False, ranges=None) if len(dims.X) + len(dims.Y) == 4: if len(dims.Y) == 3 and _varName not in metrics: fac = -1 + elif _varName == "CS": + fac = -1 arct = fac * ds[_varName].isel(**da_arg) Mask = mask2.isel(**mask_arg) if opt: @@ -493,6 +572,8 @@ def arct_connect(ds, varName, faces=None, masking=False, opt=False, ranges=None) if len(dims.X) + len(dims.Y) == 4: if len(dims.Y) == 1 and _varName not in metrics: fac = -1 + elif _varName == "SN": + fac = -1 da_arg = {"face": arc_cap, dims.X: xslice, dims.Y: yslice} mask_arg = {dims.X: xslice, dims.Y: yslice} arct = fac * ds[_varName].isel(**da_arg) @@ -521,6 +602,9 @@ def arct_connect(ds, varName, faces=None, masking=False, opt=False, ranges=None) def mates(ds): + """Defines, when needed, the variable pair and stores the name of the pair (mate) + variable as an attribute. This is needed to accurately rotate a vector field. + """ vars_mates = [ "ADVx_SLT", "ADVy_SLT", @@ -546,6 +630,8 @@ def mates(ds): "HFacS", "rAw", "rAs", + "CS", + "SN", ] for k in range(int(len(vars_mates) / 2)): nk = 2 * k @@ -556,7 +642,7 @@ def mates(ds): def rotate_vars(_ds): - """using the attribures `mates`, when this function is called it swaps the + """Using the attribures `mates`, when this function is called it swaps the variables names. This issue is only applicable to llc grid in which the grid topology makes it so that u on a rotated face transforms to `+- v` on a lat lon grid. @@ -575,7 +661,7 @@ def rotate_vars(_ds): def shift_dataset(_ds, dims_c, dims_g): - """shifts a dataset along a dimension, setting its first element to zero. Need + """Shifts a dataset along a dimension, setting its first element to zero. Need to provide the dimensions in the form of [center, corner] points. This rotation is only used in the horizontal, and so dims_c is either one of `i` or `j`, and dims_g is either one of `i_g` or `j_g`. The pair most correspond to the same @@ -605,7 +691,7 @@ def shift_dataset(_ds, dims_c, dims_g): def reverse_dataset(_ds, dims_c, dims_g, transpose=False): - """reverses the dataset along a dimension. Need to provide the dimensions in the + """Reverses the dataset along a dimension. Need to provide the dimensions in the form of [center, corner] points. This rotation is only used in the horizontal, and so dims_c is either one of `i` or `j`, and dims_g is either one of `i_g` or `j_g`. The pair most correspond to the same dimension.""" @@ -684,8 +770,9 @@ def rotate_dataset( def shift_list_ds(_DS, dims_c, dims_g, Ni, facet=1): - """given a list of n-datasets, each element of the list gets shifted along the - dimensions provided (dims_c and dims_g) so that there is no overlap between them. + """Given a list of n-datasets with matching dimensions, each element of the list + gets shifted along the dimensions provided (by dims_c and dims_g) so that there + is no overlap of values between them. """ _DS = _copy.deepcopy(_DS) fac = 1 @@ -725,7 +812,10 @@ def shift_list_ds(_DS, dims_c, dims_g, Ni, facet=1): def combine_list_ds(_DSlist): - """combines a list of n-datasets""" + """Combines a list of N-xarray.datasets along a dimension. Datasets must have + matching dimensions. See `xr.combine_first()` + + """ if len(_DSlist) == 0: _DSFacet = 0 # No dataset to combine. Return empty elif len(_DSlist) == 1: # a single face @@ -750,7 +840,7 @@ def combine_list_ds(_DSlist): def flip_v(_ds, co_list=metrics, dims=True, _len=3): - """reverses the sign of the vector fields by default along the corner coordinate + """Reverses the sign of the vector fields by default along the corner coordinate (Xp1 or Yp1). If dims is True, for each variable we infer the dimensions. Otherwise, dims is given @@ -763,6 +853,10 @@ def flip_v(_ds, co_list=metrics, dims=True, _len=3): if "mate" in _ds[_varName].attrs: if _varName not in co_list and len(_dims.X) == _len: _ds[_varName] = -_ds[_varName] + elif _varName == "SN": + _ds[_varName] = -_ds[_varName] + # elif _varName == "CS": + # _ds[_varName] = -_ds[_varName] return _ds @@ -976,7 +1070,9 @@ def _LLC_check_sizes(_DS): YG = _DS["YG"].dropna("Yp1", "all") y0 = int(YG["Yp1"][0]) y1 = int(YG["Yp1"][-1]) + 1 + _DS = _copy.deepcopy(_DS.isel(Yp1=slice(y0, y1))) + _DS = _copy.deepcopy(_DS.isel(Y=slice(y0, y1 - 1))) DIMS = [dim for dim in _DS["XC"].dims] dims_c = Dims(DIMS[::-1]) @@ -1012,7 +1108,7 @@ def _LLC_check_sizes(_DS): def _reorder_ds(_ds, dims_c, dims_g): - """Only needed when Pacific -centered data. Corrects the ordering + """Only needed when Pacific-centered data. Corrects the ordering of y-dim and transposes the data, all lazily.""" _DS = _copy.deepcopy(_ds) @@ -1035,7 +1131,68 @@ def _reorder_ds(_ds, dims_c, dims_g): return _DS +def llc_local_to_lat_lon(ds, co_list=metrics): + """ + Takes all vector fields and rotates them to orient them along geographical + coordinates. + """ + _ds = mates(_copy.deepcopy(ds)) + + grid_coords = { + "Y": {"center": "Y", "outer": "Yp1"}, + "X": {"center": "X", "outer": "Xp1"}, + "Z": {"center": "Z", "outer": "Zp1", "right": "Zu", "left": "Zl"}, + "time": {"center": "time_midp", "left": "time"}, + } + + # create grid object to interpolate + grid = Grid(_ds, coords=grid_coords, periodic=[]) + + CS = _ds["CS"] # cosine of angle between logical and geo axis. At tracer points + SN = _ds["SN"] # sine of angle between logical and geo axis. At tracer points + + CSU = grid.interp(CS, axis="X", boundary="extend") # cos at u-point + CSV = grid.interp(CS, axis="Y", boundary="extend") # cos at v-point + + SNU = grid.interp(SN, axis="X", boundary="extend") # sin at u-point + SNV = grid.interp(SN, axis="Y", boundary="extend") # sin at v-point + + data_vars = [var for var in _ds.data_vars if len(ds[var].dims) > 1] + + for var in data_vars: + DIMS = [dim for dim in _ds[var].dims] + dims = Dims(DIMS[::-1]) + if len(dims.X) + len(dims.Y) == 4: # vector field (metric) + if len(dims.Y) == 1 and var not in co_list: # u vector + _da = _copy.deepcopy(_ds[var]) + if "mate" in _ds[var].attrs: + mate = _ds[var].mate + _ds = _ds.drop_vars([var]) + VU = grid.interp( + grid.interp(_ds[mate], axis="Y", boundary="extend"), + axis="X", + boundary="extend", + ) + _ds[var] = _da * CSU - VU * SNU + elif len(dims.Y) == 3 and var not in co_list: # v vector + _da = _copy.deepcopy(_ds[var]) + if "mate" in _ds[var].attrs: + mate = _ds[var].mate + _ds = _ds.drop_vars([var]) + UV = grid.interp( + grid.interp(_ds[mate], axis="X", boundary="extend"), + axis="Y", + boundary="extend", + ) + _ds[var] = UV * SNV + _da * CSV + + return _ds + + class Dims: + """Creates a shortcut for dimension`s names associated with an arbitrary + variable.""" + axes = "XYZT" # shortcut axis names def __init__(self, vars): diff --git a/oceanspy/subsample.py b/oceanspy/subsample.py index 8ce601cb..20d7622d 100644 --- a/oceanspy/subsample.py +++ b/oceanspy/subsample.py @@ -353,7 +353,6 @@ def cutout( "XRange": XRange, "YRange": YRange, "centered": centered, - "drop": True, # required to calculate U-V grid points "chunks": chunks, } dsnew = _llc_trans.arctic_crown(**arg) diff --git a/oceanspy/tests/test_llc_rearrange.py b/oceanspy/tests/test_llc_rearrange.py index be15985b..bde4b0fa 100644 --- a/oceanspy/tests/test_llc_rearrange.py +++ b/oceanspy/tests/test_llc_rearrange.py @@ -2262,18 +2262,18 @@ def test_arc_connect( @pytest.mark.parametrize( "od, faces, varList, XRange, YRange, X0, X1, Y0, Y1", [ - (od, None, varList, None, None, 0, 359, 0, 314), - (od, [2, 5, 6, 7, 10], varList, None, None, 0, 359, 0, 134), - (od, [2, 5, 7, 10], varList, None, None, 0, 359, 0, 89), - (od, [1, 4, 8, 11], varList, None, None, 0, 359, 0, 89), - (od, [0, 3, 9, 12], varList, None, None, 0, 359, 0, 89), - (od, [6, 7, 8, 9], varList, None, None, 0, 89, 0, 314), - (od, [6, 10, 11, 12], varList, None, None, 0, 89, 0, 314), - (od, [0, 1, 2, 6], varList, None, None, 0, 89, 0, 314), - (od, [3, 4, 5, 6], varList, None, None, 0, 89, 0, 314), - (od, [2, 6, 10], varList, None, None, 0, 179, 0, 134), - (od, [6, 7, 10], varList, None, None, 0, 179, 0, 134), - (od, [6, 7, 8, 10, 11], varList, None, None, 0, 179, 0, 224), + (od, None, varList, None, None, 0, 358, 0, 313), + (od, [2, 5, 6, 7, 10], varList, None, None, 0, 358, 0, 133), + (od, [2, 5, 7, 10], varList, None, None, 0, 358, 0, 89), + (od, [1, 4, 8, 11], varList, None, None, 0, 358, 0, 89), + (od, [0, 3, 9, 12], varList, None, None, 0, 358, 0, 89), + (od, [6, 7, 8, 9], varList, None, None, 0, 88, 0, 312), + (od, [6, 10, 11, 12], varList, None, None, 0, 88, 0, 312), + (od, [0, 1, 2, 6], varList, None, None, 0, 88, 0, 313), + (od, [3, 4, 5, 6], varList, None, None, 0, 88, 0, 313), + (od, [2, 6, 10], varList, None, None, 0, 178, 0, 133), + (od, [6, 7, 10], varList, None, None, 0, 178, 0, 132), + (od, [6, 7, 8, 10, 11], varList, None, None, 0, 178, 0, 222), ( od, [4, 5, 6, 7, 8, 10, 11], @@ -2281,54 +2281,54 @@ def test_arc_connect( None, None, 0, - 269, + 268, 0, - 224, + 223, ), - (od, [7, 8, 10, 11], varList, None, None, 0, 179, 0, 179), - (od, [1, 2, 10, 11], varList, None, None, 0, 179, 0, 179), - (od, [8, 9, 11, 12], varList, None, None, 0, 179, 0, 179), - (od, [0, 1, 11, 12], varList, None, None, 0, 179, 0, 179), - (od, [0, 1, 3, 4], varList, None, None, 0, 179, 0, 179), - (od, [9, 12], varList, None, None, 0, 179, 0, 89), - (od, [0, 12], varList, None, None, 0, 179, 0, 89), - (od, [0, 3], varList, None, None, 0, 179, 0, 89), - (od, [0, 9, 12], varList, None, None, 0, 269, 0, 89), - (od, [0, 3, 12], varList, None, None, 0, 269, 0, 89), - (od, [0], varList[0], None, None, 0, 89, 0, 89), - (od, [1], varList[0], None, None, 0, 89, 0, 89), - (od, [2], varList[0], None, None, 0, 89, 0, 89), - (od, [3], varList[0], None, None, 0, 89, 0, 89), - (od, [4], varList[0], None, None, 0, 89, 0, 89), - (od, [5], varList[0], None, None, 0, 89, 0, 89), - (od, [7], varList[0], None, None, 0, 89, 0, 89), - (od, [8], varList[0], None, None, 0, 89, 0, 89), - (od, [9], varList[0], None, None, 0, 89, 0, 89), - (od, [10], varList[0], None, None, 0, 89, 0, 89), - (od, [11], varList[0], None, None, 0, 89, 0, 89), - (od, [12], varList[0], None, None, 0, 89, 0, 89), - (od, None, varList, _XRanges_[0], _YRanges_[0], 0, 44, 0, 24), - (od, None, varList, _XRanges_[1], _YRanges_[1], 0, 68, 0, 8), - (od, None, varList, _XRanges_[2], _YRanges_[2], 0, 67, 0, 6), - (od, None, varList, _XRanges_[3], _YRanges_[3], 0, 66, 0, 7), - (od, None, varList, _XRanges_[4], _YRanges_[4], 0, 71, 0, 72), - (od, None, varList, _XRanges_[5], _YRanges_[5], 0, 21, 0, 179), - (od, None, varList, _XRanges_[6], _YRanges_[6], 0, 39, 0, 55), - (od, None, varList, _XRanges_[7], _YRanges_[7], 0, 72, 0, 6), - (od, None, varList, _XRanges_[8], _YRanges_[8], 0, 112, 0, 8), - (od, None, varList, _XRanges_[9], _YRanges_[9], 0, 151, 0, 6), - (od, None, varList, _XRanges_[10], _YRanges_[10], 0, 109, 0, 5), - (od, None, varList, _XRanges_[11], _YRanges_[11], 0, 12, 0, 56), - (od, None, varList, _XRanges_[12], _YRanges_[12], 0, 16, 0, 84), - (od, None, varList, _XRanges_[13], _YRanges_[13], 0, 13, 0, 127), - (od, None, varList, _XRanges_[14], _YRanges_[14], 0, 7, 0, 191), - (od, None, varList, _XRanges_[15], _YRanges_[15], 0, 87, 0, 9), - (od, None, varList, _XRanges_[16], _YRanges_[16], 0, 6, 0, 68), - (od, None, varList, _XRanges_[17], _YRanges_[17], 0, 20, 0, 118), - (od, None, varList, [-31, 25], [58, 85.2], 0, 60, 0, 71), - (od, None, varList, [-120, -60], [58, 85.2], 0, 64, 0, 70), - (od, None, varList, [160, -150], [58, 85.2], 0, 54, 0, 70), - (od, None, varList, [60, 130], [58, 85.2], 0, 74, 0, 70), + (od, [7, 8, 10, 11], varList, None, None, 0, 178, 0, 178), + (od, [1, 2, 10, 11], varList, None, None, 0, 178, 0, 179), + (od, [8, 9, 11, 12], varList, None, None, 0, 178, 0, 178), + (od, [0, 1, 11, 12], varList, None, None, 0, 178, 0, 179), + (od, [0, 1, 3, 4], varList, None, None, 0, 178, 0, 178), + (od, [9, 12], varList, None, None, 0, 178, 0, 88), + (od, [0, 12], varList, None, None, 0, 178, 0, 89), # 19 + (od, [0, 3], varList, None, None, 0, 178, 0, 88), + (od, [0, 9, 12], varList, None, None, 0, 268, 0, 89), + (od, [0, 3, 12], varList, None, None, 0, 268, 0, 89), + (od, [0], varList[0], None, None, 0, 88, 0, 88), + (od, [1], varList[0], None, None, 0, 88, 0, 88), + (od, [2], varList[0], None, None, 0, 88, 0, 88), + (od, [3], varList[0], None, None, 0, 88, 0, 88), + (od, [4], varList[0], None, None, 0, 88, 0, 88), + (od, [5], varList[0], None, None, 0, 88, 0, 88), + (od, [7], varList[0], None, None, 0, 88, 0, 88), + (od, [8], varList[0], None, None, 0, 88, 0, 88), # 30 + (od, [9], varList[0], None, None, 0, 88, 0, 88), + (od, [10], varList[0], None, None, 0, 88, 0, 88), + (od, [11], varList[0], None, None, 0, 88, 0, 88), + (od, [12], varList[0], None, None, 0, 88, 0, 88), + (od, None, varList, _XRanges_[0], _YRanges_[0], 0, 43, 0, 23), + (od, None, varList, _XRanges_[1], _YRanges_[1], 0, 67, 0, 7), + (od, None, varList, _XRanges_[2], _YRanges_[2], 0, 66, 0, 5), + (od, None, varList, _XRanges_[3], _YRanges_[3], 0, 65, 0, 6), + (od, None, varList, _XRanges_[4], _YRanges_[4], 0, 70, 0, 71), + (od, None, varList, _XRanges_[5], _YRanges_[5], 0, 20, 0, 178), # 40 + (od, None, varList, _XRanges_[6], _YRanges_[6], 0, 38, 0, 54), + (od, None, varList, _XRanges_[7], _YRanges_[7], 0, 71, 0, 5), + (od, None, varList, _XRanges_[8], _YRanges_[8], 0, 111, 0, 7), + (od, None, varList, _XRanges_[9], _YRanges_[9], 0, 150, 0, 5), + (od, None, varList, _XRanges_[10], _YRanges_[10], 0, 108, 0, 4), + (od, None, varList, _XRanges_[11], _YRanges_[11], 0, 11, 0, 55), + (od, None, varList, _XRanges_[12], _YRanges_[12], 0, 15, 0, 83), + (od, None, varList, _XRanges_[13], _YRanges_[13], 0, 12, 0, 126), + (od, None, varList, _XRanges_[14], _YRanges_[14], 0, 6, 0, 190), + (od, None, varList, _XRanges_[15], _YRanges_[15], 0, 86, 0, 8), # 50 + (od, None, varList, _XRanges_[16], _YRanges_[16], 0, 5, 0, 67), + (od, None, varList, _XRanges_[17], _YRanges_[17], 0, 19, 0, 117), + (od, None, varList, [-31, 25], [58, 85.2], 0, 59, 0, 70), + (od, None, varList, [-120, -60], [58, 85.2], 0, 63, 0, 69), + (od, None, varList, [160, -150], [58, 85.2], 0, 53, 0, 69), + (od, None, varList, [60, 130], [58, 85.2], 0, 73, 0, 69), ], ) def test_transformation(od, faces, varList, XRange, YRange, X0, X1, Y0, Y1): diff --git a/oceanspy/tests/test_utils.py b/oceanspy/tests/test_utils.py index ca5877d4..f913f19f 100644 --- a/oceanspy/tests/test_utils.py +++ b/oceanspy/tests/test_utils.py @@ -18,6 +18,7 @@ def test_RNone(): def test_error_viewer_to_range(): with pytest.raises(TypeError): + viewer_to_range("does not eval to a list") viewer_to_range(0) viewer_to_range(["not from viewer"]) viewer_to_range([{"type": "other"}]) @@ -41,6 +42,7 @@ def test_error_path(): coords1 = [[[0, 0], [1, 1], [2, 2], [3, 3], [4, 4], [5, 5]]] coords2 = [[[5, 0], [4, 1], [3, 2], [2, 3], [1, 4], [0, 5]]] coords3 = [[[0, 6], [0, 7], [0, 8], [0, 9], [0, 10], [0, 11]]] +coords4 = '[{"type":"Point","coordinates":[-169.23960833202577,22.865677261831266]}]' @pytest.mark.parametrize( @@ -52,10 +54,14 @@ def test_error_path(): (coords1, "LineString", [[0, 0]], [[1, 1]]), (coords2, "LineString", [[5, 0]], [[4, 1]]), (coords3, "LineString", [[0, 6]], [[0, 7]]), + (coords4, "Point", [-169.23960833202577], [22.865677261831266]), ], ) def test_viewer_to_range(coords, types, lon, lat): - p = [{"type": types, "coordinates": list(coords)}] + if type(coords) == list: + p = [{"type": types, "coordinates": list(coords)}] + elif type(coords) == str: + p = coords x, y = viewer_to_range(p) assert x == lon assert y == lat diff --git a/oceanspy/utils.py b/oceanspy/utils.py index 8c7c73f3..b2deafe0 100644 --- a/oceanspy/utils.py +++ b/oceanspy/utils.py @@ -2,6 +2,7 @@ OceanSpy utilities that don't need OceanDataset objects. """ +import ast as _ast import copy as _copy import numpy as _np @@ -37,6 +38,7 @@ def viewer_to_range(p): Takes the output from the poseidon viewer `p` and returns the coordinate trajectories in X and Y in a way that is compatible with oceanspy.subsample functions. + Parameters ---------- p: list. @@ -46,6 +48,12 @@ def viewer_to_range(p): lon: list. lat: list. """ + + if type(p) == str: + if p[0] == "[" and p[-1] == "]": + p = _ast.literal_eval(p) # turn string into list + else: + TypeError("not a type extracted by poseidon viewer") _check_instance({"p": p}, {"p": ["list"]}) _check_instance({"p[0]": p[0]}, {"p[0]": ["dict"]}) _check_instance({"type": p[0]["type"]}, {"type": "str"}) @@ -61,15 +69,21 @@ def viewer_to_range(p): if p_type == "Polygon": coords = p[0]["coordinates"][0] - elif p_type in ["LineString", "Point"]: + elif p_type == "Point": + coords = p[0]["coordinates"] + elif p_type == "LineString": coords = p[0]["coordinates"] lon = [] lat = [] - for i in range(len(coords)): - lon.append(coords[i][0]) - lat.append(coords[i][1]) + if p_type == "Point": + lon.append(coords[0]) + lat.append(coords[1]) + else: + for i in range(len(coords)): + lon.append(coords[i][0]) + lat.append(coords[i][1]) return lon, lat @@ -826,33 +840,6 @@ def to_180(x): return x + (-1) * (x // 180) * 360 -def local_to_latlon(u, v, cs, sn): - """ - convert local vector to north-east - - Parameters - ---------- - u: float, numpy.array - the x component of the vector - v: float, numpy.array - the y component of the vector - cs: float, numpy.array - the cosine of the angle between grid and compass - sn: float, numpy.array - the sine of the angle between grid and compass - - Returns - ------- - uu: float, numpy.array - the eastward component of the vector - vv: float, numpy.array - the northward component of the vector - """ - uu = u * cs - v * sn - vv = u * sn + v * cs - return uu, vv - - def get_combination(lst, select): """ Iteratively find all the combination that @@ -883,42 +870,3 @@ def get_combination(lst, select): # print(sub_lst) the_lst += sub_lst return the_lst - - -def find_cs_sn(thetaA, phiA, thetaB, phiB): - """ - theta is the angle - between the meridian crossing point A - and the geodesic connecting A and B - - this function return cos and sin of theta - - Parameters - ---------- - thetaA: float,numpy.array - the latitude of the vertex of the angle - phiA: float,numpy.array - the longitude of the vertex of the angle - thetaB: float,numpy.array - the latitude of the point B - phiB: float,numpy.array - the longitude of the point B - - Returns - ------- - CS: float,numpy.array - the cosine of the angle - SN: float, numpy.array - the sine of the angle - """ - # O being north pole - AO = _np.pi / 2 - thetaA - BO = _np.pi / 2 - thetaB - dphi = phiB - phiA - # Spherical law of cosine on AOB - cos_AB = _np.cos(BO) * _np.cos(AO) + _np.sin(BO) * _np.sin(AO) * _np.cos(dphi) - sin_AB = _np.sqrt(1 - cos_AB**2) - # spherical law of sine on triangle AOB - SN = _np.sin(BO) * _np.sin(dphi) / sin_AB - CS = _np.sign(thetaB - thetaA) * _np.sqrt(1 - SN**2) - return CS, SN