LonLatGrid.hh

Go to the documentation of this file.

/* Copyright (C) 2025 PISM Authors
 *
 * This file is part of PISM.
 *
 * PISM is free software; you can redistribute it and/or modify it under the
 * terms of the GNU General Public License as published by the Free Software
 * Foundation; either version 3 of the License, or (at your option) any later
 * version.
 *
 * PISM is distributed in the hope that it will be useful, but WITHOUT ANY
 * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
 * FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
 * details.
 *
 * You should have received a copy of the GNU General Public License
 * along with PISM; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
 */
 
#ifndef PISM_LONLATGRID_H
#define PISM_LONLATGRID_H
 
#include <vector>
#include <cmath>                // M_PI
#include <cstddef>              // size_t
 
#include "pism/util/LonLatCalculator.hh"
 
namespace pism {
/*!
 * Grid definition using coordinates in radians.
 */
struct LonLatGrid {
  std::vector<double> lon;
  std::vector<double> lat;
 
  /*!
   *
   * Converts a Cartesian grid in a `projection` that uses coordinates
   * `x` and `y` in meters into the form that can be used to define a
   * curvilinear grid in YAC.
   *
   * The `projection` string has to use the format compatible with PROJ.
   */
  LonLatGrid(const std::vector<double> &x, const std::vector<double> &y,
             const std::string &projection) {
 
    size_t nrow = y.size();
    size_t ncol = x.size();
    size_t N    = nrow * ncol;
 
    lon.resize(N);
    lat.resize(N);
 
    // convert from (row, col) to the linear index in "cell" arrays
    auto C = [ncol](size_t row, size_t col) { return row * ncol + col; };
 
    // convert from degrees to radians
    auto deg2rad = [](double degree) { return degree * M_PI / 180; };
 
    pism::LonLatCalculator mapping(projection);
 
    for (size_t row = 0; row < nrow; ++row) {
      for (size_t col = 0; col < ncol; ++col) {
        auto coords = mapping.lonlat(x[col], y[row]);
 
        lon[C(row, col)] = deg2rad(coords[0]);
        lat[C(row, col)] = deg2rad(coords[1]);
      }
    }
  }
};
 
} // namespace pism
 
#endif /* PISM_LONLATGRID_H */