Skip to content

GitLab

Commit 48e123ea authored by Purves, Murray's avatar Purves, Murray
Browse files

Adding height as upper variable; radixsnd2arl now creates files that can be used by HYSPLIT

parent 05630fbd
Pipeline #16762 failed with stages
in 8 minutes and 39 seconds
Hide whitespace changes
Inline Side-by-side
Showing with 56 additions and 20 deletions
radixio/arldatastream.cc
View file @ 48e123ea
......@@ -117,11 +117,12 @@ bool ARLDataStream::write_record_header(const ARLRecordHeader& rheader)
// Remove century from year
int year = rheader.year % 100;
radix_line("Write index header: year = " << year);
radix_line("Write record header: year = " << year);
radix_line(" month = " << rheader.month);
radix_line(" day = " << rheader.day);
radix_line(" hour = " << rheader.hour);
radix_line(" variable = " << rheader.kvar);
radix_line(" level = " << rheader.il);
radix_line(" scaling exponent = " << rheader.nexp);
radix_line(" initial value = " << rheader.var1);
......@@ -144,7 +145,7 @@ bool ARLDataStream::read_next_record_header(ARLRecordHeader& rheader)
{
int bytesToSkip = roundUpInt(p->stream->bytesRead(), p->recordSize) -
p->stream->bytesRead();
radix_line("Skipping " << bytesToSkip << " bytes to the next index header");
radix_line("Skipping " << bytesToSkip << " bytes to the next record header");
p->stream->skipBytes(bytesToSkip);
......@@ -164,7 +165,7 @@ bool ARLDataStream::read_index_header(const ARLRecordHeader& rheader,
p->recordSize =
(iheader.nx * iheader.ny) + ARLDataStream::PImpl::recordHeaderLength;
radix_line(" Read record header: model = " << iheader.model_id);
radix_line(" Read index header: model = " << iheader.model_id);
radix_line(" nx = " << iheader.nx);
radix_line(" ny = " << iheader.ny);
radix_line(" nz = " << iheader.nz);
......@@ -184,7 +185,7 @@ bool ARLDataStream::write_index_header(const ARLRecordHeader& rheader,
p->recordSize =
(iheader.nx * iheader.ny) + ARLDataStream::PImpl::recordHeaderLength;
radix_line(" Write record header: model = " << iheader.model_id);
radix_line(" Write index header: model = " << iheader.model_id);
radix_line(" nx = " << iheader.nx);
radix_line(" ny = " << iheader.ny);
radix_line(" nz = " << iheader.nz);
......
radixio/examples/radixsnd2arl.cc
View file @ 48e123ea
......@@ -199,9 +199,10 @@ int main(int argc, char **argv)
// Header strings to denote fields
std::string pressureString = "pres", tempString = "temp", relHumString = "rh",
dewPtString = "tdew", wDirString = "wdir", wSpdString = "wspd";
dewPtString = "tdew", wDirString = "wdir", wSpdString = "wspd",
heightString = "zhgt";
int pressureIndex = -1, tempIndex = -1, relHumIndex = -1, wDirIndex = -1,
wSpdIndex = -1;
wSpdIndex = -1, heightIndex = -1;
// Open and read the csv
CSVFile inputCsv(inputCsvPath);
......@@ -255,6 +256,11 @@ int main(int argc, char **argv)
std::cout << " Found wind speed at index " << entry << std::endl;
wSpdIndex = entry;
}
else if (inputData[1][entry] == heightString)
{
std::cout << " Found height at index " << entry << std::endl;
heightIndex = entry;
}
}
// If we are missing one, we can't continue
if (pressureIndex == -1 || tempIndex == -1 || relHumIndex == -1 ||
......@@ -286,13 +292,17 @@ int main(int argc, char **argv)
std::cout << " Missing wind speed field (denoted by " << wSpdString
<< ")";
}
if (heightIndex == -1)
{
std::cout << " Missing height field (denoted by " << heightString << ")";
}
std::cout << "Please ensure these fields are present and rerun.";
return -4;
}
// Read the data
std::vector<float> inputPressures, inputTemps, inputRelHums, inputWDirs,
inputWSpds;
inputWSpds, inputHeights;
std::cout << "Data fields found - reading data..." << std::endl;
float missingValue = -9999.f;
for (int row = 4; row < inputData.size(); ++row)
......@@ -305,7 +315,7 @@ int main(int argc, char **argv)
}
bool foundPressure = false, foundTemp = false, foundRelHum = false,
foundWDir = false, foundWSpd = false;
foundWDir = false, foundWSpd = false, foundHeight = false;
for (int entry = 0; entry < inputData[row].size(); ++entry)
{
float thisValue = from_string(inputData[row][entry], missingValue);
......@@ -335,6 +345,11 @@ int main(int argc, char **argv)
foundWDir = true;
inputWDirs.push_back(thisValue);
}
else if (entry == heightIndex)
{
foundHeight = true;
inputHeights.push_back(thisValue);
}
}
// Add missing data if any of the above weren't found
if (!foundPressure)
......@@ -367,6 +382,12 @@ int main(int argc, char **argv)
<< std::endl;
inputWDirs.push_back(missingValue);
}
if (!foundHeight)
{
std::cout << " Warning: couldn't find height in row " << row
<< std::endl;
inputHeights.push_back(missingValue);
}
}
std::cout << "Data read complete: " << inputPressures.size()
<< " entries read." << std::endl;
......@@ -381,6 +402,8 @@ int main(int argc, char **argv)
interpolateValues(inputPressures, inputWSpds);
std::cout << "Wind direction:" << std::endl;
interpolateValues(inputPressures, inputWDirs);
std::cout << "Height:" << std::endl;
interpolateValues(inputPressures, inputHeights);
std::cout << "Interpolation complete." << std::endl;
// Convert the data into ARL format
......@@ -411,7 +434,7 @@ int main(int argc, char **argv)
thisIndexHeader.pole_lon = centreLon;
thisIndexHeader.ref_lat = centreLat;
thisIndexHeader.ref_lon = centreLon;
thisIndexHeader.size = 000.f;
thisIndexHeader.size = 10.f;
thisIndexHeader.orient = 0.f;
thisIndexHeader.tang_lat = centreLat;
thisIndexHeader.sync_xp = (numberGridCells + 1) / 2;
......@@ -427,7 +450,7 @@ int main(int argc, char **argv)
thisIndexHeader.levels = inputPressures;
thisIndexHeader.nVarsAtLevels = std::vector<int>(inputPressures.size(), 5);
thisIndexHeader.surfaceVarNames = {"PRSS", "TEMP", "RELH", "UWND", "VWND"};
thisIndexHeader.varNames = {"PRES", "TEMP", "RELH", "UWND", "VWND"};
thisIndexHeader.varNames = {"HGTS", "TEMP", "RELH", "UWND", "VWND"};
thisIndexHeader.checkSums = std::vector<int>(5, 0);
// Write the headers
......@@ -441,24 +464,36 @@ int main(int argc, char **argv)
std::cout << " Writing meteorological data for timestep " << timestep
<< ", level " << level << "..." << std::endl;
// Write pressure variables
thisRecordHeader.il = level;
// Write pressure/height variables
{
if (level == 0)
{
thisRecordHeader.kvar = "PRSS";
thisRecordHeader.var1 = inputPressures[level];
std::vector<std::vector<float>> thisData(
numberGridCells,
std::vector<float>(numberGridCells, inputPressures[level]));
std::cout << " pressure...";
outputStream.write_record_header(thisRecordHeader);
outputStream.write_record(thisRecordHeader, thisIndexHeader,
thisData);
std::cout << "written" << std::endl;
}
else
{
thisRecordHeader.kvar = "PRES";
thisRecordHeader.kvar = "HGTS";
thisRecordHeader.var1 = inputHeights[level];
std::vector<std::vector<float>> thisData(
numberGridCells,
std::vector<float>(numberGridCells, inputHeights[level]));
std::cout << " height...";
outputStream.write_record_header(thisRecordHeader);
outputStream.write_record(thisRecordHeader, thisIndexHeader,
thisData);
std::cout << "written" << std::endl;
}
thisRecordHeader.var1 = inputPressures[level];
std::vector<std::vector<float>> thisData(
numberGridCells,
std::vector<float>(numberGridCells, inputPressures[level]));
std::cout << " pressure...";
outputStream.write_record_header(thisRecordHeader);
outputStream.write_record(thisRecordHeader, thisIndexHeader, thisData);
std::cout << "written" << std::endl;
}
// Write temperature levels
......
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