hysplitcdump.i.hh 10.1 KB
Newer Older
1
2
#include <cstdio>
#include <cmath>
3
4
5
6
7
8
9
10
#include <string>
#include "radixbug/bug.hh"
#include "radixio/eafstream.hh"
#include "radixio/hysplitcdump.hh"
namespace radix
{

template<typename data_type>
11
HysplitCDumpStream<data_type>::HysplitCDumpStream(data_type container)
12
13
14
15
16
{
    mData = container;
}

template<typename data_type>
17
bool HysplitCDumpStream<data_type>::read_from(const std::string& file)
18
19
20
{
    bool result = false;
    eafstream fstr(file.c_str()
21
                   , std::ios::in | std::ios::binary);
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
    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();
41
    radix_ensure(fstr.header() == fstr.footer());
42
43
    mData->setDateTime(year, month, day, hour);
    mData->setForecastHour(forecastHour);
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60

    float lat, lon, z;
    int minutes;
    // for the number of locations
    for( int i = 0; i < numLocations; ++i)
    {
        fstr.readHeader();
        fstr >> year
                >> month
                >> day
                >> hour
                >> lat
                >> lon
                >> z
                >> minutes;
        mData->addLocation(year, month, day, hour, lat, lon, z, minutes);
        fstr.readFooter();
61
        radix_ensure(fstr.header() == fstr.footer());
62
63
64
65
66
67
68
69
70
71
72
73
    }

    int numLat, numLon;
    float deltaLat, deltaLon, cLat, cLon;
    fstr.readHeader();
    fstr >> numLat
            >> numLon
            >> deltaLat
            >> deltaLon
            >> cLat
            >> cLon;
    fstr.readFooter();
74
    radix_ensure(fstr.header() == fstr.footer());
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
    mData->setNumLat(numLat);
    mData->setNumLon(numLon);
    mData->setDeltaLat(deltaLat);
    mData->setDeltaLon(deltaLon);
    mData->setCornerLat(cLat);
    mData->setCornerLon(cLon);

    fstr.readHeader();
    int numLevels;
    fstr >> numLevels;
    for(int i = 0; i < numLevels; ++i)
    {
        int level;
        fstr >> level;
        mData->addLevel(level);
    }
    fstr.readFooter();
92
    radix_ensure(fstr.header() == fstr.footer());
93
94
95
96
97
98
99
100
101
102
103

    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();
104
    radix_ensure(fstr.header() == fstr.footer());
105
106
107


    // infinite loop over observation times
108
    int timeIndex = 0;
109
110
111
112
113
114
115
116
117
118
    for(; fstr.good();)
    {
        fstr.readHeader();
        fstr >> year
                >> month
                >> day
                >> hour
                >> minutes
                >> forecastHour;
        fstr.readFooter();
119
        radix_ensure(fstr.header() == fstr.footer());
120
121
122
123
124
125
126
127
128
129
        mData->addStartTime(year, month, day, hour, minutes, forecastHour);

        fstr.readHeader();
        fstr >> year
                >> month
                >> day
                >> hour
                >> minutes
                >> forecastHour;
        fstr.readFooter();
130
        radix_ensure(fstr.header() == fstr.footer());
131
132
133
134
135
136
137
138
139
140
141
142
143
144
        mData->addEndTime(year, month, day, hour, minutes, forecastHour);

        for(int i = 0; i < numPollutants; ++i)
        {
            for(int j = 0; j < numLevels; ++j)
            {
                fstr.readHeader();
                std::string pollutant;
                int z;
                int nxyp;
                if(packing)
                {
                    pollutant = fstr.readString(4);
                    fstr >> z >> nxyp;
145
                    mData->setNumNonZeroPoints(timeIndex,i,j,nxyp);
146
147
148
149
150
                    for(int np = 0; np < nxyp; ++np)
                    {
                        short xi = fstr.readShort();
                        short yi = fstr.readShort();
                        float c = fstr.readFloat();
151
                        mData->setPoint(timeIndex,i,j,np,xi,yi,c);
152
153
154
                    }
                } // if packing is enabled
                fstr.readFooter();
155
                radix_ensure(fstr.header() == fstr.footer());
156
157
158
                fstr.peek(); // peek ahead to check file integrity
            } // loop over num levels
        } // loop num pollutants
159
        timeIndex++;
160
    } // infinite loop over output times
161
    result = true;
162
163
164
    return result;
} // HysplitCDump::read_from

165
template<typename data_type>
166
bool HysplitCDumpStream<data_type>::write_to(const std::string& file) const
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
{
    bool result = false;
    eafstream fstr(file.c_str()
                   , std::ios::out | std::ios::binary);
    if(fstr.is_open() == false)
    {
        return result;
    }
    //big endian
    fstr.setReverseBytes(true);
    int record_length;

    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;
189
190
    mData->dateTime(year, month, day, hour);
    forecastHour = mData->forecastHour();
191
192
    numLocations = mData->numLocations();

193
    radix_tagged_line("record_length=" << record_length);
194
    fstr.writeHeader(record_length);
195
    fstr.writeString(id);
196
    fstr << year << month << day << hour << forecastHour << numLocations << packing;
197
    radix_tagged_line("bytes_written=" << fstr.bytesWritten());
198
199
200
201
202
203
204
    fstr.writeFooter(record_length);

    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;
205
        radix_tagged_line("record_length=" << record_length);
206
207
        fstr.writeHeader(record_length);
        fstr << year << month << day << hour << lat << lon << z << minutes;
208
        radix_tagged_line("bytes_written=" << fstr.bytesWritten());
209
210
211
212
213
214
215
216
217
218
        fstr.writeFooter(record_length);
    }
    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();
219

220
    record_length = sizeof(int)*2+sizeof(float)*4;
221
    radix_tagged_line("record_length=" << record_length);
222
223
    fstr.writeHeader(record_length);
    fstr << numLat << numLon << deltaLat << deltaLon << cLat << cLon;
224
    radix_tagged_line("bytes_written=" << fstr.bytesWritten());
225
226
227
    fstr.writeFooter(record_length);

    int numLevels = mData->numLevels();
228
229
    record_length = int(sizeof(int))+int(sizeof(int))*numLevels;
    radix_tagged_line("record_length=" << record_length);
230
    fstr.writeHeader(record_length);
231
    fstr << numLevels;
232
233
234
235
    for(int i = 0; i < numLevels; ++i)
    {
        fstr << mData->level(i);
    }
236
    radix_tagged_line("bytes_written=" << fstr.bytesWritten());
237
238
239
    fstr.writeFooter(record_length);

    int numPollutants = mData->numPollutants();
240
241
    record_length = int(sizeof(int))+numPollutants*4;
    radix_tagged_line("record_length=" << record_length);
242
    fstr.writeHeader(record_length);
243
    fstr << numPollutants;
244
245
246
247
    std::vector<std::string> local_pols;
    for(int i = 0; i < numPollutants; ++i)
    {
        std::string pol = mData->pollutant(i);
248
        radix_tagged_line("pollutant=" << pol);
249
250
251
        // force string to be only 4 characters
        char cid[5];
        sprintf(cid, "%4s", pol.c_str());
252
253
254
255
        cid[4] = '\0';
        std::string tmp = cid;
        fstr.writeString(tmp);
        local_pols.push_back(tmp);
256
    }
257
    radix_tagged_line("bytes_written=" << fstr.bytesWritten());
258
259
260
261
262
263
264
265
266
267
    fstr.writeFooter(record_length);

    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;
268
        radix_tagged_line("record_length=" << record_length);
269
270
        fstr.writeHeader(record_length);
        fstr << year << month << day << hour << minutes << forecastHour;
271
        radix_tagged_line("bytes_written=" << fstr.bytesWritten());
272
273
274
275
        fstr.writeFooter(record_length);

        mData->endTime(timeIndex, year, month, day, hour, minutes, forecastHour);
        record_length = sizeof(int)*6;
276
        radix_tagged_line("record_length=" << record_length);
277
278
        fstr.writeHeader(record_length);
        fstr << year << month << day << hour << minutes << forecastHour;
279
        radix_tagged_line("bytes_written=" << fstr.bytesWritten());
280
281
282
283
284
285
286
287
288
289
290
        fstr.writeFooter(record_length);

        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;
291
292
293
294
                radix_tagged_line("record_length=" << record_length
                        << " pol=" << local_pols[size_t(i)]
                        << " z=" << level
                        << " nxyp=" << nxyp);
295
296
297
298
299
300
301
302
303
304
                fstr.writeHeader(record_length);
                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;
                }
305
                radix_tagged_line("bytes_written=" << fstr.bytesWritten());
306
307
308
309
                fstr.writeFooter(record_length);
            }
        }
    }
310
    fstr.close();
311
312
313
314
    result = true;
    return result;
}

315
} // namespace radix