arldatastream.cc 13.7 KB
Newer Older
1
2
3
4
#include "radixio/arldatastream.hh"

#include "radixbug/bug.hh"
#include "radixcore/stringfunctions.hh"
5
6
#include "radixio/eafstream.hh"
#include "radixmath/util.hh"
7
8

#include <fstream>
9

10
#ifdef _MSC_VER
11
12
13
14
15
#if _MSC_VER < 1900
#define snprintf _snprintf
#endif
#endif

16
17
18
19
20
21
namespace radix
{
class ARLDataStream::PImpl
{
 public:
  std::string file;
22
  std::shared_ptr<radix::eafstream> stream;
23
24
  int recordSize                      = 0;
  static const int recordHeaderLength = 50, indexHeaderLength = 108;
25
26
};

27
void ARLDataStream::expand(const std::string& val, ARLRecordHeader& rheader)
28
{
29
  if (val.size() < ARLDataStream::PImpl::recordHeaderLength)
30
31
  {
    throw std::runtime_error(
32
        "Incorrect size for string expansion to ARLRecordHeader.");
33
  }
34
35
36
37
38
39
40
41
42
43
44
  rheader.year  = std::atoi(val.substr(0, 2).c_str());
  rheader.month = std::atoi(val.substr(2, 2).c_str());
  rheader.day   = std::atoi(val.substr(4, 2).c_str());
  rheader.hour  = std::atoi(val.substr(6, 2).c_str());
  rheader.ic    = std::atoi(val.substr(8, 2).c_str());
  rheader.il    = std::atoi(val.substr(10, 2).c_str());
  rheader.cgrid = val.substr(12, 2).c_str();
  rheader.kvar  = val.substr(14, 4);
  rheader.nexp  = std::atoi(val.substr(18, 4).c_str());
  rheader.prec  = float(std::atof(val.substr(22, 14).c_str()));
  rheader.var1  = float(std::atof(val.substr(36, 14).c_str()));
45
}
46
47
48
void ARLDataStream::expand(const std::string& val,
                           const ARLRecordHeader& rheader,
                           ARLIndexHeader& iheader)
49
{
50
  if (val.size() < ARLDataStream::PImpl::indexHeaderLength)
51
52
53
54
  {
    throw std::runtime_error(
        "Incorrect size for string expansion to ARLHeader.");
  }
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
  iheader.model_id = val.substr(0, 4);
  iheader.icx      = std::atoi(val.substr(4, 3).c_str());
  iheader.mn       = std::atoi(val.substr(7, 2).c_str());
  iheader.pole_lat = float(std::atof(val.substr(9, 7).c_str()));
  iheader.pole_lon = float(std::atof(val.substr(16, 7).c_str()));
  iheader.ref_lat  = float(std::atof(val.substr(23, 7).c_str()));
  iheader.ref_lon  = float(std::atof(val.substr(30, 7).c_str()));
  iheader.size     = float(std::atof(val.substr(37, 7).c_str()));
  iheader.orient   = float(std::atof(val.substr(44, 7).c_str()));
  iheader.tang_lat = float(std::atof(val.substr(51, 7).c_str()));
  iheader.sync_xp  = float(std::atof(val.substr(58, 7).c_str()));
  iheader.sync_yp  = float(std::atof(val.substr(65, 7).c_str()));
  iheader.sync_lat = float(std::atof(val.substr(72, 7).c_str()));
  iheader.sync_lon = float(std::atof(val.substr(79, 7).c_str()));
  iheader.dummy    = float(std::atof(val.substr(86, 7).c_str()));
  iheader.nx       = std::atoi(val.substr(93, 3).c_str());
  iheader.ny       = std::atoi(val.substr(96, 3).c_str());
  iheader.nz       = std::atoi(val.substr(99, 3).c_str());
  iheader.z_flag   = std::atoi(val.substr(102, 2).c_str());
  iheader.lenh     = std::atoi(val.substr(104, 4).c_str());

  int knx = ordinal(rheader.cgrid[0]);
  int kny = ordinal(rheader.cgrid[1]);
78
79
80
  // Check for the grid domain extending beyond 3 digits
  if (knx >= 64 || kny >= 64)
  {
81
82
    iheader.nx = (knx - 64) * 1000 + iheader.nx;
    iheader.ny = (kny - 64) * 1000 + iheader.ny;
83
84
85
  }
}

86
87
ARLDataStream::ARLDataStream(const std::string& file,
                             std::ios_base::openmode mode)
Norby, Tom's avatar
Norby, Tom committed
88
    : p(new PImpl())
89
{
90
  p->file   = file;
91
92
  p->stream = std::make_shared<radix::eafstream>(p->file.c_str(),
                                                 mode | std::ios::binary);
93
94
  if (!p->stream->is_open())
  {
95
    throw std::runtime_error("Error opening file '" + p->file + "'");
96
97
98
  }
}

Norby, Tom's avatar
Norby, Tom committed
99
100
ARLDataStream::~ARLDataStream() { delete p; }

101
bool ARLDataStream::read_record_header(ARLRecordHeader& rheader)
102
{
103
104
  bool result = false;

105
  std::string headerString =
106
107
108
109
110
111
112
113
114
115
      p->stream->readString(ARLDataStream::PImpl::recordHeaderLength);
  expand(headerString, rheader);

  radix_line("Read record header: year = " << rheader.year);
  radix_line("                   month = " << rheader.month);
  radix_line("                     day = " << rheader.day);
  radix_line("                    hour = " << rheader.hour);
  radix_line("                variable = " << rheader.kvar);
  radix_line("        scaling exponent = " << rheader.nexp);
  radix_line("           initial value = " << rheader.var1);
116
  radix_line("  header as string:\n  " << headerString);
117
118
119

  result = true;
  return result;
120
121
}

122
bool ARLDataStream::write_record_header(const ARLRecordHeader& rheader)
123
{
124
125
  bool result = false;

126
127
128
  // Remove century from year
  int year = rheader.year % 100;

129
  radix_line("Write record header: year = " << year);
130
131
132
133
  radix_line("                   month = " << rheader.month);
  radix_line("                     day = " << rheader.day);
  radix_line("                    hour = " << rheader.hour);
  radix_line("                variable = " << rheader.kvar);
134
  radix_line("                   level = " << rheader.il);
135
136
  radix_line("        scaling exponent = " << rheader.nexp);
  radix_line("           initial value = " << rheader.var1);
137
138

  // Construct index header string
139
  char recordHeader[51];
Norby, Tom's avatar
Norby, Tom committed
140
141
142
143
144
145
  snprintf(recordHeader, sizeof(recordHeader),
           "%2d%2d%2d%2d%2d%2d%2d%4s%4d%14.7E%14.7E", year, rheader.month,
           rheader.day, rheader.hour, rheader.ic, rheader.il,
           std::stoi(rheader.cgrid.c_str()), rheader.kvar.c_str(), rheader.nexp,
           rheader.prec, rheader.var1);
  recordHeader[50] = '\0';  // null-terminate
146

147
148
  p->stream->writeString(std::string(recordHeader),
                         ARLDataStream::PImpl::recordHeaderLength);
149

150
  radix_line("  header as string:\n  " << recordHeader);
151

152
153
  result = true;
  return result;
154
}
155

156
bool ARLDataStream::read_next_record_header(ARLRecordHeader& rheader)
157
158
159
{
  int bytesToSkip = roundUpInt(p->stream->bytesRead(), p->recordSize) -
                    p->stream->bytesRead();
160
  radix_line("Skipping " << bytesToSkip << " bytes to the next record header");
161
162
163

  p->stream->skipBytes(bytesToSkip);

164
  return read_record_header(rheader);
165
166
}

167
168
bool ARLDataStream::read_index_header(const ARLRecordHeader& rheader,
                                      ARLIndexHeader& iheader)
169
{
170
171
  bool result = false;

172
  // Read main part of index header
173
  radix_line("  Reading main part of index header...");
174
  std::string headerString =
175
176
      p->stream->readString(ARLDataStream::PImpl::indexHeaderLength);
  expand(headerString, rheader, iheader);
177
178

  // Calculate and save record size
179
  p->recordSize =
180
      (iheader.nx * iheader.ny) + ARLDataStream::PImpl::recordHeaderLength;
181

182
  // Read variable description part of index header
183
184
  radix_line("  Reading variable part of index header... " << iheader.nz
                                                           << " levels");
185
  for (size_t level = 0; level < size_t(iheader.nz); ++level)
186
187
188
189
190
  {
    // Read levels information
    std::string levelString = p->stream->readString(6);
    iheader.levels.push_back(std::stof(levelString));
    std::string nVarsString = p->stream->readString(2);
191
    iheader.num_vars_at_levels.push_back(std::stoi(nVarsString));
192
193

    // Read names and checksums
194
195
196
197
198
    iheader.var_names.push_back(std::vector<std::string>());
    iheader.check_sums.push_back(std::vector<int>());
    radix_line("    Level " << level << "... "
                            << iheader.num_vars_at_levels[level] << " vars");
    for (int var = 0; var < iheader.num_vars_at_levels[level]; ++var)
199
200
    {
      std::string nameString = p->stream->readString(4);
201
      iheader.var_names[level].push_back(nameString);
202
      std::string sumString = p->stream->readString(3);
203
      iheader.check_sums[level].push_back(std::stoi(sumString));
204
205
206
207
      p->stream->skipBytes(1);
    }
  }

208
  radix_line("  Read index header: model = " << iheader.model_id);
209
210
211
  radix_line("                         nx = " << iheader.nx);
  radix_line("                         ny = " << iheader.ny);
  radix_line("                         nz = " << iheader.nz);
212
  radix_line("        Size of each record = " << p->recordSize);
213
  radix_line("    header as string:\n    " << headerString);
214
215
216

  result = true;
  return result;
217
218
}

219
220
bool ARLDataStream::write_index_header(const ARLRecordHeader& rheader,
                                       const ARLIndexHeader& iheader)
221
{
222
223
  bool result = false;

224
  // Calculate and save record size
225
  p->recordSize =
226
      (iheader.nx * iheader.ny) + ARLDataStream::PImpl::recordHeaderLength;
227

228
  radix_line("  Write index header: model = " << iheader.model_id);
229
230
231
  radix_line("                          nx = " << iheader.nx);
  radix_line("                          ny = " << iheader.ny);
  radix_line("                          nz = " << iheader.nz);
232
233
  radix_line("         Size of each record = " << p->recordSize);

Purves, Murray's avatar
Purves, Murray committed
234
  const size_t MAX_HEADER_LENGTH = 10000;
Norby, Tom's avatar
Norby, Tom committed
235
  char indexHeaderMain[ARLDataStream::PImpl::indexHeaderLength + 1],
236
      indexHeaderVars[MAX_HEADER_LENGTH];
237
238
239
  int pos = 0;

  // Write the variable description part of the header
240
  for (size_t level = 0; level < size_t(iheader.nz); ++level)
241
  {
242
    pos += sprintf(indexHeaderVars + pos, "%6.1f", iheader.levels[level]);
243
244
245
    pos += sprintf(indexHeaderVars + pos, "%2d",
                   iheader.num_vars_at_levels[level]);
    for (size_t var = 0; var < size_t(iheader.num_vars_at_levels[level]); ++var)
246
    {
247
      pos += sprintf(indexHeaderVars + pos, "%4s",
248
249
250
                     iheader.var_names[level][var].c_str());
      pos += sprintf(indexHeaderVars + pos, "%3d ",
                     iheader.check_sums[level][var]);
251
252
253
254
    }
  }

  // Construct the main part of the header (inc. size calculation)
255
  int headerLength = ARLDataStream::PImpl::indexHeaderLength + pos;
256
  sprintf(indexHeaderMain,
257
258
          "%4s%3d%2d%7.2f%7.2f%7.2f%7.2f%7.2f%7.2f%7.2f%7.2f%7.2f%7.2f%7.2f%7."
          "2f%3d%3d%3d%2d%4d",
259
260
261
262
          iheader.model_id.c_str(), iheader.icx, iheader.mn, iheader.pole_lat,
          iheader.pole_lon, iheader.ref_lat, iheader.ref_lon, iheader.size,
          iheader.orient, iheader.tang_lat, iheader.sync_xp, iheader.sync_yp,
          iheader.sync_lat, iheader.sync_lon, iheader.dummy, iheader.nx,
263
          iheader.ny, iheader.nz, iheader.z_flag, headerLength);
Norby, Tom's avatar
Norby, Tom committed
264
265
  indexHeaderMain[ARLDataStream::PImpl::indexHeaderLength] =
      '\0';  // null-terminate
266

267
  // Write the two elements of the header
268
269
  p->stream->writeString(std::string(indexHeaderMain),
                         ARLDataStream::PImpl::indexHeaderLength);
270

271
  p->stream->writeString(std::string(indexHeaderVars), pos);
272
273
274
275
276
277
278

  // Skip enough bytes to get to the next index header
  int bytesToSkip = roundUpInt(p->stream->bytesWritten(), p->recordSize) -
                    p->stream->bytesWritten();
  radix_line("    Write complete: skipping "
             << bytesToSkip << " bytes to the next index header");

279
  p->stream->writeString("", size_t(bytesToSkip));
280

281
282
  radix_line("    header as string:\n    " << indexHeaderMain
                                           << indexHeaderVars);
283

284
285
  result = true;
  return result;
286
287
}

288
289
bool ARLDataStream::read_record(const ARLRecordHeader& rheader,
                                const ARLIndexHeader& iheader,
290
                                std::vector<std::vector<float> >& record)
291
{
292
293
  bool result = false;

294
  // Set up the vector size
Purves, Murray's avatar
Purves, Murray committed
295
  record.clear();
296
  record.resize(iheader.nx);
297
  for (std::vector<float>& vec : record)
298
  {
299
    vec.resize(size_t(iheader.ny));
300
301
302
  }

  // Calculate the scaling factor
303
304
  float scaleFactor = std::pow(2.0f, 7.0f - float(rheader.nexp)),
        lastValue   = 0.f;
305
306

  radix_line("    Reading record data:");
307
  for (int y = 0; y < iheader.ny; ++y)
308
  {
309
    for (int x = 0; x < iheader.nx; ++x)
310
311
    {
      // Read the raw value
312
      unsigned char ch = (unsigned char)(p->stream->readChar());
313
314
315
      int packedValue  = ch - 127;

      // Calculate the unpacked value
316
      float unpackedValue = 0.0;
317
318
319
320
321
      // Get the correct 'last value' if we are at a 0 index
      if (x == 0)
      {
        if (y == 0)
        {
322
          lastValue = rheader.var1;
323
324
325
326
327
328
329
330
331
332
333
        }
        else
        {
          lastValue = record[x][y - 1];
        }
      }

      unpackedValue = (packedValue / scaleFactor) + lastValue;
      record[x][y]  = unpackedValue;
      lastValue     = unpackedValue;

334
      if ((x < 2 && y < 2) || ((iheader.nx - x < 3) && (iheader.ny - y < 3)))
335
336
337
338
339
340
341
      {
        radix_line("      [" << x << "," << y << "]: packed = " << packedValue
                             << ", unpacked = " << unpackedValue);
      }
    }
  }

342
343
  result = true;
  return result;
344
345
}

346
347
bool ARLDataStream::write_record(const ARLRecordHeader& rheader,
                                 const ARLIndexHeader& iheader,
348
                                 const std::vector<std::vector<float> >& record)
349
{
350
351
  bool result = false;

352
353
354
355
356
357
358
  if (record.size() == 0)
  {
    radix_line("No data here - returning false");
    return false;
  }

  // Calculate the scaling factor
359
360
  float scaleFactor = std::pow(2.0f, 7.0f - float(rheader.nexp)),
        lastValue   = 0.f;
361
362

  radix_line("    Writing record data:");
363
  for (int y = 0; y < iheader.ny; ++y)
364
  {
365
    for (int x = 0; x < iheader.nx; ++x)
366
367
368
369
370
371
    {
      // Get the correct 'last value' if we are at a 0 index
      if (x == 0)
      {
        if (y == 0)
        {
372
          lastValue = rheader.var1;
373
374
375
        }
        else
        {
376
          lastValue = record[size_t(x)][size_t(y - 1)];
377
378
379
380
        }
      }

      // Calculate the packed value
381
      float unpackedValue = record[x][y];
382
      int packedInt       = int((unpackedValue - lastValue) * scaleFactor);
383
384
385
386
387
388
389

      unsigned char packedValue = (unsigned char)(packedInt + 127);

      p->stream->writeChar(packedValue);

      lastValue = unpackedValue;

390
      if ((x < 2 && y < 2) || ((iheader.nx - x < 3) && (iheader.ny - y < 3)))
391
392
393
394
395
396
397
      {
        radix_line("      [" << x << "," << y << "]: packed = " << packedInt
                             << ", unpacked = " << unpackedValue);
      }
    }
  }

398
399
  result = true;
  return result;
400
401
}

402
403
void ARLDataStream::close_stream() { p->stream->close(); }

404
}  // namespace radix