Added Selregion.cc.

src/Selregion.cc

+#include <cdi.h>
+
+#include "process_int.h"
+#include "param_conversion.h"
+#include "pmlist.h"
+#include "grid_point_search.h"
+#include "mpim_grid.h"
+
+int gengridcell(const int gridID1, const size_t gridsize2, const std::vector<long> &cellidx);
+void window_cell(const Field &field1, Field &field2, const std::vector<long> &cellidx);
+double radiusDegToKm(const double radiusInDeg);
+
+struct CirclePoint
+{
+  double radius = 1.0;
+  double lon = 0.0;
+  double lat = 0.0;
+  size_t maxpoints = SIZE_MAX;
+};
+
+struct RegionInfo
+{
+  int gridtype = -1;
+  int gridID1 = -1, gridID2 = -1;
+  std::vector<long> cellidx;
+  long nvals = 0;
+};
+
+static int
+generate_circle_grid(int gridID1, long &gridsize2, std::vector<long> &cellidx, const CirclePoint &cpoint)
+{
+  const auto gridID0 = gridID1;
+
+  gridID1 = generate_full_point_grid(gridID1);
+
+  if (!gridHasCoordinates(gridID1)) cdoAbort("Cell center coordinates missing!");
+
+  {
+    const auto gridsize1 = gridInqSize(gridID1);
+
+    Varray<double> xvals(gridsize1), yvals(gridsize1);
+    gridInqXvals(gridID1, xvals.data());
+    gridInqYvals(gridID1, yvals.data());
+
+    // Convert lat/lon units if required
+    cdo_grid_to_radian(gridID1, CDI_XAXIS, gridsize1, xvals.data(), "grid center lon");
+    cdo_grid_to_radian(gridID1, CDI_YAXIS, gridsize1, yvals.data(), "grid center lat");
+
+    GridPointSearch gps;
+    gridPointSearchCreate(gps, xvals, yvals, PointSearchMethod::spherepart);
+
+    gps.searchRadius = cpoint.radius;
+
+    auto numNeighbors = cpoint.maxpoints;
+    if (numNeighbors > gridsize1) numNeighbors = gridsize1;
+
+    knnWeightsType knnWeights(numNeighbors);
+    gridSearchPointSmooth(gps, cpoint.lon, cpoint.lat, knnWeights);
+
+    const auto nvals = knnWeights.m_numNeighbors;
+    cellidx.resize(nvals);
+  
+    for (size_t i = 0; i < nvals; ++i) cellidx[i] = knnWeights.m_addr[i];
+
+    gridPointSearchDelete(gps);
+
+    if (nvals == 0) cdoAbort("No grid points found!");
+
+    gridsize2 = nvals;
+  }
+
+  const auto gridID2 = gengridcell(gridID1, gridsize2, cellidx);
+
+  if (gridID0 != gridID1) gridDestroy(gridID1);
+
+  return gridID2;
+}
+
+static void
+selcircleGetParameter(CirclePoint &cpoint)
+{
+  const auto pargc = operatorArgc();
+  if (pargc)
+    {
+      const auto pargv = cdoGetOperArgv();
+
+      KVList kvlist;
+      kvlist.name = "SELCIRCLE";
+      if (kvlist.parseArguments(pargc, pargv) != 0) cdoAbort("Parse error!");
+      if (Options::cdoVerbose) kvlist.print();
+
+      for (const auto &kv : kvlist)
+        {
+          const auto &key = kv.key;
+          if (kv.nvalues > 1) cdoAbort("Too many values for parameter key >%s<!", key.c_str());
+          if (kv.nvalues < 1) cdoAbort("Missing value for parameter key >%s<!", key.c_str());
+          const auto &value = kv.values[0];
+
+          // clang-format off
+          if      (key == "maxpoints") cpoint.maxpoints = parameter2sizet(value);
+          else if (key == "lon")       cpoint.lon = parameter2double(value);
+          else if (key == "lat")       cpoint.lat = parameter2double(value);
+          else if (key == "radius")    cpoint.radius = radius_str_to_deg(value.c_str());
+          else cdoAbort("Invalid parameter key >%s<!", key.c_str());
+          // clang-format on
+        }
+    }
+}
+
+void *
+Selregion(void *process)
+{
+  int index;
+
+  cdoInitialize(process);
+
+  // clang-format off
+  const auto SELCIRCLE  = cdoOperatorAdd("selcircle",   0,   0, nullptr);
+  // clang-format on
+
+  const auto operatorID = cdoOperatorID();
+
+  CirclePoint cpoint;
+  if (operatorID == SELCIRCLE) selcircleGetParameter(cpoint);
+
+  if (cpoint.radius < 0 || cpoint.radius > 180) cdoAbort("%s=%g out of bounds (0-180 deg)!", "radius", cpoint.radius);
+
+  const auto streamID1 = cdoOpenRead(0);
+
+  const auto vlistID1 = cdoStreamInqVlist(streamID1);
+  const auto vlistID2 = vlistDuplicate(vlistID1);
+
+  const auto taxisID1 = vlistInqTaxis(vlistID1);
+  const auto taxisID2 = taxisDuplicate(taxisID1);
+  vlistDefTaxis(vlistID2, taxisID2);
+
+  const auto nvars = vlistNvars(vlistID1);
+  std::vector<bool> varIDs(nvars, false);
+
+  const auto ngrids = vlistNgrids(vlistID1);
+  std::vector<RegionInfo> regions(ngrids);
+
+  const auto gridsizemax = vlistGridsizeMax(vlistID1);
+  if (gridsizemax < cpoint.maxpoints) cpoint.maxpoints = gridsizemax;
+  if (Options::cdoVerbose && operatorID == SELCIRCLE)
+    cdoPrint("lon = %g, lat = %g, radius = %gdeg(%gkm)",
+             cpoint.lon, cpoint.lat, cpoint.radius, radiusDegToKm(cpoint.radius));
+
+  cpoint.radius *= DEG2RAD;
+  cpoint.lon *= DEG2RAD;
+  cpoint.lat *= DEG2RAD;
+
+  for (index = 0; index < ngrids; index++)
+    {
+      auto &region = regions[index];
+      const auto gridID1 = vlistGrid(vlistID1, index);
+      const auto gridtype = gridInqType(gridID1);
+      if (is_point_grid(gridID1))
+        {
+          int gridID2 = -1;
+          if (operatorID == SELCIRCLE)
+            {
+              const auto gridsize = gridInqSize(gridID1);
+              if (gridsize == 1) continue;
+
+              gridID2 = generate_circle_grid(gridID1, region.nvals, region.cellidx, cpoint);
+            }
+
+          region.gridtype = gridtype;
+          region.gridID1 = gridID1;
+          region.gridID2 = gridID2;
+
+          vlistChangeGridIndex(vlistID2, index, gridID2);
+
+          for (int varID = 0; varID < nvars; varID++)
+            if (gridID1 == vlistInqVarGrid(vlistID1, varID)) varIDs[varID] = true;
+        }
+      else if (gridtype == GRID_GENERIC && gridInqXsize(gridID1) <= 1 && gridInqYsize(gridID1) <= 1)
+        {
+        }
+      else
+        {
+          cdoPrint("Unsupported grid type: %s", gridNamePtr(gridtype));
+          if (gridtype == GRID_GAUSSIAN_REDUCED) cdoPrint("Use option -R to convert Gaussian reduced grid to a regular grid!");
+          cdoAbort("Unsupported grid type!");
+        }
+    }
+
+  VarList varList1, varList2;
+  varListInit(varList1, vlistID1);
+  varListInit(varList2, vlistID2);
+
+  Field field1, field2;
+
+  const auto streamID2 = cdoOpenWrite(1);
+  cdoDefVlist(streamID2, vlistID2);
+
+  int tsID = 0;
+  while (true)
+    {
+      const auto nrecs = cdoStreamInqTimestep(streamID1, tsID);
+      if (nrecs == 0) break;
+
+      taxisCopyTimestep(taxisID2, taxisID1);
+      cdoDefTimestep(streamID2, tsID);
+
+      for (int recID = 0; recID < nrecs; recID++)
+        {
+          int varID, levelID;
+          cdoInqRecord(streamID1, &varID, &levelID);
+          field1.init(varList1[varID]);
+          cdoReadRecord(streamID1, field1);
+
+          cdoDefRecord(streamID2, varID, levelID);
+
+          if (varIDs[varID])
+            {
+              const auto gridID1 = varList1[varID].gridID;
+              for (index = 0; index < ngrids; index++)
+                if (gridID1 == regions[index].gridID1) break;
+              if (index == ngrids) cdoAbort("Internal problem, grid not found!");
+
+              field2.init(varList2[varID]);
+              window_cell(field1, field2, regions[index].cellidx);
+
+              if (field1.nmiss) fieldNumMV(field2);
+
+              cdoWriteRecord(streamID2, field2);
+            }
+          else
+            {
+              cdoWriteRecord(streamID2, field1);
+            }
+        }
+
+      tsID++;
+    }
+
+  cdoStreamClose(streamID2);
+  cdoStreamClose(streamID1);
+
+  cdoFinish();
+
+  return nullptr;
+}