gfsfile.cc 30.6 KB
Newer Older
1
2
#include "radixio/gfsfile.hh"

3
#include "radixbug/bug.hh"
4
#include "radixio/eafstream.hh"
5
#include "radixmath/util.hh"
6

7
8
9
10
namespace radix
{
float total_seconds(int year, int month, int day, int hour)
{
11
12
13
14
15
16
17
18
19
  // 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,
20
21
                                       float prec, int nexp, float var1)
{
22
23
24
25
26
27
28
29
  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
30
31
  for (int j = 0; j < ny; ++j)
  {
32
33
34
35
    k    = j * nx;  // position at column 1
    jj   = j - ny1;
    rnew = (float(ord(cvar[k]) - 127) * scexp) + rold;
    rold = rnew;
36
37
    if (jj >= 0 && jj <= ly)
    {
38
39
40
      rvar[0][jj] = rnew;
    }
  }  // 1st for j < ny
41
42
  for (int j = ny1; j < (ny1 + ly); ++j)
  {
43
44
    jj   = j - ny1;  // sub-grid array (1 to ly)
    rold = rvar[0][jj];
45
46
    for (int i = 1; i < (nx1 + lx); ++i)
    {
47
48
49
50
51
      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;
52
53
      if (ii >= 0 && ii <= lx)
      {
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
        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};
79
80
GFSFile::GFSFile(std::string file)
    : mFile(file)
81
82
    , mStrcmp(9, 0.0f)
{
83
84
85
86
87
88
89
90
  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);
LEFEBVREJP email's avatar
LEFEBVREJP email committed
91
  mHeader.expand(header, mLabel);
92
93
94
95
96
97
98
99
100
101
102
103
104
105
  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
106
107
  for (int i = 0; i < nndx; ++i)
  {
108
109
110
    std::string recl = rstream->readString(mLrec);
    label            = recl.substr(0, 50);
    header           = recl.substr(50);
111
    mLabel.expand(label);
LEFEBVREJP email's avatar
LEFEBVREJP email committed
112
    mHeader.expand(header, mLabel);
113
114
115
116
    radix_tagged_line("Found grid: " << mHeader.toString());
  }
  //
  // determine if this is a lat-lon grid
117
118
  if (mHeader.size == 0.f)
  {
119
120
    mHeader.latlon = true;
  }
121
122
  if (mHeader.model_id.compare("RAMS") == 0)
  {
123
124
    mHeader.tang_lat = mHeader.pole_lat;
  }
125
126
  if (!mHeader.latlon)
  {
127
128
129
    //
    // initialize grid conversion variable (into gbase)
    stlmbr(mHeader.tang_lat, mHeader.ref_lon);
130
    //
131
132
133
134
135
136
137
138
139
    // 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);
  }
LEFEBVREJP email's avatar
LEFEBVREJP email committed
140
141
142
143
  int kol     = 108;
  int nrec    = nndx;
  size_t nlvl = size_t(mHeader.nz);
  radix_line("Number of levels (" << nlvl << ")");
144
145
146
147
148
149
  mNumVarb.clear();
  mNumVarb.resize(nlvl);
  mVarbId.clear();
  mVarbId.resize(nlvl);
  mHeight.clear();
  mHeight.resize(nlvl);
LEFEBVREJP email's avatar
LEFEBVREJP email committed
150
151
  std::vector<std::vector<int>> chk_sum(size_t(mHeader.nz));
  for (size_t l = 0; l < nlvl; ++l)
152
  {
LEFEBVREJP email's avatar
LEFEBVREJP email committed
153
    mHeight[l]  = float(std::atof(header.substr(kol, 6).c_str()));
154
    mNumVarb[l] = std::atoi(header.substr(kol + 6, 2).c_str());
LEFEBVREJP email's avatar
LEFEBVREJP email committed
155
    radix_line("Height(" << mHeight[l] << ") NumVarb(" << mNumVarb[l] << ")");
156
157
158

    kol += 8;

LEFEBVREJP email's avatar
LEFEBVREJP email committed
159
160
    mVarbId[l].resize(size_t(mNumVarb[l]));
    chk_sum[l].resize(size_t(mNumVarb[l]));
161
162
    for (int k = 0; k < mNumVarb[l]; ++k)
    {
163
164
      mVarbId[l][k] = header.substr(kol, 4);
      chk_sum[l][k] = std::atoi(header.substr(kol + 4, 3).c_str());
LEFEBVREJP email's avatar
LEFEBVREJP email committed
165
166
      radix_line("\tVarbId(" << mVarbId[l][k] << ") chk_sum(" << chk_sum[l][k]
                             << ")");
167
168
169

      kol += 8;
      nrec++;
170
    }
171
172
173
174
175
176
177
  }
  // 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();
178
179
  while (rstream->good())
  {
180
    std::string recl = rstream->readString(mLrec);
181
182
    if (recl.empty())
    {
183
      break;
184
    }
185
186
187
    label = recl.substr(0, 50);
    mLabel.expand(label);
    mRecordTimes.push_back(mLabel.totalSeconds());
188
189
    if (!first_date_loaded)
    {
190
191
192
193
194
195
      first_date_loaded = true;
      std::stringstream ss;
      ss << mLabel.month << "/" << mLabel.day << "/" << mLabel.year << " "
         << mLabel.hour;
      mStartTime = ss.str();
      mProfiles.push_back(ss.str());
196
197
198
    }
    else
    {
199
      // if we changed time then record profile time.
200
201
      if (mRecordTimes[mRecordTimes.size() - 2] != mLabel.totalSeconds())
      {
202
203
204
205
206
        std::stringstream ss;
        ss << mLabel.month << "/" << mLabel.day << "/" << mLabel.year << " "
           << mLabel.hour;
        mProfiles.push_back(ss.str());
      }
207
    }
208
209
    // peek ahead to check for eof etc...
    rstream->peek();
210
211
    if (!rstream->good())
    {
212
213
214
215
216
      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();
217
    }
218
219
220
221
222
223
    //        radix_tagged_line("Profile time: " <<
    //        mProfiles[mProfiles.size()-1]
    //                << mLabel.toString());
  }
  rstream->close();
  delete rstream;
224
}
225
226
std::pair<float, float> GFSFile::gbl2xy(float clat, float clon, float sync_lat,
                                        float ref_lat, float sync_lon,
227
228
                                        float ref_lon) const
{
229
230
231
232
233
234
235
236
237
238
239
  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;
240
241
  if (!mHeader.prime)
  {
242
243
244
245
246
247
248
249
    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;
250
}
251

252
253
std::pair<float, float> GFSFile::gbl2ll(int x, int y, float sync_lat,
                                        float ref_lat, float sync_lon,
254
255
                                        float ref_lon) const
{
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
  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;
271
272
}

273
274
std::pair<float, float> GFSFile::cnxyll(float xi, float eta) const
{
275
276
277
278
279
280
281
282
  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)
283
284
  if (arg1 >= 1.f)
  {
285
286
    xlat  = std::copysign(90., mStrcmp[0]);
    xlong = 90. + xlat;
287
    return std::make_pair(xlat, xlong);
288
  }
289
290
  if (std::abs(arg1) < 0.01f)
  {
291
292
293
294
295
    // 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)));
296
297
298
  }
  else
  {
299
300
301
302
303
304
305
306
    // 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
307
308
  if (std::abs(gxi) < 0.01f * cgeta)
  {
309
310
311
312
313
    // 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)));
314
315
316
  }
  else
  {
317
318
319
320
321
    along = std::atan2(gxi, cgeta) / gamma;
  }
  xlong = mStrcmp[1] + PI_BELOW_180 * along;
  xlat  = xlat * PI_BELOW_180;
  return std::make_pair(xlat, xlong);
322
}
323
324
std::pair<float, float> GFSFile::cnllxy(float clat, float clon) const
{
325
326
327
328
329
330
331
332
  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;
333
334
  if (std::abs(gdlong) < 0.01)
  {
335
336
337
338
339
340
341
342
343
344
345
    // 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)));
346
347
348
  }
  else
  {
349
350
351
352
    sndgam = std::sin(gdlong) / gamma;
    csdgam = (1.f - std::cos(gdlong)) / gamma / gamma;
  }
  float slat = std::sin(dlat * PI_ON_180);
353
354
  if ((slat >= almost1) || (slat <= -almost1))
  {
355
356
    result.first  = 0.0f;
    result.second = 1.f / gamma;
357
    return result;
358
359
360
361
  }
  float mercy  = .5f * std::log((1.f + slat) / (1.f - slat));
  float gmercy = gamma * mercy;
  float rhog1  = 0.f;
362
363
  if (std::abs(gmercy) < .001f)
  {
364
365
366
367
368
    // 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)));
369
370
371
  }
  else
  {
372
373
374
375
376
377
    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;
378
}
379

380
381
std::pair<float, float> GFSFile::cll2xy(float clat, float clon) const
{
382
383
384
385
386
387
388
389
390
391
392
393
394
  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);
395
396
}

397
398
std::pair<float, float> GFSFile::cxy2ll(float x, float y) const
{
399
400
401
402
403
404
405
406
407
408
  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);
409
410
}

411
std::pair<float, float> GFSFile::cg2cll(float xlat, float xlong, float ug,
412
413
                                        float vg) const
{
414
415
416
417
418
419
420
421
422
423
424
425
  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);
426
427
}

428
429
std::pair<float, float> GFSFile::cg2cxy(float x, float y, float ug, float vg)
{
430
431
432
433
434
435
436
437
  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;
438
439
  if (temp <= 1e-3)
  {
440
441
    std::pair<float, float> ll = cxy2ll(x, y);
    xy                         = cg2cll(ll.first, ll.second, ug, vg);
442
443
444
  }
  else
  {
445
446
447
448
449
450
451
    // 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;
452
453
}

454
455
void GFSFile::stlmbr(float tnglat, float xlong)
{
456
457
458
459
460
461
462
463
464
465
466
467
  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;
468
469
}

470
void GFSFile::stcm1p(float x1, float y1, float xlat1, float xlong1, float xlatg,
471
472
                     float xlongg, float gridsz, float orient)
{
473
474
475
  radix_tagged_line("stcm1p(" << x1 << "," << y1 << "," << xlat1 << ","
                              << xlong1 << "," << xlatg << "," << xlongg << ","
                              << gridsz << "," << orient << ")");
476
477
  for (size_t i = 2; i < 4; ++i)
  {
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
    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]);
498
499
}

500
501
float GFSFile::cgszll(float xlat, float xlong) const
{
502
503
  radix_tagged_line("cgszll(" << xlat << "," << xlong);
  float slat = 0.f, ymerc = 0.f, efact = 0.f;
504
505
  if (xlat > 89.995f)
  {
506
    // close to north pole
507
508
    if (mStrcmp[0] > 0.9999f)
    {
509
510
      // and to gamma == 1
      return 2.f * mStrcmp[6];
511
    }
512
    efact = std::cos(PI_ON_180 * xlat);
513
514
    if (efact <= 0.f)
    {
515
      return 0.f;
516
517
518
    }
    else
    {
519
      ymerc = -std::log(efact / (1.f + std::sin(PI_ON_180 * xlat)));
520
    }
521
522
523
  }
  else if (xlat < -89.995f)
  {
524
    // close to south pole
525
526
    if (mStrcmp[0] < -0.9999f)
    {
527
528
      // and to gamma == -1.0
      return 2.f * mStrcmp[6];
529
    }
530
    efact = std::cos(PI_ON_180 * xlat);
531
532
    if (efact <= 0.f)
    {
533
      return 0.f;
534
535
536
    }
    else
    {
537
      ymerc = std::log(efact / (1.f - std::sin(PI_ON_180 * xlat)));
538
    }
539
540
541
  }
  else
  {
542
543
544
545
    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);
546
}
547
548
std::pair<int, int> GFSFile::nearestPoint(float lat, float lon) const
{
549
550
551
  std::pair<float, float> point;
  radix_tagged_line("nearstPoint(" << lat << "," << lon << ")");
  radix_tagged_line("\tlatlon=" << std::boolalpha << mHeader.latlon);
552
553
  if (mHeader.latlon)
  {
554
555
    point = gbl2xy(lat, lon, mHeader.sync_lat, mHeader.ref_lat,
                   mHeader.sync_lon, mHeader.ref_lon);
556
557
558
  }
  else
  {
559
560
561
562
563
564
    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;
565
}
566

567
568
std::pair<float, float> GFSFile::latlon(int y, int x) const
{
569
  std::pair<float, float> result;
570
571
  if (mHeader.latlon)
  {
572
573
    result = gbl2ll(x, y, mHeader.sync_lat, mHeader.ref_lat, mHeader.sync_lon,
                    mHeader.ref_lon);
574
575
576
  }
  else
  {
577
578
579
    result = cxy2ll(x, y);
  }
  return result;
580
}
581
582
583
std::string GFSFile::startTime() const { return mStartTime; }
std::string GFSFile::endTime() const { return mEndTime; }
std::string GFSFile::profileTime() const { return mProfileTime; }
584
585
const std::vector<std::string> &GFSFile::profileTimes() const
{
586
  return mProfiles;
587
}
588
589
590
591
std::vector<std::vector<float>> GFSFile::query(float lat, float lon, int month,
                                               int day, int year, int hour,
                                               std::vector<std::string> columns)
{
592
593
594
595
596
597
598
599
600
  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);
601
602
  if (x < 0 || x >= mHeader.nx || y < 0 || y >= mHeader.ny)
  {
603
604
605
606
607
    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;
608
609
  for (size_t i = 0; i < mRecordTimes.size(); ++i)
  {
610
    double delta = mRecordTimes.at(i) - searchTime;
611
612
613
614
    if (delta > 0)
    {
      if (minDelta > delta)
      {
615
616
617
        minDelta = delta;
        minIndex = i;
      }
618
619
620
621
622
    }
    else
    {
      if (std::abs(minDelta) > std::abs(delta))
      {
623
624
625
        minDelta = delta;
        minIndex = i;
      }
626
    }
627
628
  }
  std::vector<size_t> matchingIndex;
629
630
  for (size_t i = minIndex; i < mRecordTimes.size(); ++i)
  {
631
632
633
634
635
    // 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;
LEFEBVREJP email's avatar
LEFEBVREJP email committed
636
  float msle = 1013.15f;
637
638
639
640
641
642
643
644
645
646
647
648
649
  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);
650
651
  for (size_t i = 0; i < matchingIndex.size(); ++i)
  {
652
653
654
655
656
657
658
659
660
    // 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);
661
662
    if (i == 0)
    {
663
664
665
666
      std::stringstream ss;
      ss << mLabel.month << "/" << mLabel.day << "/" << mLabel.year << " "
         << mLabel.hour;
      mProfileTime = ss.str();
667
    }
668
669
670
671
672
673
674
675
676
    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
677
678
679
680
    if (ll != lp || irec == (matchingIndex.size() - 1))
    {
      if (lp != 0 && !mHeader.latlon)
      {
681
682
683
684
685
        std::pair<float, float> xy = cg2cxy(x - 1, y - 1, utw[lp], vtw[lp]);
        utw[lp]                    = xy.first;
        vtw[lp]                    = xy.second;
      }
      lp = ll;
686
687
    }

688
689
690
691
    // 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;
692
693
    for (int kk = 0; kk < nvar; ++kk)
    {
694
695
696
697
698
699
700
701
702
703
704
705
      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];
LEFEBVREJP email's avatar
LEFEBVREJP email committed
706
    if (varb.compare("MSLE") == 0) msle = vdata[kvar][ll];
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726

    // 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);
727
728
  for (int ll = 0; ll < mHeader.nz; ++ll)
  {
729
730
731
    int nvar = mNumVarb[ll];
    // default vertical motion units in mb/s
    for (size_t nn = 0; nn < mUnits.size(); ++nn)
LEFEBVREJP email's avatar
LEFEBVREJP email committed
732
      if (0 == mVarb[nn].compare("WWND")) mUnits[nn] = "mb/h";
733

734
735
    if (mHeader.z_flag == 1)
    {
736
737
738
      // pressure sigma levels
      offset = mHeader.dummy;
      plevel = offset + (sfcp - offset) * mHeight[ll];
739
740
741
    }
    else if (mHeader.z_flag == 2)
    {
742
743
      plevel = mHeight[ll];
      if (ll == 0) plevel = sfcp;
744
745
746
    }
    else if (mHeader.z_flag == 3)
    {
747
748
749
750
751
      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
752
753
      for (size_t nn = 0; nn < mUnits.size(); ++nn)
      {
754
755
        if (0 == mVarb[ll].compare("WWND")) mUnits[nn] = " m/h";
      }
756
757
758
    }
    else if (mHeader.z_flag == 4)
    {
759
760
761
762
763
764
765
766
767
768
769
770
771
      // 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);
772
773
774
775
    for (int kk = 0; kk < nvar; ++kk)
    {
      for (size_t nn = 0; nn < mUnits.size(); ++nn)
      {
776
777
778
        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) &&
779
780
            (mVarb[nn].compare("U10M") == 0))
        {
781
          sfcwnd = true;
782
783
784
785
        }
        else if ((ll == 0) && (mVarbId[ll][kk].compare("V10M") == 0) &&
                 (mVarb[nn].compare("V10M") == 0))
        {
786
          sfcwnd = true;
787
        }
788
789
790
791
792
793
      }
    }
    //
    // 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
794
795
    for (int kk = 0; kk < nvar; ++kk)
    {
796
797
798
799
800
801
802
      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");
803
804
    if (hIt != columns.end())
    {
805
      results[ll][hIt - columns.begin()] =
LEFEBVREJP email's avatar
LEFEBVREJP email committed
806
          hpaToAltitude(plevel, temp, msle) - surfaceAltitude;
807
808
809
    }
    // check for time
    auto timeIt = std::find(columns.begin(), columns.end(), "TIME");
810
811
812
813
    if (timeIt != columns.end())
    {
      if (minIndex >= mRecordTimes.size())
      {
814
815
816
817
818
819
820
        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");
821
822
    if (presIt != columns.end())
    {
823
824
      results[ll][presIt - columns.begin()] = plevel;
    }
825

826
827
    for (int kk = 0; kk < nvar; ++kk)
    {
828
829
      std::string varb = mVarbId[ll][kk];
      if (varb.compare("PRES") == 0) plevel = vdata[kk][ll];
830
831
      if (varb.compare("THET") == 0)
      {
832
833
834
835
836
837
838
839
840
841
        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";
842
843
      if (varb.compare("PRSS") == 0)
      {
LEFEBVREJP email's avatar
LEFEBVREJP email committed
844
        surfaceAltitude = hpaToAltitude(vdata[kk][ll], temp, msle) - 2.f;
845
      }
846
847
848
849
850
851
852
853
854
855
856
857
      // check whether relative humidity has been requested, and convert
      // specific humidity if this is present instead
      auto relHumIt = std::find(columns.begin(), columns.end(), "RELH");
      if (relHumIt != columns.end() && varb.compare("SPHU") == 0)
      {
        varb                    = "RELH";
        double specificHumidity = vdata[kk][ll];
        // Convert specific humidity to relative humidity
        vdata[kk][ll] =
            specificHumidityToRelativeHumidity(specificHumidity, temp, plevel);
      }

858
      auto it = std::find(columns.begin(), columns.end(), varb);
859
860
      if (it != columns.end())
      {
861
862
863
864
        results[ll][it - columns.begin()] = vdata[kk][ll];
      }
    }
    // check surface data (2 meters)
865
866
    if (ll == 0)
    {
867
868
869
      // check for surface height
      {
        auto it = std::find(columns.begin(), columns.end(), "HGTS");
870
871
        if (it != columns.end())
        {
872
          results[ll][it - columns.begin()] = 2.0f;
873
        }
874
875
876
877
878
      }
    }
    bool hwind = false;  // have wind?
    float wd   = 0.f;
    float ws   = 0.f;
879
880
    if (ll > 1)
    {
881
882
      // potential temperature not defined, then compute
      if (tpot == 0.0f) tpot = temp * std::pow(1000.0f / plevel, 0.286f);
883
884
885
886
      if (kwnd)
      {
        if (utw[ll] != 0.0f || vtw[ll] != 0.0f)
        {
887
888
889
890
891
          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;
892
        }
893
894
895
      }
      else
      {
896
897
898
899
        wd    = utw[ll];
        ws    = vtw[ll];
        hwind = true;
      }
900
901
902
903
904
905
906
    }
    else
    {
      if (kwnd && sfcwnd)
      {
        if (utw[ll] != 0.0f || vtw[ll] != 0.0f)
        {
907
908
909
910
911
          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;
912
        }
913
914
      }
    }
915
916
    if (hwind)
    {
917
918
919
      // check for WD
      {
        auto it = std::find(columns.begin(), columns.end(), "WD");
920
921
        if (it != columns.end())
        {
922
          results[ll][it - columns.begin()] = wd;
923
        }
924
925
926
      }  // check for WS
      {
        auto it = std::find(columns.begin(), columns.end(), "WS");
927
928
        if (it != columns.end())
        {
929
          results[ll][it - columns.begin()] = ws;
930
        }
931
932
933
      }
    }
  }  // for ll < nz
934

935
  return results;
936
}
937
}  // namespace radix