Commit 7c4e68e7 authored by Lefebvre, Jordan's avatar Lefebvre, Jordan
Browse files

WIP #13. Initial crack at writing hysplit cdump file. Need to add testing.

parent af6c865e
......@@ -29,8 +29,8 @@ public:
* setDateTime(int year, int month, int day, int hour, int minute)
* setForcastHour(int)
* addLocation(int year, int month, int day, int hour, float lat, float lon, float z, int minutes)
* setNumLat(float)
* setNumLon(float)
* setNumLat(int)
* setNumLon(int)
* setDeltaLat(float)
* setDeltaLon(float)
* setCornerLat(float)
......@@ -39,8 +39,8 @@ public:
* addPollutant(const std::string&)
* addStartTime(int year, int month, int day, int hour, int minute, int forecast)
* addEndTime(int year, int month, int day, int hour, int minute, int forecast)
* setNumNonZeroPoints(int pollutantIndex, int levelIndex, int numPoints)
* addPoint(int pollutantIndex, int levelIndex, int nonZeroIndex, short xIndex, short yIndex, float concentration)
* setNumNonZeroPoints(int timeIndex, int pollutantIndex, int levelIndex, int numPoints)
* addPoint(int timeIndex, int pollutantIndex, int levelIndex, int nonZeroIndex, short xIndex, short yIndex, float concentration)
* @return true on success, false otherwise
*/
bool read_from(const std::string& file);
......@@ -50,16 +50,29 @@ public:
*
* Requires data_type to have the following methods:
* id() const
* int year() const
* int month() const
* int day() const
* int hour() const
* int minute() const
* int dateTime(int& year, int& month, int& day, int& hour, int& minute) const
* int forecastHour() const
*
* int numLocations() const
* void location(int i, int& year, int& month, int& day, int& hour, float& lat, float& lon, float& z, int& minutes) const
* int numLat() const
* int numLon() const
* float deltaLat() const
* float deltaLon() const
* float cornerLat() const
* float cornerLon() const
* int numLevels() const
* int level(int) const
* int numPollutants() const
* std::string pollutant(int i) const
* int numStartTimes() const
* void startTime(int i, int& year, int& month, int& day, int& hour, int& minute, int& forecast) const
* int numEndTimes() const
* void endTime(int i, int& year, int& month, int& day, int& hour, int& minute, int& forecast) const
* int numNonZeroPoints(int timeIndex,int pollutantIndex, int levelIndex) const
* void point(int timeIndex,int pollutantIndex, int levelIndex, int nonZeroIndex, short& xIndex, short& yIndex, float& concentration) const
* @return true on success, false otherwise
*/
bool write_to(const std::string& file);
bool write_to(const std::string& file) const;
}; // class
} // namespace radix
......
#include <cstdio>
#include <cmath>
#include <string>
#include "radixbug/bug.hh"
#include "radixio/eafstream.hh"
......@@ -16,7 +18,7 @@ bool HysplitCDump<data_type>::read_from(const std::string& file)
{
bool result = false;
eafstream fstr(file.c_str()
, std::ios::in | std::ios::binary);
, std::ios::in | std::ios::binary);
if(fstr.is_open() == false)
{
return result;
......@@ -97,6 +99,7 @@ bool HysplitCDump<data_type>::read_from(const std::string& file)
// infinite loop over observation times
int timeIndex = 0;
for(; fstr.good();)
{
fstr.readHeader();
......@@ -131,21 +134,155 @@ bool HysplitCDump<data_type>::read_from(const std::string& file)
{
pollutant = fstr.readString(4);
fstr >> z >> nxyp;
mData->setNumNonZeroPoints(i,j,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->addPoint(i,j,np,xi,yi,c);
mData->addPoint(timeIndex,i,j,np,xi,yi,c);
}
} // if packing is enabled
fstr.readFooter();
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
template<typename data_type>
bool HysplitCDump<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)
{
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;
mData->dateTime(year, month, day, hour, forecastHour);
numLocations = mData->numLocations();
fstr.writeHeader(record_length);
fstr << year << month << day << hour << forecastHour << numLocations << packing;
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();
{
// 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 = 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();
record_length = int(sizeof(int))*numLevels;
fstr.writeHeader(record_length);
for(int i = 0; i < numLevels; ++i)
{
fstr << mData->level(i);
}
fstr.writeFooter(record_length);
int numPollutants = mData->numPollutants();
record_length = numPollutants*4;
fstr.writeHeader(record_length);
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);
}
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;
fstr.writeHeader(record_length);
fstr << year << month << day << hour << minutes << forecastHour;
fstr.writeFooter(record_length);
mData->endTime(timeIndex, year, month, day, hour, minutes, forecastHour);
record_length = sizeof(int)*6;
fstr.writeHeader(record_length);
fstr << year << month << day << hour << minutes << forecastHour;
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;
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;
}
fstr.writeFooter(record_length);
}
}
}
result = true;
return result;
}
} // namespace radix
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment