arldatastream.cc 10.8 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
9
10
11
12
13
14

#include <fstream>
namespace radix
{
class ARLDataStream::PImpl
{
 public:
  std::string file;
15
  std::shared_ptr<radix::eafstream> stream;
16
17
  int recordSize                     = 0;
  static const int indexHeaderLength = 50, recordHeaderLength = 108;
18
19
20
21
};

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

  int knx = ordinal(index.cgrid[0]);
  int kny = ordinal(index.cgrid[1]);
  // Check for the grid domain extending beyond 3 digits
  if (knx >= 64 || kny >= 64)
  {
    header.nx = (knx - 64) * 1000 + header.nx;
    header.ny = (kny - 64) * 1000 + header.ny;
  }
}

78
79
ARLDataStream::ARLDataStream(const std::string& file,
                             std::ios_base::openmode mode)
80
81
    : p(new PImpl(), [](PImpl* impl) { delete impl; })
{
82
  p->file   = file;
83
84
  p->stream = std::make_shared<radix::eafstream>(p->file.c_str(),
                                                 mode | std::ios::binary);
85
86
87
88
89
90
91
92
  if (!p->stream->is_open())
  {
    radix_line("Failed to open ARL file: " << file);
  }
}

bool ARLDataStream::read_index_header(ARLIndexHeader& iheader)
{
93
94
  bool result = false;

95
96
  std::string headerString =
      p->stream->readString(ARLDataStream::PImpl::indexHeaderLength);
97
98
99
100
101
102
103
104
105
  expand(headerString, iheader);

  radix_line("Read index header: year = " << iheader.year);
  radix_line("                  month = " << iheader.month);
  radix_line("                    day = " << iheader.day);
  radix_line("                   hour = " << iheader.hour);
  radix_line("               variable = " << iheader.kvar);
  radix_line("       scaling exponent = " << iheader.nexp);
  radix_line("          initial value = " << iheader.var1);
106
  radix_line("  header as string:\n  " << headerString);
107
108
109

  result = true;
  return result;
110
111
}

112
113
bool ARLDataStream::write_index_header(const ARLIndexHeader& iheader)
{
114
115
  bool result = false;

116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
  radix_line("Write index header: year = " << iheader.year);
  radix_line("                   month = " << iheader.month);
  radix_line("                     day = " << iheader.day);
  radix_line("                    hour = " << iheader.hour);
  radix_line("                variable = " << iheader.kvar);
  radix_line("        scaling exponent = " << iheader.nexp);
  radix_line("           initial value = " << iheader.var1);

  // Construct index header string
  char indexHeader[50];
  sprintf(indexHeader, "%2d%2d%2d%2d%2d%2d%2d%4s%4d%14.7E%14.7E", iheader.year,
          iheader.month, iheader.day, iheader.hour, iheader.ic, iheader.il,
          std::stoi(iheader.cgrid.c_str()), iheader.kvar.c_str(), iheader.nexp,
          iheader.prec, iheader.var1);

  p->stream->writeString(std::string(indexHeader), 50);

  radix_line("  header as string:\n  " << indexHeader);

135
136
  result = true;
  return result;
137
}
138

139
140
141
142
bool ARLDataStream::read_next_index_header(ARLIndexHeader& iheader)
{
  int bytesToSkip = roundUpInt(p->stream->bytesRead(), p->recordSize) -
                    p->stream->bytesRead();
Purves, Murray's avatar
Purves, Murray committed
143
  radix_line("Skipping " << bytesToSkip << " bytes to the next index header");
144
145
146
147
148
149
150
151
152

  p->stream->skipBytes(bytesToSkip);

  return read_index_header(iheader);
}

bool ARLDataStream::read_record_header(const ARLIndexHeader& iheader,
                                       ARLRecordHeader& rheader)
{
153
154
  bool result = false;

155
156
  std::string headerString =
      p->stream->readString(ARLDataStream::PImpl::recordHeaderLength);
157
158
159
  expand(headerString, iheader, rheader);

  // Calculate and save record size
160
161
  p->recordSize =
      (rheader.nx * rheader.ny) + ARLDataStream::PImpl::indexHeaderLength;
162
163
164
165
166
167

  radix_line("  Read record header: model = " << rheader.model_id);
  radix_line("                         nx = " << rheader.nx);
  radix_line("                         ny = " << rheader.ny);
  radix_line("                         nz = " << rheader.nz);
  radix_line("        Size of each record = " << p->recordSize);
168
  radix_line("    header as string:\n    " << headerString);
169
170
171

  result = true;
  return result;
172
173
}

174
175
176
bool ARLDataStream::write_record_header(const ARLIndexHeader& iheader,
                                        const ARLRecordHeader& rheader)
{
177
178
  bool result = false;

179
  // Calculate and save record size
180
181
  p->recordSize =
      (rheader.nx * rheader.ny) + ARLDataStream::PImpl::indexHeaderLength;
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210

  radix_line("  Write record header: model = " << rheader.model_id);
  radix_line("                          nx = " << rheader.nx);
  radix_line("                          ny = " << rheader.ny);
  radix_line("                          nz = " << rheader.nz);
  radix_line("         Size of each record = " << p->recordSize);

  char recordHeader[108];
  sprintf(recordHeader,
          "%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",
          rheader.model_id.c_str(), rheader.icx, rheader.mn, rheader.pole_lat,
          rheader.pole_lon, rheader.ref_lat, rheader.ref_lon, rheader.size,
          rheader.orient, rheader.tang_lat, rheader.sync_xp, rheader.sync_yp,
          rheader.sync_lat, rheader.sync_lon, rheader.dummy, rheader.nx,
          rheader.ny, rheader.nz, rheader.z_flag, rheader.lenh);

  p->stream->writeString(std::string(recordHeader), 108);

  // 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");

  p->stream->writeString("", bytesToSkip);

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

211
212
  result = true;
  return result;
213
214
}

215
216
bool ARLDataStream::read_record(const ARLIndexHeader& iheader,
                                const ARLRecordHeader& rheader,
217
                                std::vector<std::vector<float> >& record)
218
{
219
220
  bool result = false;

221
  // Set up the vector size
Purves, Murray's avatar
Purves, Murray committed
222
  record.clear();
223
  record.resize(rheader.nx);
224
  for (std::vector<float>& vec : record)
225
226
227
228
229
  {
    vec.resize(rheader.ny);
  }

  // Calculate the scaling factor
230
  float scaleFactor = pow(2.0, 7.0 - (float)iheader.nexp), lastValue = 0.0;
231
232
233
234
235
236
237
238
239
240
241

  radix_line("    Reading record data:");
  for (int y = 0; y < rheader.ny; ++y)
  {
    for (int x = 0; x < rheader.nx; ++x)
    {
      // Read the raw value
      unsigned char ch = (unsigned char)p->stream->readChar();
      int packedValue  = ch - 127;

      // Calculate the unpacked value
242
      float unpackedValue = 0.0;
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
      // Get the correct 'last value' if we are at a 0 index
      if (x == 0)
      {
        if (y == 0)
        {
          lastValue = iheader.var1;
        }
        else
        {
          lastValue = record[x][y - 1];
        }
      }

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

      if ((x < 2 && y < 2) || ((rheader.nx - x < 3) && (rheader.ny - y < 3)))
      {
        radix_line("      [" << x << "," << y << "]: packed = " << packedValue
                             << ", unpacked = " << unpackedValue);
      }
    }
  }

268
269
  result = true;
  return result;
270
271
}

272
273
274
bool ARLDataStream::write_record(const ARLIndexHeader& iheader,
                                 const ARLRecordHeader& rheader,
                                 const std::vector<std::vector<float> >& record)
275
{
276
277
  bool result = false;

278
279
280
281
282
283
284
  if (record.size() == 0)
  {
    radix_line("No data here - returning false");
    return false;
  }

  // Calculate the scaling factor
285
  float scaleFactor = pow(2.0, 7.0 - (float)iheader.nexp), lastValue = 0.0;
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305

  radix_line("    Writing record data:");
  for (int y = 0; y < rheader.ny; ++y)
  {
    for (int x = 0; x < rheader.nx; ++x)
    {
      // Get the correct 'last value' if we are at a 0 index
      if (x == 0)
      {
        if (y == 0)
        {
          lastValue = iheader.var1;
        }
        else
        {
          lastValue = record[x][y - 1];
        }
      }

      // Calculate the packed value
306
307
      float unpackedValue = record[x][y];
      int packedInt       = (int)((unpackedValue - lastValue) * scaleFactor);
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322

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

      p->stream->writeChar(packedValue);

      lastValue = unpackedValue;

      if ((x < 2 && y < 2) || ((rheader.nx - x < 3) && (rheader.ny - y < 3)))
      {
        radix_line("      [" << x << "," << y << "]: packed = " << packedInt
                             << ", unpacked = " << unpackedValue);
      }
    }
  }

323
324
  result = true;
  return result;
325
326
}

327
328
void ARLDataStream::close_stream() { p->stream->close(); }

329
}  // namespace radix