Commit 190e9909 authored by LEFEBVREJP email's avatar LEFEBVREJP email
Browse files

Adding hectorpascal conversion utility. Updating tstUtil after earth radius...

Adding hectorpascal conversion utility. Updating tstUtil after earth radius update. Updating gfsfile to convert pressure levels to altitude m-agl.
parent 52d0a190
...@@ -230,8 +230,8 @@ GFSFile::GFSFile(std::string file) ...@@ -230,8 +230,8 @@ GFSFile::GFSFile(std::string file)
<< " " << mLabel.hour; << " " << mLabel.hour;
mEndTime = ss.str(); mEndTime = ss.str();
} }
// radix_tagged_line("Profile time: " << mProfiles[mProfiles.size()-1] // radix_tagged_line("Profile time: " << mProfiles[mProfiles.size()-1]
// << mLabel.toString()); // << mLabel.toString());
} }
rstream->close(); rstream->close();
delete rstream; delete rstream;
...@@ -708,6 +708,7 @@ std::vector<std::vector<float>> GFSFile::query(float lat ...@@ -708,6 +708,7 @@ std::vector<std::vector<float>> GFSFile::query(float lat
bool sfcwnd = false; bool sfcwnd = false;
float offset = 0.0f; float offset = 0.0f;
float plevel = 0.0f; float plevel = 0.0f;
float surfaceAltitude = 0.f;
std::vector<std::vector<float>> results(mHeader.nz); std::vector<std::vector<float>> results(mHeader.nz);
for(int ll = 0; ll < mHeader.nz; ++ll) for(int ll = 0; ll < mHeader.nz; ++ll)
...@@ -747,7 +748,6 @@ std::vector<std::vector<float>> GFSFile::query(float lat ...@@ -747,7 +748,6 @@ std::vector<std::vector<float>> GFSFile::query(float lat
} }
// by default assume level = pressure unless PRES variable appears // by default assume level = pressure unless PRES variable appears
// (i.e. terrain data (type=3) will have local pressure variable // (i.e. terrain data (type=3) will have local pressure variable
int level = static_cast<int>(plevel);
// match variables defined in file's index record with those variables // match variables defined in file's index record with those variables
// that have been defined in this subroutine and create a variable number // that have been defined in this subroutine and create a variable number
...@@ -782,6 +782,13 @@ std::vector<std::vector<float>> GFSFile::query(float lat ...@@ -782,6 +782,13 @@ std::vector<std::vector<float>> GFSFile::query(float lat
} }
// initialize space for results vector // initialize space for results vector
results[ll] = std::vector<float>(columns.size(), 0.0f); results[ll] = std::vector<float>(columns.size(), 0.0f);
// if "HGTS" has been requested add level as the default HGTS
auto hIt = std::find(columns.begin(), columns.end(), "HGTS");
if(hIt != columns.end())
{
results[ll][hIt-columns.begin()] = hpaToAltitude(plevel) - surfaceAltitude;
}
// check for time // check for time
auto timeIt = std::find(columns.begin(), columns.end(), "TIME"); auto timeIt = std::find(columns.begin(), columns.end(), "TIME");
if(timeIt != columns.end()) if(timeIt != columns.end())
...@@ -793,25 +800,29 @@ std::vector<std::vector<float>> GFSFile::query(float lat ...@@ -793,25 +800,29 @@ std::vector<std::vector<float>> GFSFile::query(float lat
auto presIt = std::find(columns.begin(), columns.end(), "PRSS"); auto presIt = std::find(columns.begin(), columns.end(), "PRSS");
if(presIt != columns.end()) if(presIt != columns.end())
{ {
results[ll][presIt-columns.begin()] = (float)level; results[ll][presIt-columns.begin()] = plevel;
} }
for(int kk = 0; kk < nvar; ++kk) for(int kk = 0; kk < nvar; ++kk)
{ {
if(mVarbId[ll][kk].compare("PRES") == 0) plevel = vdata[kk][ll]; std::string varb = mVarbId[ll][kk];
if(mVarbId[ll][kk].compare("THET") == 0) if(varb.compare("PRES") == 0) plevel = vdata[kk][ll];
if(varb.compare("THET") == 0)
{ {
tpot = vdata[kk][ll]; tpot = vdata[kk][ll];
// potential temperature defined; replace with ambient // potential temperature defined; replace with ambient
vdata[kk][ll] = (tpot*std::pow(plevel/1000.0f,0.286f))-273.16f; vdata[kk][ll] = (tpot*std::pow(plevel/1000.0f,0.286f))-273.16f;
} }
if(mVarbId[ll][kk].compare("TEMP") == 0) temp = vdata[kk][ll]+273.16f; if(varb.compare("TEMP") == 0) temp = vdata[kk][ll]+273.16f;
std::string varb = mVarbId[ll][kk];
//map certain varb //map certain surface temperature(T02M) or surface relative humidity(RH2M)
if(varb.compare("T02M") == 0) varb = "TEMP"; if(varb.compare("T02M") == 0) varb = "TEMP";
if(varb.compare("RH2M") == 0) varb = "RELH"; if(varb.compare("RH2M") == 0) varb = "RELH";
if(varb.compare("PRSS") == 0)
{
surfaceAltitude = hpaToAltitude(vdata[kk][ll]) - 2.f;
}
auto it = std::find(columns.begin(), columns.end(), varb); auto it = std::find(columns.begin(), columns.end(), varb);
if(it != columns.end()) if(it != columns.end())
{ {
......
...@@ -4,6 +4,15 @@ ...@@ -4,6 +4,15 @@
#include "radixmath/util.hh" #include "radixmath/util.hh"
using namespace radix; using namespace radix;
TEST(radix, hpaToAltitude)
{
Real altitude = hpaToAltitude(835);
EXPECT_FLOAT_EQ(1601.6289, altitude);
altitude = hpaToAltitude(67.7);
EXPECT_FLOAT_EQ(17830.553, altitude);
}
TEST(radix, cspanf) TEST(radix, cspanf)
{ {
Real lon = -190; Real lon = -190;
...@@ -62,7 +71,7 @@ TEST(radix, greatCircleDistance){ ...@@ -62,7 +71,7 @@ TEST(radix, greatCircleDistance){
, lat1 , lat1
, lon2 , lon2
, EARTH_RADIUS_MEAN); , EARTH_RADIUS_MEAN);
EXPECT_FLOAT_EQ(4430.4688, distance); EXPECT_FLOAT_EQ(4430.6079, distance);
// north // north
distance = greatCircleDistance(lat2 distance = greatCircleDistance(lat2
...@@ -70,7 +79,7 @@ TEST(radix, greatCircleDistance){ ...@@ -70,7 +79,7 @@ TEST(radix, greatCircleDistance){
, lat2 , lat2
, lon2 , lon2
, EARTH_RADIUS_MEAN); , EARTH_RADIUS_MEAN);
EXPECT_FLOAT_EQ(4427.5361, distance); EXPECT_FLOAT_EQ(4427.6748, distance);
// west // west
double west = greatCircleDistance(lat1 double west = greatCircleDistance(lat1
...@@ -78,7 +87,7 @@ TEST(radix, greatCircleDistance){ ...@@ -78,7 +87,7 @@ TEST(radix, greatCircleDistance){
, lat2 , lat2
, lon1 , lon1
, EARTH_RADIUS_MEAN); , EARTH_RADIUS_MEAN);
EXPECT_FLOAT_EQ(5559.7461, west); EXPECT_FLOAT_EQ(5559.9209, west);
// east (should equal west) // east (should equal west)
distance = greatCircleDistance(lat1 distance = greatCircleDistance(lat1
...@@ -102,7 +111,7 @@ TEST(radix, greatCircleDistance){ ...@@ -102,7 +111,7 @@ TEST(radix, greatCircleDistance){
, lat1 , lat1
, lon2 , lon2
, EARTH_RADIUS_MEAN); , EARTH_RADIUS_MEAN);
EXPECT_FLOAT_EQ(4430.4688, distance); EXPECT_FLOAT_EQ(4430.6079, distance);
// north // north
distance = greatCircleDistance(lat2 distance = greatCircleDistance(lat2
...@@ -110,7 +119,7 @@ TEST(radix, greatCircleDistance){ ...@@ -110,7 +119,7 @@ TEST(radix, greatCircleDistance){
, lat2 , lat2
, lon2 , lon2
, EARTH_RADIUS_MEAN); , EARTH_RADIUS_MEAN);
EXPECT_FLOAT_EQ(4427.5361, distance); EXPECT_FLOAT_EQ(4427.6748, distance);
// west // west
double west = greatCircleDistance(lat1 double west = greatCircleDistance(lat1
...@@ -118,7 +127,7 @@ TEST(radix, greatCircleDistance){ ...@@ -118,7 +127,7 @@ TEST(radix, greatCircleDistance){
, lat2 , lat2
, lon1 , lon1
, EARTH_RADIUS_MEAN); , EARTH_RADIUS_MEAN);
EXPECT_FLOAT_EQ(5559.7461, west); EXPECT_FLOAT_EQ(5559.9209, west);
// east (should equal west) // east (should equal west)
distance = greatCircleDistance(lat1 distance = greatCircleDistance(lat1
...@@ -142,7 +151,7 @@ TEST(radix, greatCircleDistance){ ...@@ -142,7 +151,7 @@ TEST(radix, greatCircleDistance){
, lat1 , lat1
, lon2 , lon2
, EARTH_RADIUS_MEAN); , EARTH_RADIUS_MEAN);
EXPECT_FLOAT_EQ(4430.4688, distance); EXPECT_FLOAT_EQ(4430.6079, distance);
// north // north
distance = greatCircleDistance(lat2 distance = greatCircleDistance(lat2
...@@ -150,7 +159,7 @@ TEST(radix, greatCircleDistance){ ...@@ -150,7 +159,7 @@ TEST(radix, greatCircleDistance){
, lat2 , lat2
, lon2 , lon2
, EARTH_RADIUS_MEAN); , EARTH_RADIUS_MEAN);
EXPECT_FLOAT_EQ(4427.5361, distance); EXPECT_FLOAT_EQ(4427.6748, distance);
// west // west
double west = greatCircleDistance(lat1 double west = greatCircleDistance(lat1
...@@ -158,7 +167,7 @@ TEST(radix, greatCircleDistance){ ...@@ -158,7 +167,7 @@ TEST(radix, greatCircleDistance){
, lat2 , lat2
, lon1 , lon1
, EARTH_RADIUS_MEAN); , EARTH_RADIUS_MEAN);
EXPECT_FLOAT_EQ(5559.7461, west); EXPECT_FLOAT_EQ(5559.9209, west);
// east (should equal west) // east (should equal west)
distance = greatCircleDistance(lat1 distance = greatCircleDistance(lat1
...@@ -182,7 +191,7 @@ TEST(radix, greatCircleDistance){ ...@@ -182,7 +191,7 @@ TEST(radix, greatCircleDistance){
, lat1 , lat1
, lon2 , lon2
, EARTH_RADIUS_MEAN); , EARTH_RADIUS_MEAN);
EXPECT_FLOAT_EQ(4430.4688, distance); EXPECT_FLOAT_EQ(4430.6079, distance);
// north // north
distance = greatCircleDistance(lat2 distance = greatCircleDistance(lat2
...@@ -190,7 +199,7 @@ TEST(radix, greatCircleDistance){ ...@@ -190,7 +199,7 @@ TEST(radix, greatCircleDistance){
, lat2 , lat2
, lon2 , lon2
, EARTH_RADIUS_MEAN); , EARTH_RADIUS_MEAN);
EXPECT_FLOAT_EQ(4427.5361, distance); EXPECT_FLOAT_EQ(4427.6748, distance);
// west // west
double west = greatCircleDistance(lat1 double west = greatCircleDistance(lat1
...@@ -198,7 +207,7 @@ TEST(radix, greatCircleDistance){ ...@@ -198,7 +207,7 @@ TEST(radix, greatCircleDistance){
, lat2 , lat2
, lon1 , lon1
, EARTH_RADIUS_MEAN); , EARTH_RADIUS_MEAN);
EXPECT_FLOAT_EQ(5559.7461, west); EXPECT_FLOAT_EQ(5559.9209, west);
// east (should equal west) // east (should equal west)
distance = greatCircleDistance(lat1 distance = greatCircleDistance(lat1
...@@ -224,5 +233,5 @@ TEST(radix, greatCircleVolume) ...@@ -224,5 +233,5 @@ TEST(radix, greatCircleVolume)
, lon2 , lon2
, EARTH_RADIUS_MEAN , EARTH_RADIUS_MEAN
, EARTH_RADIUS_MEAN+1); , EARTH_RADIUS_MEAN+1);
EXPECT_FLOAT_EQ(24624134, volume); EXPECT_FLOAT_EQ(24625680, volume);
} }
...@@ -167,4 +167,9 @@ Real cspanf(Real value, Real begin, Real end) ...@@ -167,4 +167,9 @@ Real cspanf(Real value, Real begin, Real end)
} }
} }
Real hpaToAltitude(Real hpa)
{
return 0.3048 * (1 - std::pow(hpa/1013.25, 0.190284)) * 145366.45;
}
} // namespace radix } // namespace radix
...@@ -16,6 +16,12 @@ ...@@ -16,6 +16,12 @@
namespace radix namespace radix
{ {
/**
* @brief hpaToAltitude converts hectorpascals or millibars to altitude (meters)
* @param hpa hectorpascals in millibars
* @return altitude in meters
*/
Real hpaToAltitude(Real hpa);
/** /**
* @brief cspanf returns a value in the interval (begin,end] * @brief cspanf returns a value in the interval (begin,end]
* This function is used to reduce periodic variables to a standard range. * This function is used to reduce periodic variables to a standard range.
......
Supports Markdown
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