Commit ce40eecf authored by PIOLLE's avatar PIOLLE
Browse files

adaptation to cerbere2

parent f11f2f0d
......@@ -16,12 +16,11 @@ try:
except:
pass
from cerbere.datamodel.field import Field
from cerbere.datamodel.variable import Variable
from cerbere.datamodel.swath import Swath
from cerbere.datamodel.grid import Grid
from cerbere.datamodel.image import Image
from cerbere.datamodel.trajectory import Trajectory
from cerbere.dataset.field import Field
from cerbere.feature.swath import Swath
from cerbere.feature.grid import Grid
from cerbere.feature.image import Image
from cerbere.feature.trajectory import Trajectory
try:
# from cerbereutils.resampling import fastregrid
......@@ -225,8 +224,7 @@ def interpolation(source,
# extrapolate the source mask
extrapolate_mask = __closest_neighbour(source_lats,
source_lons,
srcmask.astype(
'int8'),
srcmask.astype('int8'),
source.get_model_name(),
target_lats,
target_lons,
......@@ -255,30 +253,26 @@ def interpolation(source,
raise Exception('Field already existing in target grid')
else:
prefix = 'interpolated_'
var = source.get_field(fieldname).variable
sourcefield = source.get_field(fieldname)
var.shortname = prefix + var.shortname
new_dims = target.get_geolocation_dimsizes()
field = Field(var,
new_dims,
values=interpolated_data,
units=sourcefield.units,
attributes=copy.copy(sourcefield.attributes),
fillvalue=sourcefield.fillvalue,
valid_min=sourcefield.valid_min,
valid_max=sourcefield.valid_max,
)
new_dims = target.geodimsizes
field = Field(
interpolated_data,
name=prefix + fieldname,
dims=new_dims,
units=sourcefield.units,
attrs=copy.copy(sourcefield.attrs),
fillvalue=sourcefield.fillvalue)
outfeature.add_field(field)
# Add geolocation
if add_geolocation:
# add time and lat/lon
for fieldname in ['time']:
field = source.get_geolocation_field(fieldname).clone()
if len(field.dimensions) == 1:
field.variable.shortname = 'source_time'
field = source.get_coord(fieldname).clone()
if len(field.dims) == 1:
field.name = 'source_time'
field.set_values(
source.get_geolocation_field(fieldname).get_values().copy()
source.get_coord(fieldname).get_values().copy()
)
outfeature.add_field(field)
return outfeature
......@@ -301,31 +295,23 @@ def __get_pyresample_def(srclat,
"""
if (srclon > 180).any():
srclon = utils.wrap_longitudes(srclon)
if srcmodel == 'Swath' or srcmodel == 'Image' or \
srcmodel == 'Trajectory' or srcmodel == 'PointCollection':
source_def = geometry.SwathDefinition(lons=srclon,
lats=srclat
)
elif srcmodel == 'Grid':
if srcmodel == 'Grid':
if len(srclon.shape) == 1:
srclon, srclat = numpy.meshgrid(srclon, srclat)
source_def = geometry.GridDefinition(lons=srclon, lats=srclat)
else:
source_def = geometry.SwathDefinition(lons=srclon, lats=srclat)
source_def = geometry.GridDefinition(lons=srclon,
lats=srclat
)
if tgtmodel == 'Swath' or tgtmodel == 'Image' or \
tgtmodel == 'Trajectory':
target_def = geometry.SwathDefinition(lons=tgtlon,
lats=tgtlat
)
elif tgtmodel == 'Grid':
if tgtmodel == 'Grid':
if len(tgtlon.shape) == 1:
tgtlon, tgtlat = numpy.meshgrid(tgtlon, tgtlat)
target_def = geometry.GridDefinition(lons=tgtlon,
lats=tgtlat
)
target_def = geometry.GridDefinition(lons=tgtlon, lats=tgtlat)
else:
target_def = geometry.SwathDefinition(lons=tgtlon, lats=tgtlat)
if radius is None:
if srcmodel == 'PointCollection':
if 'Collection' in srcmodel:
# resolution = abs(srclat[0] - srclat[1])
# resolution = numpy.amin(abs(numpy.diff(srclat)))
resolution = 1
......@@ -449,7 +435,7 @@ def closest_neighbour(source,
"""
if not type(fields) is list:
fields = [fields]
if fields == []:
if len(fields) == 0:
logging.warning("No fields requested? slipping resampling")
return
......@@ -509,10 +495,10 @@ def closest_neighbour(source,
if len(values.shape) > 2:
field = source.get_field(fieldname)
dims = field.dimensions
dims = field.dims
# get non spatial dimensions of source
dims = field.dimensions
dims = field.dims
if source.__class__.__name__ == 'Swath':
spatial_dim_indices = [
dims.keys().index('row'),
......@@ -535,9 +521,8 @@ def closest_neighbour(source,
# calculate dimensions of resampled field
new_dims = copy.copy(nonspatial_dims)
new_dims.update(target.get_geolocation_dimsizes())
res = numpy.ma.masked_all(new_dims.values(),
dtype=values.dtype)
new_dims.update(target.geodimsizes)
res = numpy.ma.masked_all(new_dims.values(), dtype=values.dtype)
def layers(dims):
indices = range(dims[0])
......@@ -575,10 +560,10 @@ def closest_neighbour(source,
index_array,
fill_value=None
)
target_shape = target.get_geolocation_dimsizes().values()
target_shape = target.geodimsizes
if isinstance(target, Grid):
res = res.reshape(target_shape)
new_dims = target.get_geolocation_dimsizes()
new_dims = target.geodimsizes
if fieldname in ['time', 'lat', 'lon']:
if not res.any():
......@@ -593,27 +578,23 @@ def closest_neighbour(source,
'Field {} already existing in target grid'.format(newfieldname))
else:
if fieldname in ['lat', 'lon', 'time']:
var = copy.copy(source._geolocation_fields[fieldname].variable)
sourcefield = source._geolocation_fields[fieldname]
sourcefield = source.get_coord(fieldname)
else:
var = source.get_field(fieldname).variable
sourcefield = source.get_field(fieldname)
var.shortname = newfieldname
field = Field(var,
new_dims,
values=res,
units=sourcefield.units,
attributes=copy.copy(sourcefield.attributes),
fillvalue=sourcefield.fillvalue,
valid_min=sourcefield.valid_min,
valid_max=sourcefield.valid_max,
)
field = Field(
res,
name=newfieldname,
dims=new_dims,
units=sourcefield.units,
attrs=copy.copy(sourcefield.attributes),
fillvalue=sourcefield.fillvalue)
outfeature.add_field(field)
if add_reference:
# Add the indices (in the source feature data) of the
# selected pixels
target_dims = target.get_geolocation_dimsizes()
target_dims = target.geodimsizes
if isinstance(source, Swath) or isinstance(source, Image) or isinstance(
source, Trajectory):
xdimname = 'cell'
......@@ -642,41 +623,37 @@ def closest_neighbour(source,
index_array = index_array.reshape(target_shape)
reference_x = numpy.ma.array(index_array % nx)
reference_y = numpy.ma.array(index_array / nx)
varx = Variable(
shortname='resampled_%s' % xdimname,
description='%s index of the selected pixel' % xdimname
)
fieldx = Field(varx,
target_dims,
values=reference_x,
datatype=numpy.dtype(numpy.int32),
fillvalue=-2363645994848887)
vary = Variable(
shortname='resampled_%s' % ydimname,
description='%s index of the selected pixel' % ydimname
)
fieldy = Field(vary,
target_dims,
values=reference_y,
datatype=numpy.dtype(numpy.int32),
fillvalue=-2363645994848887)
fieldx = Field(
reference_x,
name='resampled_%{}'.format(xdimname),
description='{} index of the selected pixel'.format(xdimname),
dims=target_dims,
dtype=numpy.dtype(numpy.int32),
fillvalue=-2363645994848887)
fieldy = Field(
reference_y,
name='resampled_%{}'.format(ydimname),
description='{} index of the selected pixel'.format(ydimname),
dims=target_dims,
dtype=numpy.dtype(numpy.int32),
fillvalue=-2363645994848887)
outfeature.add_field(fieldx)
outfeature.add_field(fieldy)
# Add the distance
vardist = Variable(
shortname='distance_to_pixel',
description='distance of the selected pixel to the pixel center'
)
fielddist = Field(
vardist,
target_dims,
values=numpy.ma.array(distance_array.reshape(target_shape),
mask=res.mask
),
numpy.ma.array(distance_array.reshape(target_shape), mask=res.mask),
name='distance_to_pixel',
description='distance of the selected pixel to the pixel center',
dims=target_dims,
units='m',
datatype=numpy.dtype(numpy.float32)
dtype=numpy.dtype(numpy.float32)
)
outfeature.add_field(fielddist)
return outfeature
......@@ -802,19 +779,20 @@ def resample_to_regular_grid(src, target, fields, data=None,
field.dimensions = geodims
field.set_values(result)
else:
var = Variable(shortname=fieldname)
field = Field(
var,
dimensions=geodims,
values=result
)
result,
name=fieldname,
dims=geodims)
field.attributes['binning_method'] = 'mean'
target.add_field(field)
# return None if there is valid remapped data
if empty:
logging.warning("There were no valid data in any remapped fields."
" No field was returned.")
return None
return
if use_src_time:
times = src.get_times()
if times.shape == (1,):
......@@ -843,16 +821,19 @@ def resample_to_regular_grid(src, target, fields, data=None,
# convert back to initial (integer) type
if result.dtype != valdtype:
times = numpy.ma.floor(times).astype(valdtype)
timefield = target.get_geolocation_field('time')
timefield = target.get_coord('time')
if timefield is None:
timefield = Field(
variable=copy.copy(
src.get_geolocation_field('time').variable),
dimensions=geodims,
times,
name='time',
standard_name=src.get_coord('time').standard_name,
description=src.get_coord('time').description,
dims=geodims,
units=src.get_time_units(),
values=times
)
target._geolocation_fields['time'] = timefield
target.add_field(timefield)
return target
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment