Commit 194d426e authored by Purves, Murray's avatar Purves, Murray
Browse files

Merge branch 'arldatastream' into python-bindings

parents 0c2ba96d 943e3531
Pipeline #28290 failed with stages
in 58 seconds
......@@ -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,94 @@ 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;
float thisPressure;
for (thisPressure = resampleStart; thisPressure > minPressure;
thisPressure -= resampleInterval)
{
resampledPressures.push_back(thisPressure);
}
// Add another entry below the min. pressure
// - but only if that pressure won't be 0
if (thisPressure > 0)
{
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 +474,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);
......@@ -451,7 +524,6 @@ int main(int argc, char **argv)
"VWND"};
thisIndexHeader.var_names.push_back(surfaceVarNames);
for (size_t level = 1; level < inputPressures.size(); ++level)
{
thisIndexHeader.var_names.push_back(varNames);
......@@ -480,11 +552,10 @@ int main(int argc, char **argv)
if (level == 0)
{
thisRecordHeader.kvar = "PRSS";
// Add ground elevation to height (default is 0)
float thisHeight = inputHeights[level] + groundElevation;
thisRecordHeader.var1 = thisHeight;
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,
......@@ -494,10 +565,10 @@ int main(int argc, char **argv)
else
{
thisRecordHeader.kvar = "HGTS";
thisRecordHeader.var1 = inputHeights[level];
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,
......
......@@ -164,21 +164,41 @@ std::vector<T> interpolateToOtherBaseValues(
}
radix_line(" Initial interpolation complete");
// Work out if the base array goes low to high or high to low
bool lowToHigh = true;
if (baseValues.front() < baseValues.back())
{
radix_line(" Base values go low to high");
}
else
{
lowToHigh = false;
radix_line(" Base values go high to low");
}
// Iterate over the new base vector to calculate values at these positions
radix_line(" Interpolating to new base values");
for (size_t i = 0; i < newBaseValues.size(); ++i)
{
// Fill in values with a new base lower than the lowest original
if (newBaseValues[i] < baseValues.front())
if ((lowToHigh && newBaseValues[i] < baseValues.front()) ||
(!lowToHigh && newBaseValues[i] < baseValues.back()))
{
radix_line(" Base value below range: filling with lowest");
finalInterpolatedValues[i] = initialInterpolatedValues.front();
radix_line(" Base value "
<< newBaseValues[i] << " below range " << baseValues.front()
<< "-" << baseValues.back() << ": filling with lowest");
finalInterpolatedValues[i] = std::min(initialInterpolatedValues.front(),
initialInterpolatedValues.back());
}
// Fill in values with a new base higher than the highest original
else if (newBaseValues[i] > baseValues.back())
else if ((lowToHigh && newBaseValues[i] > baseValues.back()) ||
(!lowToHigh && newBaseValues[i] > baseValues.front()))
{
radix_line(" Base value above range: filling with highest");
finalInterpolatedValues[i] = initialInterpolatedValues.back();
radix_line(" Base value "
<< newBaseValues[i] << " above range " << baseValues.front()
<< "-" << baseValues.back() << ": filling with highest");
finalInterpolatedValues[i] = std::max(initialInterpolatedValues.front(),
initialInterpolatedValues.back());
}
// Fill in central values
else
......@@ -188,7 +208,8 @@ std::vector<T> interpolateToOtherBaseValues(
size_t lowerBaseIndex = 0, higherBaseIndex = 0;
for (size_t j = 0; j < baseValues.size(); ++j)
{
if (baseValues[j] > newBaseValues[i])
if ((lowToHigh && (baseValues[j] >= newBaseValues[i])) ||
(!lowToHigh && (baseValues[j] <= newBaseValues[i])))
{
higherBaseIndex = j;
break;
......@@ -218,19 +239,44 @@ std::vector<T> interpolateToOtherBaseValues(
}
}
finalInterpolatedValues[i] =
lastGoodValue +
((nextGoodValue - lastGoodValue) *
((newBaseValues[i] - baseValues[lowerBaseIndex]) /
(baseValues[higherBaseIndex] - baseValues[lowerBaseIndex])));
// If the two values are the same, we are on an existing data point
// - don't interpolate
if (higherBaseIndex == lowerBaseIndex)
{
finalInterpolatedValues[i] = lastGoodValue;
}
else
{
finalInterpolatedValues[i] =
lastGoodValue +
((nextGoodValue - lastGoodValue) *
((newBaseValues[i] - baseValues[lowerBaseIndex]) /
(baseValues[higherBaseIndex] - baseValues[lowerBaseIndex])));
}
if (circular)
{
finalInterpolatedValues[i] = fmod(finalInterpolatedValues[i], 360.0);
}
// radix_line("interp: previous " << lastGoodValue << ", next "
// << nextGoodValue << ", result "
// << finalInterpolatedValues[i]);
}
}
radix_line("Interpolation to new base values complete");
radix_block(
radix(" Initial base values: ");
for (float f
: baseValues) { radix(" " << f); } radix_line("");
radix(" Initial interpolated values: ");
for (float f
: initialInterpolatedValues) { radix(" " << f); } radix_line("");
radix(" New base values: ");
for (float f
: newBaseValues) { radix(" " << f); } radix_line("");
radix(" Final interpolated values: ");
for (float f
: finalInterpolatedValues) { radix(" " << f); } radix_line(""););
return finalInterpolatedValues;
}
......
......@@ -150,20 +150,18 @@ TEST(Interpolate, NewBaseValues)
expectValues1[i] * tolerance);
}
// Uncomment if negative gradient base values required
// std::vector<float> baseValues2{3000.f, 2000.f, 1000.f};
// std::vector<float> newBaseValues2{3500.f, 2500.f, 1500.f, 500.f};
// std::vector<float> interpValues2{30.5, 20.5, 10.5};
// std::vector<float> expectValues2{30.5, 25.5, 15.5, 10.5};
// std::vector<float> testValues2 =
// interpolateToOtherBaseValues(baseValues2, newBaseValues2,
// interpValues2);
// ASSERT_EQ(expectValues2.size(), testValues2.size());
// for (int i = 0; i < testValues2.size(); ++i)
// {
// EXPECT_NEAR(expectValues2[i], testValues2[i], expectValues2[i] *
// tolerance);
// }
// High-to-low base values
std::vector<float> baseValues2{3000.f, 2000.f, 1000.f};
std::vector<float> newBaseValues2{3500.f, 2500.f, 1500.f, 500.f};
std::vector<float> interpValues2{30.5, 20.5, 10.5};
std::vector<float> expectValues2{30.5, 25.5, 15.5, 10.5};
std::vector<float> testValues2 =
interpolateToOtherBaseValues(baseValues2, newBaseValues2, interpValues2);
ASSERT_EQ(expectValues2.size(), testValues2.size());
for (int i = 0; i < testValues2.size(); ++i)
{
EXPECT_NEAR(expectValues2[i], testValues2[i], expectValues2[i] * tolerance);
}
}
TEST(Interpolate, NewBaseValuesDoubles)
......
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