cdo-2.0.0/src/Verifygrid.cc

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

  Author: Asela Rajapakse
          Uwe Schulzweida

*/

/*
   This module contains the following operators:
*/

#include <cdi.h>

#include <utility>

#include "cdo_options.h"
#include "process_int.h"
#include <mpim_grid.h>
#include "gridreference.h"
#include "verifygrid.h"

static constexpr double eps = 0.000000001;

static void
printIndex2D(const size_t cell_no, const size_t nx)
{
  const auto iy = cell_no / nx;
  const auto ix = cell_no - iy * nx;
  fprintf(stdout, " [i=%zu j=%zu]", ix + 1, iy + 1);
}

static void
quick_sort(double *array, size_t array_length)
{
  if (array_length < 2) return;
  const auto p = array[array_length / 2];

  size_t i, j;
  for (i = 0, j = array_length - 1;; i++, j--)
    {
      while (array[i] < p) i++;
      while (p < array[j]) j--;
      if (i >= j) break;
      std::swap(array[i], array[j]);
    }
  quick_sort(array, i);
  quick_sort(array + i, array_length - i);
}

/* Quicksort is called with a pointer to the array of center points to be sorted and an integer indicating its length.
 * It sorts the array by its longitude coordinates */
static void
quick_sort_by_lon(double *array, size_t array_length)
{
  if (array_length < 4) return;

  const auto p = ((array_length / 2) % 2) ? array[(array_length / 2) + 1] : array[array_length / 2];

  size_t i, j;
  for (i = 0, j = array_length - 2;; i += 2, j -= 2)
    {
      while (array[i] < p) i += 2;

      while (p < array[j]) j -= 2;

      if (i >= j) break;

      std::swap(array[i], array[j]);
      std::swap(array[i + 1], array[j + 1]);
    }

  quick_sort_by_lon(array, i);
  quick_sort_by_lon(array + i, array_length - i);
}

// This uses quicksort to sort the latitude coordinates in a subarray of all coordinates.
static void
quick_sort_of_subarray_by_lat(double *array, size_t subarray_start, size_t subarray_end)
{
  const auto subarray_length = (subarray_end - subarray_start) / 2 + 1;
  Varray<double> subarray(subarray_length);
  size_t subarray_index = 0;

  for (size_t index = subarray_start + 1; index <= subarray_end + 1; index += 2)
    {
      subarray[subarray_index] = array[index];
      subarray_index += 1;
    }

  quick_sort(subarray.data(), subarray_length);

  subarray_index = 0;

  for (size_t index = subarray_start + 1; index <= subarray_end + 1; index += 2)
    {
      array[index] = subarray[subarray_index];
      subarray_index += 1;
    }
}

static double
determinant(const double (&matrix)[3][3])
{
  // Calculates the determinant for a 3 x 3 matrix.

  return matrix[0][0] * matrix[1][1] * matrix[2][2] + matrix[0][1] * matrix[1][2] * matrix[2][0]
         + matrix[0][2] * matrix[1][0] * matrix[2][1] - matrix[0][2] * matrix[1][1] * matrix[2][0]
         - matrix[0][1] * matrix[1][0] * matrix[2][2] - matrix[0][0] * matrix[1][2] * matrix[2][1];
}

static void
find_unit_normal(const Point3D &a, const Point3D &b, const Point3D &c, Point3D &unit_normal)
{
  // Calculates the unit normal for a plane defined on three points a, b, c in Euclidean space.

  const double matrix_for_x[3][3] = { { 1, a.Y, a.Z }, { 1, b.Y, b.Z }, { 1, c.Y, c.Z } };
  const double matrix_for_y[3][3] = { { a.X, 1, a.Z }, { b.X, 1, b.Z }, { c.X, 1, c.Z } };
  const double matrix_for_z[3][3] = { { a.X, a.Y, 1 }, { b.X, b.Y, 1 }, { c.X, c.Y, 1 } };

  const auto x = determinant(matrix_for_x);
  const auto y = determinant(matrix_for_y);
  const auto z = determinant(matrix_for_z);

  const auto magnitude = std::sqrt(x * x + y * y + z * z);

  unit_normal.X = x / magnitude;
  unit_normal.Y = y / magnitude;
  unit_normal.Z = z / magnitude;
}

int
find_coordinate_to_ignore(const Varray<Point3D> &cell_corners_xyz)
{
  // Takes the first three corners/vertices of the cell and calculates the unit normal via determinants.

  const auto &pcorner_coordinates1 = cell_corners_xyz[0];
  const auto &pcorner_coordinates2 = cell_corners_xyz[1];
  const auto &pcorner_coordinates3 = cell_corners_xyz[2];

  Point3D surface_normal_of_the_cell;
  find_unit_normal(pcorner_coordinates1, pcorner_coordinates2, pcorner_coordinates3, surface_normal_of_the_cell);

  // The surface normal is used to choose the coordinate to ignore.

  const auto abs_x = std::fabs(surface_normal_of_the_cell.X);
  const auto abs_y = std::fabs(surface_normal_of_the_cell.Y);
  const auto abs_z = std::fabs(surface_normal_of_the_cell.Z);

  int coordinate_to_ignore = 3;

  if (abs_x > abs_y)
    {
      if (abs_x > abs_z) coordinate_to_ignore = 1;
    }
  else
    {
      if (abs_y > abs_z) coordinate_to_ignore = 2;
    }

  return coordinate_to_ignore;
}

static inline double
is_point_left_of_edge(const Point &point1, const Point &point2, const Point &point)
{
  /*
     Computes whether a point is left of the line through point1 and point2.
     This is part of the solution to the point in polygon problem.

     Returns:  0 if the point is on the line through point1 and point2
             > 0 if the point is left of the line
             < 0 if the point is right of the line

     This algorithm is by Dan Sunday (geomalgorithms.com) and is completely free for use and modification.
  */

  return ((point2.x - point1.x) * (point.y - point1.y) - (point.x - point1.x) * (point2.y - point1.y));
}

int
winding_numbers_algorithm(const Varray<Point> &cell_corners, int number_corners, const Point &point)
{
  /*
     Computes whether a point is inside the bounds of a cell. This is the solution to the point in polygon problem.
     Returns 0 if the point is outside, returns 1 if the point is inside the cell. Based on an algorithm by Dan Sunday
     (geomalgorithms.com). His algorithm is completely free for use and modification.
  */

  int winding_number = 0;

  for (int k = 0; k < number_corners - 1; k++)
    {
      if (cell_corners[k].y <= point.y)
        {
          if (cell_corners[k + 1].y > point.y)
            {
              if (is_point_left_of_edge(cell_corners[k], cell_corners[k + 1], point) > 0) winding_number++;
            }
        }
      else
        {
          if (cell_corners[k + 1].y <= point.y)
            {
              if (is_point_left_of_edge(cell_corners[k], cell_corners[k + 1], point) < 0) winding_number--;
            }
        }
    }

  return winding_number;
}

static double
sign(double x)
{
  // Is +1 if x is positive, -1 if x is negative and 0 if x is zero.

  return (x > 0.0) - (x < 0.0);
}

static bool
is_simple_polygon_convex(const Varray<Point> &cell_corners, int number_corners)
{
  // Tests in which direction the polygon winds when walking along its edges. Does so for all edges of the polygon.

  double direction = 0.0;

  for (int i = 0; i < number_corners - 2; i++)
    {
      auto turns_to = (cell_corners[i].x - cell_corners[i + 1].x) * (cell_corners[i + 1].y - cell_corners[i + 2].y)
                      - (cell_corners[i].y - cell_corners[i + 1].y) * (cell_corners[i + 1].x - cell_corners[i + 2].x);

      // In the first iteration the direction of winding of the entire polygon is set. Better not be 0.

      if (i == 1) direction = turns_to;

      if (IS_NOT_EQUAL(sign(direction), sign(turns_to)))
        {
          if (IS_NOT_EQUAL(direction, 0.0)) return false;
        }
      else
        {
          direction = turns_to;
        }
    }

  return true;
}

double
calculate_the_polygon_area(const Varray<Point> &cell_corners, int number_corners)
{
  /* This algorithm is based on the calculation from Wolfram Mathworld Polygon Area. It results in the area of planar
   * non-self-intersecting polygon. */

  double twice_the_polygon_area = 0.0;

  for (int i = 0; i < number_corners - 1; i++)
    twice_the_polygon_area += (cell_corners[i].x * cell_corners[i + 1].y) - (cell_corners[i + 1].x * cell_corners[i].y);

  return twice_the_polygon_area / 2.0;
}

bool
are_polygon_vertices_arranged_in_clockwise_order(double cell_area)
{
  /* A negative area indicates a clockwise arrangement of vertices, a positive area a counterclockwise arrangement.
   * There should be an area to begin with.
   */
  return cell_area < 0.0;
}

static void
create_sorted_point_array(size_t gridsize, const Varray<double> &grid_center_lon, const Varray<double> &grid_center_lat,
                          Varray<double> &center_point_array)
{
  for (size_t cell_no = 0; cell_no < gridsize; cell_no++)
    {
      center_point_array[cell_no * 2 + 0] = grid_center_lon[cell_no];
      center_point_array[cell_no * 2 + 1] = grid_center_lat[cell_no];
    }

  // The cell center points are sorted by their first coordinate (lon) with quicksort.

  quick_sort_by_lon(center_point_array.data(), gridsize * 2);

  // Now the lat coordinates in subarrays that reflect equal lon coordinates are sorted with quicksort.

  int subarray_start = 0;
  int subarray_end = 0;

  for (size_t cell_no = 0; cell_no < gridsize - 1; cell_no++)
    {
      if (cell_no == gridsize - 2)
        {
          subarray_end = gridsize * 2 - 2;
          quick_sort_of_subarray_by_lat(center_point_array.data(), subarray_start, subarray_end);
        }

      if (std::fabs(center_point_array[cell_no * 2 + 0] - center_point_array[(cell_no + 1) * 2 + 0]) > eps)
        {
          subarray_end = cell_no * 2;
          if ((subarray_end - subarray_start) > 1)
            quick_sort_of_subarray_by_lat(center_point_array.data(), subarray_start, subarray_end);

          subarray_start = subarray_end + 2;
        }
    }
}

static void
print_header(int gridtype, size_t gridsize, size_t nx, int gridno, int ngrids)
{
  if (ngrids == 1)
    {
      if (nx)
        cdo_print(Blue("Grid consists of %zu (%zux%zu) cells (type: %s), of which"), gridsize, nx, gridsize / nx,
                  gridNamePtr(gridtype));
      else
        cdo_print(Blue("Grid consists of %zu cells (type: %s), of which"), gridsize, gridNamePtr(gridtype));
    }
  else
    {
      if (nx)
        cdo_print(Blue("Grid no %u (of %u) consists of %zu (%zux%zu) cells (type: %s), of which"), gridno + 1, ngrids, gridsize, nx,
                  gridsize / nx, gridNamePtr(gridtype));
      else
        cdo_print(Blue("Grid no %u (of %u) consists of %zu cells (type: %s), of which"), gridno + 1, ngrids, gridsize,
                  gridNamePtr(gridtype));
    }
}

static size_t
get_no_unique_center_points(size_t gridsize, size_t nx, const Varray<double> &grid_center_lon,
                            const Varray<double> &grid_center_lat)
{
  (void)nx;
  // For performing the first test, an array of all center point coordinates is built.
  Varray<double> center_points(gridsize * 2);
  create_sorted_point_array(gridsize, grid_center_lon, grid_center_lat, center_points);

  size_t no_unique_center_points = 1;
  for (size_t cell_no = 0; cell_no < gridsize - 1; cell_no++)
    {
      if (std::fabs(center_points[cell_no * 2 + 0] - center_points[(cell_no + 1) * 2 + 0]) < eps
          && std::fabs(center_points[cell_no * 2 + 1] - center_points[(cell_no + 1) * 2 + 1]) < eps)
        {
          if (Options::cdoVerbose)
            {
              fprintf(stdout, "Duplicate point [lon=%.5g lat=%.5g] was found", center_points[cell_no * 2 + 0],
                      center_points[cell_no * 2 + 1]);
              /*
              fprintf(stdout, "Duplicate point [lon=%.5g lat=%.5g] was found in cell no %zu",
                      center_points[cell_no * 2 + 0], center_points[cell_no * 2 + 1], cell_no + 1);
              if (nx) printIndex2D(cell_no, nx);
              */
              fprintf(stdout, "\n");
            }
        }
      else
        {
          no_unique_center_points++;
        }
    }

  return no_unique_center_points;
}

void
set_cell_corners_xyz(int ncorner, const double *cell_corners_lon, const double *cell_corners_lat, Varray<Point3D> &cell_corners_xyz)
{
  double corner_coordinates[3];
  for (int k = 0; k < ncorner; k++)
    {
      // Conversion of corner spherical coordinates to Cartesian coordinates.
      gcLLtoXYZ(cell_corners_lon[k], cell_corners_lat[k], corner_coordinates);

      // The components of the result vector are appended to the list of cell corner coordinates.
      cell_corners_xyz[k].X = corner_coordinates[0];
      cell_corners_xyz[k].Y = corner_coordinates[1];
      cell_corners_xyz[k].Z = corner_coordinates[2];
    }
}

void
set_center_point_plane_projection(int coordinate_to_ignore, const Point3D &center_point_xyz, Point &center_point_plane_projection)
{
  // clang-format off
  switch (coordinate_to_ignore)
    {
    case 1: center_point_plane_projection = Point {center_point_xyz.Y, center_point_xyz.Z}; break;
    case 2: center_point_plane_projection = Point {center_point_xyz.Z, center_point_xyz.X}; break;
    case 3: center_point_plane_projection = Point {center_point_xyz.X, center_point_xyz.Y}; break;
    }
  // clang-format on
}

void
set_cell_corners_plane_projection(int coordinate_to_ignore, int ncorner, const Varray<Point3D> &cell_corners_xyz,
                                  Varray<Point> &cell_corners_plane)
{
  switch (coordinate_to_ignore)
    {
    case 1:
      for (int k = 0; k <= ncorner; k++) cell_corners_plane[k] = Point{ cell_corners_xyz[k].Y, cell_corners_xyz[k].Z };
      break;
    case 2:
      for (int k = 0; k <= ncorner; k++) cell_corners_plane[k] = Point{ cell_corners_xyz[k].Z, cell_corners_xyz[k].X };
      break;
    case 3:
      for (int k = 0; k <= ncorner; k++) cell_corners_plane[k] = Point{ cell_corners_xyz[k].X, cell_corners_xyz[k].Y };
      break;
    }
}

int
get_actual_number_of_corners(int ncorner, const Varray<Point3D> &cell_corners_xyz_open_cell)
{
  auto actual_number_of_corners = ncorner;

  for (int k = ncorner - 1; k > 0; k--)
    {
      if (IS_EQUAL(cell_corners_xyz_open_cell[k].X, cell_corners_xyz_open_cell[k - 1].X)
          && IS_EQUAL(cell_corners_xyz_open_cell[k].Y, cell_corners_xyz_open_cell[k - 1].Y)
          && IS_EQUAL(cell_corners_xyz_open_cell[k].Z, cell_corners_xyz_open_cell[k - 1].Z))
        actual_number_of_corners--;
      else
        break;
    }

  return actual_number_of_corners;
}

int
get_no_duplicates(int actual_number_of_corners, const Varray<Point3D> &cell_corners_xyz_open_cell,
                  std::vector<bool> &marked_duplicate_indices)
{
  for (int i = 0; i < actual_number_of_corners; i++) marked_duplicate_indices[i] = false;

  int no_duplicates = 0;

  for (int i = 0; i < actual_number_of_corners; i++)
    for (int j = i + 1; j < actual_number_of_corners; j++)
      if (std::fabs(cell_corners_xyz_open_cell[i].X - cell_corners_xyz_open_cell[j].X) < eps
          && std::fabs(cell_corners_xyz_open_cell[i].Y - cell_corners_xyz_open_cell[j].Y) < eps
          && std::fabs(cell_corners_xyz_open_cell[i].Z - cell_corners_xyz_open_cell[j].Z) < eps)
        {
          no_duplicates++;
          marked_duplicate_indices[j] = true;
        }

  return no_duplicates;
}

void
copy_unique_corners(int actual_number_of_corners, const Varray<Point3D> &cell_corners_xyz_open_cell,
                    const std::vector<bool> &marked_duplicate_indices, Varray<Point3D> &cell_corners_xyz_without_duplicates)
{
  int unique_corner_number = 0;

  for (int k = 0; k < actual_number_of_corners; k++)
    {
      if (marked_duplicate_indices[k] == false)
        {
          cell_corners_xyz_without_duplicates[unique_corner_number] = cell_corners_xyz_open_cell[k];
          unique_corner_number++;
        }
    }
}

static void
verify_grid(size_t gridsize, size_t nx, int ncorner, const Varray<double> &grid_center_lon, const Varray<double> &grid_center_lat,
            const Varray<double> &grid_corner_lon, const Varray<double> &grid_corner_lat)
{
  /*
     First, this function performs the following test:

     1) it tests whether there are duplicate cells in the given grid by comparing their center point

     Additionally, on each cell of a given grid:

     2) it tests whether all cells are convex and all cell bounds have the same orientation,
        i.e. the corners of the cell are in clockwise or counterclockwise order

     3) it tests whether the center point is within the bounds of the cell

     The results of the tests are printed on stdout.
  */
  Point3D center_point_xyz;
  Varray<Point3D> cell_corners_xyz_open_cell(ncorner);

  Point center_point_plane_projection;

  size_t no_of_cells_with_duplicates = 0;
  size_t no_usable_cells = 0;
  size_t no_convex_cells = 0;
  size_t no_clockwise_cells = 0;
  size_t no_counterclockwise_cells = 0;
  size_t no_of_cells_with_center_points_out_of_bounds = 0;

  std::vector<int> no_cells_with_a_specific_no_of_corners(ncorner, 0);

  // Checking for the number of unique center point coordinates.
  const auto no_unique_center_points = get_no_unique_center_points(gridsize, nx, grid_center_lon, grid_center_lat);

  // used only actual_number_of_corners
  std::vector<bool> marked_duplicate_indices(ncorner);
  Varray<Point3D> cell_corners_xyz(ncorner + 1);
  Varray<Point> cell_corners_plane_projection(ncorner + 1);
  Varray<size_t> cells_with_duplicates;
  cells_with_duplicates.reserve(gridsize);

  /*
     Latitude and longitude are spherical coordinates on a unit circle. Each such coordinate tuple is transformed into a
     triple of Cartesian coordinates in Euclidean space. This is first done for the presumed center point of the cell
     and then for all the corners of the cell.
  */

  for (size_t cell_no = 0; cell_no < gridsize; cell_no++)
    {
      // Conversion of center point spherical coordinates to Cartesian coordinates.
      double center_coordinates[3];
      gcLLtoXYZdeg(grid_center_lon[cell_no], grid_center_lat[cell_no], center_coordinates);
      center_point_xyz.X = center_coordinates[0];
      center_point_xyz.Y = center_coordinates[1];
      center_point_xyz.Z = center_coordinates[2];

      double corner_coordinates[3];
      for (int k = 0; k < ncorner; k++)
        {
          // Conversion of corner spherical coordinates to Cartesian coordinates.
          gcLLtoXYZdeg(grid_corner_lon[cell_no * ncorner + k], grid_corner_lat[cell_no * ncorner + k], corner_coordinates);

          // The components of the result vector are appended to the list of cell corner coordinates.
          cell_corners_xyz_open_cell[k].X = corner_coordinates[0];
          cell_corners_xyz_open_cell[k].Y = corner_coordinates[1];
          cell_corners_xyz_open_cell[k].Z = corner_coordinates[2];
        }

      /*
         Not all cells have the same number of corners. The array, however, has ncorner * 3  values for each cell, where
         ncorner is the maximum number of corners. Unused values have been filled with the values of the final cell. The
         following identifies the surplus corners and gives the correct length of the cell.
      */

      auto actual_number_of_corners = get_actual_number_of_corners(ncorner, cell_corners_xyz_open_cell);

      no_cells_with_a_specific_no_of_corners[actual_number_of_corners - 1] += 1;

      // If there are less than three corners in the cell, it is unusable and considered degenerate. No area can be computed.

      if (actual_number_of_corners < 3)
        {
          if (Options::cdoVerbose)
            {
              fprintf(stdout, "Less than three vertices found in cell no %zu", cell_no + 1);
              if (nx) printIndex2D(cell_no, nx);
              fprintf(stdout, ", omitted!");
              fprintf(stdout, "\n");
            }
          continue;
        }

      no_usable_cells++;

      // Checks if there are any duplicate vertices in the list of corners. Note that the last (additional) corner has not been set
      // yet.

      auto no_duplicates = get_no_duplicates(actual_number_of_corners, cell_corners_xyz_open_cell, marked_duplicate_indices);

      if (Options::cdoVerbose)
        {
          for (int i = 0; i < actual_number_of_corners; i++)
            {
              if (marked_duplicate_indices[i])
                {
                  fprintf(stdout, "Duplicate vertex [lon=%.5g lat=%.5g] was found in cell no %zu",
                          grid_corner_lon[cell_no * ncorner + i], grid_corner_lat[cell_no * ncorner + i], cell_no + 1);
                  if (nx) printIndex2D(cell_no, nx);
                  fprintf(stdout, "\n");
                }
            }
        }

      // Writes the unique corner vertices in a new array.

      copy_unique_corners(actual_number_of_corners, cell_corners_xyz_open_cell, marked_duplicate_indices, cell_corners_xyz);

      actual_number_of_corners -= no_duplicates;

      // We are creating a closed polygon/cell by setting the additional last corner to be the same as the first one.

      cell_corners_xyz[actual_number_of_corners] = cell_corners_xyz[0];

      if (no_duplicates != 0)
        {
          no_of_cells_with_duplicates += 1;
          cells_with_duplicates.push_back(cell_no);
        }

      /* If there are less than three corners in the cell left after removing duplicates, it is unusable and considered
       * degenerate. No area can be computed. */

      if (actual_number_of_corners < 3)
        {
          if (Options::cdoVerbose)
            fprintf(stdout,
                    "Less than three vertices found in cell no %zu. This cell is considered degenerate and "
                    "will be omitted from further computation!\n",
                    cell_no + 1);

          continue;
        }

      const auto coordinate_to_ignore = find_coordinate_to_ignore(cell_corners_xyz);

      /* The remaining two-dimensional coordinates are extracted into one array for all the cell's corners and into one
       * array for the center point. */

      /* The following projection on the plane that two coordinate axes lie on changes the arrangement of the polygon
         vertices if the coordinate to be ignored along the third axis is smaller than 0. In this case, the result of
         the computation of the orientation of vertices needs to be  inverted. Clockwise becomes counterclockwise and
         vice versa. */

      bool invert_result = false;
      const auto cval = (coordinate_to_ignore == 1) ? cell_corners_xyz[0].X
                                                    : ((coordinate_to_ignore == 2) ? cell_corners_xyz[0].Y : cell_corners_xyz[0].Z);
      if (cval < 0.0) invert_result = true;

      set_center_point_plane_projection(coordinate_to_ignore, center_point_xyz, center_point_plane_projection);
      set_cell_corners_plane_projection(coordinate_to_ignore, actual_number_of_corners, cell_corners_xyz,
                                        cell_corners_plane_projection);

      // Checking for convexity of the cell.

      if (is_simple_polygon_convex(cell_corners_plane_projection, actual_number_of_corners + 1))
        {
          no_convex_cells += 1;
        }
      else if (Options::cdoVerbose)
        {
          fprintf(stdout, "Vertices are not convex in cell no %zu", cell_no + 1);
          if (nx) printIndex2D(cell_no, nx);
          fprintf(stdout, "\n");
        }

      // Checking the arrangement or direction of cell vertices.

      const auto polygon_area = calculate_the_polygon_area(cell_corners_plane_projection, actual_number_of_corners + 1);
      auto is_clockwise = are_polygon_vertices_arranged_in_clockwise_order(polygon_area);

      /* If the direction of the vertices was flipped during the projection onto the two-dimensional plane, the previous
       * result needs to be inverted now. */

      if (invert_result) is_clockwise = !is_clockwise;

      // The overall counter of (counter)clockwise cells is increased by one.
      if (is_clockwise)
        {
          no_clockwise_cells += 1;
          if (Options::cdoVerbose)
            {
              fprintf(stdout, "Vertices arranged in a clockwise order in cell no %zu", cell_no + 1);
              if (nx) printIndex2D(cell_no, nx);
              fprintf(stdout, "\n");
            }
        }
      else
        {
          no_counterclockwise_cells += 1;
        }

      // The winding numbers algorithm is used to test whether the presumed center point is within the bounds of the cell.
      auto winding_number
          = winding_numbers_algorithm(cell_corners_plane_projection, actual_number_of_corners + 1, center_point_plane_projection);

      // if ( winding_number == 0 ) printf("%d,", cell_no+1);
      if (winding_number == 0) no_of_cells_with_center_points_out_of_bounds += 1;

      if (Options::cdoVerbose && winding_number == 0)
        {
          fprintf(stdout, "Center point [lon=%.5g lat=%.5g] outside bounds [", grid_center_lon[cell_no], grid_center_lat[cell_no]);
          for (int k = 0; k < ncorner; k++)
            fprintf(stdout, " %.5g/%.5g", grid_corner_lon[cell_no * ncorner + k], grid_corner_lat[cell_no * ncorner + k]);
          fprintf(stdout, "] in cell no %zu", cell_no + 1);
          if (nx) printIndex2D(cell_no, nx);
          fprintf(stdout, "\n");
        }
    }

  const auto no_nonunique_cells = gridsize - no_unique_center_points;
  const auto no_nonconvex_cells = gridsize - no_convex_cells;
  const auto no_nonusable_cells = gridsize - no_usable_cells;

  for (int i = 2; i < ncorner; i++)
    if (no_cells_with_a_specific_no_of_corners[i])
      cdo_print(Blue("%9d cells have %d vertices"), no_cells_with_a_specific_no_of_corners[i], i + 1);

  if (no_of_cells_with_duplicates) cdo_print(Blue("%9zu cells have duplicate vertices"), no_of_cells_with_duplicates);

  if (no_nonusable_cells) cdo_print(Red("%9zu cells have unusable vertices"), no_nonusable_cells);

  if (no_nonunique_cells) cdo_print(Red("%9zu cells are not unique"), no_nonunique_cells);

  if (no_nonconvex_cells) cdo_print(Red("%9zu cells are non-convex"), no_nonconvex_cells);

  if (no_clockwise_cells) cdo_print(Red("%9zu cells have their vertices arranged in a clockwise order"), no_clockwise_cells);

  if (no_of_cells_with_center_points_out_of_bounds)
    cdo_print(Red("%9zu cells have their center point located outside their boundaries"),
              no_of_cells_with_center_points_out_of_bounds);

  // cdo_print("");
}

static void
print_lonlat(int dig, int gridID, size_t gridsize, const Varray<double> &lon, const Varray<double> &lat)
{
  size_t num_lons_nan = 0;
  for (size_t i = 0; i < gridsize; ++i)
    if (std::isnan(lon[i])) num_lons_nan++;
  if (num_lons_nan) cdo_print(Red("%9zu longitudes are NaN"), num_lons_nan);

  size_t num_lats_nan = 0;
  for (size_t i = 0; i < gridsize; ++i)
    if (std::isnan(lat[i])) num_lats_nan++;
  if (num_lats_nan) cdo_print(Red("%9zu latitudes are NaN"), num_lats_nan);

  char name[CDI_MAX_NAME];
  int length = CDI_MAX_NAME;
  cdiInqKeyString(gridID, CDI_XAXIS, CDI_KEY_NAME, name, &length);
  const auto xmm = varray_min_max(gridsize, lon);
  cdo_print("%10s : %.*g to %.*g degrees", name, dig, xmm.min, dig, xmm.max);
  if (xmm.min < -720.0 || xmm.max > 720.0) cdo_warning("Grid cell center longitudes out of range!");

  length = CDI_MAX_NAME;
  cdiInqKeyString(gridID, CDI_YAXIS, CDI_KEY_NAME, name, &length);
  const auto ymm = varray_min_max(gridsize, lat);
  cdo_print("%10s : %.*g to %.*g degrees", name, dig, ymm.min, dig, ymm.max);
  if (ymm.min < -90.00001 || ymm.max > 90.00001) cdo_warning("Grid cell center latitudes out of range!");
}

void *
Verifygrid(void *argument)
{
  cdo_initialize(argument);

  cdo_operator_add("verifygrid", 0, 0, nullptr);

  operator_check_argc(0);

  const auto streamID = cdo_open_read(0);
  const auto vlistID = cdo_stream_inq_vlist(streamID);

  const auto ngrids = vlistNgrids(vlistID);
  for (int gridno = 0; gridno < ngrids; ++gridno)
    {
      bool lgeo = true;

      auto gridID = vlistGrid(vlistID, gridno);
      const auto gridtype = gridInqType(gridID);

      if (gridtype == GRID_GENERIC || gridtype == GRID_SPECTRAL) lgeo = false;

      if (gridtype == GRID_GME) gridID = gridToUnstructured(gridID, 1);

      if (lgeo && gridtype != GRID_UNSTRUCTURED && gridtype != GRID_CURVILINEAR) gridID = gridToCurvilinear(gridID, 1);

      if (gridtype == GRID_UNSTRUCTURED && !gridHasCoordinates(gridID))
        {
          const auto reference = dereferenceGrid(gridID);
          if (reference.isValid) gridID = reference.gridID;
          if (reference.notFound)
            {
              cdo_print("Reference to source grid not found!");
              lgeo = false;
            }
        }

      const auto gridsize = gridInqSize(gridID);
      auto nx = gridInqXsize(gridID);
      auto ny = gridInqYsize(gridID);
      nx = (nx * ny == gridsize) ? nx : 0;

      if (lgeo)
        {
          const auto ncorner = (gridInqType(gridID) == GRID_UNSTRUCTURED) ? gridInqNvertex(gridID) : 4;

          if (!gridHasCoordinates(gridID)) cdo_abort("Cell center coordinates missing!");

          Varray<double> grid_center_lat(gridsize), grid_center_lon(gridsize);

          gridInqYvals(gridID, grid_center_lat.data());
          gridInqXvals(gridID, grid_center_lon.data());

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

          print_header(gridtype, gridsize, nx, gridno, ngrids);

          if (gridHasBounds(gridID))
            {
              if (ncorner == 0) cdo_abort("Number of cell corners undefined!");
              const auto nalloc = ncorner * gridsize;
              Varray<double> grid_corner_lat(nalloc), grid_corner_lon(nalloc);

              gridInqYbounds(gridID, grid_corner_lat.data());
              gridInqXbounds(gridID, grid_corner_lon.data());

              cdo_grid_to_degree(gridID, CDI_XAXIS, ncorner * gridsize, grid_corner_lon.data(), "grid corner lon");
              cdo_grid_to_degree(gridID, CDI_YAXIS, ncorner * gridsize, grid_corner_lat.data(), "grid corner lat");

              verify_grid(gridsize, nx, ncorner, grid_center_lon, grid_center_lat, grid_corner_lon, grid_corner_lat);
            }
          else
            {
              cdo_warning("Grid cell corner coordinates missing!");
            }

          const int dig = Options::CDO_flt_digits;
          print_lonlat(dig, gridID, gridsize, grid_center_lon, grid_center_lat);
        }
      else
        {
          if (ngrids == 1)
            cdo_print(Blue("Grid consists of %zu points (type: %s)"), gridsize, gridNamePtr(gridtype));
          else
            cdo_print(Blue("Grid no %u (of %u) consists of %zu points (type: %s)"), gridno + 1, ngrids, gridsize,
                      gridNamePtr(gridtype));
          // cdo_print("");
        }
    }

  cdo_stream_close(streamID);

  cdo_finish();

  return nullptr;
}