Commit 5ed17114 authored by Jeff Piollé's avatar Jeff Piollé
Browse files

added resampling of multi-dimensional fields

parent 4f6bb624
......@@ -499,7 +499,68 @@ class Resampler(object):
values = source_times
else:
values = source.get_values(fieldname)
res = kd_tree.get_sample_from_neighbour_info(
# multi-dimensional array : interpolate all layers
if len(values.shape) > 2:
field = source.get_field(fieldname)
dims = field.dimensions
# get non spatial dimensions of source
dims = field.dimensions
if source.__class__.__name__ == 'Swath':
spatial_dim_indices = [
dims.keys().index('row'),
dims.keys().index('cell')
]
nonspatial_dims = copy.copy(dims)
del nonspatial_dims['row']
del nonspatial_dims['cell']
elif source.__class__.__name__ == 'Grid':
spatial_dim_indices = [
dims.keys().index('y'),
dims.keys().index('x')
]
nonspatial_dims = copy.copy(dims)
del nonspatial_dims['y']
del nonspatial_dims['x']
else:
raise NotImplementedError(
"Class: ".format(source.__class__.__name__))
# 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)
def layers(dims):
indices = range(dims[0])
if len(dims) == 1:
return [[i] for i in indices]
layerlist = []
for i in indices:
for j in layers(dims[1:]):
layerlist.append([i] + j)
return layerlist
for layer in layers(nonspatial_dims.values()):
subset = tuple(
[slice(None, None, None) if i in spatial_dim_indices else layer[i]
for i in range(len(dims))])
res[..., :, :] = kd_tree.get_sample_from_neighbour_info(
'nn',
target_def.shape,
values[subset],
valid_input_index,
valid_output_index,
index_array,
fill_value=None
)
# two-dimensional array
else:
res = kd_tree.get_sample_from_neighbour_info(
'nn',
target_def.shape,
values,
......@@ -508,9 +569,11 @@ class Resampler(object):
index_array,
fill_value=None
)
target_shape = target.get_geolocation_dimsizes().values()
if isinstance(target, Grid):
res = res.reshape(target_shape)
target_shape = target.get_geolocation_dimsizes().values()
if isinstance(target, Grid):
res = res.reshape(target_shape)
new_dims = target.get_geolocation_dimsizes()
if fieldname in ['time', 'lat', 'lon']:
if not res.any():
# return None if no valid source data can be resampled on
......@@ -529,7 +592,6 @@ class Resampler(object):
var = source.get_field(fieldname).variable
sourcefield = source.get_field(fieldname)
var.shortname = newfieldname
new_dims = target.get_geolocation_dimsizes()
field = Field(var,
new_dims,
values=res,
......
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