cdo-2.0.0/src/Selbox.cc

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

  Author: Uwe Schulzweida

*/

/*
   This module contains the following operators:

      Selbox     sellonlatbox    Select lon/lat box
      Selbox     selindexbox     Select index box
*/

#include <cdi.h>

#include <mpim_grid.h>
#include "grid_define.h"
#include "cdo_options.h"
#include "process_int.h"
#include "cdo_vlist.h"
#include "param_conversion.h"
#include "gridreference.h"
#include "selboxinfo.h"

static void
correct_xvals(const long nlon, const long inc, double *xvals)
{
  if (nlon > 1 && IS_EQUAL(xvals[0], xvals[(nlon - 1) * inc])) xvals[(nlon - 1) * inc] += 360;

  if (xvals[0] > xvals[(nlon - 1) * inc])
    for (long i = 0; i < nlon; i++)
      if (xvals[i * inc] >= 180) xvals[i * inc] -= 360;

  for (long i = 0; i < nlon; i++)
    {
      if (xvals[i * inc] < -180) xvals[i * inc] += 360;
      if (xvals[i * inc] > 360) xvals[i * inc] -= 360;
    }

  if (xvals[0] > xvals[(nlon - 1) * inc])
    for (long i = 1; i < nlon; i++)
      if (xvals[i * inc] < xvals[(i - 1) * inc]) xvals[i * inc] += 360;
}

static void
gengridxyvals(const int gridtype, const int gridID1, const int gridID2, const long nlon, const long nlat, const long nlon2,
              const long nlat2, const SelboxInfo &sbox, const bool unitsIsDegree)
{
  Varray<double> xvals1, yvals1, xvals2, yvals2;

  const bool lxvals = gridInqXvals(gridID1, nullptr) > 0;
  const bool lyvals = gridInqYvals(gridID1, nullptr) > 0;

  if (gridtype == GRID_CURVILINEAR)
    {
      if (lxvals && lyvals)
        {
          xvals1.resize(nlon * nlat);
          yvals1.resize(nlon * nlat);
          xvals2.resize(nlon2 * nlat2);
          yvals2.resize(nlon2 * nlat2);
        }
    }
  else
    {
      if (lxvals) xvals1.resize(nlon);
      if (lyvals) yvals1.resize(nlat);
      if (lxvals) xvals2.resize(nlon2);
      if (lyvals) yvals2.resize(nlat2);
    }

  auto pxvals2 = xvals2.data();
  auto pyvals2 = yvals2.data();

  if (lxvals) gridInqXvals(gridID1, xvals1.data());
  if (lyvals) gridInqYvals(gridID1, yvals1.data());

  const auto &lat1 = sbox.lat1;
  const auto &lat2 = sbox.lat2;
  const auto &lon11 = sbox.lon11;
  const auto &lon12 = sbox.lon12;
  const auto &lon21 = sbox.lon21;
  const auto &lon22 = sbox.lon22;

  if (gridtype == GRID_CURVILINEAR)
    {
      if (lxvals && lyvals)
        for (long ilat = lat1; ilat <= lat2; ilat++)
          {
            for (long ilon = lon21; ilon <= lon22; ilon++)
              {
                *pxvals2++ = xvals1[ilat * nlon + ilon];
                *pyvals2++ = yvals1[ilat * nlon + ilon];
              }
            for (long ilon = lon11; ilon <= lon12; ilon++)
              {
                *pxvals2++ = xvals1[ilat * nlon + ilon];
                *pyvals2++ = yvals1[ilat * nlon + ilon];
              }
          }
    }
  else
    {
      if (lxvals)
        {
          for (long i = lon21; i <= lon22; i++) *pxvals2++ = xvals1[i];
          for (long i = lon11; i <= lon12; i++) *pxvals2++ = xvals1[i];
          if (unitsIsDegree) correct_xvals(nlon2, 1, xvals2.data());
        }

      if (lyvals)
        for (long i = lat1; i <= lat2; i++) *pyvals2++ = yvals1[i];
    }
  // for ( int i = 0; i < nlat2; i++ ) printf("lat : %d %g\n", i+1, yvals2[i]);
  // for ( int i = 0; i < nlon2; i++ ) printf("lon : %d %g\n", i+1, xvals2[i]);
  if (lxvals) gridDefXvals(gridID2, xvals2.data());
  if (lyvals) gridDefYvals(gridID2, yvals2.data());
}

static void
gengridXboundsCurvi(const int gridID1, const long nlon, const long nlat, const long nlon2, const long nlat2, const SelboxInfo &sbox,
                    Varray<double> &xbounds1, Varray<double> &xbounds2)
{
  xbounds1.resize(4 * nlon * nlat);
  xbounds2.resize(4 * nlon2 * nlat2);

  auto pxbounds2 = xbounds2.data();

  gridInqXbounds(gridID1, xbounds1.data());

  const auto &lat1 = sbox.lat1;
  const auto &lat2 = sbox.lat2;
  const auto &lon11 = sbox.lon11;
  const auto &lon12 = sbox.lon12;
  const auto &lon21 = sbox.lon21;
  const auto &lon22 = sbox.lon22;
  for (long ilat = lat1; ilat <= lat2; ilat++)
    {
      for (long ilon = 4 * lon21; ilon < 4 * (lon22 + 1); ilon++) *pxbounds2++ = xbounds1[4 * ilat * nlon + ilon];

      for (long ilon = 4 * lon11; ilon < 4 * (lon12 + 1); ilon++) *pxbounds2++ = xbounds1[4 * ilat * nlon + ilon];
    }
}

static void
gengridYboundsCurvi(const int gridID1, const long nlon, const long nlat, const long nlon2, const long nlat2, const SelboxInfo &sbox,
                    Varray<double> &ybounds1, Varray<double> &ybounds2)
{
  ybounds1.resize(4 * nlon * nlat);
  ybounds2.resize(4 * nlon2 * nlat2);

  auto pybounds2 = ybounds2.data();

  gridInqYbounds(gridID1, ybounds1.data());

  const auto &lat1 = sbox.lat1;
  const auto &lat2 = sbox.lat2;
  const auto &lon11 = sbox.lon11;
  const auto &lon12 = sbox.lon12;
  const auto &lon21 = sbox.lon21;
  const auto &lon22 = sbox.lon22;
  for (long ilat = lat1; ilat <= lat2; ilat++)
    {
      for (long ilon = 4 * lon21; ilon < 4 * (lon22 + 1); ilon++) *pybounds2++ = ybounds1[4 * ilat * nlon + ilon];

      for (long ilon = 4 * lon11; ilon < 4 * (lon12 + 1); ilon++) *pybounds2++ = ybounds1[4 * ilat * nlon + ilon];
    }
}

static void
gengridXboundsRect2D(const int gridID1, const long nlon, const long nlon2, const SelboxInfo &sbox, const bool unitsIsDegree,
                     Varray<double> &xbounds1, Varray<double> &xbounds2)
{
  xbounds1.resize(2 * nlon);
  xbounds2.resize(2 * nlon2);

  auto pxbounds2 = xbounds2.data();

  gridInqXbounds(gridID1, xbounds1.data());

  const auto &lon11 = sbox.lon11;
  const auto &lon12 = sbox.lon12;
  const auto &lon21 = sbox.lon21;
  const auto &lon22 = sbox.lon22;
  for (long i = 2 * lon21; i < 2 * (lon22 + 1); i++) *pxbounds2++ = xbounds1[i];
  for (long i = 2 * lon11; i < 2 * (lon12 + 1); i++) *pxbounds2++ = xbounds1[i];

  if (unitsIsDegree)
    {
      correct_xvals(nlon2, 2, xbounds2.data());
      correct_xvals(nlon2, 2, xbounds2.data() + 1);
      // make sure that the first bound is less than the second
      long nx = 0;
      for (long i = 0; i < nlon2; ++i)
        if (xbounds2[2 * i] > xbounds2[2 * i + 1]) nx++;
      if (nx == nlon2 && xbounds2[0] > -180.)
        for (long i = 0; i < nlon2; ++i) xbounds2[2 * i] -= 360.;
    }
}

static void
gengridYboundsRect2D(const int gridID1, const long nlat, const long nlat2, const SelboxInfo &sbox, Varray<double> &ybounds1,
                     Varray<double> &ybounds2)
{
  ybounds1.resize(2 * nlat);
  ybounds2.resize(2 * nlat2);

  auto pybounds2 = ybounds2.data();

  gridInqYbounds(gridID1, ybounds1.data());

  const auto &lat1 = sbox.lat1;
  const auto &lat2 = sbox.lat2;
  for (long i = 2 * lat1; i < 2 * (lat2 + 1); i++) *pybounds2++ = ybounds1[i];
}

static void
gengridXbounds(const int gridID1, const int gridID2, const long nlon, const long nlat, const long nlon2, const long nlat2,
               const SelboxInfo &sbox, const bool unitsIsDegree)
{
  Varray<double> xbounds1, xbounds2;
  const auto gridtype = gridInqType(gridID1);

  if (gridtype == GRID_CURVILINEAR)
    {
      gridDefNvertex(gridID2, 4);
      gengridXboundsCurvi(gridID1, nlon, nlat, nlon2, nlat2, sbox, xbounds1, xbounds2);
    }
  else
    {
      gridDefNvertex(gridID2, 2);
      gengridXboundsRect2D(gridID1, nlon, nlon2, sbox, unitsIsDegree, xbounds1, xbounds2);
    }

  gridDefXbounds(gridID2, xbounds2.data());
}

static void
gengridYbounds(const int gridID1, const int gridID2, const long nlon, const long nlat, const long nlon2, const long nlat2,
               const SelboxInfo &sbox)
{
  Varray<double> ybounds1, ybounds2;
  const auto gridtype = gridInqType(gridID1);

  if (gridtype == GRID_CURVILINEAR)
    {
      gridDefNvertex(gridID2, 4);
      gengridYboundsCurvi(gridID1, nlon, nlat, nlon2, nlat2, sbox, ybounds1, ybounds2);
    }
  else
    {
      gridDefNvertex(gridID2, 2);
      gengridYboundsRect2D(gridID1, nlat, nlat2, sbox, ybounds1, ybounds2);
    }

  gridDefYbounds(gridID2, ybounds2.data());
}

static int
gengrid(const int gridID1, const SelboxInfo &sbox)
{
  const auto &lat1 = sbox.lat1;
  const auto &lat2 = sbox.lat2;
  const auto &lon11 = sbox.lon11;
  const auto &lon12 = sbox.lon12;
  const auto &lon21 = sbox.lon21;
  const auto &lon22 = sbox.lon22;
  const long nlon = gridInqXsize(gridID1);
  const long nlat = gridInqYsize(gridID1);

  const long nlon21 = sbox.lon12 - sbox.lon11 + 1;
  const long nlon22 = sbox.lon22 - sbox.lon21 + 1;
  const long nlon2 = nlon21 + nlon22;
  const long nlat2 = sbox.lat2 - sbox.lat1 + 1;

  if (Options::cdoVerbose) cdo_print("nlon1=%ld  nlat1=%ld", nlon, nlat);
  if (Options::cdoVerbose) cdo_print("nlon2=%ld  nlat2=%ld", nlon2, nlat2);

  const auto gridtype = gridInqType(gridID1);

  const auto gridID2 = gridCreate(gridtype, nlon2 * nlat2);
  if (nlon > 0) gridDefXsize(gridID2, nlon2);
  if (nlat > 0) gridDefYsize(gridID2, nlat2);
  if (nlat > 0) gridDefNP(gridID2, gridInqNP(gridID1));

  int datatype = CDI_UNDEFID;
  cdiInqKeyInt(gridID1, CDI_GLOBAL, CDI_KEY_DATATYPE, &datatype);
  cdiDefKeyInt(gridID2, CDI_GLOBAL, CDI_KEY_DATATYPE, datatype);

  grid_copy_keys(gridID1, gridID2);

  if (gridtype == GRID_PROJECTION) grid_copy_mapping(gridID1, gridID2);

  char xunits[CDI_MAX_NAME];
  int length = CDI_MAX_NAME;
  cdiInqKeyString(gridID1, CDI_XAXIS, CDI_KEY_UNITS, xunits, &length);

  const bool unitsIsDegree = strncmp(xunits, "degree", 6) == 0;

  gengridxyvals(gridtype, gridID1, gridID2, nlon, nlat, nlon2, nlat2, sbox, unitsIsDegree);

  if (gridInqXbounds(gridID1, nullptr)) gengridXbounds(gridID1, gridID2, nlon, nlat, nlon2, nlat2, sbox, unitsIsDegree);
  if (gridInqYbounds(gridID1, nullptr)) gengridYbounds(gridID1, gridID2, nlon, nlat, nlon2, nlat2, sbox);

  if (gridtype == GRID_CURVILINEAR && gridHasArea(gridID1))
    {
      Varray<double> areaIn(nlon * nlat), areaOut(nlon2 * nlat2);
      gridInqArea(gridID1, areaIn.data());
      auto pareaOut = areaOut.data();

      for (long ilat = lat1; ilat <= lat2; ilat++)
        {
          for (long ilon = lon21; ilon < (lon22 + 1); ilon++) *pareaOut++ = areaIn[ilat * nlon + ilon];
          for (long ilon = lon11; ilon < (lon12 + 1); ilon++) *pareaOut++ = areaIn[ilat * nlon + ilon];
        }

      gridDefArea(gridID2, areaOut.data());
    }

  const auto projID1 = gridInqProj(gridID1);
  if (projID1 != CDI_UNDEFID && gridInqType(projID1) == GRID_PROJECTION)
    {
      const auto projID2 = gridCreate(GRID_PROJECTION, nlon2 * nlat2);
      gridDefXsize(projID2, nlon2);
      gridDefYsize(projID2, nlat2);

      grid_copy_keys(projID1, projID2);
      grid_copy_mapping(projID1, projID2);

      gengridxyvals(GRID_PROJECTION, projID1, projID2, nlon, nlat, nlon2, nlat2, sbox, false);

      gridDefProj(gridID2, projID2);
    }

  return gridID2;
}

int
gengridcell(const int gridID1, const size_t gridsize2, const std::vector<long> &cellidx)
{
  int gridID2 = -1;
  auto gridtype = gridInqType(gridID1);
  const auto gridsize1 = gridInqSize(gridID1);

  const auto setDimName = (gridtype == GRID_CURVILINEAR);
  if (gridtype == GRID_CURVILINEAR) gridtype = GRID_UNSTRUCTURED;

  if (gridtype == GRID_UNSTRUCTURED || gridtype == GRID_GENERIC)
    gridID2 = gridCreate(gridtype, gridsize2);
  else
    return gridID2;

  int datatype = CDI_UNDEFID;
  cdiInqKeyInt(gridID1, CDI_GLOBAL, CDI_KEY_DATATYPE, &datatype);
  cdiDefKeyInt(gridID2, CDI_GLOBAL, CDI_KEY_DATATYPE, datatype);

  grid_copy_keys(gridID1, gridID2);

  if (setDimName) cdiDefKeyString(gridID2, CDI_XAXIS, CDI_KEY_DIMNAME, "ncells");

  if (gridHasCoordinates(gridID1))
    {
      Varray<double> xvals1(gridsize1), yvals1(gridsize1);
      Varray<double> xvals2(gridsize2), yvals2(gridsize2);

      gridInqXvals(gridID1, xvals1.data());
      gridInqYvals(gridID1, yvals1.data());

      for (size_t i = 0; i < gridsize2; ++i)
        {
          xvals2[i] = xvals1[cellidx[i]];
          yvals2[i] = yvals1[cellidx[i]];
        }

      gridDefXvals(gridID2, xvals2.data());
      gridDefYvals(gridID2, yvals2.data());
    }

  if (gridHasBounds(gridID1))
    {
      const auto nv = gridInqNvertex(gridID1);

      Varray<double> xbounds1(nv * gridsize1), ybounds1(nv * gridsize1);
      Varray<double> xbounds2(nv * gridsize2), ybounds2(nv * gridsize2);

      gridInqXbounds(gridID1, xbounds1.data());
      gridInqYbounds(gridID1, ybounds1.data());

      gridDefNvertex(gridID2, nv);

      for (size_t i = 0; i < gridsize2; ++i)
        {
          for (int k = 0; k < nv; ++k)
            {
              xbounds2[i * nv + k] = xbounds1[cellidx[i] * nv + k];
              ybounds2[i * nv + k] = ybounds1[cellidx[i] * nv + k];
            }
        }

      gridDefXbounds(gridID2, xbounds2.data());
      gridDefYbounds(gridID2, ybounds2.data());
    }

  if (gridHasArea(gridID1))
    {
      Varray<double> areaIn(gridsize1), areaOut(gridsize2);
      gridInqArea(gridID1, areaIn.data());

      for (size_t i = 0; i < gridsize2; ++i) areaOut[i] = areaIn[cellidx[i]];

      gridDefArea(gridID2, areaOut.data());
    }

  return gridID2;
}

void
genlonlatboxReg2d(int gridID, double xlon1, double xlon2, double xlat1, double xlat2, SelboxInfo &sbox)
{
  auto &lat1 = sbox.lat1;
  auto &lat2 = sbox.lat2;
  auto &lon11 = sbox.lon11;
  auto &lon12 = sbox.lon12;
  auto &lon21 = sbox.lon21;
  auto &lon22 = sbox.lon22;
  lat1 = 0;
  lat2 = 0;
  lon11 = 1;
  lon12 = 0;
  lon21 = 0;
  lon22 = 0;

  if (IS_NOT_EQUAL(xlon1, xlon2))
    {
      xlon2 -= 360 * std::floor((xlon2 - xlon1) / 360);
      if (IS_EQUAL(xlon1, xlon2)) xlon2 += 360;
    }
  else
    {
      xlon2 += 0.00001;
    }

  const long nlon = gridInqXsize(gridID);
  const long nlat = gridInqYsize(gridID);

  Varray<double> xvals(nlon), yvals(nlat);

  if (nlon > 0)
    {
      gridInqXvals(gridID, xvals.data());

      cdo_grid_to_degree(gridID, CDI_XAXIS, nlon, xvals.data(), "grid center lon");

      xlon2 -= 360 * std::floor((xlon1 - xvals[0]) / 360);
      xlon1 -= 360 * std::floor((xlon1 - xvals[0]) / 360);
    }

  if (nlat > 0)
    {
      gridInqYvals(gridID, yvals.data());

      cdo_grid_to_degree(gridID, CDI_YAXIS, nlat, yvals.data(), "grid center lat");
    }

  if (nlon > 0)
    {
      // while ( nlon == 1 || (xvals[nlon-1] - xvals[0]) >= 360 ) nlon--;

      for (lon21 = 0; lon21 < nlon && xvals[lon21] < xlon1; lon21++)
        ;
      for (lon22 = lon21; lon22 < nlon && xvals[lon22] < xlon2; lon22++)
        ;

      if (lon22 >= nlon || xvals[lon22] > xlon2) lon22--;

      xlon1 -= 360;
      xlon2 -= 360;

      for (lon11 = 0; xvals[lon11] < xlon1; lon11++)
        ;
      for (lon12 = lon11; lon12 < nlon && xvals[lon12] < xlon2; lon12++)
        ;

      // lon12--;
      if (lon12 >= nlon || xvals[lon12] > xlon2) lon12--;
      if (lon12 >= 0 && IS_EQUAL(xvals[lon12], xvals[lon21])) lon12--;

      if (lon12 - lon11 + 1 + lon22 - lon21 + 1 <= 0) cdo_abort("Longitudinal dimension is too small!");
    }

  if (nlat > 0)
    {
      if (yvals[0] > yvals[nlat - 1])
        {
          if (xlat1 > xlat2)
            {
              for (lat1 = 0; lat1 < nlat && yvals[lat1] > xlat1; lat1++)
                ;
              for (lat2 = nlat - 1; lat2 && yvals[lat2] < xlat2; lat2--)
                ;
            }
          else
            {
              for (lat1 = 0; lat1 < nlat && yvals[lat1] > xlat2; lat1++)
                ;
              for (lat2 = nlat - 1; lat2 && yvals[lat2] < xlat1; lat2--)
                ;
            }
        }
      else
        {
          if (xlat1 < xlat2)
            {
              for (lat1 = 0; lat1 < nlat && yvals[lat1] < xlat1; lat1++)
                ;
              for (lat2 = nlat - 1; lat2 && yvals[lat2] > xlat2; lat2--)
                ;
            }
          else
            {
              for (lat1 = 0; lat1 < nlat && yvals[lat1] < xlat2; lat1++)
                ;
              for (lat2 = nlat - 1; lat2 && yvals[lat2] > xlat1; lat2--)
                ;
            }
        }

      if (lat2 - lat1 + 1 <= 0) cdo_abort("Latitudinal dimension is too small!");
    }
}

static void
genlonlatboxCurv(const int gridID, double xlon1, double xlon2, double xlat1, double xlat2, SelboxInfo &sbox)
{
  auto &lat1 = sbox.lat1;
  auto &lat2 = sbox.lat2;
  auto &lon11 = sbox.lon11;
  auto &lon12 = sbox.lon12;
  auto &lon21 = sbox.lon21;
  auto &lon22 = sbox.lon22;
  const long nlon = gridInqXsize(gridID);
  const long nlat = gridInqYsize(gridID);
  const size_t gridsize = nlon * nlat;

  const auto grid_is_circular = gridIsCircular(gridID);

  Varray<double> xvals(gridsize), yvals(gridsize);

  gridInqXvals(gridID, xvals.data());
  gridInqYvals(gridID, yvals.data());

  // Convert lat/lon units if required
  cdo_grid_to_degree(gridID, CDI_XAXIS, gridsize, xvals.data(), "grid center lon");
  cdo_grid_to_degree(gridID, CDI_YAXIS, gridsize, yvals.data(), "grid center lat");

  if (xlon1 > xlon2) cdo_abort("The second longitude have to be greater than the first one!");
  /*
  if ( xlon2 < 180. )
    {
      xlon1 += 360;
      xlon2 += 360;
    }
  */
  if (xlat1 > xlat2) std::swap(xlat1, xlat2);

  lat1 = nlat - 1;
  lat2 = 0;
  lon11 = 0;
  lon12 = -1;
  // lon11 = nlon-1;
  // lon12 = 0;
  lon21 = nlon - 1;
  lon22 = 0;

  bool lp2 = false;
  double xfirst, xlast, ylast;

  if (grid_is_circular)
    {
      for (long ilat = 0; ilat < nlat; ilat++)
        {
          xlast = xvals[ilat * nlon + nlon - 1];
          ylast = yvals[ilat * nlon + nlon - 1];
          if (ylast >= xlat1 && ylast <= xlat2)
            if (xlon1 <= xlast && xlon2 > xlast && (xlon2 - xlon1) < 360)
              {
                lon11 = nlon - 1;
                lon12 = 0;
                lp2 = true;
                break;
              }
        }
    }

  for (long ilat = 0; ilat < nlat; ilat++)
    for (long ilon = 0; ilon < nlon; ilon++)
      {
        const auto xval = xvals[ilat * nlon + ilon];
        const auto yval = yvals[ilat * nlon + ilon];
        if (yval >= xlat1 && yval <= xlat2)
          {
            if (lp2)
              {
                xfirst = xvals[ilat * nlon];
                if (xfirst < xlon1) xfirst = xlon1;

                xlast = xvals[ilat * nlon + nlon - 1];
                if (xlast > xlon2) xlast = xlon2;

                // printf("%g %g %g %g %g %g\n", yval, xval, xlon1, xlast,  xfirst, xlon2);
                if (xval >= xlon1 && xval <= xlast)
                  {
                    if (ilon < lon21) lon21 = ilon;
                    if (ilon > lon22) lon22 = ilon;
                    if (ilat < lat1) lat1 = ilat;
                    if (ilat > lat2) lat2 = ilat;
                  }
                else if (xval >= xfirst && xval <= xlon2)
                  {
                    if (ilon < lon11) lon11 = ilon;
                    if (ilon > lon12) lon12 = ilon;
                    if (ilat < lat1) lat1 = ilat;
                    if (ilat > lat2) lat2 = ilat;
                  }
              }
            else
              {
                if (((xval >= xlon1 && xval <= xlon2) || (xval - 360 >= xlon1 && xval - 360 <= xlon2)
                     || (xval + 360 >= xlon1 && xval + 360 <= xlon2)))
                  {
                    if (ilon < lon21) lon21 = ilon;
                    if (ilon > lon22) lon22 = ilon;
                    if (ilat < lat1) lat1 = ilat;
                    if (ilat > lat2) lat2 = ilat;
                  }
              }
            /*
            if ( xval >= xlon1 && xval <= xlon2 )
              {
                if ( ilon < lon21 ) lon21 = ilon;
                if ( ilon > lon22 ) lon22 = ilon;
                if ( ilat < lat1 ) lat1 = ilat;
                if ( ilat > lat2 ) lat2 = ilat;
              }
            else if ( xval >= xlon1-360 && xval <= xlon2-360 )
              {
                if ( ilon < lon11 ) lon11 = ilon;
                if ( ilon > lon12 ) lon12 = ilon;
                if ( ilat < lat1 ) lat1 = ilat;
                if ( ilat > lat2 ) lat2 = ilat;
              }
            */
          }
      }

  // while ( lon12 >= lon21 ) lon12--;
  // if ( lon12 <= lon11 ) { lon11 = nlon-1; lon12 = 0; }
  if (lon12 == 0 && lon11 > 0)
    {
      lon11 = 0;
      lon12 = -1;
    }
  /*
  printf("lon21, lon22, lon11, lon12  idx: %ld %ld %ld %ld  lon: %g %g %g %g\n",
         lon21, lon22, lon11, lon12, xvals[lon21], xvals[lon22], xvals[lon11], xvals[lon12]);
  */
  if (lat2 - lat1 + 1 <= 0) cdo_abort("Latitudinal dimension is too small!");
}

void
getlonlatparams(const int argc_offset, double &xlon1, double &xlon2, double &xlat1, double &xlat2)
{
  bool lset = false;

  const auto nargc = cdo_operator_argc() - argc_offset;
  if (nargc == 1)
    {
      const auto gridname = cdo_operator_argv(argc_offset + 0).c_str();
      if (strcmp(gridname, "europe") == 0)
        {
          xlon1 = -30;
          xlon2 = 60;
          xlat1 = 30;
          xlat2 = 80;
          lset = true;
        }
    }

  if (!lset)
    {
      operator_check_argc(argc_offset + 4);

      xlon1 = parameter_to_double(cdo_operator_argv(argc_offset + 0));
      xlon2 = parameter_to_double(cdo_operator_argv(argc_offset + 1));
      xlat1 = parameter_to_double(cdo_operator_argv(argc_offset + 2));
      xlat2 = parameter_to_double(cdo_operator_argv(argc_offset + 3));
    }
}

void
genlonlatbox(const int argc_offset, const int gridID, SelboxInfo &sbox)
{
  double xlon1 = 0, xlon2 = 0, xlat1 = 0, xlat2 = 0;
  getlonlatparams(argc_offset, xlon1, xlon2, xlat1, xlat2);

  const auto gridtype = gridInqType(gridID);

  if (gridtype == GRID_CURVILINEAR)
    genlonlatboxCurv(gridID, xlon1, xlon2, xlat1, xlat2, sbox);
  else
    genlonlatboxReg2d(gridID, xlon1, xlon2, xlat1, xlat2, sbox);
}

static int
genlonlatgrid(const int gridID1, SelboxInfo &sbox)
{
  genlonlatbox(0, gridID1, sbox);

  return gengrid(gridID1, sbox);
}

static int
gencellgrid(int gridID1, long &gridsize2, std::vector<long> &cellidx)
{
  const int argc_offset = 0;
  operator_check_argc(argc_offset + 4);

  auto xlon1 = parameter_to_double(cdo_operator_argv(argc_offset + 0));
  auto xlon2 = parameter_to_double(cdo_operator_argv(argc_offset + 1));
  auto xlat1 = parameter_to_double(cdo_operator_argv(argc_offset + 2));
  auto xlat2 = parameter_to_double(cdo_operator_argv(argc_offset + 3));

  if (xlon1 >= xlon2) std::swap(xlon1, xlon2);
  if (xlat1 >= xlat2) std::swap(xlat1, xlat2);

  const auto gridtype = gridInqType(gridID1);
  const auto gridsize1 = gridInqSize(gridID1);

  if (gridtype != GRID_UNSTRUCTURED) cdo_abort("Internal problem, wrong grid type!");

  if (!gridHasCoordinates(gridID1))
    {
      const auto reference = dereferenceGrid(gridID1);
      if (reference.isValid) gridID1 = reference.gridID;
      if (reference.notFound) cdo_abort("Reference to source grid not found!");
    }

  if (!gridHasCoordinates(gridID1)) cdo_abort("Grid cell coordinates missing!");

  Varray<double> xvals(gridsize1), yvals(gridsize1);
  gridInqXvals(gridID1, xvals.data());
  gridInqYvals(gridID1, yvals.data());

  // Convert lat/lon units if required
  cdo_grid_to_degree(gridID1, CDI_XAXIS, gridsize1, xvals.data(), "grid center lon");
  cdo_grid_to_degree(gridID1, CDI_YAXIS, gridsize1, yvals.data(), "grid center lat");

  // find gridsize2
  long maxcell = 0;
  long nvals = 0;
  for (size_t i = 0; i < gridsize1; ++i)
    {
      const auto xval = xvals[i];
      const auto yval = yvals[i];
      if (yval >= xlat1 && yval <= xlat2)
        if ((xval >= xlon1 && xval <= xlon2) || (xval + 360 >= xlon1 && xval + 360 <= xlon2)
            || (xval - 360 >= xlon1 && xval - 360 <= xlon2))
          {
            nvals++;
            if (nvals > maxcell)
              {
                constexpr long cellinc = 4096;
                maxcell += cellinc;
                cellidx.resize(maxcell);
              }
            cellidx[nvals - 1] = i;
          }
    }

  if (nvals == 0) cdo_abort("No grid points found!");

  gridsize2 = nvals;

  return gengridcell(gridID1, gridsize2, cellidx);
}

static void
checkBounds(const char *txt, const long minVal, const long maxVal, long &lval)
{
  if (lval < minVal)
    {
      cdo_warning("%s index out of range, set to %ld!", txt, minVal);
      lval = minVal;
    }

  if (lval > maxVal)
    {
      cdo_warning("%s index out of range, set to %ld!", txt, maxVal);
      lval = maxVal;
    }
}

void
genindexbox(const int argc_offset, const int gridID1, SelboxInfo &sbox)
{
  auto &lat1 = sbox.lat1;
  auto &lat2 = sbox.lat2;
  auto &lon11 = sbox.lon11;
  auto &lon12 = sbox.lon12;
  auto &lon21 = sbox.lon21;
  auto &lon22 = sbox.lon22;

  operator_check_argc(argc_offset + 4);

  lon11 = parameter_to_int(cdo_operator_argv(argc_offset + 0));
  lon12 = parameter_to_int(cdo_operator_argv(argc_offset + 1));
  lat1 = parameter_to_int(cdo_operator_argv(argc_offset + 2));
  lat2 = parameter_to_int(cdo_operator_argv(argc_offset + 3));

  if (lat1 > lat2) std::swap(lat1, lat2);

  const long nlon = gridInqXsize(gridID1);
  const long nlat = gridInqYsize(gridID1);

  checkBounds("First latitude", 1, nlat, lat1);
  checkBounds("Last latitude", 1, nlat, lat2);
  checkBounds("First longitude", 1, nlon, lon11);
  checkBounds("Last longitude", 1, nlon, lon12);

  lon11--;
  lon12--;
  lat1--;
  lat2--;

  if (lon11 > lon12)
    {
      lon21 = lon11;
      lon22 = nlon - 1;
      lon11 = 0;
    }
  else
    {
      if (lon12 > nlon - 1)
        {
          lon21 = lon11;
          lon22 = nlon - 1;
          lon11 = 0;
          lon12 = 0;
        }
      else
        {
          lon21 = 0;
          lon22 = -1;
        }
    }
}

static int
genindexgrid(const int gridID1, SelboxInfo &sbox)
{
  genindexbox(0, gridID1, sbox);

  if (gridInqType(gridID1) == GRID_PROJECTION && gridInqProjType(gridID1) == CDI_PROJ_LCC)
    return cdo_define_subgrid_grid(gridID1, sbox.lon11, sbox.lon12, sbox.lat1, sbox.lat2);
  else
    return gengrid(gridID1, sbox);
}

template <typename T>
static void
window_box(const int nwpv, const T *array1, const int gridID1, T *array2, const long lat1, const long lat2, const long lon11,
           const long lon12, const long lon21, const long lon22)
{
  const long nlon = gridInqXsize(gridID1);

  if (nwpv == 2)
    {
      for (long ilat = lat1; ilat <= lat2; ilat++)
        {
          for (long ilon = lon21; ilon <= lon22; ilon++)
            {
              *array2++ = array1[ilat * nlon * 2 + ilon * 2];
              *array2++ = array1[ilat * nlon * 2 + ilon * 2 + 1];
            }
          for (long ilon = lon11; ilon <= lon12; ilon++)
            {
              *array2++ = array1[ilat * nlon * 2 + ilon * 2];
              *array2++ = array1[ilat * nlon * 2 + ilon * 2 + 1];
            }
        }
    }
  else
    {
      for (long ilat = lat1; ilat <= lat2; ilat++)
        {
          for (long ilon = lon21; ilon <= lon22; ilon++) *array2++ = array1[ilat * nlon + ilon];
          for (long ilon = lon11; ilon <= lon12; ilon++) *array2++ = array1[ilat * nlon + ilon];
        }
    }
}

static void
window_box(const Field &field1, Field &field2, const long lat1, const long lat2, const long lon11, const long lon12,
           const long lon21, const long lon22)
{
  if (field1.memType == MemType::Float)
    window_box(field1.nwpv, field1.vec_f.data(), field1.grid, field2.vec_f.data(), lat1, lat2, lon11, lon12, lon21, lon22);
  else
    window_box(field1.nwpv, field1.vec_d.data(), field1.grid, field2.vec_d.data(), lat1, lat2, lon11, lon12, lon21, lon22);
}

template <typename T>
static void
window_cell(const int nwpv, const Varray<T> &array1, Varray<T> &array2, const long gridsize2, const std::vector<long> &cellidx)
{
  if (nwpv == 2)
    {
      for (long i = 0; i < gridsize2; ++i)
        {
          array2[i * 2] = array1[cellidx[i] * 2];
          array2[i * 2 + 1] = array1[cellidx[i] * 2 + 1];
        }
    }
  else
    {
      for (long i = 0; i < gridsize2; ++i) array2[i] = array1[cellidx[i]];
    }
}

void
window_cell(const Field &field1, Field &field2, const std::vector<long> &cellidx)
{
  if (field1.memType == MemType::Float)
    window_cell(field1.nwpv, field1.vec_f, field2.vec_f, field2.gridsize, cellidx);
  else
    window_cell(field1.nwpv, field1.vec_d, field2.vec_d, field2.gridsize, cellidx);
}

void *
Selbox(void *process)
{
  int varID, levelID;
  int index;

  cdo_initialize(process);

  // clang-format off
  const auto SELLONLATBOX = cdo_operator_add("sellonlatbox", 0, 0, "western and eastern longitude and southern and northern latitude");
  const auto SELINDEXBOX  = cdo_operator_add("selindexbox",  0, 0, "index of first and last longitude and index of first and last latitude");
  // clang-format on

  const auto operatorID = cdo_operator_id();

  operator_input_arg(cdo_operator_enter(operatorID));

  const auto streamID1 = cdo_open_read(0);

  const auto vlistID1 = cdo_stream_inq_vlist(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> vars(nvars, false);

  const auto ngrids = vlistNgrids(vlistID1);
  std::vector<SelboxInfo> sbox(ngrids);

  for (index = 0; index < ngrids; index++)
    {
      auto &sb = sbox[index];
      const auto gridID1 = vlistGrid(vlistID1, index);
      const auto gridtype = gridInqType(gridID1);
      const auto projtype = gridInqProjType(gridID1);

      const auto lprojection = (gridtype == GRID_PROJECTION && projtype == CDI_PROJ_RLL);
      //                    || (gridtype == GRID_PROJECTION && projtype == CDI_PROJ_LCC)
      //                    || (gridtype == GRID_PROJECTION && projtype == CDI_PROJ_STERE);
      if (gridtype == GRID_LONLAT || gridtype == GRID_GAUSSIAN || gridtype == GRID_CURVILINEAR || lprojection
          || (operatorID == SELINDEXBOX && (gridtype == GRID_GENERIC || gridtype == GRID_PROJECTION) && gridInqXsize(gridID1) > 0
              && gridInqYsize(gridID1) > 0)
          || (operatorID == SELLONLATBOX && gridtype == GRID_UNSTRUCTURED))
        {
          int gridID2;
          if (operatorID == SELLONLATBOX)
            {
              const auto gridsize = gridInqSize(gridID1);
              if (gridsize == 1) continue;

              if (gridtype == GRID_UNSTRUCTURED)
                gridID2 = gencellgrid(gridID1, sb.nvals, sb.cellidx);
              else
                gridID2 = genlonlatgrid(gridID1, sb);
            }
          else
            gridID2 = genindexgrid(gridID1, sb);

          sb.gridtype = gridtype;
          sb.gridID1 = gridID1;
          sb.gridID2 = gridID2;

          vlistChangeGridIndex(vlistID2, index, gridID2);

          for (varID = 0; varID < nvars; varID++)
            if (gridID1 == vlistInqVarGrid(vlistID1, varID)) vars[varID] = true;
        }
      else if (gridtype == GRID_GENERIC && gridInqXsize(gridID1) <= 1 && gridInqYsize(gridID1) <= 1)
        {
        }
      else
        {
          cdo_print("Unsupported grid type: %s", gridNamePtr(gridtype));
          if (gridtype == GRID_GAUSSIAN_REDUCED) cdo_print("Use option -R to convert Gaussian reduced grid to a regular grid!");
          cdo_abort("Unsupported grid type!");
        }
    }

  if (Options::cdoVerbose)
    {
      for (const auto &sb : sbox)
        if (sb.gridtype != GRID_UNSTRUCTURED)
          {
            cdo_print("box1 - idx1,idx2,idy1,idy2: %ld,%ld,%ld,%ld", sb.lon21 + 1, sb.lon22 + 1, sb.lat1 + 1, sb.lat2 + 1);
            cdo_print("box2 - idx1,idx2,idy1,idy2: %ld,%ld,%ld,%ld", sb.lon11 + 1, sb.lon12 + 1, sb.lat1 + 1, sb.lat2 + 1);
          }
    }

  for (varID = 0; varID < nvars; varID++)
    if (vars[varID]) break;

  if (varID >= nvars) cdo_warning("No variables selected!");

  const auto streamID2 = cdo_open_write(1);
  cdo_def_vlist(streamID2, vlistID2);

  VarList varList1, varList2;
  varListInit(varList1, vlistID1);
  varListInit(varList2, vlistID2);

  Field field1, field2;

  int tsID = 0;
  while (true)
    {
      const auto nrecs = cdo_stream_inq_timestep(streamID1, tsID);
      if (nrecs == 0) break;

      cdo_taxis_copy_timestep(taxisID2, taxisID1);
      cdo_def_timestep(streamID2, tsID);

      for (int recID = 0; recID < nrecs; recID++)
        {
          cdo_inq_record(streamID1, &varID, &levelID);
          field1.init(varList1[varID]);
          cdo_read_record(streamID1, field1);

          cdo_def_record(streamID2, varID, levelID);

          if (vars[varID])
            {
              field2.init(varList2[varID]);

              const auto gridID1 = varList1[varID].gridID;
              for (index = 0; index < ngrids; index++)
                if (gridID1 == sbox[index].gridID1) break;
              if (index == ngrids) cdo_abort("Internal problem, grid not found!");

              const auto &sb = sbox[index];

              if (operatorID == SELLONLATBOX && sb.gridtype == GRID_UNSTRUCTURED)
                window_cell(field1, field2, sb.cellidx);
              else
                window_box(field1, field2, sb.lat1, sb.lat2, sb.lon11, sb.lon12, sb.lon21, sb.lon22);

              if (field1.nmiss) field_num_mv(field2);

              cdo_write_record(streamID2, field2);
            }
          else
            {
              cdo_write_record(streamID2, field1);
            }
        }

      tsID++;
    }

  cdo_stream_close(streamID2);
  cdo_stream_close(streamID1);

  vlistDestroy(vlistID2);

  cdo_finish();

  return nullptr;
}