gfsfile.cc

#include "radixio/gfsfile.hh"

#include "radixbug/bug.hh"
#include "radixio/eafstream.hh"
#include "radixmath/util.hh"

namespace radix
{
float total_seconds(int year, int month, int day, int hour)
{
  // assume 1900 as start of time
  float start = 365 * 86400;
  // years = hours*days*years
  return year * 365.f * 86400.f + (month - 1.f) * (365.f / 12.f) * 86400.f +
         (day - 1.f) * 86400.f + hour * 3600.f - start;
}  // total_seconds
int ord(char c) { return (unsigned char)c; }
std::vector<std::vector<float>> pakinp(const std::string &cvar, int nx, int ny,
                                       int nx1, int ny1, int lx, int ly,
                                       float prec, int nexp, float var1)
{
  int k, jj, ii;
  float rnew;
  float rold  = var1;
  float scexp = 1.0f / std::pow(2.0f, float(7 - nexp));  // scaling exponent
  std::vector<std::vector<float>> rvar(nx);
  for (size_t i = 0; i < rvar.size(); ++i)
    rvar[i] = std::vector<float>(ny, 0.0);
  // initialize column 1 data
  for (int j = 0; j < ny; ++j)
  {
    k    = j * nx;  // position at column 1
    jj   = j - ny1;
    rnew = (float(ord(cvar[k]) - 127) * scexp) + rold;
    rold = rnew;
    if (jj >= 0 && jj <= ly)
    {
      rvar[0][jj] = rnew;
    }
  }  // 1st for j < ny
  for (int j = ny1; j < (ny1 + ly); ++j)
  {
    jj   = j - ny1;  // sub-grid array (1 to ly)
    rold = rvar[0][jj];
    for (int i = 1; i < (nx1 + lx); ++i)
    {
      k    = j * nx + i;
      rnew = (float(ord(cvar[k]) - 127) * scexp) + rold;
      rold = rnew;
      ii   = i - nx1;
      if (std::abs(rnew) < prec) rnew = 0.0f;
      if (ii >= 0 && ii <= lx)
      {
        rvar[ii][jj] = rnew;
      }
    }  // for i < (ny1+ly)
  }    // 2nd for j < ny
  return rvar;
}  // pakinp

std::vector<std::string> GFSFile::mVarb = {
    "    ", "PRSS", "TPPA", "TPPT", "TPP6", "PRT6", "TPP1", "CPP1", "TPP3",
    "CPP3", "MSLP", "SHGT", "U10M", "V10M", "RH2M", "DP2M", "MXHT", "VSBY",
    "T02M", "LHTF", "SHTF", "USTR", "RGHS", "DSWF", "UWND", "VWND", "WWND",
    "SPHU", "TEMP", "RELH", "HGTS", "TKEN", "TMPS", "SOLT", "SOLW", "P10M",
    "LCLD", "MCLD", "HCLD", "TCLD", "PBLH", "THET", "DZDT", "PRT3"};
std::vector<std::string> GFSFile::mUnits = {
    "    ", " hPa", "  mm", "  mm", "  mm", "mm/h", "  mm", "  mm", "  mm",
    "  mm", " hPa", "   m", " m/s", " m/s", "   %", "  oC", "   m", "  km",
    "  oC", "W/m2", "W/m2", "cm/s", "   m", "W/m2", " m/s", " m/s", "mb/h",
    "g/kg", "  oC", "   %", "   m", "Joul", "  oC", "  oK", "kgm2", "  oK",
    "   %", "   %", "   %", "   %", "   m", "  oC", " m/h", "mm/h"};
std::vector<float> GFSFile::mFact = {
    1.0f,   1.0f, 1000.f, 1000.f, 1000.f, 60000.f, 1000.f, 1000.f, 1000.f,
    1000.f, 1.0f, 1.0f,   1.0f,   1.0f,   1.0f,    1.0f,   1.0f,   0.001f,
    1.0f,   1.0f, 1.0f,   100.f,  1.0f,   1.0f,    1.0f,   1.0f,   3600.f,
    1000.f, 1.0f, 1.0f,   1.0f,   1.0f,   1.0f,    1.0f,   1.0f,   1.0f,
    1.0f,   1.0f, 1.0f,   1.0f,   1.0f,   1.0f,    3600.f, 60000.f};
GFSFile::GFSFile(std::string file)
    : mFile(file)
    , mStrcmp(9, 0.0f)
{
  radix_tagged_line("GFSFile(" << file << ")");
  // initialize data structure
  eafstream *rstream =
      new eafstream(file.c_str(), std::ifstream::in | std::ifstream::binary);
  std::string label  = rstream->readString(50);
  std::string header = rstream->readString(108);
  // initialize
  mLabel.expand(label);
  mHeader.expand(header, mLabel);
  mHeader.latlon = false;
  mHeader.global = false;
  mHeader.gbldat = false;
  mHeader.prime  = false;

  // calculate length of records
  int ldat    = mHeader.nx * mHeader.ny;
  int rec_len = ldat + 50;
  mLrec       = rec_len;
  int nndx    = mHeader.lenh / ldat + 1;
  // rewind to beginning of file
  rstream->seekg(0, rstream->beg);

  // loop over remaining index records
  for (int i = 0; i < nndx; ++i)
  {
    std::string recl = rstream->readString(mLrec);
    label            = recl.substr(0, 50);
    header           = recl.substr(50);
    mLabel.expand(label);
    mHeader.expand(header, mLabel);
    radix_tagged_line("Found grid: " << mHeader.toString());
  }
  //
  // determine if this is a lat-lon grid
  if (mHeader.size == 0.f)
  {
    mHeader.latlon = true;
  }
  if (mHeader.model_id.compare("RAMS") == 0)
  {
    mHeader.tang_lat = mHeader.pole_lat;
  }
  if (!mHeader.latlon)
  {
    //
    // initialize grid conversion variable (into gbase)
    stlmbr(mHeader.tang_lat, mHeader.ref_lon);
    //
    // use single point grid definition
    radix_line("sync_xp=" << mHeader.sync_xp << " sync_yp=" << mHeader.sync_yp
                          << " sync_lat=" << mHeader.sync_lat << " sync_lon="
                          << mHeader.sync_lon << " ref_lat=" << mHeader.ref_lat
                          << " ref_lon=" << mHeader.ref_lon << " size="
                          << mHeader.size << " orient=" << mHeader.orient);
    stcm1p(mHeader.sync_xp, mHeader.sync_yp, mHeader.sync_lat, mHeader.sync_lon,
           mHeader.ref_lat, mHeader.ref_lon, mHeader.size, mHeader.orient);
  }
  int kol     = 108;
  int nrec    = nndx;
  size_t nlvl = size_t(mHeader.nz);
  radix_line("Number of levels (" << nlvl << ")");
  mNumVarb.clear();
  mNumVarb.resize(nlvl);
  mVarbId.clear();
  mVarbId.resize(nlvl);
  mHeight.clear();
  mHeight.resize(nlvl);
  std::vector<std::vector<int>> chk_sum(size_t(mHeader.nz));
  for (size_t l = 0; l < nlvl; ++l)
  {
    mHeight[l]  = float(std::atof(header.substr(kol, 6).c_str()));
    mNumVarb[l] = std::atoi(header.substr(kol + 6, 2).c_str());
    radix_line("Height(" << mHeight[l] << ") NumVarb(" << mNumVarb[l] << ")");

    kol += 8;

    mVarbId[l].resize(size_t(mNumVarb[l]));
    chk_sum[l].resize(size_t(mNumVarb[l]));
    for (int k = 0; k < mNumVarb[l]; ++k)
    {
      mVarbId[l][k] = header.substr(kol, 4);
      chk_sum[l][k] = std::atoi(header.substr(kol + 4, 3).c_str());
      radix_line("\tVarbId(" << mVarbId[l][k] << ") chk_sum(" << chk_sum[l][k]
                             << ")");

      kol += 8;
      nrec++;
    }
  }
  // skip to the next time period index record to find the time interval
  // between date periods (minutes)
  nrec++;

  bool first_date_loaded = false;
  mRecordTimes.clear();
  while (rstream->good())
  {
    std::string recl = rstream->readString(mLrec);
    if (recl.empty())
    {
      break;
    }
    label = recl.substr(0, 50);
    mLabel.expand(label);
    mRecordTimes.push_back(mLabel.totalSeconds());
    if (!first_date_loaded)
    {
      first_date_loaded = true;
      std::stringstream ss;
      ss << mLabel.month << "/" << mLabel.day << "/" << mLabel.year << " "
         << mLabel.hour;
      mStartTime = ss.str();
      mProfiles.push_back(ss.str());
    }
    else
    {
      // if we changed time then record profile time.
      if (mRecordTimes[mRecordTimes.size() - 2] != mLabel.totalSeconds())
      {
        std::stringstream ss;
        ss << mLabel.month << "/" << mLabel.day << "/" << mLabel.year << " "
           << mLabel.hour;
        mProfiles.push_back(ss.str());
      }
    }
    // peek ahead to check for eof etc...
    rstream->peek();
    if (!rstream->good())
    {
      std::stringstream ss;
      // if we are at the end of the file dump the ending time
      ss << mLabel.month << "/" << mLabel.day << "/" << mLabel.year << " "
         << mLabel.hour;
      mEndTime = ss.str();
    }
    //        radix_tagged_line("Profile time: " <<
    //        mProfiles[mProfiles.size()-1]
    //                << mLabel.toString());
  }
  rstream->close();
  delete rstream;
}
std::pair<float, float> GFSFile::gbl2xy(float clat, float clon, float sync_lat,
                                        float ref_lat, float sync_lon,
                                        float ref_lon) const
{
  radix_tagged_line("gbl2xy(" << clat << "," << clon << "," << sync_lat << ","
                              << ref_lat << "," << sync_lon << "," << ref_lon
                              << ")");
  float tlat = clat;
  std::pair<float, float> result;
  if (tlat > 90.0f) tlat = 180.0f - tlat;
  if (tlat < -90.0f) tlat = -180.0f - tlat;
  result.second = 1.0f + (tlat - sync_lat) / ref_lat;
  radix_tagged_line("\tcomputed y =" << result.second);

  float tlon = clon;
  if (!mHeader.prime)
  {
    if (tlon < 0.0f) tlon = 360.0f + tlon;
    if (tlon > 360.0) tlon = tlon - 360.0f;
  }
  tlon = tlon - sync_lon;
  if (tlon < 0.0f) tlon = tlon + 360.0f;
  result.first = 1.0f + tlon / ref_lon;
  radix_tagged_line("\tcomputed x =" << result.first);
  return result;
}

std::pair<float, float> GFSFile::gbl2ll(int x, int y, float sync_lat,
                                        float ref_lat, float sync_lon,
                                        float ref_lon) const
{
  std::pair<float, float> result;
  radix_tagged_line("gbl2ll(" << x << "," << y << "," << sync_lat << ","
                              << ref_lat << "," << sync_lon << "," << ref_lon
                              << ")");
  if (!mHeader.latlon) return result;
  float clat = sync_lat + (y - 1.0f) * ref_lat;
  if (clat > 90.0f) clat = 180.0f - clat;
  if (clat < -90.0f) clat = -180.0f - clat;

  float clon = sync_lon + (x - 1.0f) * ref_lon;
  clon       = std::fmod(clon, 360.f);
  if (clon > 180.0f) clon = clon - 360.f;
  result.first  = clat;
  result.second = clon;
  return result;
}

std::pair<float, float> GFSFile::cnxyll(float xi, float eta) const
{
  float gamma = mStrcmp[0];
  float cgeta = 1.f - gamma * eta;
  float gxi   = gamma * xi;
  // calculate equivalent mercator coordinate
  float arg2 = eta + (eta * cgeta - gxi * xi);
  float arg1 = gamma * arg2;
  float xlat = 0, xlong = 0, temp = 0, ymerc = 0, along = 0;
  // distance to north (or south) pole is zero (or imaginary)
  if (arg1 >= 1.f)
  {
    xlat  = std::copysign(90., mStrcmp[0]);
    xlong = 90. + xlat;
    return std::make_pair(xlat, xlong);
  }
  if (std::abs(arg1) < 0.01f)
  {
    // this avoids round-off error or divide-by zero in case of mercator
    // projects
    temp  = std::pow(arg1 / (2.f - arg1), 2);
    ymerc = arg2 / (2.f - arg1) *
            (1.f + temp * (1.f / 3.f + temp * (1.f / 5.f + temp * 1.f / 7.f)));
  }
  else
  {
    // code for moderate values of gamma
    ymerc = -std::log(1.f - arg1) / 2.f / gamma;
  }
  // convert ymerc to latitude
  temp = std::exp(-std::abs(ymerc));
  xlat =
      std::copysign(std::atan2((1.f - temp) * (1.f + temp), 2.f * temp), ymerc);
  // find longitudes
  if (std::abs(gxi) < 0.01f * cgeta)
  {
    // this avoids round-off error or divide-by zero in case of mercator
    // projects
    temp  = std::pow(gxi / cgeta, 2);
    along = xi / cgeta *
            (1.f - temp * (1.f / 3.f - temp * (1.f / 5.f - temp * 1.f / 7.f)));
  }
  else
  {
    along = std::atan2(gxi, cgeta) / gamma;
  }
  xlong = mStrcmp[1] + PI_BELOW_180 * along;
  xlat  = xlat * PI_BELOW_180;
  return std::make_pair(xlat, xlong);
}
std::pair<float, float> GFSFile::cnllxy(float clat, float clon) const
{
  std::pair<float, float> result;
  float almost1 = .99999;
  float gamma   = mStrcmp[0];
  float dlat    = clat;
  float dlong   = cspanf(clon - mStrcmp[1], -180, 180);
  dlong         = dlong * PI_ON_180;
  float gdlong  = gamma * dlong;
  float csdgam = 0.0, sndgam = 0.0;
  if (std::abs(gdlong) < 0.01)
  {
    // For gamma small or zero. avoids round-off error or division
    // by zero in the case of mercator or near-mercator projections.
    gdlong = gdlong * gdlong;
    sndgam =
        dlong *
        (1.f - 1.f / 6.f * gdlong *
                   (1.f - 1.f / 20.f * gdlong * (1.f - 1.f / 42.f * gdlong)));
    csdgam =
        dlong * dlong * .5f *
        (1.f - 1.f / 12.f * gdlong *
                   (1.f - 1.f / 30.f * gdlong * (1.f - 1.f / 56.f * gdlong)));
  }
  else
  {
    sndgam = std::sin(gdlong) / gamma;
    csdgam = (1.f - std::cos(gdlong)) / gamma / gamma;
  }
  float slat = std::sin(dlat * PI_ON_180);
  if ((slat >= almost1) || (slat <= -almost1))
  {
    result.first  = 0.0f;
    result.second = 1.f / gamma;
    return result;
  }
  float mercy  = .5f * std::log((1.f + slat) / (1.f - slat));
  float gmercy = gamma * mercy;
  float rhog1  = 0.f;
  if (std::abs(gmercy) < .001f)
  {
    // For gamma small or zero. avoids round-off error or division
    // by zero in the case of mercator or near-mercator projections.
    rhog1 = mercy *
            (1.f - .5f * gmercy *
                       (1.f - 1.f / 3.f * gmercy * (1.f - 1.f / 4.f * gmercy)));
  }
  else
  {
    rhog1 = (1.f - std::exp(-gmercy)) / gamma;
  }
  result.first  = (1.f - gamma * rhog1) * sndgam;
  result.second = rhog1 + (1.f - gamma * rhog1) * gamma * csdgam;

  return result;
}

std::pair<float, float> GFSFile::cll2xy(float clat, float clon) const
{
  radix_tagged_line("cll2xy(" << clat << "," << clon << ")");
  std::pair<float, float> xi_eta = cnllxy(clat, clon);
  radix_tagged_line("\txi=" << xi_eta.first << " eta=" << xi_eta.second);

  float radius = EARTH_RADIUS_MEAN / 1000.f;
  float x =
      mStrcmp[2] + radius / mStrcmp[6] *
                       (xi_eta.first * mStrcmp[4] + xi_eta.second * mStrcmp[5]);
  float y =
      mStrcmp[3] + radius / mStrcmp[6] *
                       (xi_eta.second * mStrcmp[4] - xi_eta.first * mStrcmp[5]);
  radix_tagged_line("\tx=" << x << " y=" << y);
  return std::make_pair(x, y);
}

std::pair<float, float> GFSFile::cxy2ll(float x, float y) const
{
  radix_tagged_line("cxy2ll(" << x << "," << y << ")");
  float radius               = EARTH_RADIUS_MEAN / 1000.f;
  float xi0                  = (x - mStrcmp[2]) * mStrcmp[6] / radius;
  float eta0                 = (y - mStrcmp[3]) * mStrcmp[6] / radius;
  float xi                   = xi0 * mStrcmp[4] - eta0 * mStrcmp[5];
  float eta                  = eta0 * mStrcmp[4] + xi0 * mStrcmp[5];
  std::pair<float, float> ll = cnxyll(xi, eta);
  radix_line("\tcnxy2ll result lat=" << ll.first << " lon=" << ll.second);
  float xlong = cspanf(ll.second, -180.f, 180.f);
  return std::make_pair(ll.first, xlong);
}

std::pair<float, float> GFSFile::cg2cll(float xlat, float xlong, float ug,
                                        float vg) const
{
  float along = cspanf(xlong - mStrcmp[1], -180.f, 180.f);
  float rot   = -mStrcmp[0] + along;
  // allow cartographic wind vector transformations everywhere
  // with rotation to nominal longitudes at the poles, to match u,v values
  // on a lat-lon grid
  float slong = std::sin(PI_ON_180 * rot);
  float clong = std::cos(PI_ON_180 * rot);
  float xpolg = slong * mStrcmp[4] + clong * mStrcmp[5];
  float ypolg = clong * mStrcmp[4] - slong * mStrcmp[5];
  float vn    = ypolg * ug + xpolg * vg;
  float ue    = ypolg * vg + xpolg * ug;
  return std::make_pair(ue, vn);
}

std::pair<float, float> GFSFile::cg2cxy(float x, float y, float ug, float vg)
{
  float xpolg = 0.f, ypolg = 0.f, temp = 0.f, xi0 = 0.f, eta0 = 0.f;
  float radius = EARTH_RADIUS_MEAN / 1000.f;
  xi0          = (x - mStrcmp[2]) * mStrcmp[6] / radius;
  eta0         = (y - mStrcmp[3]) * mStrcmp[6] / radius;
  xpolg        = mStrcmp[5] - mStrcmp[0] * xi0;
  ypolg        = mStrcmp[4] - mStrcmp[0] * eta0;
  temp         = std::sqrt(std::pow(xpolg, 2.f) + std::pow(ypolg, 2.f));
  std::pair<float, float> xy;
  if (temp <= 1e-3)
  {
    std::pair<float, float> ll = cxy2ll(x, y);
    xy                         = cg2cll(ll.first, ll.second, ug, vg);
  }
  else
  {
    // use vector alegbra instread of time consuming trig
    xpolg     = xpolg / temp;
    ypolg     = ypolg / temp;
    xy.first  = ypolg * ug - xpolg * vg;
    xy.second = ypolg * vg + xpolg * ug;
  }
  return xy;
}

void GFSFile::stlmbr(float tnglat, float xlong)
{
  float radius                   = EARTH_RADIUS_MEAN / 1000.f;
  mStrcmp[0]                     = std::sin(PI_ON_180 * tnglat);
  mStrcmp[1]                     = cspanf(xlong, -180., 180.);
  mStrcmp[2]                     = 0.;
  mStrcmp[3]                     = 0.;
  mStrcmp[4]                     = 1.;
  mStrcmp[5]                     = 0.;
  mStrcmp[6]                     = radius;
  std::pair<float, float> xi_eta = cnllxy(89., xlong);
  mStrcmp[7] = 2. * xi_eta.second - mStrcmp[0] * xi_eta.second * xi_eta.second;
  xi_eta     = cnllxy(-89., xlong);
  mStrcmp[8] = 2. * xi_eta.second - mStrcmp[0] * xi_eta.second * xi_eta.second;
}

void GFSFile::stcm1p(float x1, float y1, float xlat1, float xlong1, float xlatg,
                     float xlongg, float gridsz, float orient)
{
  radix_tagged_line("stcm1p(" << x1 << "," << y1 << "," << xlat1 << ","
                              << xlong1 << "," << xlatg << "," << xlongg << ","
                              << gridsz << "," << orient << ")");
  for (size_t i = 2; i < 4; ++i)
  {
    mStrcmp[i] = 0.f;
  }
  float turn = PI_ON_180 *
               (orient - mStrcmp[0] * cspanf(xlongg - mStrcmp[1], -180., 180.));
  radix_line("turn=" << turn);
  mStrcmp[4]         = std::cos(turn);
  mStrcmp[5]         = -std::sin(turn);
  mStrcmp[6]         = 1.f;
  float cgszllResult = cgszll(xlatg, mStrcmp[1]);
  radix_line("cgszll=" << cgszllResult);
  mStrcmp[6] = gridsz * mStrcmp[6] / cgszllResult;
  radix_line("mStrcmp[7]=" << mStrcmp[6]);
  std::pair<float, float> a1 = cll2xy(xlat1, xlong1);
  radix_line("x1a=" << a1.first << " y1a=" << a1.second);
  mStrcmp[2] = mStrcmp[2] + x1 - a1.first;
  mStrcmp[3] = mStrcmp[3] + y1 - a1.second;
  radix_line("1=" << mStrcmp[0] << ", 2=" << mStrcmp[1] << ", 3=" << mStrcmp[2]
                  << ", 4=" << mStrcmp[3] << ", 5=" << mStrcmp[4]
                  << ", 6=" << mStrcmp[5] << ", 7=" << mStrcmp[6]
                  << ", 8=" << mStrcmp[7]);
}

float GFSFile::cgszll(float xlat, float xlong) const
{
  radix_tagged_line("cgszll(" << xlat << "," << xlong);
  float slat = 0.f, ymerc = 0.f, efact = 0.f;
  if (xlat > 89.995f)
  {
    // close to north pole
    if (mStrcmp[0] > 0.9999f)
    {
      // and to gamma == 1
      return 2.f * mStrcmp[6];
    }
    efact = std::cos(PI_ON_180 * xlat);
    if (efact <= 0.f)
    {
      return 0.f;
    }
    else
    {
      ymerc = -std::log(efact / (1.f + std::sin(PI_ON_180 * xlat)));
    }
  }
  else if (xlat < -89.995f)
  {
    // close to south pole
    if (mStrcmp[0] < -0.9999f)
    {
      // and to gamma == -1.0
      return 2.f * mStrcmp[6];
    }
    efact = std::cos(PI_ON_180 * xlat);
    if (efact <= 0.f)
    {
      return 0.f;
    }
    else
    {
      ymerc = std::log(efact / (1.f - std::sin(PI_ON_180 * xlat)));
    }
  }
  else
  {
    slat  = std::sin(PI_ON_180 * xlat);
    ymerc = std::log((1.f + slat) / (1.f - slat)) / 2.f;
  }
  return mStrcmp[6] * std::cos(PI_ON_180 * xlat) * std::exp(mStrcmp[0] * ymerc);
}
std::pair<int, int> GFSFile::nearestPoint(float lat, float lon) const
{
  std::pair<float, float> point;
  radix_tagged_line("nearstPoint(" << lat << "," << lon << ")");
  radix_tagged_line("\tlatlon=" << std::boolalpha << mHeader.latlon);
  if (mHeader.latlon)
  {
    point = gbl2xy(lat, lon, mHeader.sync_lat, mHeader.ref_lat,
                   mHeader.sync_lon, mHeader.ref_lon);
  }
  else
  {
    point = cll2xy(lat, lon);
  }
  std::pair<int, int> ipoint;
  ipoint.first  = (int)std::round(point.first);
  ipoint.second = (int)std::round(point.second);
  return ipoint;
}

std::pair<float, float> GFSFile::latlon(int y, int x) const
{
  std::pair<float, float> result;
  if (mHeader.latlon)
  {
    result = gbl2ll(x, y, mHeader.sync_lat, mHeader.ref_lat, mHeader.sync_lon,
                    mHeader.ref_lon);
  }
  else
  {
    result = cxy2ll(x, y);
  }
  return result;
}
std::string GFSFile::startTime() const { return mStartTime; }
std::string GFSFile::endTime() const { return mEndTime; }
std::string GFSFile::profileTime() const { return mProfileTime; }
const std::vector<std::string> &GFSFile::profileTimes() const
{
  return mProfiles;
}
std::vector<std::vector<float>> GFSFile::query(float lat, float lon, int month,
                                               int day, int year, int hour,
                                               std::vector<std::string> columns)
{
  float searchTime = total_seconds(year, month, day, hour);
  // assume class was correctly initialized
  // get the grid points for the lon, lat in the met file
  std::pair<int, int> point = nearestPoint(lat, lon);  // gbl2xy(lat, lon
  //, mHeader.sync_lat, mHeader.ref_lat
  //, mHeader.sync_lon, mHeader.ref_lon);

  int x = point.first;   //(int)std::round(point.first);
  int y = point.second;  //(int)std::round(point.second);
  if (x < 0 || x >= mHeader.nx || y < 0 || y >= mHeader.ny)
  {
    std::cerr << "Selected location is outside of file boundary." << std::endl;
    return std::vector<std::vector<float>>();
  }
  double minDelta = 99999999.f;
  size_t minIndex = 999999999;
  for (size_t i = 0; i < mRecordTimes.size(); ++i)
  {
    double delta = mRecordTimes.at(i) - searchTime;
    if (delta > 0)
    {
      if (minDelta > delta)
      {
        minDelta = delta;
        minIndex = i;
      }
    }
    else
    {
      if (std::abs(minDelta) > std::abs(delta))
      {
        minDelta = delta;
        minIndex = i;
      }
    }
  }
  std::vector<size_t> matchingIndex;
  for (size_t i = minIndex; i < mRecordTimes.size(); ++i)
  {
    // if time has changed then lets break out of loop
    if (mRecordTimes.at(i) != mRecordTimes.at(minIndex)) break;
    matchingIndex.push_back(i);
  }
  float sfcp = 1013.0f;
  float msle = 1013.15f;
  float sfct = 0.0f;
  int lp     = 0;
  std::vector<std::vector<float>> vdata(mvar);
  for (size_t i = 0; i < vdata.size(); ++i)
    vdata[i] = std::vector<float>(mlvl, 0.0f);
  std::vector<float> utw(mlvl, 0.0f);
  std::vector<float> vtw(mlvl, 0.0f);

  std::string label, header;
  std::vector<std::vector<float>> rdata;
  // open the file for reading
  radix::eafstream *rstream = new radix::eafstream(
      mFile.c_str(), std::ifstream::in | std::ifstream::binary);
  for (size_t i = 0; i < matchingIndex.size(); ++i)
  {
    // get the fortran record index
    size_t irec = matchingIndex.at(i) + 1;
    // calculate the file offset
    size_t foffset = mLrec * irec;
    // seek to the position in the file
    rstream->seekg(foffset, rstream->beg);
    std::string recl = rstream->readString(mLrec);
    label            = recl.substr(0, 256);
    mLabel.expand(label);
    if (i == 0)
    {
      std::stringstream ss;
      ss << mLabel.month << "/" << mLabel.day << "/" << mLabel.year << " "
         << mLabel.hour;
      mProfileTime = ss.str();
    }
    header           = recl.substr(50);
    std::string varb = mLabel.kvar;
    if (varb.compare("INDX") == 0) continue;
    rdata = pakinp(header, mHeader.nx, mHeader.ny, 0, 0, mHeader.nx, mHeader.ny,
                   mLabel.prec, mLabel.nexp, mLabel.var1);

    int ll = mLabel.il;
    // convert level number to array index because input data
    // level index starts at 0 for the surface
    if (ll != lp || irec == (matchingIndex.size() - 1))
    {
      if (lp != 0 && !mHeader.latlon)
      {
        std::pair<float, float> xy = cg2cxy(x - 1, y - 1, utw[lp], vtw[lp]);
        utw[lp]                    = xy.first;
        vtw[lp]                    = xy.second;
      }
      lp = ll;
    }

    // find the variable array element number - match the input
    // variable with its position as indicated in the index record
    int nvar = mNumVarb.at(ll);
    int kvar = 0;
    for (int kk = 0; kk < nvar; ++kk)
    {
      if (varb.compare(mVarbId[ll][kk]) == 0) kvar = kk;
    }
    vdata[kvar][ll] = rdata[x - 1][y - 1];
    // convert unit of temperature to oC
    if (varb.compare("TEMP") == 0 || varb.compare("T02M") == 0 ||
        varb.compare("TMPS") == 0 || varb.compare("DP2M") == 0)
      vdata[kvar][ll] = vdata[kvar][ll] - 273.16f;

    // save the surface pressure and terrain mHeight for scaling
    // of the vertical coordinate system (mHeight = signma*scaling)
    if (varb.compare("PRSS") == 0) sfcp = vdata[kvar][ll];
    if (varb.compare("SHGT") == 0) sfct = vdata[kvar][ll];
    if (varb.compare("MSLE") == 0) msle = vdata[kvar][ll];

    // load the winds for subsequent rotation to true
    if (varb.compare("UWND") == 0 || varb.compare("U10M") == 0)
      utw[ll] = vdata[kvar][ll];
    if (varb.compare("VWND") == 0 || varb.compare("V10M") == 0)
      vtw[ll] = vdata[kvar][ll];
  }  // for matching records
  // close the file
  rstream->close();
  delete rstream;

  // SOUND section of Fortran
  float tpot            = 0.0f;
  float temp            = 0.0f;
  bool sfcwnd           = false;
  float offset          = 0.0f;
  float plevel          = 0.0f;
  float surfaceAltitude = 0.f;

  std::vector<std::vector<float>> results(mHeader.nz);
  for (int ll = 0; ll < mHeader.nz; ++ll)
  {
    int nvar = mNumVarb[ll];
    // default vertical motion units in mb/s
    for (size_t nn = 0; nn < mUnits.size(); ++nn)
      if (0 == mVarb[nn].compare("WWND")) mUnits[nn] = "mb/h";

    if (mHeader.z_flag == 1)
    {
      // pressure sigma levels
      offset = mHeader.dummy;
      plevel = offset + (sfcp - offset) * mHeight[ll];
    }
    else if (mHeader.z_flag == 2)
    {
      plevel = mHeight[ll];
      if (ll == 0) plevel = sfcp;
    }
    else if (mHeader.z_flag == 3)
    {
      float ztop = 20000.0f;
      if (mHeight[mHeader.nz - 1] > ztop) ztop = 34800.0f;
      float factor = 1.0f - sfct / ztop;
      plevel       = factor * mHeight[ll];
      // terrain follow Z system units in m/s
      for (size_t nn = 0; nn < mUnits.size(); ++nn)
      {
        if (0 == mVarb[ll].compare("WWND")) mUnits[nn] = " m/h";
      }
    }
    else if (mHeader.z_flag == 4)
    {
      // ecmwf hubrid coordinate system
      offset       = static_cast<int>(mHeight[ll]);
      float psigma = mHeight[ll] - offset;
      plevel       = sfcp * psigma + offset;
      if (ll == 0) plevel = sfcp;
    }
    // by default assume level = pressure unless PRES variable appears
    // (i.e. terrain data (type=3) will have local pressure variable

    // match variables defined in file's index record with those variables
    // that have been defined in this subroutine and create a variable number
    // for simple table lookup
    std::vector<int> nt(nvar, 0);
    for (int kk = 0; kk < nvar; ++kk)
    {
      for (size_t nn = 0; nn < mUnits.size(); ++nn)
      {
        if (mVarbId[ll][kk].compare(mVarb[nn]) == 0) nt[kk] = (int)(nn);
        // check for 10 meter winds
        if ((ll == 0) && (mVarbId[ll][kk].compare("U10M") == 0) &&
            (mVarb[nn].compare("U10M") == 0))
        {
          sfcwnd = true;
        }
        else if ((ll == 0) && (mVarbId[ll][kk].compare("V10M") == 0) &&
                 (mVarb[nn].compare("V10M") == 0))
        {
          sfcwnd = true;
        }
      }
    }
    //
    // convert each variable at that level to standard units as defined
    // from the table lookup. Variales not found are not converted and
    // have no specific units label
    for (int kk = 0; kk < nvar; ++kk)
    {
      vdata[kk][ll] = vdata[kk][ll] * mFact[nt[kk]];
    }
    // initialize space for results vector
    results[ll] = std::vector<float>(columns.size(), 0.0f);

    // if "HGTS" has been requested add level as the default HGTS
    auto hIt = std::find(columns.begin(), columns.end(), "HGTS");
    if (hIt != columns.end())
    {
      results[ll][hIt - columns.begin()] =
          hpaToAltitude(plevel, temp, msle) - surfaceAltitude;
    }
    // check for time
    auto timeIt = std::find(columns.begin(), columns.end(), "TIME");
    if (timeIt != columns.end())
    {
      if (minIndex >= mRecordTimes.size())
      {
        minIndex = *mRecordTimes.end();
      }
      int hour = (mRecordTimes[minIndex] - searchTime) / 3600.0;
      results[ll][timeIt - columns.begin()] = (float)hour;
    }
    // check for pressure
    auto presIt = std::find(columns.begin(), columns.end(), "PRSS");
    if (presIt != columns.end())
    {
      results[ll][presIt - columns.begin()] = plevel;
    }

    for (int kk = 0; kk < nvar; ++kk)
    {
      std::string varb = mVarbId[ll][kk];
      if (varb.compare("PRES") == 0) plevel = vdata[kk][ll];
      if (varb.compare("THET") == 0)
      {
        tpot = vdata[kk][ll];
        // potential temperature defined; replace with ambient
        vdata[kk][ll] = (tpot * std::pow(plevel / 1000.0f, 0.286f)) - 273.16f;
      }
      if (varb.compare("TEMP") == 0) temp = vdata[kk][ll] + 273.16f;

      // map certain surface temperature(T02M) or surface relative
      // humidity(RH2M)
      if (varb.compare("T02M") == 0) varb = "TEMP";
      if (varb.compare("RH2M") == 0) varb = "RELH";
      if (varb.compare("PRSS") == 0)
      {
        surfaceAltitude = hpaToAltitude(vdata[kk][ll], temp, msle) - 2.f;
      }
      auto it = std::find(columns.begin(), columns.end(), varb);
      if (it != columns.end())
      {
        results[ll][it - columns.begin()] = vdata[kk][ll];
      }
    }
    // check surface data (2 meters)
    if (ll == 0)
    {
      // check for surface height
      {
        auto it = std::find(columns.begin(), columns.end(), "HGTS");
        if (it != columns.end())
        {
          results[ll][it - columns.begin()] = 2.0f;
        }
      }
    }
    bool hwind = false;  // have wind?
    float wd   = 0.f;
    float ws   = 0.f;
    if (ll > 1)
    {
      // potential temperature not defined, then compute
      if (tpot == 0.0f) tpot = temp * std::pow(1000.0f / plevel, 0.286f);
      if (kwnd)
      {
        if (utw[ll] != 0.0f || vtw[ll] != 0.0f)
        {
          wd    = 57.29578f * std::atan2(utw[ll], vtw[ll]) + 360.0f;
          wd    = std::fmod(wd, 360.0f);
          wd    = std::fmod((wd + 180.0f), 360.0f);
          ws    = std::sqrt(utw[ll] * utw[ll] + vtw[ll] * vtw[ll]);
          hwind = true;
        }
      }
      else
      {
        wd    = utw[ll];
        ws    = vtw[ll];
        hwind = true;
      }
    }
    else
    {
      if (kwnd && sfcwnd)
      {
        if (utw[ll] != 0.0f || vtw[ll] != 0.0f)
        {
          wd    = 57.295778f * std::atan2(utw[ll], vtw[ll]) + 360.0f;
          wd    = std::fmod(wd, 360.0f);
          wd    = std::fmod((wd + 180.0f), 360.0f);
          ws    = std::sqrt(utw[ll] * utw[ll] + vtw[ll] * vtw[ll]);
          hwind = true;
        }
      }
    }
    if (hwind)
    {
      // check for WD
      {
        auto it = std::find(columns.begin(), columns.end(), "WD");
        if (it != columns.end())
        {
          results[ll][it - columns.begin()] = wd;
        }
      }  // check for WS
      {
        auto it = std::find(columns.begin(), columns.end(), "WS");
        if (it != columns.end())
        {
          results[ll][it - columns.begin()] = ws;
        }
      }
    }
  }  // for ll < nz

  return results;
}
}  // namespace radix