Skip to content

Commit

Permalink
Merge pull request #575 from LBL-EESA/cf_restripe_regrid
Browse files Browse the repository at this point in the history 
Cf restripe regrid
  • Loading branch information
@burlen
burlen authored Mar 12, 2021
2 parents 02f6b12 + d91e005 commit 5329148
Show file tree
Hide file tree
Showing 16 changed files with 912 additions and 167 deletions.
2 changes: 1 addition & 1 deletion .travis.yml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,7 +17,7 @@ env:
- BUILD_TYPE=Debug
- TECA_DIR=/travis_teca_dir
- TECA_PYTHON_VERSION=3
- TECA_DATA_REVISION=107
- TECA_DATA_REVISION=108
jobs:
- DOCKER_IMAGE=ubuntu IMAGE_VERSION=20.04 IMAGE_NAME=ubuntu_20_04 REQUIRE_NETCDF_MPI=TRUE
- DOCKER_IMAGE=ubuntu IMAGE_VERSION=20.04 IMAGE_NAME=ubuntu_20_04 REQUIRE_NETCDF_MPI=FALSE
Expand Down
173 changes: 141 additions & 32 deletions alg/teca_cartesian_mesh_regrid.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,6 +18,28 @@ using std::endl;

//#define TECA_DEBUG

// always use nearest neighbor interpolation for integers
// to avoid truncation errors. an alternative would be to
// implement rounding in the interpolator for integer types
template <typename data_t>
int get_interpolation_mode(int desired_mode,
typename std::enable_if<std::is_integral<data_t>::value>::type* = 0)
{
(void)desired_mode;
return teca_cartesian_mesh_regrid::nearest;
}

// use the requested interpolation mode for floating point
// data
template <typename data_t>
int get_interpolation_mode(int desired_mode,
typename std::enable_if<std::is_floating_point<data_t>::value>::type* = 0)
{
return desired_mode;
}


// 3D
template<typename NT1, typename NT2, typename NT3, class interp_t>
int interpolate(unsigned long target_nx, unsigned long target_ny,
unsigned long target_nz, const NT1 *p_target_xc, const NT1 *p_target_yc,
Expand Down Expand Up @@ -53,31 +75,100 @@ int interpolate(unsigned long target_nx, unsigned long target_ny,
return 0;
}

// 2D - x-y
template<typename NT1, typename NT2, typename NT3, class interp_t>
int interpolate(unsigned long target_nx, unsigned long target_ny,
const NT1 *p_target_xc, const NT1 *p_target_yc,
NT3 *p_target_a, const NT2 *p_source_xc,
const NT2 *p_source_yc, const NT3 *p_source_a,
unsigned long source_ihi, unsigned long source_jhi,
unsigned long source_nx)
{
interp_t f;
unsigned long q = 0;
for (unsigned long j = 0; j < target_ny; ++j)
{
NT2 ty = static_cast<NT2>(p_target_yc[j]);
for (unsigned long i = 0; i < target_nx; ++i, ++q)
{
NT2 tx = static_cast<NT2>(p_target_xc[i]);
if (f(tx,ty,
p_source_xc, p_source_yc,
p_source_a, source_ihi, source_jhi,
source_nx, p_target_a[q]))
{
TECA_ERROR("failed to interpolate i=(" << i << ", " << j
<< ") x=(" << tx << ", " << ty << ", " << ")")
return -1;
}
}
}
return 0;
}

template<typename taget_coord_t, typename source_coord_t, typename array_t>
int interpolate(int mode, unsigned long target_nx, unsigned long target_ny,
unsigned long target_nz, const taget_coord_t *p_target_xc, const taget_coord_t *p_target_yc,
const taget_coord_t *p_target_zc, array_t *p_target_a, const source_coord_t *p_source_xc,
const source_coord_t *p_source_yc, const source_coord_t *p_source_zc, const array_t *p_source_a,
unsigned long source_ihi, unsigned long source_jhi, unsigned long source_khi,
unsigned long source_nx, unsigned long source_nxy)
unsigned long target_nz, const taget_coord_t *p_target_xc,
const taget_coord_t *p_target_yc, const taget_coord_t *p_target_zc,
array_t *p_target_a, const source_coord_t *p_source_xc,
const source_coord_t *p_source_yc, const source_coord_t *p_source_zc,
const array_t *p_source_a, unsigned long source_ihi, unsigned long source_jhi,
unsigned long source_khi, unsigned long source_nx, unsigned long source_ny,
unsigned long source_nz)
{
using nearest_interp_t = teca_coordinate_util::interpolate_t<0>;
using linear_interp_t = teca_coordinate_util::interpolate_t<1>;

switch (mode)
unsigned long source_nxy = source_nx*source_ny;

switch (get_interpolation_mode<array_t>(mode))
{
case teca_cartesian_mesh_regrid::nearest:
return interpolate<taget_coord_t,source_coord_t,array_t,nearest_interp_t>(
target_nx, target_ny, target_nz, p_target_xc, p_target_yc, p_target_zc,
p_target_a, p_source_xc, p_source_yc, p_source_zc, p_source_a,
source_ihi, source_jhi, source_khi, source_nx, source_nxy);
{
if ((target_nz == 1) && (source_nz == 1))
{
// 2D in the x-y plane
return interpolate<taget_coord_t,
source_coord_t, array_t, nearest_interp_t>(
target_nx, target_ny, p_target_xc, p_target_yc,
p_target_a, p_source_xc, p_source_yc, p_source_a,
source_ihi, source_jhi, source_nx);
}
else
{
// 3D
return interpolate<taget_coord_t,
source_coord_t, array_t, nearest_interp_t>(
target_nx, target_ny, target_nz, p_target_xc,
p_target_yc, p_target_zc, p_target_a, p_source_xc,
p_source_yc, p_source_zc, p_source_a, source_ihi,
source_jhi, source_khi, source_nx, source_nxy);
}
break;
}
case teca_cartesian_mesh_regrid::linear:
return interpolate<taget_coord_t,source_coord_t,array_t,linear_interp_t>(
target_nx, target_ny, target_nz, p_target_xc, p_target_yc, p_target_zc,
p_target_a, p_source_xc, p_source_yc, p_source_zc, p_source_a,
source_ihi, source_jhi, source_khi, source_nx, source_nxy);
{
if ((target_nz == 1) && (source_nz == 1))
{
// 2D in the x-y plane
return interpolate<taget_coord_t,
source_coord_t, array_t, linear_interp_t>(
target_nx, target_ny, p_target_xc, p_target_yc,
p_target_a, p_source_xc, p_source_yc, p_source_a,
source_ihi, source_jhi, source_nx);
}
else
{
// 3D
return interpolate<taget_coord_t,
source_coord_t, array_t, linear_interp_t>(
target_nx, target_ny, target_nz, p_target_xc,
p_target_yc, p_target_zc, p_target_a, p_source_xc,
p_source_yc, p_source_zc, p_source_a, source_ihi,
source_jhi, source_khi, source_nx, source_nxy);
}
break;
}
}

TECA_ERROR("invalid interpolation mode \"" << mode << "\"")
Expand All @@ -86,6 +177,7 @@ int interpolate(int mode, unsigned long target_nx, unsigned long target_ny,




// --------------------------------------------------------------------------
teca_cartesian_mesh_regrid::teca_cartesian_mesh_regrid()
: target_input(0), interpolation_mode(nearest)
Expand Down Expand Up @@ -308,7 +400,6 @@ std::vector<teca_metadata> teca_cartesian_mesh_regrid::get_upstream_request(
return up_reqs;
}


// --------------------------------------------------------------------------
const_p_teca_dataset teca_cartesian_mesh_regrid::execute(
unsigned int port, const std::vector<const_p_teca_dataset> &input_data,
Expand Down Expand Up @@ -351,7 +442,7 @@ const_p_teca_dataset teca_cartesian_mesh_regrid::execute(

// get the list of arrays to move
std::vector<std::string> req_arrays;
request.get("regrid_arrays", req_arrays);
request.get("arrays", req_arrays);

// add any explicitly named
std::copy(this->arrays.begin(), this->arrays.end(),
Expand All @@ -369,6 +460,13 @@ const_p_teca_dataset teca_cartesian_mesh_regrid::execute(
source_arrays.push_back(*it);
}

// catch a user error
if (!source_arrays.size() &&
teca_mpi_util::mpi_rank_0(this->get_communicator()))
{
TECA_WARNING("No arrays will be interpolated")
}

// move the arrays
const_p_teca_variant_array target_xc = target->get_x_coordinates();
const_p_teca_variant_array target_yc = target->get_y_coordinates();
Expand All @@ -388,28 +486,27 @@ const_p_teca_dataset teca_cartesian_mesh_regrid::execute(
unsigned long source_nx = source_xc->size();
unsigned long source_ny = source_yc->size();
unsigned long source_nz = source_zc->size();
unsigned long source_nxy = source_nx*source_ny;
unsigned long source_ihi = source_nx - 1;
unsigned long source_jhi = source_ny - 1;
unsigned long source_khi = source_nz - 1;

NESTED_TEMPLATE_DISPATCH_FP(
const teca_variant_array_impl,
target_xc.get(),
1,
_TGT,

const NT1 *p_target_xc = std::dynamic_pointer_cast<TT1>(target_xc)->get();
const NT1 *p_target_yc = std::dynamic_pointer_cast<TT1>(target_yc)->get();
const NT1 *p_target_zc = std::dynamic_pointer_cast<TT1>(target_zc)->get();
const NT_TGT *p_target_xc = std::dynamic_pointer_cast<TT_TGT>(target_xc)->get();
const NT_TGT *p_target_yc = std::dynamic_pointer_cast<TT_TGT>(target_yc)->get();
const NT_TGT *p_target_zc = std::dynamic_pointer_cast<TT_TGT>(target_zc)->get();

NESTED_TEMPLATE_DISPATCH_FP(
const teca_variant_array_impl,
source_xc.get(),
2,
_SRC,

const NT2 *p_source_xc = std::dynamic_pointer_cast<TT2>(source_xc)->get();
const NT2 *p_source_yc = std::dynamic_pointer_cast<TT2>(source_yc)->get();
const NT2 *p_source_zc = std::dynamic_pointer_cast<TT2>(source_zc)->get();
const NT_SRC *p_source_xc = std::dynamic_pointer_cast<TT_SRC>(source_xc)->get();
const NT_SRC *p_source_yc = std::dynamic_pointer_cast<TT_SRC>(source_yc)->get();
const NT_SRC *p_source_zc = std::dynamic_pointer_cast<TT_SRC>(source_zc)->get();

size_t n_arrays = source_arrays.size();
for (size_t i = 0; i < n_arrays; ++i)
Expand All @@ -421,25 +518,37 @@ const_p_teca_dataset teca_cartesian_mesh_regrid::execute(
NESTED_TEMPLATE_DISPATCH(
teca_variant_array_impl,
target_a.get(),
3,
_DATA,

const NT3 *p_source_a = std::static_pointer_cast<const TT3>(source_a)->get();
NT3 *p_target_a = std::static_pointer_cast<TT3>(target_a)->get();
const NT_DATA *p_source_a = std::static_pointer_cast<const TT_DATA>(source_a)->get();
NT_DATA *p_target_a = std::static_pointer_cast<TT_DATA>(target_a)->get();

if (interpolate(this->interpolation_mode, target_nx, target_ny, target_nz,
p_target_xc, p_target_yc, p_target_zc, p_target_a, p_source_xc,
p_source_yc, p_source_zc, p_source_a, source_ihi, source_jhi,
source_khi, source_nx, source_nxy))
source_khi, source_nx, source_ny, source_nz))
{
TECA_ERROR("Failed to move \"" << source_arrays[i] << "\"")
return nullptr;
}

target_ac->set(source_arrays[i], target_a);
)
else
{
TECA_ERROR("Unsupported array type " << source_a->get_class_name())
}
)

target_ac->set(source_arrays[i], target_a);
}
)
else
{
TECA_ERROR("Unupported coordinate type " << source_xc->get_class_name())
}
)
else
{
TECA_ERROR("Unupported coordinate type " << target_xc->get_class_name())
}

return target;
}
2 changes: 1 addition & 1 deletion alg/teca_cartesian_mesh_regrid.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,7 +21,7 @@ By default the first input is the target mesh. the second input is the source
mesh. This can be changed by setting the target_input property.
the arrays to move from source to target can be selected using add_array api or
in the request key regrid_source_arrays. this is a spatial regriding operation
in the request key "arrays". this is a spatial regriding operation
for temporal regriding see teca_mesh_temporal_regrid.
*/
class teca_cartesian_mesh_regrid : public teca_algorithm
Expand Down
66 changes: 65 additions & 1 deletion alg/teca_cartesian_mesh_source.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -94,7 +94,7 @@ teca_cartesian_mesh_source::teca_cartesian_mesh_source() :
field_type_code(teca_variant_array_code<double>::get()),
x_axis_variable("lon"), y_axis_variable("lat"), z_axis_variable("plev"),
t_axis_variable("time"), x_axis_units("degrees_east"),
y_axis_units("degrees_north"), z_axis_units("pascals"),
y_axis_units("degrees_north"), z_axis_units("Pa"),
calendar("Gregorian"), time_units("seconds since 1970-01-01 00:00:00"),
whole_extents{0l, 359l, 0l, 179l, 0l, 0l, 0l, 0l},
bounds{0., 360, -90., 90., 0., 0., 0., 0.},
Expand Down Expand Up @@ -129,6 +129,70 @@ void teca_cartesian_mesh_source::set_properties(const std::string &prefix,
}
#endif

// --------------------------------------------------------------------------
int teca_cartesian_mesh_source::set_spatial_bounds(const teca_metadata &md)
{
teca_metadata coords;
if (md.get("coordinates", coords))
return -1;

// get the bounds in the x direction
p_teca_variant_array x = coords.get("x");

if (!x)
return -1;

x->get(0lu, this->bounds[0]);
x->get(x->size() - 1lu, this->bounds[1]);

// get the bounds in the y direction
p_teca_variant_array y = coords.get("y");

if (!y)
return -1;

y->get(0lu, this->bounds[2]);
y->get(y->size() - 1lu, this->bounds[3]);

// get the bounds in the z direction
p_teca_variant_array z = coords.get("z");

if (!z)
return -1;

z->get(0lu, this->bounds[4]);
z->get(z->size() - 1lu, this->bounds[5]);

// set the coordinate type
this->set_coordinate_type_code(x->type_code());

return 0;
}

// --------------------------------------------------------------------------
int teca_cartesian_mesh_source::set_calendar(const teca_metadata &md)
{
teca_metadata atts;
if (md.get("attributes", atts))
return -1;

teca_metadata time_atts;
if (atts.get("time", time_atts))
return -1;

std::string calendar;
if (time_atts.get("calendar", calendar))
return -1;

std::string units;
if (time_atts.get("units", units))
return -1;

this->calendar = calendar;
this->time_units = units;

return 0;
}

// --------------------------------------------------------------------------
void teca_cartesian_mesh_source::set_modified()
Expand Down
Loading

0 comments on commit 5329148

Please sign in to comment.