Commit 19080530 authored by PIOLLE's avatar PIOLLE
Browse files
parents 335bdc30 79655a61
...@@ -1034,8 +1034,8 @@ class Dataset(ABC): ...@@ -1034,8 +1034,8 @@ class Dataset(ABC):
# probably not all of them... # probably not all of them...
warnings.warn( warnings.warn(
"Entering a special case where not xarray variable " "Entering a special case where not xarray variable "
"assignment was possible, probably because of fill " "assignment was possible, probably because of fill or duplicate"
"values in a coordinate. Field: {}".format(field.name)) " values in a coordinate. Field: {}".format(field.name))
self._std_dataset[field.name] = xr.DataArray( self._std_dataset[field.name] = xr.DataArray(
coords=[(d, self._std_dataset.coords[d]) for d in field.dims], coords=[(d, self._std_dataset.coords[d]) for d in field.dims],
data=field.get_values(), data=field.get_values(),
...@@ -2171,6 +2171,8 @@ class Dataset(ABC): ...@@ -2171,6 +2171,8 @@ class Dataset(ABC):
return attrval.strftime(cls._get_time_format()) return attrval.strftime(cls._get_time_format())
elif attrval is None: elif attrval is None:
return '' return ''
elif isinstance(attrval, shapely.geometry.base.BaseGeometry):
return attrval.wkt
return attrval return attrval
def get_collection_id(self) -> str: def get_collection_id(self) -> str:
......
...@@ -649,6 +649,10 @@ class Field(object): ...@@ -649,6 +649,10 @@ class Field(object):
pad with fill values the ``subset`` array extracted from ``array`` pad with fill values the ``subset`` array extracted from ``array``
where ``index`` is beyond the limits of ``array``. where ``index`` is beyond the limits of ``array``.
""" """
if len(array.shape) == 0:
# dimensionless field
return subset.to_masked_array()
pad_edges = [] pad_edges = []
for dim in list(array.dims): for dim in list(array.dims):
if dim in index: if dim in index:
......
...@@ -56,7 +56,8 @@ class CRACollection(Feature): ...@@ -56,7 +56,8 @@ class CRACollection(Feature):
@property @property
def _feature_geodimnames(self): def _feature_geodimnames(self):
return self.feature_class._instance_dimname, 'z', return (self.feature_class._instance_dimname,) + \
self.feature_class._feature_geodimnames
def get_geocoord_dimnames( def get_geocoord_dimnames(
self, fieldname: str, self, fieldname: str,
......
...@@ -206,13 +206,24 @@ class Feature(Dataset, ABC): ...@@ -206,13 +206,24 @@ class Feature(Dataset, ABC):
return dim_validity return dim_validity
def append(self, feature, prefix, fields=None): def append(self,
feature,
prefix: str = None,
add_coords: bool = False,
as_new_dims: bool = False,
fields: str = None):
"""Append the fields from another feature """Append the fields from another feature
It is assumed the two features do not share any dimension. The appended The fields can share the same coordinates and dimensions (if the added
feature (child) is considered as ancillary fields and looses its feature is an overlay of the current feature). In such case the
model properties. The feature model is the one of the receiving shared coordinates of the added feature are not copied, unless
(parent) feature. `add_coords` is set to True (in which case they are first prefixed
using `prefix`).
If the dimensions of the added field have unrelated to the current
feature (even if they have similar names), `as_new_dims` must be set
to True. The dimensions (and corresponding coordinates) are prefixed
with `prefix` and added to with the fields to the current feature.
Args: Args:
feature (AbstractFeature derived class): the feature to append. feature (AbstractFeature derived class): the feature to append.
...@@ -221,49 +232,48 @@ class Feature(Dataset, ABC): ...@@ -221,49 +232,48 @@ class Feature(Dataset, ABC):
feature (to avoid conflicts with the existing fields of the feature (to avoid conflicts with the existing fields of the
current feature). current feature).
add_coords: if True, add the feature coordinates as variables. If
false, they are not added.
as_new_dims: rename the field and coordinate dimensions with prefix
fields (list of str): a list of fieldnames specifying which fields fields (list of str): a list of fieldnames specifying which fields
are to be appended. By default, all fields of the feature are are to be appended. By default, all fields of the feature are
appended. appended.
""" """
child_fields = {} if as_new_dims and prefix is None:
# create new fields raise ValueError(
for fieldname in feature.get_geolocation_fields(): "a valid prefix must be provided to rename the dimensions")
newname = ''.join([prefix, fieldname])
if newname in self.fieldnames: if add_coords and prefix is None:
logging.debug("Field already existing: {}".format(newname)) raise ValueError(
continue "a valid prefix must be provided to add the coordinates")
childfield = feature.get_geolocation_field(fieldname)
if feature.get_geolocation_field(fieldname) is not None: added_fields = fields
child_fields[newname] = childfield.clone() if added_fields is None:
# additional _ ensures no field and dimensions have the same added_fields = feature.fieldnames
# name if add_coords:
dims = OrderedDict([('_'.join([prefix, k]), v) added_fields.extend(feature.coordnames)
for k, v in childfield.dimensions.items()])
child_fields[newname].dimensions = dims
child_fields[newname].variable.shortname = newname
child_fields[newname].set_values(childfield.get_values())
fieldlist = fields
if fieldlist is None:
fieldlist = feature.fieldnames
for fieldname in fieldlist:
newname = prefix + fieldname
childfield = feature.get_field(fieldname)
child_fields[newname] = childfield.clone()
# additional _ ensures no field and dimensions have the same
# name
dims = OrderedDict([('_'.join([prefix, k]), v)
for k, v in childfield.dimensions.items()])
child_fields[newname].dimensions = dims
child_fields[newname].variable.shortname = newname
child_fields[newname].set_values(
feature.get_values(fieldname)
)
# append fields # append fields
for fieldname, field in child_fields.items(): for fieldname in added_fields:
if fieldname in feature.fieldnames:
field = feature.get_field(fieldname).clone()
elif fieldname in feature.coords:
field = feature.get_coord(fieldname).clone()
else:
raise ValueError("field {} not found in feature".format(
fieldname))
if prefix is not None:
field.rename('{}{}'.format(prefix, field.name))
if as_new_dims:
field._array = field._array.rename(
{_: '{}{}'.format(prefix, _) for _ in field.dims})
self.add_field(field) self.add_field(field)
def get_values(self, *args, **kwargs): def get_values(self, *args, **kwargs):
""" """
See: See:
......
...@@ -81,6 +81,8 @@ class Grid(Feature): ...@@ -81,6 +81,8 @@ class Grid(Feature):
values of the existing attributes with the same name in the file). values of the existing attributes with the same name in the file).
Use with caution in this context. Use with caution in this context.
""" """
_feature_geodimnames = 'y', 'x',
def __init__(self, def __init__(self,
*args, *args,
projection: 'Projection' = Projection(), projection: 'Projection' = Projection(),
...@@ -103,10 +105,6 @@ class Grid(Feature): ...@@ -103,10 +105,6 @@ class Grid(Feature):
self._std_dataset = self._std_dataset.squeeze(dim='time') self._std_dataset = self._std_dataset.squeeze(dim='time')
self._std_dataset.coords['time'] = ctime self._std_dataset.coords['time'] = ctime
@property
def _feature_geodimnames(self) -> Tuple[str, ...]:
return 'y', 'x'
def get_geocoord_dimnames( def get_geocoord_dimnames(
self, fieldname: str, self, fieldname: str,
values: 'xr.DataArray') -> Tuple[str, ...]: values: 'xr.DataArray') -> Tuple[str, ...]:
......
...@@ -52,6 +52,8 @@ class GridTimeSeries(Feature): ...@@ -52,6 +52,8 @@ class GridTimeSeries(Feature):
values of the existing attributes with the same name in the file). values of the existing attributes with the same name in the file).
Use with caution in this context. Use with caution in this context.
""" """
_feature_geodimnames = 'time', 'y', 'x',
def __init__(self, def __init__(self,
*args, *args,
projection=Projection(), projection=Projection(),
...@@ -70,10 +72,6 @@ class GridTimeSeries(Feature): ...@@ -70,10 +72,6 @@ class GridTimeSeries(Feature):
"""Returns the iterator""" """Returns the iterator"""
return GridIterator(self) return GridIterator(self)
@property
def _feature_geodimnames(self) -> Tuple[str, ...]:
return 'time', 'y', 'x',
@property @property
def _grid_class(self): def _grid_class(self):
return Grid return Grid
......
...@@ -25,6 +25,8 @@ class Image(Feature): ...@@ -25,6 +25,8 @@ class Image(Feature):
that we don't need a time value for each pixel or even scan line. A single that we don't need a time value for each pixel or even scan line. A single
time value is enough to describe it, as if it was an instant snapshot. time value is enough to describe it, as if it was an instant snapshot.
""" """
_feature_geodimnames = 'row', 'cell',
def __init__(self, *args, **kwargs): def __init__(self, *args, **kwargs):
# create feature # create feature
super(Image, self).__init__( super(Image, self).__init__(
...@@ -38,10 +40,6 @@ class Image(Feature): ...@@ -38,10 +40,6 @@ class Image(Feature):
self._std_dataset = self._std_dataset.squeeze(dim='time') self._std_dataset = self._std_dataset.squeeze(dim='time')
self._std_dataset.coords['time'] = ctime self._std_dataset.coords['time'] = ctime
@property
def _feature_geodimnames(self) -> Tuple[str, ...]:
return 'row', 'cell',
def get_geocoord_dimnames( def get_geocoord_dimnames(
self, fieldname: str, self, fieldname: str,
values: 'xr.DataArray') -> Tuple[str, ...]: values: 'xr.DataArray') -> Tuple[str, ...]:
......
...@@ -54,13 +54,15 @@ class IMDCollection(Feature): ...@@ -54,13 +54,15 @@ class IMDCollection(Feature):
@property @property
def _feature_geodimnames(self): def _feature_geodimnames(self):
return self.feature_class._instance_dimname, 'z', return (self.feature_class._instance_dimname,) + \
self.feature_class._feature_geodimnames
def get_geocoord_dimnames( def get_geocoord_dimnames(
self, fieldname: str, self, fieldname: str,
values: 'xr.DataArray') -> Tuple[str, ...]: values: 'xr.DataArray') -> Tuple[str, ...]:
if fieldname in ['depth', 'height']: if fieldname in ['depth', 'height']:
return self.feature_class._instance_dimname, 'z', return ((self.feature_class._instance_dimname,) +
self.feature_class._feature_geodimnames)
else: else:
return self.feature_class._instance_dimname return self.feature_class._instance_dimname
......
...@@ -54,12 +54,13 @@ class OMDCollection(Feature): ...@@ -54,12 +54,13 @@ class OMDCollection(Feature):
@property @property
def _feature_geodimnames(self): def _feature_geodimnames(self):
return self.feature_class._instance_dimname, 'z', return (self.feature_class._instance_dimname,) + \
self.feature_class._feature_geodimnames
def get_geocoord_dimnames( def get_geocoord_dimnames(
self, fieldname: str, self, fieldname: str,
values: 'xr.DataArray') -> Tuple[str, ...]: values: 'xr.DataArray') -> Tuple[str, ...]:
if fieldname in ['depth', 'height']: if fieldname in ['depth', 'height']:
return 'z', return self.feature_class._feature_geodimnames
else: else:
return self.feature_class._instance_dimname return self.feature_class._instance_dimname
...@@ -23,15 +23,12 @@ class Point(Feature): ...@@ -23,15 +23,12 @@ class Point(Feature):
dimension: ``obs``. dimension: ``obs``.
""" """
_instance_dimname = 'obs' _instance_dimname = 'obs'
_feature_geodimnames = ()
def __init__(self, *args, **kwargs): def __init__(self, *args, **kwargs):
super(Point, self).__init__(*args, **kwargs) super(Point, self).__init__(*args, **kwargs)
self._std_dataset.attrs['featureType'] = 'point' self._std_dataset.attrs['featureType'] = 'point'
@property
def _feature_geodimnames(self):
return 'obs',
def get_geocoord_dimnames(self, *args, **kwargs): def get_geocoord_dimnames(self, *args, **kwargs):
return 'obs', return ()
...@@ -26,6 +26,7 @@ class Profile(Feature): ...@@ -26,6 +26,7 @@ class Profile(Feature):
dimension: ``profile``. dimension: ``profile``.
""" """
_instance_dimname = 'profile' _instance_dimname = 'profile'
_feature_geodimnames = 'z',
def __init__(self, *args, **kwargs): def __init__(self, *args, **kwargs):
super(Feature, self).__init__(*args, **kwargs) super(Feature, self).__init__(*args, **kwargs)
...@@ -39,10 +40,6 @@ class Profile(Feature): ...@@ -39,10 +40,6 @@ class Profile(Feature):
for _, v in self._std_dataset.data_vars.items(): for _, v in self._std_dataset.data_vars.items():
v.attrs['coordinates'] = 'time lat lon z' v.attrs['coordinates'] = 'time lat lon z'
@property
def _feature_geodimnames(self):
return 'z',
def get_geocoord_dimnames( def get_geocoord_dimnames(
self, fieldname: str, self, fieldname: str,
values: 'xr.DataArray') -> Tuple[str, ...]: values: 'xr.DataArray') -> Tuple[str, ...]:
......
...@@ -20,6 +20,8 @@ class Swath(Feature): ...@@ -20,6 +20,8 @@ class Swath(Feature):
Feature class for representing a swath, a two-dimensional irregular grid Feature class for representing a swath, a two-dimensional irregular grid
along the satellite track. along the satellite track.
""" """
_feature_geodimnames = 'row', 'cell',
def __init__(self, *args, **kwargs): def __init__(self, *args, **kwargs):
# create feature # create feature
super(Swath, self).__init__( super(Swath, self).__init__(
...@@ -27,10 +29,6 @@ class Swath(Feature): ...@@ -27,10 +29,6 @@ class Swath(Feature):
**kwargs **kwargs
) )
@property
def _feature_geodimnames(self):
return 'row', 'cell',
def get_geocoord_dimnames(self, fieldname, shape=None): def get_geocoord_dimnames(self, fieldname, shape=None):
if fieldname == 'depth': if fieldname == 'depth':
return 'depth', return 'depth',
......
...@@ -24,6 +24,7 @@ class TimeSeries(Feature): ...@@ -24,6 +24,7 @@ class TimeSeries(Feature):
(except ``time`` for T axis coordinate). (except ``time`` for T axis coordinate).
""" """
_instance_dimname = 'station' _instance_dimname = 'station'
_feature_geodimnames = 'time',
def __init__(self, *args, **kwargs): def __init__(self, *args, **kwargs):
super(TimeSeries, self).__init__(*args, **kwargs) super(TimeSeries, self).__init__(*args, **kwargs)
...@@ -35,10 +36,6 @@ class TimeSeries(Feature): ...@@ -35,10 +36,6 @@ class TimeSeries(Feature):
for _, v in self._std_dataset.data_vars.items(): for _, v in self._std_dataset.data_vars.items():
v.attrs['coordinates'] = 'time lat lon z' v.attrs['coordinates'] = 'time lat lon z'
@property
def _feature_geodimnames(self) -> Tuple[str, ...]:
return 'time'
def get_geocoord_dimnames( def get_geocoord_dimnames(
self, fieldname: str, *args, **kwargs) -> Tuple[str, ...]: self, fieldname: str, *args, **kwargs) -> Tuple[str, ...]:
if (fieldname in STANDARD_Z_FIELD + ['lat', 'lon'] or if (fieldname in STANDARD_Z_FIELD + ['lat', 'lon'] or
......
...@@ -24,6 +24,8 @@ class TimeSeriesProfile(Feature): ...@@ -24,6 +24,8 @@ class TimeSeriesProfile(Feature):
* ``time(profile)`` * ``time(profile)``
* ``depth/alt (profile, z)`` * ``depth/alt (profile, z)``
""" """
_feature_geodimnames = 'profile', 'z',
def __init__(self, *args, **kwargs): def __init__(self, *args, **kwargs):
super(TimeSeriesProfile, self).__init__(*args, **kwargs) super(TimeSeriesProfile, self).__init__(*args, **kwargs)
...@@ -34,10 +36,6 @@ class TimeSeriesProfile(Feature): ...@@ -34,10 +36,6 @@ class TimeSeriesProfile(Feature):
for _, v in self._std_dataset.data_vars.items(): for _, v in self._std_dataset.data_vars.items():
v.attrs['coordinates'] = 'time lat lon z' v.attrs['coordinates'] = 'time lat lon z'
@property
def _feature_geodimnames(self) -> Tuple[str, ...]:
return 'profile', 'z',
def get_geocoord_dimnames( def get_geocoord_dimnames(
self, fieldname: str, *args, **kwargs) -> Tuple[str, ...]: self, fieldname: str, *args, **kwargs) -> Tuple[str, ...]:
if fieldname in STANDARD_Z_FIELD or self.coords[fieldname].is_axis('Z'): if fieldname in STANDARD_Z_FIELD or self.coords[fieldname].is_axis('Z'):
...@@ -59,10 +57,6 @@ class UniZTimeSeriesProfile(Feature): ...@@ -59,10 +57,6 @@ class UniZTimeSeriesProfile(Feature):
* ``time(profile)`` * ``time(profile)``
* ``depth/alt (z)`` * ``depth/alt (z)``
""" """
@property
def _feature_geodimnames(self) -> Tuple[str, ...]:
return 'profile', 'z',
def get_geocoord_dimnames( def get_geocoord_dimnames(
self, fieldname: str, *args, **kwargs) -> Tuple[str, ...]: self, fieldname: str, *args, **kwargs) -> Tuple[str, ...]:
if fieldname in STANDARD_Z_FIELD or self.coords[fieldname].is_axis('Z'): if fieldname in STANDARD_Z_FIELD or self.coords[fieldname].is_axis('Z'):
......
...@@ -22,9 +22,8 @@ class Trajectory(Feature): ...@@ -22,9 +22,8 @@ class Trajectory(Feature):
``time(time)``, optionally ``z(time)``, all having a single geolocation ``time(time)``, optionally ``z(time)``, all having a single geolocation
dimension: ``time``. dimension: ``time``.
""" """
@property _instance_dimname = 'trajectory'
def _feature_geodimnames(self): _feature_geodimnames = 'time',
return 'time',
def get_geocoord_dimnames(self, *args, **kwargs): def get_geocoord_dimnames(self, *args, **kwargs):
return 'time', return 'time',
......
...@@ -22,6 +22,8 @@ class TrajectoryProfile(Feature): ...@@ -22,6 +22,8 @@ class TrajectoryProfile(Feature):
``time(profile)``, ``depth/alt (profile, z)``, all having a single ``time(profile)``, ``depth/alt (profile, z)``, all having a single
geolocation dimension: ``profile`` (except ``z`` for Z axis coordinate). geolocation dimension: ``profile`` (except ``z`` for Z axis coordinate).
""" """
_feature_geodimnames = 'profile', 'z',
def __init__(self, *args, **kwargs): def __init__(self, *args, **kwargs):
super(TrajectoryProfile, self).__init__(*args, **kwargs) super(TrajectoryProfile, self).__init__(*args, **kwargs)
...@@ -32,10 +34,6 @@ class TrajectoryProfile(Feature): ...@@ -32,10 +34,6 @@ class TrajectoryProfile(Feature):
for _, v in self._std_dataset.data_vars.items(): for _, v in self._std_dataset.data_vars.items():
v.attrs['coordinates'] = 'time lat lon z' v.attrs['coordinates'] = 'time lat lon z'
@property
def _feature_geodimnames(self) -> Tuple[str, ...]:
return 'profile', 'z',
def get_geocoord_dimnames( def get_geocoord_dimnames(
self, fieldname: str, *args, **kwargs) -> Tuple[str, ...]: self, fieldname: str, *args, **kwargs) -> Tuple[str, ...]:
if fieldname in ['lat', 'lon', 'time']: if fieldname in ['lat', 'lon', 'time']:
......
...@@ -166,3 +166,21 @@ class TestFeature(unittest.TestCase): ...@@ -166,3 +166,21 @@ class TestFeature(unittest.TestCase):
ncf = NCDataset(dataset=TEST_FILE) ncf = NCDataset(dataset=TEST_FILE)
feat = self.get_feature(ncf) feat = self.get_feature(ncf)
self.assertRaises(ValueError, feat.get_field, 'unknown_field') self.assertRaises(ValueError, feat.get_field, 'unknown_field')
def test_append(self):
feat = self.get_feature(self.define_base_feature())
# another feature of same dimensions & coordinates
feat2 = self.get_feature(self.define_base_feature())
# append with shared dimension
feat.append(feat2, prefix="v2_")
for fieldname in feat2.fieldnames:
self.assertIn("v2_"+fieldname, feat.fieldnames)
np.allclose(feat.get_values(fieldname), feat2.get_values(fieldname))
# test append with unshared dimensions
feat = self.get_feature(self.define_base_feature())
feat.append(feat2, prefix="v2_", as_new_dims=True)
self.assertIn("v2_time", feat.dims)
self.assertIn("v2_time", feat.coordnames)
import numpy as np
import xarray as xr
from cerbere.feature.point import Point
from cerbere.feature.imdcollection import IMDCollection
from .test_feature import TestFeature
class TestIMDCollectionProfileFeature(TestFeature):
"""Test class for incomplete multidimensional collection of profile
features"""
def __init__(self, methodName="runTest"):
super().__init__(methodName)
def get_feature_class(self):
return IMDCollection
def get_feature(self, *args, **kwargs):
return self.get_feature_class()(
*args, feature_class=Point, **kwargs
)
def define_base_feature(self):
# creates a test xarray object
lon = xr.DataArray(data=[100., 100], dims=['obs'])
lat = xr.DataArray(data=[50., 50], dims=['obs'])
time = xr.DataArray(
[np.datetime64('2018-01-01 00:00:00'),
np.datetime64('2018-01-01 00:00:00')],
dims=['obs']
)
var = xr.DataArray(
data=np.ones(shape=(2,)),
dims=['obs',],
attrs={'myattr': 'test_attr_val'}
)