hysplitcdump.i.hh 9.41 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
42
    mData->setDateTime(year, month, day, hour);
    mData->setForecastHour(forecastHour);
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102

    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();
    }

    int numLat, numLon;
    float deltaLat, deltaLon, cLat, cLon;
    fstr.readHeader();
    fstr >> numLat
            >> numLon
            >> deltaLat
            >> deltaLon
            >> cLat
            >> cLon;
    fstr.readFooter();
    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();

    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();


    // infinite loop over observation times
103
    int timeIndex = 0;
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
    for(; fstr.good();)
    {
        fstr.readHeader();
        fstr >> year
                >> month
                >> day
                >> hour
                >> minutes
                >> forecastHour;
        fstr.readFooter();
        mData->addStartTime(year, month, day, hour, minutes, forecastHour);

        fstr.readHeader();
        fstr >> year
                >> month
                >> day
                >> hour
                >> minutes
                >> forecastHour;
        fstr.readFooter();
        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;
138
                    mData->setNumNonZeroPoints(timeIndex,i,j,nxyp);
139
140
141
142
143
                    for(int np = 0; np < nxyp; ++np)
                    {
                        short xi = fstr.readShort();
                        short yi = fstr.readShort();
                        float c = fstr.readFloat();
144
                        mData->setPoint(timeIndex,i,j,np,xi,yi,c);
145
146
147
148
149
150
                    }
                } // if packing is enabled
                fstr.readFooter();
                fstr.peek(); // peek ahead to check file integrity
            } // loop over num levels
        } // loop num pollutants
151
        timeIndex++;
152
    } // infinite loop over output times
153
    result = true;
154
155
156
    return result;
} // HysplitCDump::read_from

157
template<typename data_type>
158
bool HysplitCDumpStream<data_type>::write_to(const std::string& file) const
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
{
    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;
181
182
    mData->dateTime(year, month, day, hour);
    forecastHour = mData->forecastHour();
183
184
    numLocations = mData->numLocations();

185
    radix_tagged_line("record_length=" << record_length);
186
    fstr.writeHeader(record_length);
187
    fstr.writeString(id);
188
    fstr << year << month << day << hour << forecastHour << numLocations << packing;
189
    radix_tagged_line("bytes_written=" << fstr.bytesWritten());
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
    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;
        fstr.writeHeader(record_length);
        fstr << year << month << day << hour << lat << lon << z << minutes;
        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();
209

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

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

    int numPollutants = mData->numPollutants();
228
229
    record_length = int(sizeof(int))+numPollutants*4;
    radix_tagged_line("record_length=" << record_length);
230
    fstr.writeHeader(record_length);
231
    fstr << numPollutants;
232
233
234
235
236
237
238
239
240
241
242
243
    std::vector<std::string> local_pols;
    for(int i = 0; i < numPollutants; ++i)
    {
        std::string pol = mData->pollutant(i);
        // force string to be only 4 characters
        char cid[5];
        sprintf(cid, "%4s", pol.c_str());
        cid[4] = '\0'; // null terminate
        pol = cid;
        fstr.writeString(pol);
        local_pols.push_back(pol);
    }
244
    radix_tagged_line("bytes_written=" << fstr.bytesWritten());
245
246
247
248
249
250
251
252
253
254
    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;
255
        radix_tagged_line("record_length=" << record_length);
256
257
        fstr.writeHeader(record_length);
        fstr << year << month << day << hour << minutes << forecastHour;
258
        radix_tagged_line("bytes_written=" << fstr.bytesWritten());
259
260
261
262
        fstr.writeFooter(record_length);

        mData->endTime(timeIndex, year, month, day, hour, minutes, forecastHour);
        record_length = sizeof(int)*6;
263
        radix_tagged_line("record_length=" << record_length);
264
265
        fstr.writeHeader(record_length);
        fstr << year << month << day << hour << minutes << forecastHour;
266
        radix_tagged_line("bytes_written=" << fstr.bytesWritten());
267
268
269
270
271
272
273
274
275
276
277
        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;
278
279
280
281
                radix_tagged_line("record_length=" << record_length
                        << " pol=" << local_pols[size_t(i)]
                        << " z=" << level
                        << " nxyp=" << nxyp);
282
283
284
285
286
287
288
289
290
291
                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;
                }
292
                radix_tagged_line("bytes_written=" << fstr.bytesWritten());
293
294
295
296
                fstr.writeFooter(record_length);
            }
        }
    }
297
    fstr.close();
298
299
300
301
    result = true;
    return result;
}

302
} // namespace radix