hysplitcdump.i.hh 9.71 KB
Newer Older
1
#include <cmath>
2
#include <cstdio>
3
4
5
6
#include <string>
#include "radixbug/bug.hh"
#include "radixio/eafstream.hh"
#include "radixio/hysplitcdump.hh"
7

8
#ifdef _MSC_VER
9
10
11
12
13
#if _MSC_VER < 1900
#define snprintf _snprintf
#endif
#endif

14
15
namespace radix
{
16
template <typename data_type>
17
18
HysplitCDumpStream<data_type>::HysplitCDumpStream(data_type container)
{
19
  mData = container;
20
21
}

22
template <typename data_type>
23
24
bool HysplitCDumpStream<data_type>::read_from(const std::string &file)
{
25
26
  bool result = false;
  eafstream fstr(file.c_str(), std::ios::in | std::ios::binary);
27
28
  if (fstr.is_open() == false)
  {
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
    return result;
  }
  // big endian
  fstr.setReverseBytes(true);
  // read fortran header
  fstr.readHeader();
  std::string id = fstr.readString(4);
  mData->setId(id);
  int packing, year, month, day, hour, forecastHour, numLocations;
  fstr >> year >> month >> day >> hour >> forecastHour >> numLocations >>
      packing;
  fstr.readFooter();
  radix_ensure(fstr.header() == fstr.footer());
  mData->setDateTime(year, month, day, hour);
  mData->setForecastHour(forecastHour);

  float lat, lon, z;
  int minutes;
  // for the number of locations
48
49
  for (int i = 0; i < numLocations; ++i)
  {
50
    fstr.readHeader();
51
52
    fstr >> year >> month >> day >> hour >> lat >> lon >> z >> minutes;
    mData->addLocation(year, month, day, hour, lat, lon, z, minutes);
53
    fstr.readFooter();
54
    radix_ensure(fstr.header() == fstr.footer());
55
  }
56

57
58
59
60
61
62
63
64
65
66
67
68
  int numLat, numLon;
  float deltaLat, deltaLon, cLat, cLon;
  fstr.readHeader();
  fstr >> numLat >> numLon >> deltaLat >> deltaLon >> cLat >> cLon;
  fstr.readFooter();
  radix_ensure(fstr.header() == fstr.footer());
  mData->setNumLat(numLat);
  mData->setNumLon(numLon);
  mData->setDeltaLat(deltaLat);
  mData->setDeltaLon(deltaLon);
  mData->setCornerLat(cLat);
  mData->setCornerLon(cLon);
69

70
71
72
  fstr.readHeader();
  int numLevels;
  fstr >> numLevels;
73
74
  for (int i = 0; i < numLevels; ++i)
  {
75
76
77
78
79
80
81
82
83
84
    int level;
    fstr >> level;
    mData->addLevel(level);
  }
  fstr.readFooter();
  radix_ensure(fstr.header() == fstr.footer());

  fstr.readHeader();
  int numPollutants;
  fstr >> numPollutants;
85
86
  for (int i = 0; i < numPollutants; ++i)
  {
87
88
89
90
91
92
    std::string pol;
    pol = fstr.readString(4);
    mData->addPollutant(pol);
  }
  fstr.readFooter();
  radix_ensure(fstr.header() == fstr.footer());
93

94
95
  // infinite loop over observation times
  int timeIndex = 0;
96
97
  for (; fstr.good();)
  {
98
    fstr.readHeader();
99
    fstr >> year >> month >> day >> hour >> minutes >> forecastHour;
100
    fstr.readFooter();
101
    radix_ensure(fstr.header() == fstr.footer());
102
    mData->addStartTime(year, month, day, hour, minutes, forecastHour);
103
104

    fstr.readHeader();
105
    fstr >> year >> month >> day >> hour >> minutes >> forecastHour;
106
    fstr.readFooter();
107
    radix_ensure(fstr.header() == fstr.footer());
108
    mData->addEndTime(year, month, day, hour, minutes, forecastHour);
109

110
111
112
113
    for (int i = 0; i < numPollutants; ++i)
    {
      for (int j = 0; j < numLevels; ++j)
      {
114
        fstr.readHeader();
115
116
117
        std::string pollutant;
        int z;
        int nxyp;
118
119
        if (packing)
        {
120
121
122
          pollutant = fstr.readString(4);
          fstr >> z >> nxyp;
          mData->setNumNonZeroPoints(timeIndex, i, j, nxyp);
123
124
          for (int np = 0; np < nxyp; ++np)
          {
125
126
127
128
129
130
            short xi = fstr.readShort();
            short yi = fstr.readShort();
            float c  = fstr.readFloat();
            mData->setPoint(timeIndex, i, j, np, xi, yi, c);
          }
        }  // if packing is enabled
131
        fstr.readFooter();
132
        radix_ensure(fstr.header() == fstr.footer());
133
134
135
136
137
138
139
140
        fstr.peek();  // peek ahead to check file integrity
      }               // loop over num levels
    }                 // loop num pollutants
    timeIndex++;
  }  // infinite loop over output times
  result = true;
  return result;
}  // HysplitCDump::read_from
141

142
template <typename data_type>
143
144
bool HysplitCDumpStream<data_type>::write_to(const std::string &file) const
{
145
146
  bool result = false;
  eafstream fstr(file.c_str(), std::ios::out | std::ios::binary);
147
148
  if (fstr.is_open() == false)
  {
149
    return result;
150
151
152
153
  }
  // big endian
  fstr.setReverseBytes(true);
  int record_length;
154

155
156
  std::string id = mData->id();
  {
Norby, Tom's avatar
Norby, Tom committed
157
    // force string to be no more than 4 characters
158
    char cid[5];
Norby, Tom's avatar
Norby, Tom committed
159
    snprintf(cid, sizeof(cid), "%4s", id.c_str());
160
161
    cid[4] = '\0';  // null terminate
    id     = cid;
162
163
164
165
  }
  record_length = id.size() + sizeof(int) * 7;
  int year, month, day, hour, forecastHour, minutes, numLocations, packing = 1;
  mData->dateTime(year, month, day, hour);
166
167
168
169
170
171
172
  //
  // Adjust the date from YYYY to YY format.
  // HYSPLIT utils always assume YY format
  if (year > 99)
  {
    year = year % 100;
  }
173
174
  forecastHour = mData->forecastHour();
  numLocations = mData->numLocations();
175

176
177
178
179
180
181
182
  radix_tagged_line("record_length=" << record_length);
  fstr.writeHeader(record_length);
  fstr.writeString(id);
  fstr << year << month << day << hour << forecastHour << numLocations
       << packing;
  radix_tagged_line("bytes_written=" << fstr.bytesWritten());
  fstr.writeFooter(record_length);
183

184
  float lat, lon, z;
185
186
  for (int i = 0; i < numLocations; ++i)
  {
187
188
189
190
191
192
193
    mData->location(i, year, month, day, hour, lat, lon, z, minutes);  //
    // Adjust the date from YYYY to YY format.
    // HYSPLIT utils always assume YY format
    if (year > 99)
    {
      year = year % 100;
    }
194
    record_length = sizeof(int) * 5 + sizeof(float) * 3;
195
    radix_tagged_line("record_length=" << record_length);
196
    fstr.writeHeader(record_length);
197
    fstr << year << month << day << hour << lat << lon << z << minutes;
198
    radix_tagged_line("bytes_written=" << fstr.bytesWritten());
199
    fstr.writeFooter(record_length);
200
201
202
203
204
205
206
207
208
  }
  int numLat, numLon;
  float deltaLat, deltaLon, cLat, cLon;
  numLat   = mData->numLat();
  numLon   = mData->numLon();
  deltaLat = mData->deltaLat();
  deltaLon = mData->deltaLon();
  cLat     = mData->cornerLat();
  cLon     = mData->cornerLon();
209

210
211
212
213
214
215
  record_length = sizeof(int) * 2 + sizeof(float) * 4;
  radix_tagged_line("record_length=" << record_length);
  fstr.writeHeader(record_length);
  fstr << numLat << numLon << deltaLat << deltaLon << cLat << cLon;
  radix_tagged_line("bytes_written=" << fstr.bytesWritten());
  fstr.writeFooter(record_length);
216

217
218
219
220
221
  int numLevels = mData->numLevels();
  record_length = int(sizeof(int)) + int(sizeof(int)) * numLevels;
  radix_tagged_line("record_length=" << record_length);
  fstr.writeHeader(record_length);
  fstr << numLevels;
222
223
  for (int i = 0; i < numLevels; ++i)
  {
224
225
226
227
228
229
230
231
232
233
234
    fstr << mData->level(i);
  }
  radix_tagged_line("bytes_written=" << fstr.bytesWritten());
  fstr.writeFooter(record_length);

  int numPollutants = mData->numPollutants();
  record_length     = int(sizeof(int)) + numPollutants * 4;
  radix_tagged_line("record_length=" << record_length);
  fstr.writeHeader(record_length);
  fstr << numPollutants;
  std::vector<std::string> local_pols;
235
236
  for (int i = 0; i < numPollutants; ++i)
  {
237
238
    std::string pol = mData->pollutant(i);
    radix_tagged_line("pollutant=" << pol);
Norby, Tom's avatar
Norby, Tom committed
239
    // force string to be no more than 4 characters
240
    char icid[5];
Norby, Tom's avatar
Norby, Tom committed
241
    snprintf(icid, sizeof(icid), "%4s", pol.c_str());
242
    icid[4]         = '\0';  // null terminate
243
    std::string tmp = icid;
244
245
246
247
248
    fstr.writeString(tmp);
    local_pols.push_back(tmp);
  }
  radix_tagged_line("bytes_written=" << fstr.bytesWritten());
  fstr.writeFooter(record_length);
249

250
251
252
253
254
  int numStartTimes = mData->numStartTimes();
  int numEndTimes   = mData->numEndTimes();
  radix_insist(
      numStartTimes == numEndTimes,
      "The number of start times did not match the number of end times.");
255
256
  for (int timeIndex = 0; timeIndex < numStartTimes; ++timeIndex)
  {
257
258
259
260
261
262
263
264
    mData->startTime(timeIndex, year, month, day, hour, minutes,
                     forecastHour);  //
    // Adjust the date from YYYY to YY format.
    // HYSPLIT utils always assume YY format
    if (year > 99)
    {
      year = year % 100;
    }
265
    record_length = sizeof(int) * 6;
266
    radix_tagged_line("record_length=" << record_length);
267
    fstr.writeHeader(record_length);
268
    fstr << year << month << day << hour << minutes << forecastHour;
269
    radix_tagged_line("bytes_written=" << fstr.bytesWritten());
270
271
    fstr.writeFooter(record_length);

272
273
274
275
276
277
278
279
    mData->endTime(timeIndex, year, month, day, hour, minutes,
                   forecastHour);  //
    // Adjust the date from YYYY to YY format.
    // HYSPLIT utils always assume YY format
    if (year > 99)
    {
      year = year % 100;
    }
280
    record_length = sizeof(int) * 6;
281
    radix_tagged_line("record_length=" << record_length);
282
    fstr.writeHeader(record_length);
283
    fstr << year << month << day << hour << minutes << forecastHour;
284
    radix_tagged_line("bytes_written=" << fstr.bytesWritten());
285
286
    fstr.writeFooter(record_length);

287
288
289
290
    for (int i = 0; i < numPollutants; ++i)
    {
      for (int j = 0; j < numLevels; ++j)
      {
291
292
293
294
295
296
297
298
        int level     = mData->level(j);
        int nxyp      = mData->numNonZeroPoints(timeIndex, i, j);
        record_length = 4 + int(sizeof(int)) * 2 +
                        int(sizeof(short)) * 2 * nxyp +
                        int(sizeof(float)) * nxyp;
        radix_tagged_line("record_length="
                          << record_length << " pol=" << local_pols[size_t(i)]
                          << " z=" << level << " nxyp=" << nxyp);
299
        fstr.writeHeader(record_length);
300
301
        fstr.writeString(local_pols[size_t(i)]);
        fstr << level << nxyp;
302
303
        for (int np = 0; np < nxyp; ++np)
        {
304
305
306
307
308
          short xi, yi;
          float c;
          mData->point(timeIndex, i, j, np, xi, yi, c);
          fstr << xi << yi << c;
        }
309
        radix_tagged_line("bytes_written=" << fstr.bytesWritten());
310
        fstr.writeFooter(record_length);
311
      }
312
    }
313
314
315
316
  }
  fstr.close();
  result = true;
  return result;
317
}  // namespace radix
318

319
}  // namespace radix