Skip to content

GitLab

Commit 4b272183 authored by Purves, Murray's avatar Purves, Murray
Browse files

WIP Adding resampling to radixsnd2arl.cc

parent 8cc3297c
Pipeline #28264 failed with stages
in 2 minutes and 42 seconds
Hide whitespace changes
Inline Side-by-side
radixio/examples/radixsnd2arl.cc
View file @ 4b272183
......@@ -52,6 +52,11 @@ int main(int argc, char **argv)
commandLine.addOption("g", "Add ground level elevation to height values",
false);
commandLine.addOption("o", "Output ARL file", false);
commandLine.addOption(
"rs", "Resample input data to uniform pressure levels (value=interval)",
false);
commandLine.addOption("rsstart", "Resampling start pressure level[950]",
false);
// Ensure required options present
std::vector<std::string> commandErrors;
......@@ -72,13 +77,15 @@ int main(int argc, char **argv)
std::string inputCsvPath = commandLine.get<std::string>("i");
std::string outputArlPath =
commandLine.get<std::string>("o", inputCsvPath + ".bin");
float extent = commandLine.get<float>("e", 500.0);
float resolution = commandLine.get<float>("r", 10.0);
float centreLat = commandLine.get<float>("clat");
float centreLon = commandLine.get<float>("clon");
int numberTimesteps = commandLine.get<int>("n", 1);
float groundElevation = commandLine.get<float>("g", 0.f);
std::string startTime = commandLine.get<std::string>("t", "1951111917");
float extent = commandLine.get<float>("e", 500.0);
float resolution = commandLine.get<float>("r", 10.0);
float centreLat = commandLine.get<float>("clat");
float centreLon = commandLine.get<float>("clon");
int numberTimesteps = commandLine.get<int>("n", 1);
float groundElevation = commandLine.get<float>("g", 0.f);
std::string startTime = commandLine.get<std::string>("t", "1951111917");
float resampleInterval = commandLine.get<float>("rs", -1.f);
float resampleStart = commandLine.get<float>("rsstart", 50.f);
// Get the grid size
int numberGridCells = (int)(extent / resolution);
if (!numberGridCells % 2 == 1)
......@@ -356,6 +363,87 @@ int main(int argc, char **argv)
std::cout << "Data read complete: " << inputPressures.size()
<< " entries read." << std::endl;
// We will need the ground level pressure later
// Assume it is the first entry in the input
float groundPressure = inputPressures[0];
// Interpolate the data to remove missing values
std::cout << "Interpolating values to remove missing data..." << std::endl;
// Resample the pressures if required
if (resampleInterval > 0)
{
std::cout << " Resampling input data using sampling interval: "
<< resampleInterval << " mb" << std::endl;
std::vector<float> resampledPressures;
float minPressure =
*std::min_element(std::begin(inputPressures), std::end(inputPressures));
// Calculate the resampled pressure values using the interval given
std::cout << " Calculating resampled pressure values..." << std::endl;
for (float thisPressure = resampleStart; thisPressure > minPressure;
thisPressure -= resampleInterval)
{
resampledPressures.push_back(thisPressure);
}
std::cout << " Temperature:" << std::endl;
inputTemps = interpolateToOtherBaseValues(inputPressures,
resampledPressures, inputTemps);
if (usingRelHum)
{
std::cout << " Relative humidity:" << std::endl;
inputRelHums = interpolateToOtherBaseValues(
inputPressures, resampledPressures, inputRelHums);
}
if (usingDewPt)
{
std::cout << " Dew point:" << std::endl;
inputDewPts = interpolateToOtherBaseValues(
inputPressures, resampledPressures, inputDewPts);
}
std::cout << " Wind speed:" << std::endl;
inputWSpds = interpolateToOtherBaseValues(inputPressures,
resampledPressures, inputWSpds);
std::cout << " Wind direction:" << std::endl;
inputWDirs = interpolateToOtherBaseValues(
inputPressures, resampledPressures, inputWDirs, true);
std::cout << " Height:" << std::endl;
inputHeights = interpolateToOtherBaseValues(
inputPressures, resampledPressures, inputHeights);
std::cout << " Replacing pressure data with resampled array" << std::endl;
radix_block(
radix(" Initial pressure array: ");
for (float p
: inputPressures) { radix(" " << p); } radix_line("");
radix(" Resampled pressure array: ");
for (float p
: resampledPressures) { radix(" " << p); } radix_line(""););
inputPressures = resampledPressures;
}
else
{
std::cout << " Not resampling data" << std::endl;
std::cout << " Temperature:" << std::endl;
inputTemps = interpolateValues(inputPressures, inputTemps);
if (usingRelHum)
{
std::cout << " Relative humidity:" << std::endl;
inputRelHums = interpolateValues(inputPressures, inputRelHums);
}
if (usingDewPt)
{
std::cout << " Dew point:" << std::endl;
inputDewPts = interpolateValues(inputPressures, inputDewPts);
}
std::cout << " Wind speed:" << std::endl;
inputWSpds = interpolateValues(inputPressures, inputWSpds);
std::cout << " Wind direction:" << std::endl;
inputWDirs = interpolateValues(inputPressures, inputWDirs, true);
std::cout << " Height:" << std::endl;
inputHeights = interpolateValues(inputPressures, inputHeights);
}
std::cout << " Interpolation complete." << std::endl;
// Check if the grid is of appropriate size - with too small grid/too many
// vertical levels, the index header can become oversized leading to errors
int indexBase = 108, levelBase = 48;
......@@ -379,28 +467,6 @@ int main(int argc, char **argv)
throw std::exception();
}
// Interpolate the data to remove missing values
std::cout << "Interpolating values to remove missing data..." << std::endl;
std::cout << "Temperature:" << std::endl;
inputTemps = interpolateValues(inputPressures, inputTemps);
if (usingRelHum)
{
std::cout << "Relative humidity:" << std::endl;
inputRelHums = interpolateValues(inputPressures, inputRelHums);
}
if (usingDewPt)
{
std::cout << "Dew point:" << std::endl;
inputDewPts = interpolateValues(inputPressures, inputDewPts);
}
std::cout << "Wind speed:" << std::endl;
inputWSpds = interpolateValues(inputPressures, inputWSpds);
std::cout << "Wind direction:" << std::endl;
inputWDirs = interpolateValues(inputPressures, inputWDirs, true);
std::cout << "Height:" << std::endl;
inputHeights = interpolateValues(inputPressures, inputHeights);
std::cout << "Interpolation complete." << std::endl;
// Convert the data into ARL format
std::cout << "Converting data to ARL format..." << std::endl;
ARLDataStream outputStream(outputArlPath, std::ios::out);
......@@ -479,10 +545,10 @@ int main(int argc, char **argv)
if (level == 0)
{
thisRecordHeader.kvar = "PRSS";
float groundPressure = inputPressures[level];
thisRecordHeader.var1 = groundPressure;
std::vector<std::vector<float>> thisData(
numberGridCells, std::vector<float>(numberGridCells, thisHeight));
numberGridCells,
std::vector<float>(numberGridCells, groundPressure));
std::cout << " pressure...";
outputStream.write_record_header(thisRecordHeader);
outputStream.write_record(thisRecordHeader, thisIndexHeader,
......@@ -492,10 +558,10 @@ int main(int argc, char **argv)
else
{
thisRecordHeader.kvar = "HGTS";
thisRecordHeader.var1 = inputHeights[level] + groundElevation;
float thisHeight = inputHeights[level] + groundElevation;
thisRecordHeader.var1 = thisHeight;
std::vector<std::vector<float>> thisData(
numberGridCells,
std::vector<float>(numberGridCells, inputHeights[level]));
numberGridCells, std::vector<float>(numberGridCells, thisHeight));
std::cout << " height...";
outputStream.write_record_header(thisRecordHeader);
outputStream.write_record(thisRecordHeader, thisIndexHeader,
......
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