Skip to content

GitLab

Commit c3452fcc authored by LEFEBVREJP email's avatar LEFEBVREJP email
Browse files

- Removed radixsnd2arl 

- Added description to CommandLine
parent b46f658b
Hide whitespace changes
Inline Side-by-side
radixcommand/commandline.cc
View file @ c3452fcc
......@@ -136,6 +136,11 @@ void CommandLine::printParsedLine(std::ostream &out) const
void CommandLine::help(std::ostream &out) const
{
// dump commandline example first
// if there is a description, add it up front
if (!mDescription.empty())
{
out << mDescription << std::endl;
}
out << mExecutable;
for (const auto &it : mData)
{
......@@ -287,6 +292,13 @@ bool CommandLine::validate(std::ostream &out) const
return valid;
}
const std::string &CommandLine::description() const { return mDescription; }
void CommandLine::setDescription(const std::string &description)
{
mDescription = description;
}
template <>
std::string CommandLine::get(const std::string &name) const
{
......
radixcommand/commandline.hh
View file @ c3452fcc
......@@ -82,9 +82,13 @@ class RADIX_PUBLIC CommandLine
*/
bool validate(std::ostream &out) const;
const std::string &description() const;
void setDescription(const std::string &description);
private:
int mArgc;
char **mArgv;
std::string mDescription;
std::string mExecutable;
std::vector<std::string> mArgs;
std::vector<std::pair<std::string, std::string>> mDeclArgs;
......
radixio/examples/CMakeLists.txt
View file @ c3452fcc
TRIBITS_ADD_EXECUTABLE(radixsnd2arl
NOEXEPREFIX
SOURCES radixsnd2arl.cc
)
#TRIBITS_ADD_EXECUTABLE(fmsscoreex
# NOEXEPREFIX
# SOURCES fmsscoreex.cc
#)
radixio/examples/radixsnd2arl.cc deleted 100644 → 0
View file @ b46f658b
/*
* Example utility to convert a vertical profile of meteorological
* data to the ARL format
*/
#include <algorithm>
#include <cmath>
#include <ctime>
#include <iostream>
#include <string>
#include <vector>
#include "radixcommand/commandline.hh"
#include "radixcore/stringfunctions.hh"
#include "radixmath/constants.hh"
#include "radixmath/interpolate.hh"
#include "radixmath/util.hh"
#include "radixio/arldatastream.hh"
#include "radixio/csvfile.hh"
using namespace radix;
void addHour(struct tm *time, int hours)
{
int seconds = hours * 60 * 60;
time_t date_seconds = mktime(time) + seconds;
*time = *localtime(&date_seconds);
}
int main(int argc, char **argv)
{
std::cout << "************************" << std::endl;
std::cout << "***** radixsnd2arl *****" << std::endl;
std::cout << "************************" << std::endl;
// Set up command line options
CommandLine commandLine(argc, argv);
commandLine.addOption("i", "Input csv file containing met data", true);
commandLine.addOption("clat", "Centre latitude of output ARL file (degrees)",
true);
commandLine.addOption("clon", "Centre longitude of output ARL file (degrees)",
true);
commandLine.addOption("e", "Extent of output ARL file (km) [500]", false);
commandLine.addOption("r", "Resolution of output ARL file (km) [10]", false);
commandLine.addOption("t", "Time of data start (YYYYMMDDHH) [1951111917]",
false);
commandLine.addOption("n", "Number of one hour timesteps to output [1]",
false);
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 - "
"Defaults to use input pressure extent",
false);
commandLine.addOption("rsend",
"Resampling end pressure level - "
"Defaults to use input pressure extent",
false);
// Ensure required options present
std::vector<std::string> commandErrors;
if (!commandLine.validate(commandErrors))
{
std::cout << "Error in arguments..." << std::endl;
for (std::string error : commandErrors)
{
std::cout << "\t" << error << std::endl;
}
std::cout << std::endl;
commandLine.printParsedLine(std::cout);
return -1;
}
// Get command line options
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 resampleInterval = commandLine.get<float>("rs", -1.f);
float resampleStart = commandLine.get<float>("rsstart", -1.f);
float resampleEnd = commandLine.get<float>("rsend", -1.f);
// Get the grid size
int numberGridCells = (int)(extent / resolution);
if (!numberGridCells % 2 == 1)
{
numberGridCells++;
}
// Parse the start time
int year = from_string(startTime.substr(0, 4), 1951);
int month = from_string(startTime.substr(4, 2), 11);
int day = from_string(startTime.substr(6, 2), 19);
int hour = from_string(startTime.substr(8, 2), 17);
struct tm metTime = {0, 0, 0};
metTime.tm_year = year - 1900;
metTime.tm_mon = month - 1;
metTime.tm_mday = day;
metTime.tm_hour = hour;
std::cout << "Creating ARL-formatted met file with parameters:" << std::endl;
std::cout << " " << extent << " by " << extent << " km grid centred on ("
<< centreLat << "," << centreLon << ") with resolution "
<< resolution << " (" << numberGridCells
<< " cells in each direction)" << std::endl;
std::cout << " " << numberTimesteps
<< " timesteps (of 1 hour each) starting from " << metTime.tm_year
<< "-" << metTime.tm_mon + 1 << "-" << metTime.tm_mday << " at "
<< metTime.tm_hour << "00" << std::endl
<< std::endl;
// If we have correct options, parse input file
std::cout << "Reading input csv file: " << inputCsvPath << std::endl;
// Header strings to denote fields
std::string pressureString = "pres", tempString = "temp", relHumString = "rh",
dewPtString = "tdew", wDirString = "wdir", wSpdString = "wspd",
heightString = "zhgt";
int pressureIndex = -1, tempIndex = -1, relHumIndex = -1, dewPtIndex = -1,
wDirIndex = -1, wSpdIndex = -1, heightIndex = -1;
bool usingRelHum = false, usingDewPt = false;
// Open and read the csv
CSVFile inputCsv(inputCsvPath);
std::vector<std::vector<std::string>> inputData;
bool inputReadSuccess = inputCsv.data(inputData);
if (!inputReadSuccess)
{
std::cout << "Failed to read input csv!" << std::endl;
std::cout
<< "Please ensure the file exists and has the appropriate permissions."
<< std::endl;
return -2;
}
// Ensure there is data in here
if (inputData.size() < 5)
{
std::cout << "Not enough data in input csv!" << std::endl;
std::cout
<< "There are only " << inputData.size()
<< " lines in this file - at least 5 are required (4 headers + data)"
<< std::endl;
return -3;
}
// Work out which fields are which data points - should be 2nd line (1)
std::cout << "Finding data fields..." << std::endl;
for (size_t entry = 0; entry < inputData[1].size(); ++entry)
{
if (inputData[1][entry] == pressureString)
{
std::cout << " Found pressure at index " << entry << std::endl;
pressureIndex = entry;
}
else if (inputData[1][entry] == tempString)
{
std::cout << " Found temperature at index " << entry << std::endl;
tempIndex = entry;
}
else if (inputData[1][entry] == relHumString)
{
std::cout << " Found relative humidity at index " << entry << std::endl;
relHumIndex = entry;
usingRelHum = true;
}
else if (inputData[1][entry] == dewPtString)
{
std::cout << " Found dew point at index " << entry << std::endl;
dewPtIndex = entry;
usingDewPt = true;
}
else if (inputData[1][entry] == wDirString)
{
std::cout << " Found wind direction at index " << entry << std::endl;
wDirIndex = entry;
}
else if (inputData[1][entry] == wSpdString)
{
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 && dewPtIndex == -1) || (wDirIndex == -1) ||
(wSpdIndex == -1))
{
std::cout << "Missing data fields in input!" << std::endl;
if (pressureIndex == -1)
{
std::cout << " Missing pressure field (denoted by " << pressureString
<< ")";
}
if (tempIndex == -1)
{
std::cout << " Missing temperature field (denoted by " << tempString
<< ")";
}
if (relHumIndex == -1 && dewPtIndex == -1)
{
std::cout << " Missing both relative humidity field (denoted by "
<< relHumString << ") and dew point field (denoted by "
<< dewPtString << ") - need at least one of these";
}
if (wDirIndex == -1)
{
std::cout << " Missing wind direction field (denoted by " << wDirString
<< ")";
}
if (wSpdIndex == -1)
{
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;
}
// Ensure we aren't using both relative humidity and dew point (duplicate
// data) Prefer use of relative humidity
if (usingRelHum)
{
usingDewPt = false;
}
// Read the data
std::vector<float> inputPressures, inputTemps, inputRelHums, inputDewPts,
inputWDirs, inputWSpds, inputHeights;
std::cout << "Data fields found - reading data..." << std::endl;
float missingValue = -9999.f;
for (size_t row = 4; row < inputData.size(); ++row)
{
if (inputData[row].size() < inputData[1].size())
{
std::cout << " Warning: this row (" << row + 1
<< ") has less entries than the field names row" << std::endl;
;
}
bool foundPressure = false, foundTemp = false, foundRelHum = false,
foundDewPt = false, foundWDir = false, foundWSpd = false,
foundHeight = false;
for (size_t entry = 0; entry < inputData[row].size(); ++entry)
{
float thisValue = from_string(inputData[row][entry], missingValue);
if (entry == pressureIndex)
{
foundPressure = true;
inputPressures.push_back(thisValue);
}
else if (entry == tempIndex)
{
foundTemp = true;
inputTemps.push_back(thisValue);
}
else if (usingRelHum && entry == relHumIndex)
{
foundRelHum = true;
inputRelHums.push_back(thisValue);
}
else if (usingDewPt && entry == dewPtIndex)
{
foundDewPt = true;
inputDewPts.push_back(thisValue);
}
else if (entry == wSpdIndex)
{
foundWSpd = true;
inputWSpds.push_back(thisValue);
}
else if (entry == wDirIndex)
{
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)
{
std::cout << " Warning: couldn't find pressure in row " << row
<< std::endl;
inputPressures.push_back(missingValue);
}
if (!foundTemp)
{
std::cout << " Warning: couldn't find temperature in row " << row
<< std::endl;
inputTemps.push_back(missingValue);
}
if (usingRelHum && !foundRelHum)
{
std::cout << " Warning: couldn't find relative humidity in row " << row
<< std::endl;
inputRelHums.push_back(missingValue);
}
if (usingDewPt && !foundDewPt)
{
std::cout << " Warning: couldn't find dew point in row " << row
<< std::endl;
inputDewPts.push_back(missingValue);
}
if (!foundWSpd)
{
std::cout << " Warning: couldn't find wind speed in row " << row
<< std::endl;
inputWSpds.push_back(missingValue);
}
if (!foundWDir)
{
std::cout << " Warning: couldn't find wind direction in row " << row
<< std::endl;
inputWDirs.push_back(missingValue);
}
if (!foundHeight)
{
std::cout << " Warning: couldn't find height in row " << row
<< std::endl;
inputHeights.push_back(missingValue);
}
// If pressure value in this row is -9999, remove the row
// - pressure entry is required for interpolation & by HYSPLIT
if (inputPressures.back() == missingValue)
{
std::cout << " Warning: pressure entry in row " << row
<< " is -9999; removing row" << std::endl;
inputPressures.pop_back();
inputTemps.pop_back();
if (usingRelHum)
{
inputRelHums.pop_back();
}
if (usingDewPt)
{
inputDewPts.pop_back();
}
inputWSpds.pop_back();
inputWDirs.pop_back();
inputHeights.pop_back();
}
}
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;
// Calculate the upper and lower bounds of the resampled pressure array
radix_line(" input start = " << resampleStart
<< "; end = " << resampleEnd);
if (resampleStart < 0)
{
resampleStart =
(float)roundUpInt((int)*std::max_element(std::begin(inputPressures),
std::end(inputPressures)),
resampleInterval);
}
if (resampleEnd < 0)
{
resampleEnd =
(float)roundDownInt((int)*std::min_element(std::begin(inputPressures),
std::end(inputPressures)),
resampleInterval);
}
std::cout << " Using start = " << resampleStart
<< "; end = " << resampleEnd << std::endl;
// Calculate the resampled pressure values using the interval given
std::cout << " Calculating resampled pressure values..." << std::endl;
float thisPressure;
for (thisPressure = resampleStart; thisPressure > resampleEnd;
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, true);
if (usingRelHum)
{
std::cout << " Relative humidity:" << std::endl;
inputRelHums = interpolateToOtherBaseValues(
inputPressures, resampledPressures, inputRelHums, true);
}
if (usingDewPt)
{
std::cout << " Dew point:" << std::endl;
inputDewPts = interpolateToOtherBaseValues(
inputPressures, resampledPressures, inputDewPts, true);
}
std::cout << " Wind speed:" << std::endl;
inputWSpds = interpolateToOtherBaseValues(
inputPressures, resampledPressures, inputWSpds, true);
std::cout << " Wind direction:" << std::endl;
inputWDirs = interpolateToOtherBaseValues(
inputPressures, resampledPressures, inputWDirs, true, true);
std::cout << " Height:" << std::endl;
inputHeights = interpolateToOtherBaseValues(
inputPressures, resampledPressures, inputHeights, true);
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, true);
if (usingRelHum)
{
std::cout << " Relative humidity:" << std::endl;
inputRelHums = interpolateValues(inputPressures, inputRelHums, true);
}
if (usingDewPt)
{
std::cout << " Dew point:" << std::endl;
inputDewPts = interpolateValues(inputPressures, inputDewPts, true);
}
std::cout << " Wind speed:" << std::endl;
inputWSpds = interpolateValues(inputPressures, inputWSpds, false);
std::cout << " Wind direction:" << std::endl;
inputWDirs = interpolateValues(inputPressures, inputWDirs, false, true);
std::cout << " Height:" << std::endl;
inputHeights = interpolateValues(inputPressures, inputHeights, true);
}
// First convert temperature from Celsius to Kelvin
std::cout << "Converting Temperature from C to Kelvin." << std::endl;
std::transform(inputTemps.begin(), inputTemps.end(), inputTemps.begin(),
[](double v) { return v - ABS_ZERO_CELSIUS; });
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;
int indexSize = indexBase + (inputPressures.size() * levelBase);
int recordSize = numberGridCells * numberGridCells;
if (indexSize > recordSize)
{
// The index header will be too big - instruct user to either reduce number
// of vertical levels or increase number of grid cells
std::cout << "Error: index header will be too large!" << std::endl
<< " The size of the index header in the ARL-formatted output ("
<< indexBase << " + number of vertical levels * " << levelBase
<< " = " << indexSize
<< ") will be larger than the size of each record (number of "
"grid cells ^ 2 = "
<< recordSize << ")." << std::endl;
std::cout << " Please either increase the number of horizontal grid cells "
"(increase extent/resolution) or reduce the number of input "
"vertical levels."
<< std::endl;
throw std::exception();
}
// Convert the data into ARL format
std::cout << "Converting data to ARL format..." << std::endl;
ARLDataStream outputStream(outputArlPath, std::ios::out);
for (int timestep = 0; timestep < numberTimesteps; ++timestep)
{
// Write index header section
std::cout << " Writing index headers for timestep " << timestep << "..."
<< std::endl;
ARLRecordHeader thisRecordHeader;
thisRecordHeader.year = metTime.tm_year;
thisRecordHeader.month = metTime.tm_mon + 1;
thisRecordHeader.day = metTime.tm_mday;
thisRecordHeader.hour = metTime.tm_hour;
thisRecordHeader.ic = 0;
thisRecordHeader.il = 0;
thisRecordHeader.cgrid = "99";
thisRecordHeader.kvar = "INDX";
thisRecordHeader.nexp = 0;
thisRecordHeader.prec = 0.f;
thisRecordHeader.var1 = 0.f;
ARLIndexHeader thisIndexHeader;
thisIndexHeader.model_id = "NFDB";
thisIndexHeader.icx = 0;
thisIndexHeader.mn = 0;
thisIndexHeader.pole_lat = centreLat;
thisIndexHeader.pole_lon = centreLon;
thisIndexHeader.ref_lat = centreLat;
thisIndexHeader.ref_lon = centreLon;
thisIndexHeader.size = resolution;
thisIndexHeader.orient = 0.f;
thisIndexHeader.tang_lat = centreLat;
thisIndexHeader.sync_xp = (numberGridCells + 1) / 2;
thisIndexHeader.sync_yp = (numberGridCells + 1) / 2;
thisIndexHeader.sync_lat = centreLat;
thisIndexHeader.sync_lon = centreLon;
thisIndexHeader.dummy = 0.f;
thisIndexHeader.nx = numberGridCells;
thisIndexHeader.ny = numberGridCells;
thisIndexHeader.nz = int(inputPressures.size());
thisIndexHeader.z_flag = 2;
thisIndexHeader.lenh = numberGridCells * numberGridCells;
thisIndexHeader.levels = inputPressures;
thisIndexHeader.num_vars_at_levels =