Commit 2990631b authored by Jeff Piollé's avatar Jeff Piollé
Browse files

Merge branch 'master' of http://git.cersat.fr/cerbere/cerinterp

parents 4b13cac0 5228085a
...@@ -195,17 +195,17 @@ class Resampler(object): ...@@ -195,17 +195,17 @@ class Resampler(object):
missing = source.get_field(fieldname).fillvalue missing = source.get_field(fieldname).fillvalue
else: else:
missing = fillvalue missing = fillvalue
logging.debug('missing value: %s',missing)
if mask is None: if mask is None:
srcmask = numpy.array(data.mask, copy=True) srcmask = numpy.array(data.mask, copy=True)
if extrapolation != 0: if extrapolation != 0:
logging.info('Read source and extrapolate over missing data') logging.debug('Read source and extrapolate over missing data')
# call extrapolation procedure (replace NaN with interpolated # call extrapolation procedure (replace NaN with interpolated
# values) # values)
logging.debug('data %s',data.shape)
data = cls.extrapolate(data, extrapolation) data = cls.extrapolate(data, extrapolation)
# masked data to missing value # masked data to missing value
data = data.filled(missing) data = data.filled(missing)
logging.debug('resampler | source lon %s lat %s z %s',source_lons.shape,source_lats.shape,data.shape)
logging.info("interpolation to target grid")
# re-interpolate source over target # re-interpolate source over target
# --------------------------------- # ---------------------------------
source_lons_new = Resampler.__monotonic(source_lons) source_lons_new = Resampler.__monotonic(source_lons)
...@@ -231,7 +231,7 @@ class Resampler(object): ...@@ -231,7 +231,7 @@ class Resampler(object):
) )
else: else:
# extrapolate the source mask # extrapolate the source mask
mask = Resampler.__closest_neighbour(source_lats, extrapolate_mask = Resampler.__closest_neighbour(source_lats,
source_lons, source_lons,
srcmask.astype('int8'), srcmask.astype('int8'),
source.get_model_name(), source.get_model_name(),
...@@ -240,7 +240,7 @@ class Resampler(object): ...@@ -240,7 +240,7 @@ class Resampler(object):
target.get_model_name(), target.get_model_name(),
) )
interpolated_data = numpy.ma.masked_where( interpolated_data = numpy.ma.masked_where(
mask == 1, extrapolate_mask == 1,
interpolated_data, interpolated_data,
copy=False copy=False
) )
...@@ -307,11 +307,10 @@ class Resampler(object): ...@@ -307,11 +307,10 @@ class Resampler(object):
`pyresample <https://pyresample.readthedocs.org/en/latest/index.html>`_ `pyresample <https://pyresample.readthedocs.org/en/latest/index.html>`_
package documentation. package documentation.
""" """
logging.info('%s %s | __get_pyresample_def',__name__,__file__)
if (srclon>180).any(): if (srclon>180).any():
srclon = utils.wrap_longitudes(srclon) srclon = utils.wrap_longitudes(srclon)
if srcmodel == 'Swath' or srcmodel == 'Image' or \ if srcmodel == 'Swath' or srcmodel == 'Image' or \
srcmodel == 'Trajectory': srcmodel == 'Trajectory' or srcmodel =='PointCollection':
source_def = geometry.SwathDefinition(lons=srclon, source_def = geometry.SwathDefinition(lons=srclon,
lats=srclat lats=srclat
) )
...@@ -334,10 +333,15 @@ class Resampler(object): ...@@ -334,10 +333,15 @@ class Resampler(object):
lats=tgtlat lats=tgtlat
) )
if radius is None: if radius is None:
y, x = srclat.shape if srcmodel =='PointCollection':
resolution = abs(srclat[y / 2, x / 2] - srclat[y / 2 + 1, x / 2]) # resolution = abs(srclat[0] - srclat[1])
# resolution = numpy.amin(abs(numpy.diff(srclat)))
resolution = 1
else:
y, x = srclat.shape
resolution = abs(srclat[y / 2, x / 2] - srclat[y / 2 + 1, x / 2])
radius = resolution * 100000. / 2. * 1.6 radius = resolution * 100000. / 2. * 1.6
logging.info("Radius for resampling : %s", radius) logging.debug("Radius for resampling : %s", radius)
return source_def, target_def, radius return source_def, target_def, radius
@classmethod @classmethod
...@@ -390,11 +394,12 @@ class Resampler(object): ...@@ -390,11 +394,12 @@ class Resampler(object):
valid_input_index, valid_output_index, index_array, distance_array\ valid_input_index, valid_output_index, index_array, distance_array\
= Resampler.__get_closest_neighbour_info(source_def, = Resampler.__get_closest_neighbour_info(source_def,
target_def, target_def,
radius_of_influence=radius) radius=radius)
if not valid_input_index.any(): if not valid_input_index.any():
# return None if no source data can be resampled on the # return None if no source data can be resampled on the
# target # target
return None return None
logging.debug('target %s srcdtaa %s',target_def.shape,srcdata)
res = kd_tree.get_sample_from_neighbour_info('nn', res = kd_tree.get_sample_from_neighbour_info('nn',
target_def.shape, target_def.shape,
srcdata, srcdata,
...@@ -452,7 +457,7 @@ class Resampler(object): ...@@ -452,7 +457,7 @@ class Resampler(object):
source_lons = source.get_lon().astype('float32') source_lons = source.get_lon().astype('float32')
source_lats = source.get_lat().astype('float32') source_lats = source.get_lat().astype('float32')
source_times = source.get_times() source_times = source.get_times()
if isinstance(source, Grid): if isinstance(source, Grid) and len(source_lons.shape) == 1:
source_lons, source_lats = numpy.meshgrid(source_lons, source_lats) source_lons, source_lats = numpy.meshgrid(source_lons, source_lats)
source_times = numpy.resize(source_times, source_lons.shape) source_times = numpy.resize(source_times, source_lons.shape)
target_lons = target.get_lon().astype('float32') target_lons = target.get_lon().astype('float32')
...@@ -467,7 +472,6 @@ class Resampler(object): ...@@ -467,7 +472,6 @@ class Resampler(object):
target_lons, target_lons,
target.__class__.__name__, target.__class__.__name__,
radius) radius)
if new: if new:
# create a clone of target without any field # create a clone of target without any field
outfeature = target.extract_subset(fields=[]) outfeature = target.extract_subset(fields=[])
...@@ -641,7 +645,8 @@ class Resampler(object): ...@@ -641,7 +645,8 @@ class Resampler(object):
fieldx = Field(varx, fieldx = Field(varx,
target_dims, target_dims,
values=reference_x, values=reference_x,
datatype=numpy.dtype(numpy.int32)) datatype=numpy.dtype(numpy.int32),
fillvalue=-2363645994848887)
vary = Variable( vary = Variable(
shortname='resampled_%s' % ydimname, shortname='resampled_%s' % ydimname,
description='%s index of the selected pixel' % ydimname description='%s index of the selected pixel' % ydimname
...@@ -649,7 +654,8 @@ class Resampler(object): ...@@ -649,7 +654,8 @@ class Resampler(object):
fieldy = Field(vary, fieldy = Field(vary,
target_dims, target_dims,
values=reference_y, values=reference_y,
datatype=numpy.dtype(numpy.int32)) datatype=numpy.dtype(numpy.int32),
fillvalue=-2363645994848887)
outfeature.add_field(fieldx) outfeature.add_field(fieldx)
outfeature.add_field(fieldy) outfeature.add_field(fieldy)
# Add the distance # Add the distance
...@@ -707,6 +713,8 @@ class Resampler(object): ...@@ -707,6 +713,8 @@ class Resampler(object):
logging.debug("Resample to regular grid") logging.debug("Resample to regular grid")
srclons = src.get_lon() srclons = src.get_lon()
srclats = src.get_lat() srclats = src.get_lat()
logging.debug('srclons: %s',srclons)
logging.debug('srclats: %s',srclats)
if not target.projection.is_cylindrical(): if not target.projection.is_cylindrical():
raise Exception('Projection is not valid for this type of regridding') raise Exception('Projection is not valid for this type of regridding')
tarlon = target.get_lon() tarlon = target.get_lon()
...@@ -730,6 +738,7 @@ class Resampler(object): ...@@ -730,6 +738,7 @@ class Resampler(object):
val = src.get_values(fieldname) val = src.get_values(fieldname)
else: else:
val = data[i] val = data[i]
logging.debug('val %s %s',val.shape,val.mean())
lons = numpy.ma.masked_where(val.mask, srclons).compressed() lons = numpy.ma.masked_where(val.mask, srclons).compressed()
lats = numpy.ma.masked_where(val.mask, srclats).compressed() lats = numpy.ma.masked_where(val.mask, srclats).compressed()
# create empty field if add_empty_fields is set when there are # create empty field if add_empty_fields is set when there are
...@@ -749,10 +758,14 @@ class Resampler(object): ...@@ -749,10 +758,14 @@ class Resampler(object):
logging.warning("Converting to double this field: %s" % logging.warning("Converting to double this field: %s" %
fieldname) fieldname)
val = val.astype('float64') val = val.astype('float64')
logging.debug('latmin %s latmax %s lonmin %s lonmax %s',latmin , latmax , lonmin , lonmax)
val_comp = val.compressed()
logging.debug('shape val_compressed: %s',val_comp.shape)
logging.debug('lats compressed :%s',lats.shape)
result, nbsampl = fastregrid.fastregrid( result, nbsampl = fastregrid.fastregrid(
lats, lats,
lons, lons,
val.compressed(), val_comp,
latmin, latmax, lonmin, lonmax, latmin, latmax, lonmin, lonmax,
lonbin, latbin lonbin, latbin
) )
......
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