Skip to content
GitLab
Commit 09680b4b authored by LEFEBVREJP email's avatar LEFEBVREJP email
Browse files

Fixing awr met data loading. Something is still wrong with heights.

parent e0464fb8
Hide whitespace changes
Inline Side-by-side
radixio/gfsfile.cc
View file @ 09680b4b
......@@ -94,6 +94,7 @@ std::vector<float> GFSFile::mFact = { 1.0f, 1.0f, 1000.f, 1000.f, 1000.f, 60000.
3600.f, 60000.f };
GFSFile::GFSFile(std::string file)
: mFile(file)
, mStrcmp(9, 0.0f)
{
radix_tagged_line("GFSFile(" << file << ")" );
// initialize data structure
......@@ -103,6 +104,10 @@ GFSFile::GFSFile(std::string file)
// initialize
mLabel.expand(label);
mHeader.expand(header);
mHeader.latlon = false;
mHeader.global = false;
mHeader.gbldat = false;
mHeader.prime = false;
// calculate length of records
int ldat = mHeader.nx*mHeader.ny;
......@@ -122,6 +127,36 @@ GFSFile::GFSFile(std::string file)
mHeader.expand(header);
radix_tagged_line("Found grid: " << mHeader.toString());
}
//
// determine if this is a lat-lon grid
if(mHeader.size == 0.f)
{
mHeader.latlon = true;
}
if(mHeader.model_id.compare("RAMS") == 0)
{
mHeader.tang_lat = mHeader.pole_lat;
}
if(!mHeader.latlon)
{
//
// initialize grid conversion variable (into gbase)
stlmbr(mHeader.tang_lat, mHeader.ref_lon);
//
// use single point grid definition
radix_line("sync_xp=" << mHeader.sync_xp
<< " sync_yp=" << mHeader.sync_yp
<< " sync_lat=" << mHeader.sync_lat
<< " sync_lon=" << mHeader.sync_lon
<< " ref_lat=" << mHeader.ref_lat
<< " ref_lon=" << mHeader.ref_lon
<< " size=" << mHeader.size
<< " orient=" << mHeader.orient);
stcm1p(mHeader.sync_xp, mHeader.sync_yp
, mHeader.sync_lat, mHeader.sync_lon
, mHeader.ref_lat, mHeader.ref_lon
, mHeader.size, mHeader.orient);
}
int kol = 108;
int nrec = nndx;
int nlvl = mHeader.nz;
......@@ -212,6 +247,7 @@ std::pair<float, float> GFSFile::gbl2xy(float clat
<< clat << ","
<< clon << ","
<< sync_lat << ","
<< ref_lat << ","
<< sync_lon << ","
<< ref_lon << ")");
float tlat = clat;
......@@ -222,9 +258,11 @@ std::pair<float, float> GFSFile::gbl2xy(float clat
radix_tagged_line("\tcomputed y =" << result.second);
float tlon = clon;
// default PRIME section
if(tlon < 0.0f) tlon = 360.0f+tlon;
if(tlon > 360.0) tlon = tlon-360.0f;
if(!mHeader.prime)
{
if(tlon < 0.0f) tlon = 360.0f+tlon;
if(tlon > 360.0) tlon = tlon-360.0f;
}
tlon = tlon-sync_lon;
if(tlon < 0.0f) tlon = tlon+360.0f;
result.first = 1.0f+tlon/ref_lon;
......@@ -232,14 +270,13 @@ std::pair<float, float> GFSFile::gbl2xy(float clat
return result;
}
std::pair<float, float> GFSFile::cll2xy(float clat, float clon) const
std::pair<float, float> GFSFile::cnllxy(float clat, float clon) const
{
std::pair<float, float> result;
float almost1 = .99999;
float gamma = std::sin(mHeader.tang_lat*PI_BELOW_180);
float gamma = mStrcmp[0];
float dlat = clat;
float strcmp2 = cspanf(mHeader.ref_lon, -180, 180);
float dlong = cspanf(clon - strcmp2, -180, 180);
float dlong = cspanf(clon - mStrcmp[1], -180, 180);
dlong = dlong * PI_ON_180;
float gdlong = gamma * dlong;
float csdgam = 0.0, sndgam = 0.0;
......@@ -248,47 +285,176 @@ std::pair<float, float> GFSFile::cll2xy(float clat, float clon) const
// For gamma small or zero. avoids round-off error or division
// by zero in the case of mercator or near-mercator projections.
gdlong = gdlong * gdlong;
sndgam = dlong * (1.-1./6. * gdlong *
(1.-1./20. * gdlong *
(1.-1./42. * gdlong)));
csdgam = dlong * dlong * .5 *
(1.-1./12. * gdlong *
(1.-1./30. * gdlong *
(1.-1./56. * gdlong)));
sndgam = dlong * (1.f-1.f/6.f * gdlong *
(1.f-1.f/20.f * gdlong *
(1.f-1.f/42.f * gdlong)));
csdgam = dlong * dlong * .5f *
(1.f-1.f/12.f * gdlong *
(1.f-1.f/30.f * gdlong *
(1.f-1.f/56.f * gdlong)));
} else
{
sndgam = std::sin(gdlong)/gamma;
csdgam = (1.-std::cos(gdlong))/gamma/gamma;
csdgam = (1.f-std::cos(gdlong))/gamma/gamma;
}
float slat = std::sin(dlat*PI_ON_180);
if((slat >= almost1) || (slat <= -almost1))
{
result.first = 0.0;
result.second = 1./gamma;
result.first = 0.0f;
result.second = 1.f/gamma;
return result;
}
float mercy = .5 * std::log((1.+slat)/(1.-slat));
float mercy = .5f * std::log((1.f+slat)/(1.f-slat));
float gmercy = gamma * mercy;
float rhog1 = 0;
if(std::abs(gmercy) < .001)
float rhog1 = 0.f;
if(std::abs(gmercy) < .001f)
{
// For gamma small or zero. avoids round-off error or division
// by zero in the case of mercator or near-mercator projections.
rhog1 = mercy * ( 1. -.5 * gmercy *
(1.-1./3. * gmercy *
(1.-1./4. * gmercy)));
rhog1 = mercy * ( 1.f -.5f * gmercy *
(1.f-1.f/3.f * gmercy *
(1.f-1.f/4.f * gmercy)));
} else
{
rhog1 = (1. - std::exp(-gmercy)) / gamma;
rhog1 = (1.f - std::exp(-gmercy)) / gamma;
}
result.first = (1.-gamma*rhog1)*sndgam;
result.second = rhog1 + (1.-gamma*rhog1)*gamma*csdgam;
result.first = (1.f-gamma*rhog1)*sndgam;
result.second = rhog1 + (1.f-gamma*rhog1)*gamma*csdgam;
return result;
}
std::pair<float, float> GFSFile::cll2xy(float clat, float clon) const
{
radix_tagged_line("cll2xy(" << clat << "," << clon << ")");
std::pair<float, float> xi_eta = cnllxy(clat, clon);
radix_tagged_line("\txi=" << xi_eta.first
<< " eta=" << xi_eta.second);
float radius = EARTH_RADIUS_MEAN/1000.f;
float x = mStrcmp[2] + radius/mStrcmp[6]
* (xi_eta.first*mStrcmp[4] + xi_eta.second * mStrcmp[5]);
float y = mStrcmp[3] + radius/mStrcmp[6]
* (xi_eta.second*mStrcmp[4] - xi_eta.first * mStrcmp[5]);
radix_tagged_line("\tx=" << x << " y=" << y);
return std::make_pair(x,y);
}
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;
float radius = EARTH_RADIUS_MEAN/1000.f;
xi0 = (x -mStrcmp[2]) * mStrcmp[6]/radius;
eta0 = (y - mStrcmp[3]) * mStrcmp[6]/radius;
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> lat_lon = cxy2ll
}
}
void GFSFile::stlmbr(float tnglat, float xlong)
{
float radius = EARTH_RADIUS_MEAN / 1000.f;
mStrcmp[0] = std::sin(PI_ON_180*tnglat);
mStrcmp[1] = cspanf(xlong, -180., 180.);
mStrcmp[2] = 0.;
mStrcmp[3] = 0.;
mStrcmp[4] = 1.;
mStrcmp[5] = 0.;
mStrcmp[6] = radius;
std::pair<float,float> xi_eta = cnllxy(89., xlong);
mStrcmp[7] = 2. * xi_eta.second - mStrcmp[0]
* xi_eta.second * xi_eta.second;
xi_eta = cnllxy(-89., xlong);
mStrcmp[8] = 2. * xi_eta.second - mStrcmp[0]
* xi_eta.second * xi_eta.second;
}
void GFSFile::stcm1p(float x1, float y1, float xlat1, float xlong1
, float xlatg, float xlongg, float gridsz, float orient)
{
radix_tagged_line("stcm1p(" << x1 << "," << y1
<< "," << xlat1 << "," << xlong1
<< "," << xlatg << "," << xlongg
<< "," << gridsz << "," << orient << ")");
for(size_t i = 2; i < 4; ++i)
{
mStrcmp[i] = 0.f;
}
float turn = PI_ON_180 * (orient - mStrcmp[0]
* cspanf(xlongg - mStrcmp[1], -180., 180.));
radix_line("turn=" << turn);
mStrcmp[4] = std::cos(turn);
mStrcmp[5] = -std::sin(turn);
mStrcmp[6] = 1.f;
float cgszllResult = cgszll(xlatg, mStrcmp[1]);
radix_line("cgszll=" << cgszllResult);
mStrcmp[6] = gridsz * mStrcmp[6] / cgszllResult;
radix_line("mStrcmp[7]=" << mStrcmp[6]);
std::pair<float, float> a1 = cll2xy(xlat1, xlong1);
radix_line("x1a=" << a1.first << " y1a=" << a1.second);
mStrcmp[2] = mStrcmp[2] + x1 - a1.first;
mStrcmp[3] = mStrcmp[3] + y1 - a1.second;
radix_line("1=" << mStrcmp[0]
<< ", 2=" << mStrcmp[1]
<< ", 3=" << mStrcmp[2]
<< ", 4=" << mStrcmp[3]
<< ", 5=" << mStrcmp[4]
<< ", 6=" << mStrcmp[5]
<< ", 7=" << mStrcmp[6]
<< ", 8=" << mStrcmp[7]);
}
float GFSFile::cgszll(float xlat, float xlong) const
{
radix_tagged_line("cgszll(" << xlat << "," << xlong);
float slat = 0.f, ymerc = 0.f, efact = 0.f;
if(xlat > 89.995f)
{
// close to north pole
if(mStrcmp[0] > 0.9999f)
{// and to gamma == 1
return 2.f*mStrcmp[6];
}
efact = std::cos(PI_ON_180*xlat);
if(efact <= 0.f)
{
return 0.f;
} else
{
ymerc = -std::log(efact / (1.f + std::sin(PI_ON_180*xlat)));
}
} else if(xlat < -89.995f)
{
// close to south pole
if(mStrcmp[0] < -0.9999f)
{// and to gamma == -1.0
return 2.f*mStrcmp[6];
}
efact = std::cos(PI_ON_180*xlat);
if(efact <= 0.f)
{
return 0.f;
} else
{
ymerc = std::log(efact / (1.f - std::sin(PI_ON_180*xlat)));
}
} else
{
slat = std::sin(PI_ON_180*xlat);
ymerc = std::log((1.f+slat)/(1.f-slat))/2.f;
}
return mStrcmp[6] * std::cos(PI_ON_180*xlat)*std::exp(mStrcmp[0]*ymerc);
}
std::pair<int,int> GFSFile::nearestPoint(float lat, float lon) const
{
std::pair<float,float> point;
if(mHeader.size == 0)
radix_tagged_line("nearstPoint("<<lat<<","<<lon<<")");
radix_tagged_line("\tlatlon=" << std::boolalpha << mHeader.latlon);
if(mHeader.latlon)
{
point = gbl2xy(lat, lon
, mHeader.sync_lat, mHeader.ref_lat
......@@ -406,7 +572,19 @@ std::vector<std::vector<float>> GFSFile::query(float lat
int ll = mLabel.il;
// convert level number to array index because input data
// level index starts at 0 for the surface
if(ll != lp || irec == (matchingIndex.size()-1)) lp = ll;
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;
// }
// }
lp = ll;
}
// find the variable array element number - match the input
// variable with its position as indicated in the index record
......
radixio/gfsfile.hh
View file @ 09680b4b
......@@ -48,6 +48,7 @@ private:
std::vector<float> mRecordTimes;
std::vector<int> mNumVarb;
std::vector<float> mHeight;
std::vector<float> mStrcmp;
std::string mStartTime;
std::string mEndTime;
std::string mProfileTime;
......@@ -100,6 +101,10 @@ protected:
}; // class Label
class Header {
public:
bool prime;
bool latlon;
bool gbldat;
bool global;
std::string model_id;
int icx;
int mn;
......@@ -180,12 +185,22 @@ protected:
, float ref_lon) const;
/**
* @brief cll2xy converts geographic coordinates lat&lon to xy coordinates for grid lookup
* @brief cnllxy transforms lat-lon to cononical(equator-centered,radian unit) coordinates
* @param clat real latitude
* @param clon real longitude
* @return std::pair<float,float>
*/
std::pair<float,float> cll2xy(float clat, float clon) const;
std::pair<float,float> cnllxy(float clat, float clon) const;
std::pair<float, float> cll2xy(float clat, float clon) const;
std::pair<float,float> cg2cxy(float x, float y, float ug, float vg);
void stlmbr(float tnglat, float xlong);
void stcm1p(float x1, float y1, float xlat1, float xlong1
, float xlatg, float xlongg, float gridsz, float orient);
float cgszll(float xlat, float xlong) const;
public:
/**
......
radixmath/constants.hh
View file @ 09680b4b
......@@ -16,7 +16,7 @@ namespace radix
*/
typedef unsigned short index_t;
typedef double Real;
const Real EARTH_RADIUS_MEAN = 6371000.0;
const Real EARTH_RADIUS_MEAN = 6371200.0;
const Real wgs84_minor = 6356752.314245;
const Real kEpsilon = 1.0E-12;
const Real lookAheadEpsilon = 1.0E-9;
......
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