Skip to content
GitLab
Commit 52d0a190 authored by LEFEBVREJP email's avatar LEFEBVREJP email
Browse files

Fixing wind speed and direction for awr met source.

parent 09680b4b
Hide whitespace changes
Inline Side-by-side
radixio/gfsfile.cc
View file @ 52d0a190
......@@ -270,6 +270,54 @@ std::pair<float, float> GFSFile::gbl2xy(float clat
return result;
}
std::pair<float, float> GFSFile::cnxyll(float xi, float eta) const
{
float gamma = mStrcmp[0];
float cgeta = 1.f - gamma * eta;
float gxi = gamma * xi;
// calculate equivalent mercator coordinate
float arg2 = eta + (eta*cgeta - gxi * xi);
float arg1 = gamma * arg2;
float xlat = 0, xlong = 0, temp = 0, ymerc = 0, along = 0;
// distance to north (or south) pole is zero (or imaginary)
if(arg1 >= 1.f)
{
xlat = std::copysign(90., mStrcmp[0]);
xlong = 90. + xlat;
return std::make_pair(xlat, xlong);
}
if(std::abs(arg1) < 0.01f)
{
// this avoids round-off error or divide-by zero in case of mercator projects
temp = std::pow(arg1 / (2.f - arg1), 2);
ymerc = arg2 / (2.f -arg1) * (1.f + temp *
(1.f/3.f + temp *
(1.f/5.f + temp *
1.f/7.f)));
} else
{
//code for moderate values of gamma
ymerc = - std::log(1.f-arg1) / 2.f / gamma;
}
// convert ymerc to latitude
temp = std::exp(- std::abs(ymerc));
xlat = std::copysign(std::atan2((1.f-temp)*(1.f+temp), 2.f*temp),ymerc);
// find longitudes
if(std::abs(gxi) < 0.01f*cgeta)
{ // this avoids round-off error or divide-by zero in case of mercator projects
temp = std::pow(gxi / cgeta, 2);
ymerc = xi / cgeta * (1.f - temp *
(1.f/3.f - temp *
(1.f/5.f - temp *
1.f/7.f)));
} else
{
along = std::atan2(gxi, cgeta) /gamma;
}
xlong = mStrcmp[1] + PI_BELOW_180 * along;
xlat = xlat * PI_BELOW_180;
return std::make_pair(xlat, xlong);
}
std::pair<float, float> GFSFile::cnllxy(float clat, float clon) const
{
std::pair<float, float> result;
......@@ -340,6 +388,38 @@ std::pair<float, float> GFSFile::cll2xy(float clat, float clon) const
return std::make_pair(x,y);
}
std::pair<float, float> GFSFile::cxy2ll(float x, float y) const
{
radix_tagged_line("cxy2ll(" << x
<< "," << y << ")");
float radius = EARTH_RADIUS_MEAN / 1000.f;
float xi0 = (x -mStrcmp[2]) * mStrcmp[6] / radius;
float eta0 = (y - mStrcmp[3]) * mStrcmp[6] / radius;
float xi = xi0 * mStrcmp[4] - eta0 * mStrcmp[5];
float eta = eta0 * mStrcmp[4] + xi0 * mStrcmp[5];
std::pair<float, float> ll = cnxyll(xi, eta);
radix_line("\tcnxy2ll result lat=" << ll.first
<< " lon=" << ll.second);
float xlong = cspanf(ll.second, -180.f, 180.f);
return std::make_pair(ll.first, xlong);
}
std::pair<float, float> GFSFile::cg2cll(float xlat, float xlong, float ug, float vg) const
{
float along = cspanf(xlong - mStrcmp[1], -180.f, 180.f);
float rot = - mStrcmp[0] + along;
// allow cartographic wind vector transformations everywhere
// with rotation to nominal longitudes at the poles, to match u,v values
// on a lat-lon grid
float slong = std::sin(PI_ON_180 * rot);
float clong = std::cos(PI_ON_180 * rot);
float xpolg = slong * mStrcmp[4] + clong * mStrcmp[5];
float ypolg = clong * mStrcmp[4] - slong * mStrcmp[5];
float vn = ypolg * ug + xpolg * vg;
float ue = ypolg * vg + xpolg * ug;
return std::make_pair(ue, vn);
}
std::pair<float, float> GFSFile::cg2cxy(float x, float y, float ug, float vg)
{
float xpolg = 0.f, ypolg = 0.f, temp = 0.f, xi0 = 0.f, eta0 = 0.f;
......@@ -349,10 +429,20 @@ std::pair<float, float> GFSFile::cg2cxy(float x, float y, float ug, float vg)
xpolg = mStrcmp[5] - mStrcmp[0] * xi0;
ypolg = mStrcmp[4] - mStrcmp[0] * eta0;
temp = std::sqrt( std::pow(xpolg, 2.f) + std::pow(ypolg, 2.f));
if(temp <= 1e03)
std::pair<float,float> xy;
if(temp <= 1e-3)
{
std::pair<float, float> ll = cxy2ll(x, y);
xy = cg2cll(ll.first, ll.second, ug, vg);
} else
{
// std::pair<float, float> lat_lon = cxy2ll
// use vector alegbra instread of time consuming trig
xpolg = xpolg / temp;
ypolg = ypolg / temp;
xy.first = ypolg * ug - xpolg * vg;
xy.second = ypolg * vg + xpolg * ug;
}
return xy;
}
void GFSFile::stlmbr(float tnglat, float xlong)
......@@ -574,15 +664,12 @@ std::vector<std::vector<float>> GFSFile::query(float lat
// level index starts at 0 for the surface
if(ll != lp || irec == (matchingIndex.size()-1))
{
// if(lp != 0)
// {
// if(!mHeader.latlon)
// {
// std::pair<float,float> result = cg2cxy();
// utw[lp] = result.first;
// vtw[lp] = result.second;
// }
// }
if(lp != 0 && !mHeader.latlon)
{
std::pair<float,float> xy = cg2cxy(x-1, y-1, utw[lp], vtw[lp]);
utw[lp] = xy.first;
vtw[lp] = xy.second;
}
lp = ll;
}
......@@ -599,7 +686,8 @@ std::vector<std::vector<float>> GFSFile::query(float lat
if( varb.compare("TEMP") == 0
|| varb.compare("T02M") == 0
|| varb.compare("TMPS") == 0
|| varb.compare("DP2M") == 0) vdata[kvar][ll] = vdata[kvar][ll]-273.16f;
|| varb.compare("DP2M") == 0)
vdata[kvar][ll] = vdata[kvar][ll]-273.16f;
// save the surface pressure and terrain mHeight for scaling
// of the vertical coordinate system (mHeight = signma*scaling)
......
radixio/gfsfile.hh
View file @ 52d0a190
......@@ -184,6 +184,14 @@ protected:
, float sync_lon
, float ref_lon) const;
/**
* @brief cnxyll transforms cononical (equator-centered, radian unit) coordinates
* @param xi
* @param eta
* @return std::pair<float,float> xlat, xlong
*/
std::pair<float,float> cnxyll(float xi, float eta) const;
/**
* @brief cnllxy transforms lat-lon to cononical(equator-centered,radian unit) coordinates
* @param clat real latitude
......@@ -193,8 +201,11 @@ protected:
std::pair<float,float> cnllxy(float clat, float clon) const;
std::pair<float, float> cll2xy(float clat, float clon) const;
std::pair<float, float> cxy2ll(float x, float y) const;
std::pair<float,float> cg2cxy(float x, float y, float ug, float vg);
std::pair<float, float> cg2cll(float xlat, float xlong
, float ug, float vg) const;
void stlmbr(float tnglat, float xlong);
......
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