hysplitcdump.i.hh 8.87 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
namespace radix {
8

9
10
11
template <typename data_type>
HysplitCDumpStream<data_type>::HysplitCDumpStream(data_type container) {
  mData = container;
12
13
}

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

46
47
48
49
50
51
52
53
54
55
56
57
  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);
58

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

  fstr.readHeader();
  int numPollutants;
  fstr >> numPollutants;
  for (int i = 0; i < numPollutants; ++i) {
    std::string pol;
    pol = fstr.readString(4);
    mData->addPollutant(pol);
  }
  fstr.readFooter();
  radix_ensure(fstr.header() == fstr.footer());
80

81
82
83
  // infinite loop over observation times
  int timeIndex = 0;
  for (; fstr.good();) {
84
    fstr.readHeader();
85
    fstr >> year >> month >> day >> hour >> minutes >> forecastHour;
86
    fstr.readFooter();
87
    radix_ensure(fstr.header() == fstr.footer());
88
    mData->addStartTime(year, month, day, hour, minutes, forecastHour);
89
90

    fstr.readHeader();
91
    fstr >> year >> month >> day >> hour >> minutes >> forecastHour;
92
    fstr.readFooter();
93
    radix_ensure(fstr.header() == fstr.footer());
94
    mData->addEndTime(year, month, day, hour, minutes, forecastHour);
95

96
97
    for (int i = 0; i < numPollutants; ++i) {
      for (int j = 0; j < numLevels; ++j) {
98
        fstr.readHeader();
99
100
101
102
103
104
105
106
107
108
109
110
111
112
        std::string pollutant;
        int z;
        int nxyp;
        if (packing) {
          pollutant = fstr.readString(4);
          fstr >> z >> nxyp;
          mData->setNumNonZeroPoints(timeIndex, i, j, nxyp);
          for (int np = 0; np < nxyp; ++np) {
            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
113
        fstr.readFooter();
114
        radix_ensure(fstr.header() == fstr.footer());
115
116
117
118
119
120
121
122
        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
123

124
125
126
127
128
template <typename data_type>
bool HysplitCDumpStream<data_type>::write_to(const std::string &file) const {
  bool result = false;
  eafstream fstr(file.c_str(), std::ios::out | std::ios::binary);
  if (fstr.is_open() == false) {
129
    return result;
130
131
132
133
  }
  // big endian
  fstr.setReverseBytes(true);
  int record_length;
134

135
136
137
138
139
140
141
142
143
144
145
146
147
  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);
  forecastHour = mData->forecastHour();
  numLocations = mData->numLocations();
148

149
150
151
152
153
154
155
  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);
156

157
158
159
160
  float lat, lon, z;
  for (int i = 0; i < numLocations; ++i) {
    mData->location(i, year, month, day, hour, lat, lon, z, minutes);
    record_length = sizeof(int) * 5 + sizeof(float) * 3;
161
    radix_tagged_line("record_length=" << record_length);
162
    fstr.writeHeader(record_length);
163
    fstr << year << month << day << hour << lat << lon << z << minutes;
164
    radix_tagged_line("bytes_written=" << fstr.bytesWritten());
165
    fstr.writeFooter(record_length);
166
167
168
169
170
171
172
173
174
  }
  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();
175

176
177
178
179
180
181
  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);
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
211
212
  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;
  for (int i = 0; i < numLevels; ++i) {
    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;
  for (int i = 0; i < numPollutants; ++i) {
    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);
213

214
215
216
217
218
219
220
221
  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.");
  for (int timeIndex = 0; timeIndex < numStartTimes; ++timeIndex) {
    mData->startTime(timeIndex, year, month, day, hour, minutes, forecastHour);
    record_length = sizeof(int) * 6;
222
    radix_tagged_line("record_length=" << record_length);
223
    fstr.writeHeader(record_length);
224
    fstr << year << month << day << hour << minutes << forecastHour;
225
    radix_tagged_line("bytes_written=" << fstr.bytesWritten());
226
227
    fstr.writeFooter(record_length);

228
229
    mData->endTime(timeIndex, year, month, day, hour, minutes, forecastHour);
    record_length = sizeof(int) * 6;
230
    radix_tagged_line("record_length=" << record_length);
231
    fstr.writeHeader(record_length);
232
    fstr << year << month << day << hour << minutes << forecastHour;
233
    radix_tagged_line("bytes_written=" << fstr.bytesWritten());
234
235
    fstr.writeFooter(record_length);

236
237
238
239
240
241
242
243
244
245
    for (int i = 0; i < numPollutants; ++i) {
      for (int j = 0; j < numLevels; ++j) {
        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);
246
        fstr.writeHeader(record_length);
247
248
249
250
251
252
253
254
        fstr.writeString(local_pols[size_t(i)]);
        fstr << level << nxyp;
        for (int np = 0; np < nxyp; ++np) {
          short xi, yi;
          float c;
          mData->point(timeIndex, i, j, np, xi, yi, c);
          fstr << xi << yi << c;
        }
255
        radix_tagged_line("bytes_written=" << fstr.bytesWritten());
256
        fstr.writeFooter(record_length);
257
      }
258
    }
259
260
261
262
  }
  fstr.close();
  result = true;
  return result;
263
264
}

265
}  // namespace radix