xref: /dports/science/cdo/cdo-2.0.0/src/remap_distwgt.cc

/*
  This file is part of CDO. CDO is a collection of Operators to manipulate and analyse Climate model Data.

  Author: Uwe Schulzweida

*/

#include <atomic>

#include "process_int.h"
#include "cdo_wtime.h"
#include "remap.h"
#include "remap_store_link.h"
#include "cdo_options.h"
#include "progress.h"
#include "cimdOmp.h"

// Interpolation using a distance-weighted average

// This routine computes the inverse-distance weights for a nearest-neighbor interpolation
void
remap_distwgt_weights(size_t numNeighbors, RemapSearch &rsearch, RemapVars &rv)
{
  auto src_grid = rsearch.srcGrid;
  auto tgt_grid = rsearch.tgtGrid;

  if (Options::cdoVerbose) cdo_print("Called %s()", __func__);

  progress::init();

  // Compute mappings from source to target grid

  auto tgt_grid_size = tgt_grid->size;

  std::vector<WeightLinks> weightLinks(tgt_grid_size);
  weight_links_alloc(numNeighbors, tgt_grid_size, weightLinks);

  std::vector<knnWeightsType> knnWeights;
  for (int i = 0; i < Threading::ompNumThreads; ++i) knnWeights.push_back(knnWeightsType(numNeighbors));

  const auto start = Options::cdoVerbose ? cdo_get_wtime() : 0.0;

  // Loop over destination grid

  std::atomic<size_t> atomicCount{0};

#ifdef _OPENMP
#pragma omp parallel for default(none) \
    shared(atomicCount, rsearch, weightLinks, numNeighbors, src_grid, tgt_grid, tgt_grid_size, knnWeights)
#endif
  for (size_t tgt_cell_add = 0; tgt_cell_add < tgt_grid_size; ++tgt_cell_add)
    {
      const auto ompthID = cdo_omp_get_thread_num();

      atomicCount++;
      if (ompthID == 0) progress::update(0, 1, (double)atomicCount / tgt_grid_size);

      weightLinks[tgt_cell_add].nlinks = 0;

      if (!tgt_grid->mask[tgt_cell_add]) continue;

      const auto llpoint = remapgrid_get_lonlat(tgt_grid, tgt_cell_add);

      // Find nearest grid points on source grid and distances to each point
      remap_search_points(rsearch, llpoint, knnWeights[ompthID]);

      // Compute weights based on inverse distance if mask is false, eliminate those points
      const auto nadds = knnWeights[ompthID].computeWeights(src_grid->mask);

      for (size_t n = 0; n < nadds; ++n)
        if (knnWeights[ompthID].m_mask[n]) tgt_grid->cell_frac[tgt_cell_add] = 1.0;

      // Store the link
      store_weightlinks(0, nadds, knnWeights[ompthID].m_addr.data(), knnWeights[ompthID].m_dist.data(), tgt_cell_add, weightLinks);
    }

  progress::update(0, 1, 1);

  grid_point_search_delete(rsearch.gps);

  weight_links_to_remap_links(0, tgt_grid_size, weightLinks, rv);

  if (numNeighbors == 1) rv.links_per_value = numNeighbors;

  if (Options::cdoVerbose) cdo_print("Point search nearest: %.2f seconds", cdo_get_wtime() - start);
}  // remap_distwgt_weights

template <typename T>
static void
remap_dist_wgt(size_t numNeighbors, RemapSearch &rsearch, const Varray<T> &src_array, Varray<T> &tgt_array, T missval)
{
  auto src_grid = rsearch.srcGrid;
  auto tgt_grid = rsearch.tgtGrid;

  if (Options::cdoVerbose) cdo_print("Called %s()", __func__);

  progress::init();

  // Compute mappings from source to target grid

  auto tgt_grid_size = tgt_grid->size;
  auto src_grid_size = src_grid->size;

  Varray<short> src_grid_mask(src_grid_size);
#ifdef _OPENMP
#pragma omp parallel for default(none) schedule(static) shared(src_grid_size, src_array, src_grid_mask, missval)
#endif
  for (size_t i = 0; i < src_grid_size; ++i) src_grid_mask[i] = !DBL_IS_EQUAL(src_array[i], missval);

  std::vector<knnWeightsType> knnWeights;
  for (int i = 0; i < Threading::ompNumThreads; ++i) knnWeights.push_back(knnWeightsType(numNeighbors));

  const auto start = Options::cdoVerbose ? cdo_get_wtime() : 0.0;

  // Loop over destination grid

  std::atomic<size_t> atomicCount{0};

#ifdef _OPENMP
#pragma omp parallel for default(none) shared(atomicCount, rsearch, numNeighbors, src_grid, tgt_grid, tgt_grid_size, src_array, \
                                              tgt_array, missval, src_grid_mask, knnWeights)
#endif
  for (size_t tgt_cell_add = 0; tgt_cell_add < tgt_grid_size; ++tgt_cell_add)
    {
      const auto ompthID = cdo_omp_get_thread_num();

      atomicCount++;
      if (ompthID == 0) progress::update(0, 1, (double)atomicCount / tgt_grid_size);

      tgt_array[tgt_cell_add] = missval;

      if (!tgt_grid->mask[tgt_cell_add]) continue;

      const auto llpoint = remapgrid_get_lonlat(tgt_grid, tgt_cell_add);

      // Find nearest grid points on source grid and distances to each point
      remap_search_points(rsearch, llpoint, knnWeights[ompthID]);

      // Compute weights based on inverse distance if mask is false, eliminate those points
      const auto nadds = knnWeights[ompthID].computeWeights(src_grid_mask);
      if (nadds) tgt_array[tgt_cell_add] = knnWeights[ompthID].arrayWeightsSum(src_array);
    }

  progress::update(0, 1, 1);

  if (Options::cdoVerbose) cdo_print("Point search nearest: %.2f seconds", cdo_get_wtime() - start);
}  // remap_dist_wgt

void
remap_dist_wgt(size_t numNeighbors, RemapSearch &rsearch, const Field &field1, Field &field2)
{
  if (field1.memType == MemType::Float)
    remap_dist_wgt(numNeighbors, rsearch, field1.vec_f, field2.vec_f, (float) field1.missval);
  else
    remap_dist_wgt(numNeighbors, rsearch, field1.vec_d, field2.vec_d, field1.missval);
}

void remapInit(RemapType &remap);

void
intgriddis(Field &field1, Field &field2, size_t numNeighbors)
{
  auto mapType = RemapMethod::DISTWGT;
  auto gridID1 = field1.grid;
  auto gridID2 = field2.grid;
  auto src_missval = field1.missval;
  auto tgt_missval = field2.missval;
  auto src_array = field1.vec_d.data();
  auto tgt_array = field2.vec_d.data();

  if (Options::cdoVerbose) cdo_print("Called %s()", __func__);

  progress::init();

  // Interpolate from source to target grid

  RemapType remap;
  remapInit(remap);

  bool remap_extrapolate = false;
  remap_init_grids(mapType, remap_extrapolate, gridID1, remap.src_grid, gridID2, remap.tgt_grid);

  auto src_grid_size = remap.src_grid.size;
  auto tgt_grid_size = remap.tgt_grid.size;

  Varray<short> src_mask(src_grid_size);
  for (size_t i = 0; i < src_grid_size; ++i) src_mask[i] = !DBL_IS_EQUAL(src_array[i], src_missval);
  // Varray<short> tgt_mask(tgt_grid_size, 1);

  std::vector<knnWeightsType> knnWeights;
  for (int i = 0; i < Threading::ompNumThreads; ++i) knnWeights.push_back(knnWeightsType(numNeighbors));

  remap_search_init(mapType, remap.search, remap.src_grid, remap.tgt_grid);

  const auto start = Options::cdoVerbose ? cdo_get_wtime() : 0.0;

  // Loop over destination grid

  size_t nmiss = 0;
  double findex = 0.0;

  for (size_t tgt_cell_add = 0; tgt_cell_add < tgt_grid_size; ++tgt_cell_add)
    {
      const auto ompthID = cdo_omp_get_thread_num();

      findex++;
      if (ompthID == 0) progress::update(0, 1, findex / tgt_grid_size);

      tgt_array[tgt_cell_add] = tgt_missval;

      // if (!tgt_mask[tgt_cell_add]) continue;

      const auto llpoint = remapgrid_get_lonlat(&remap.tgt_grid, tgt_cell_add);

      // Find nearest grid points on source grid and distances to each point
      remap_search_points(remap.search, llpoint, knnWeights[ompthID]);

      // Compute weights based on inverse distance if mask is false, eliminate those points
      const auto nadds = knnWeights[ompthID].computeWeights(src_mask);
      if (nadds)
        tgt_array[tgt_cell_add] = knnWeights[ompthID].arrayWeightsSum(src_array);
      else
        nmiss++;
    }

  progress::update(0, 1, 1);

  field2.nmiss = nmiss;

  remap_grid_free(remap.src_grid);
  remap_grid_free(remap.tgt_grid);
  remap_search_free(remap.search);

  if (Options::cdoVerbose) cdo_print("Point search nearest: %.2f seconds", cdo_get_wtime() - start);
}  // intgriddis