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