hysplitcdump.i.hh 9.55 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
namespace radix
{
9
template <typename data_type>
10
11
HysplitCDumpStream<data_type>::HysplitCDumpStream(data_type container)
{
12
  mData = container;
13
14
}

15
template <typename data_type>
16
17
bool HysplitCDumpStream<data_type>::read_from(const std::string &file)
{
18
19
  bool result = false;
  eafstream fstr(file.c_str(), std::ios::in | std::ios::binary);
20
21
  if (fstr.is_open() == false)
  {
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
    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
41
42
  for (int i = 0; i < numLocations; ++i)
  {
43
    fstr.readHeader();
44
45
    fstr >> year >> month >> day >> hour >> lat >> lon >> z >> minutes;
    mData->addLocation(year, month, day, hour, lat, lon, z, minutes);
46
    fstr.readFooter();
47
    radix_ensure(fstr.header() == fstr.footer());
48
  }
49

50
51
52
53
54
55
56
57
58
59
60
61
  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);
62

63
64
65
  fstr.readHeader();
  int numLevels;
  fstr >> numLevels;
66
67
  for (int i = 0; i < numLevels; ++i)
  {
68
69
70
71
72
73
74
75
76
77
    int level;
    fstr >> level;
    mData->addLevel(level);
  }
  fstr.readFooter();
  radix_ensure(fstr.header() == fstr.footer());

  fstr.readHeader();
  int numPollutants;
  fstr >> numPollutants;
78
79
  for (int i = 0; i < numPollutants; ++i)
  {
80
81
82
83
84
85
    std::string pol;
    pol = fstr.readString(4);
    mData->addPollutant(pol);
  }
  fstr.readFooter();
  radix_ensure(fstr.header() == fstr.footer());
86

87
88
  // infinite loop over observation times
  int timeIndex = 0;
89
90
  for (; fstr.good();)
  {
91
    fstr.readHeader();
92
    fstr >> year >> month >> day >> hour >> minutes >> forecastHour;
93
    fstr.readFooter();
94
    radix_ensure(fstr.header() == fstr.footer());
95
    mData->addStartTime(year, month, day, hour, minutes, forecastHour);
96
97

    fstr.readHeader();
98
    fstr >> year >> month >> day >> hour >> minutes >> forecastHour;
99
    fstr.readFooter();
100
    radix_ensure(fstr.header() == fstr.footer());
101
    mData->addEndTime(year, month, day, hour, minutes, forecastHour);
102

103
104
105
106
    for (int i = 0; i < numPollutants; ++i)
    {
      for (int j = 0; j < numLevels; ++j)
      {
107
        fstr.readHeader();
108
109
110
        std::string pollutant;
        int z;
        int nxyp;
111
112
        if (packing)
        {
113
114
115
          pollutant = fstr.readString(4);
          fstr >> z >> nxyp;
          mData->setNumNonZeroPoints(timeIndex, i, j, nxyp);
116
117
          for (int np = 0; np < nxyp; ++np)
          {
118
119
120
121
122
123
            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
124
        fstr.readFooter();
125
        radix_ensure(fstr.header() == fstr.footer());
126
127
128
129
130
131
132
133
        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
134

135
template <typename data_type>
136
137
bool HysplitCDumpStream<data_type>::write_to(const std::string &file) const
{
138
139
  bool result = false;
  eafstream fstr(file.c_str(), std::ios::out | std::ios::binary);
140
141
  if (fstr.is_open() == false)
  {
142
    return result;
143
144
145
146
  }
  // big endian
  fstr.setReverseBytes(true);
  int record_length;
147

148
149
150
151
152
153
154
155
156
157
158
  std::string id = mData->id();
  {
    // force string to be only 4 characters
    char cid[5];
    sprintf(cid, "%4s", id.c_str());
    cid[4] = '\0';  // null terminate
    id     = cid;
  }
  record_length = id.size() + sizeof(int) * 7;
  int year, month, day, hour, forecastHour, minutes, numLocations, packing = 1;
  mData->dateTime(year, month, day, hour);
159
160
161
162
163
164
165
  //
  // Adjust the date from YYYY to YY format.
  // HYSPLIT utils always assume YY format
  if (year > 99)
  {
    year = year % 100;
  }
166
167
  forecastHour = mData->forecastHour();
  numLocations = mData->numLocations();
168

169
170
171
172
173
174
175
  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);
176

177
  float lat, lon, z;
178
179
  for (int i = 0; i < numLocations; ++i)
  {
180
181
182
183
184
185
186
    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;
    }
187
    record_length = sizeof(int) * 5 + sizeof(float) * 3;
188
    radix_tagged_line("record_length=" << record_length);
189
    fstr.writeHeader(record_length);
190
    fstr << year << month << day << hour << lat << lon << z << minutes;
191
    radix_tagged_line("bytes_written=" << fstr.bytesWritten());
192
    fstr.writeFooter(record_length);
193
194
195
196
197
198
199
200
201
  }
  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();
202

203
204
205
206
207
208
  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);
209

210
211
212
213
214
  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;
215
216
  for (int i = 0; i < numLevels; ++i)
  {
217
218
219
220
221
222
223
224
225
226
227
    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;
228
229
  for (int i = 0; i < numPollutants; ++i)
  {
230
231
232
233
234
235
236
237
238
239
240
241
    std::string pol = mData->pollutant(i);
    radix_tagged_line("pollutant=" << pol);
    // force string to be only 4 characters
    char cid[5];
    sprintf(cid, "%4s", pol.c_str());
    cid[4]          = '\0';
    std::string tmp = cid;
    fstr.writeString(tmp);
    local_pols.push_back(tmp);
  }
  radix_tagged_line("bytes_written=" << fstr.bytesWritten());
  fstr.writeFooter(record_length);
242

243
244
245
246
247
  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.");
248
249
  for (int timeIndex = 0; timeIndex < numStartTimes; ++timeIndex)
  {
250
251
252
253
254
255
256
257
    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;
    }
258
    record_length = sizeof(int) * 6;
259
    radix_tagged_line("record_length=" << record_length);
260
    fstr.writeHeader(record_length);
261
    fstr << year << month << day << hour << minutes << forecastHour;
262
    radix_tagged_line("bytes_written=" << fstr.bytesWritten());
263
264
    fstr.writeFooter(record_length);

265
266
267
268
269
270
271
272
    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;
    }
273
    record_length = sizeof(int) * 6;
274
    radix_tagged_line("record_length=" << record_length);
275
    fstr.writeHeader(record_length);
276
    fstr << year << month << day << hour << minutes << forecastHour;
277
    radix_tagged_line("bytes_written=" << fstr.bytesWritten());
278
279
    fstr.writeFooter(record_length);

280
281
282
283
    for (int i = 0; i < numPollutants; ++i)
    {
      for (int j = 0; j < numLevels; ++j)
      {
284
285
286
287
288
289
290
291
        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);
292
        fstr.writeHeader(record_length);
293
294
        fstr.writeString(local_pols[size_t(i)]);
        fstr << level << nxyp;
295
296
        for (int np = 0; np < nxyp; ++np)
        {
297
298
299
300
301
          short xi, yi;
          float c;
          mData->point(timeIndex, i, j, np, xi, yi, c);
          fstr << xi << yi << c;
        }
302
        radix_tagged_line("bytes_written=" << fstr.bytesWritten());
303
        fstr.writeFooter(record_length);
304
      }
305
    }
306
307
308
309
  }
  fstr.close();
  result = true;
  return result;
310
311
}

312
}  // namespace radix