WSL/SLF GitLab Repository

Commit f3d55dbe authored by Mathias Bavay's avatar Mathias Bavay
Browse files

More unsigned int / size_t cleanup...

parent 2e1a1e2a
......@@ -168,7 +168,7 @@ void A3DIO::readStationData(const Date& timestamp, std::vector<StationData>& vec
read2DStations(timestamp, tmpvecS);
cleanup();
for(unsigned int i=0; i<tmpvecS.size(); i++) {
for(size_t i=0; i<tmpvecS.size(); i++) {
vecStation.push_back(tmpvecS[i]);
}
}
......@@ -337,7 +337,7 @@ bool A3DIO::readMeteoDataLine(std::string& line, MeteoData& tmpdata, std::string
//throw InvalidFormatException("Premature End of Line or no data for date " + date_in + " found in File " + filename, AT);
}
for (unsigned int ii=0; ii<4; ii++) {
for (size_t ii=0; ii<4; ii++) {
if (!IOUtils::convertString(tmp_ymdh[ii], tmpvec.at(ii), std::dec))
throw InvalidFormatException(filename + ": " + line, AT);
}
......@@ -346,7 +346,7 @@ bool A3DIO::readMeteoDataLine(std::string& line, MeteoData& tmpdata, std::string
//Read rest of line with values ta, iswr, vw, rh, ea, hnw
for (unsigned int ii=0; ii<6; ii++) { //go through the columns
for (size_t ii=0; ii<6; ii++) { //go through the columns
if (!IOUtils::convertString(tmp_values[ii], tmpvec.at(ii+4), std::dec)) {
throw InvalidFormatException(filename + ": " + line, AT);
}
......@@ -366,16 +366,15 @@ bool A3DIO::readMeteoDataLine(std::string& line, MeteoData& tmpdata, std::string
void A3DIO::read2DStations(const Date& timestamp, std::vector<StationData>& vecStation)
{
unsigned int stations=0;
std::vector<std::string> filenames = std::vector<std::string>();
std::map<std::string, unsigned int> hashStations = std::map<std::string, unsigned int>();
std::map<std::string, size_t> hashStations = std::map<std::string, size_t>();
constructMeteo2DFilenames(timestamp, timestamp, filenames);
stations = getNrOfStations(filenames, hashStations);
const size_t stations = getNrOfStations(filenames, hashStations);
vecStation.resize(stations); //we receive an empty vector, so we need to allocate the room it needs
try {
for (unsigned int ii=0; ii<filenames.size(); ii++){
for (size_t ii=0; ii<filenames.size(); ii++){
read2DMeteoHeader(filenames[ii], hashStations, vecStation);
}
} catch(...) {
......@@ -400,8 +399,8 @@ void A3DIO::read2DStations(const Date& timestamp, std::vector<StationData>& vecS
void A3DIO::read2DMeteo(std::vector< std::vector<MeteoData> >& vecMeteo)
{
unsigned int stations=0, bufferindex=0;
std::map<std::string, unsigned int> hashStations = std::map<std::string, unsigned int>();
size_t bufferindex=0;
std::map<std::string, size_t> hashStations = std::map<std::string, size_t>();
std::vector<std::string> filenames = std::vector<std::string>();
//Requirement: meteo1D data must exist:
......@@ -415,7 +414,7 @@ void A3DIO::read2DMeteo(std::vector< std::vector<MeteoData> >& vecMeteo)
Date endDate(vecMeteo[0][vecMeteo[0].size()-1].date.getJulianDate(), in_tz, false);
constructMeteo2DFilenames(startDate, endDate, filenames);//get all files for all years
stations = getNrOfStations(filenames, hashStations);
const size_t stations = getNrOfStations(filenames, hashStations);
constructMeteo2DFilenames(startDate, startDate, filenames);//get filenames for current year
std::cerr << "[I] Number of 2D meteo stations: " << stations << std::endl;
......@@ -427,24 +426,24 @@ void A3DIO::read2DMeteo(std::vector< std::vector<MeteoData> >& vecMeteo)
std::vector<StationData> tmpvecS = std::vector<StationData>(stations); //stores unique stations
try {
for (unsigned int ii=0; ii<filenames.size(); ii++){
for (size_t ii=0; ii<filenames.size(); ii++){
read2DMeteoHeader(filenames[ii], hashStations, tmpvecS);
}
//init vecStation with proper StationData, vecMeteo with nodata
for (unsigned int jj=0; jj<tmpvecS.size(); jj++){
for (size_t jj=0; jj<tmpvecS.size(); jj++){
vecMeteo.push_back( std::vector<MeteoData>() );
vecMeteo[jj+1].reserve(buffer_reserve);
MeteoData tmd; //create an empty data set
tmd.meta = tmpvecS[jj]; //add meta data to that data set
for (unsigned int ii=0; ii<vecMeteo[0].size(); ii++){
for (size_t ii=0; ii<vecMeteo[0].size(); ii++){
//NOTE: there needs to be the same amount of 1D and 2D data
vecMeteo[jj+1].push_back(tmd);
}
}
do {
unsigned int currentindex = bufferindex;
const size_t currentindex = bufferindex;
read2DMeteoData(filenames[0], "nswc", hashStations, vecMeteo, bufferindex);
bufferindex = currentindex;
read2DMeteoData(filenames[1], "rh", hashStations, vecMeteo, bufferindex);
......@@ -474,8 +473,8 @@ void A3DIO::read2DMeteo(std::vector< std::vector<MeteoData> >& vecMeteo)
}
//clean data and convert the units
for (unsigned int ii=1; ii<vecMeteo.size(); ii++) { //loop over all stations except 1D Meteo
for (unsigned int jj=0; jj<vecMeteo[ii].size(); jj++){ //Meteo1D data already cleaned
for (size_t ii=1; ii<vecMeteo.size(); ii++) { //loop over all stations except 1D Meteo
for (size_t jj=0; jj<vecMeteo[ii].size(); jj++){ //Meteo1D data already cleaned
//vecMeteo[ii][jj].cleanData();
convertUnits(vecMeteo[ii][jj]);
}
......@@ -512,7 +511,7 @@ void A3DIO::constructMeteo2DFilenames(const Date& startDate, const Date& endDate
filenames.push_back(wdirFilename);
}
for (unsigned int ii=0; ii<filenames.size(); ii++) {
for (size_t ii=0; ii<filenames.size(); ii++) {
if (!IOUtils::fileExists(filenames[ii])) {
throw FileNotFoundException(filenames[ii], AT);
}
......@@ -520,12 +519,12 @@ void A3DIO::constructMeteo2DFilenames(const Date& startDate, const Date& endDate
}
unsigned int A3DIO::getNrOfStations(std::vector<std::string>& filenames, std::map<std::string, unsigned int>& hashStations)
size_t A3DIO::getNrOfStations(std::vector<std::string>& filenames, std::map<std::string, size_t>& hashStations)
{
std::vector<std::string> tmpvec;
std::string line_in="";
for (unsigned int ii=0; ii<filenames.size(); ii++) {
for (size_t ii=0; ii<filenames.size(); ii++) {
//cout << *it << endl;
std::string filename = filenames[ii];
......@@ -537,11 +536,11 @@ unsigned int A3DIO::getNrOfStations(std::vector<std::string>& filenames, std::ma
IOUtils::skipLines(fin, 4, eoln);
getline(fin, line_in, eoln); //5th line holds the names of the stations
unsigned int cols = IOUtils::readLineToVec(line_in, tmpvec);
size_t cols = IOUtils::readLineToVec(line_in, tmpvec);
if ( cols > 4) { // if there are any stations
//check each station name and whether it's already hashed, otherwise: hash!
for (unsigned int ii=4; ii<cols; ii++) {
unsigned int tmp_int = hashStations.count(tmpvec.at(ii));
for (size_t ii=4; ii<cols; ii++) {
size_t tmp_int = hashStations.count(tmpvec.at(ii));
if (tmp_int == 0) {
//cout << "Found station: " << tmpvec.at(ii) << endl;
hashStations[tmpvec.at(ii)] = hashStations.size();
......@@ -555,12 +554,11 @@ unsigned int A3DIO::getNrOfStations(std::vector<std::string>& filenames, std::ma
}
void A3DIO::read2DMeteoData(const std::string& filename, const std::string& parameter,
std::map<std::string,unsigned int>& hashStations,
std::vector< std::vector<MeteoData> >& vecM, unsigned int& bufferindex)
std::map<std::string,size_t>& hashStations,
std::vector< std::vector<MeteoData> >& vecM, size_t& bufferindex)
{
std::string line_in = "";
unsigned int columns;
std::vector<std::string> tmpvec, vec_names;
Date tmp_date;
int tmp_ymdh[4];
......@@ -571,11 +569,11 @@ void A3DIO::read2DMeteoData(const std::string& filename, const std::string& para
throw FileAccessException(filename, AT);
}
char eoln = IOUtils::getEoln(fin); //get the end of line character for the file
const char eoln = IOUtils::getEoln(fin); //get the end of line character for the file
IOUtils::skipLines(fin, 4, eoln); //skip first 4 lines
getline(fin, line_in, eoln); //line containing UNIQUE station names
columns = IOUtils::readLineToVec(line_in, vec_names);
const size_t columns = IOUtils::readLineToVec(line_in, vec_names);
if (columns < 4) {
cleanup();
throw InvalidFormatException("[E] Premature end of line in file " + filename + " (line does not even contain full timestamp)", AT);
......@@ -592,7 +590,7 @@ void A3DIO::read2DMeteoData(const std::string& filename, const std::string& para
break;
}
const unsigned int cols_found = IOUtils::readLineToVec(line_in, tmpvec);
const size_t cols_found = IOUtils::readLineToVec(line_in, tmpvec);
if (cols_found!=columns) { //Every station has to have its own column
cleanup();
std::stringstream ss;
......@@ -602,7 +600,7 @@ void A3DIO::read2DMeteoData(const std::string& filename, const std::string& para
throw InvalidFormatException(ss.str(), AT);
}
for (unsigned int ii=0; ii<4; ii++) {
for (size_t ii=0; ii<4; ii++) {
if (!IOUtils::convertString(tmp_ymdh[ii], tmpvec[ii], std::dec)) {
cleanup();
throw InvalidFormatException("[E] Check date columns in " + filename + " at line " + line_in, AT);
......@@ -613,8 +611,8 @@ void A3DIO::read2DMeteoData(const std::string& filename, const std::string& para
MeteoData& currentMeteoData = vecM[0][bufferindex]; //1D Element to synchronize date
if (tmp_date == currentMeteoData.date) {
//Read in data
for (unsigned int ii=4; ii<columns; ii++) {
unsigned int stationnr = hashStations[vec_names.at(ii)];
for (size_t ii=4; ii<columns; ii++) {
size_t stationnr = hashStations[vec_names.at(ii)];
MeteoData& tmpmd = vecM[stationnr][bufferindex];
tmpmd.date = tmp_date;
......@@ -656,11 +654,10 @@ void A3DIO::read2DMeteoData(const std::string& filename, const std::string& para
cleanup();
}
void A3DIO::read2DMeteoHeader(const std::string& filename, std::map<std::string,unsigned int>& hashStations,
void A3DIO::read2DMeteoHeader(const std::string& filename, std::map<std::string,size_t>& hashStations,
std::vector<StationData>& vecS)
{
std::string line_in = "";
unsigned int columns = 0;
std::vector<std::string> vec_altitude, vec_xcoord, vec_ycoord, vec_names;
fin.clear();
......@@ -669,13 +666,13 @@ void A3DIO::read2DMeteoHeader(const std::string& filename, std::map<std::string,
throw FileAccessException(filename, AT);
}
char eoln = IOUtils::getEoln(fin); //get the end of line character for the file
const char eoln = IOUtils::getEoln(fin); //get the end of line character for the file
IOUtils::skipLines(fin, 1, eoln);
//Read all relevant lines in
getline(fin, line_in, eoln); //Altitude
columns = IOUtils::readLineToVec(line_in, vec_altitude);
const size_t columns = IOUtils::readLineToVec(line_in, vec_altitude);
getline(fin, line_in, eoln); //xcoord
if (IOUtils::readLineToVec(line_in, vec_xcoord) != columns) {
......@@ -699,9 +696,9 @@ void A3DIO::read2DMeteoHeader(const std::string& filename, std::map<std::string,
//Check for duplicate station names within one file ... station names need to be unique!
vector<string> vec_dup = vec_names;
for (unsigned int ii=0; ii<vec_names.size(); ii++) {
for (size_t ii=0; ii<vec_names.size(); ii++) {
const string& tmp = vec_names[ii];
for (unsigned int jj=0; jj<vec_dup.size(); jj++) {
for (size_t jj=0; jj<vec_dup.size(); jj++) {
if (jj != ii) {
if (vec_dup[jj] == tmp) {
cleanup();
......@@ -714,8 +711,8 @@ void A3DIO::read2DMeteoHeader(const std::string& filename, std::map<std::string,
//Build Coords object to convert easting/northing values to lat/long in WGS84
Coords coordinate(coordin, coordinparam);
for (unsigned int ii=4; ii<columns; ii++) {
unsigned int stationnr = hashStations[vec_names.at(ii)];
for (size_t ii=4; ii<columns; ii++) {
const size_t stationnr = hashStations[vec_names.at(ii)];
double altitude, easting, northing;
std::string stationName;
if ((!IOUtils::convertString(altitude, vec_altitude.at(ii), std::dec))
......@@ -747,7 +744,7 @@ void A3DIO::readSpecialPoints(std::vector<Coords>& pts)
throw FileAccessException(filename,AT);
}
char eoln = IOUtils::getEoln(fin); //get the end of line character for the file
const char eoln = IOUtils::getEoln(fin); //get the end of line character for the file
while (!fin.eof()) {
getline(fin, line_in, eoln);
......@@ -772,7 +769,7 @@ void A3DIO::readSpecialPoints(std::vector<Coords>& pts)
//Now put everything into the output vector TODO: don't do any intermediate steps... copy directly into vector!
Coords tmp_pts;
for (unsigned int jj=0; jj<mypts.size(); jj++) {
for (size_t jj=0; jj<mypts.size(); jj++) {
tmp_pts.setGridIndex(mypts.at(jj).first, mypts.at(jj).second, IOUtils::inodata, false);
pts.push_back(tmp_pts);
}
......@@ -782,11 +779,11 @@ int A3DIO::create1DFile(const std::vector< std::vector<MeteoData> >& data)
{//TODO: add check for stations' positions
std::string tmp_path;
cfg.getValue("METEOPATH", "Output", tmp_path);
const unsigned int sta_nr = data.size();
const size_t sta_nr = data.size();
if(sta_nr==0) return EXIT_FAILURE;
for(unsigned int ii=0; ii<sta_nr; ii++) {
const unsigned int size = data[ii].size();
for(size_t ii=0; ii<sta_nr; ii++) {
const size_t size = data[ii].size();
if(size>0) {
const std::string filename = tmp_path+"/meteo1D_"+data[ii][0].meta.getStationID()+".txt";
std::ofstream file(filename.c_str(), std::ios::out | std::ios::trunc);
......@@ -803,7 +800,7 @@ int A3DIO::create1DFile(const std::vector< std::vector<MeteoData> >& data)
file << "YYYY MM DD HH ta iswr vw rh ea nswc" << endl;
file.flags ( std::ios::fixed );
for(unsigned int j=0; j<size; j++) {
for(size_t j=0; j<size; j++) {
int yyyy, mm, dd, hh;
Date tmp_date(data[ii][j].date.getJulianDate(true), out_tz);
tmp_date.getDate(yyyy, mm, dd, hh);
......@@ -846,11 +843,11 @@ int A3DIO::writeHeader(std::ofstream &file, const std::vector< std::vector<Meteo
std::ostringstream str_altitudes;
std::ostringstream str_eastings;
std::ostringstream str_northings;
const unsigned int sta_nr = data.size();
const size_t sta_nr = data.size();
if(sta_nr==0) return EXIT_FAILURE;
file << "X:\\filepath " << parameter_name <<endl;
for(unsigned int ii=0; ii<sta_nr; ii++) {
for(size_t ii=0; ii<sta_nr; ii++) {
if(data[ii].size() > 0) {
str_altitudes << data[ii][0].meta.position.getAltitude() << " ";
str_eastings << data[ii][0].meta.position.getEasting() << " ";
......@@ -861,7 +858,7 @@ int A3DIO::writeHeader(std::ofstream &file, const std::vector< std::vector<Meteo
file << "YY MM DD HH " << str_eastings.str() << endl; //easting
file << "YY MM DD HH " << str_northings.str() << endl; //northing
file << "YYYY MM DD HH";
for(unsigned int ii=0; ii<sta_nr; ii++) {
for(size_t ii=0; ii<sta_nr; ii++) {
if(data[ii].size() > 0) {
file << " " << data[ii][0].meta.getStationID();
}
......@@ -892,9 +889,9 @@ int A3DIO::write2DmeteoFile(const std::vector< std::vector<MeteoData> >& data,
const unsigned int& parindex, const std::string& fileprefix,
const std::string& label)
{//HACK: we assume that all stations have data that is time synchronized...
const unsigned int sta_nr = data.size();
const size_t sta_nr = data.size();
if(sta_nr==0) return EXIT_FAILURE;
const unsigned int nb_timesteps = data[0].size();
const size_t nb_timesteps = data[0].size();
if(nb_timesteps==0) return EXIT_FAILURE;
std::ofstream file;
......@@ -905,7 +902,7 @@ int A3DIO::write2DmeteoFile(const std::vector< std::vector<MeteoData> >& data,
open2DFile(data, fileprefix, label, startyear, file);
file.flags ( ios::fixed );
for(unsigned int ii=0; ii<nb_timesteps; ii++) {
for(size_t ii=0; ii<nb_timesteps; ii++) {
tmp_date.setDate(data[0][ii].date.getJulianDate(true), out_tz);
tmp_date.getDate(year, month, day, hour);
if(year!=startyear) {
......@@ -919,7 +916,7 @@ int A3DIO::write2DmeteoFile(const std::vector< std::vector<MeteoData> >& data,
file.fill('0');
file << setw(4) << year << " " << setw(2) << month << " " << setw(2) << day << " " << setw(2) << hour;
file.fill(' ');
for(unsigned int j=0; j<sta_nr; j++) {
for(size_t j=0; j<sta_nr; j++) {
double value = data[j][ii].param(parindex);
if(value==IOUtils::nodata) {
file << " " << setw(7) << setprecision(0) << IOUtils::nodata;
......
......@@ -46,11 +46,11 @@ class A3DIO : public IOInterface {
virtual void readStationData(const Date& date, std::vector<StationData>& vecStation);
virtual void readMeteoData(const Date& dateStart, const Date& dateEnd,
virtual void readMeteoData(const Date& dateStart, const Date& dateEnd,
std::vector< std::vector<MeteoData> >& vecMeteo,
const unsigned int& stationindex=IOUtils::npos);
virtual void writeMeteoData(const std::vector< std::vector<MeteoData> >& vecMeteo,
virtual void writeMeteoData(const std::vector< std::vector<MeteoData> >& vecMeteo,
const std::string& name="");
virtual void readSpecialPoints(std::vector<Coords>& pts);
......@@ -69,19 +69,19 @@ class A3DIO : public IOInterface {
bool readMeteoDataLine(std::string& line, MeteoData& tmpdata, std::string filename);
void convertUnits(MeteoData& meteo);
void cleanup() throw();
void read2DMeteoData(const std::string&, const std::string&, std::map<std::string,unsigned int>& hashStations,
std::vector< std::vector<MeteoData> >&, unsigned int& bufferindex);
void read2DMeteoHeader(const std::string& filename, std::map<std::string, unsigned int>& hashStations,
void read2DMeteoData(const std::string&, const std::string&, std::map<std::string,size_t>& hashStations,
std::vector< std::vector<MeteoData> >&, size_t& bufferindex);
void read2DMeteoHeader(const std::string& filename, std::map<std::string, size_t>& hashStations,
std::vector<StationData>&);
unsigned int getNrOfStations(std::vector<std::string>& filenames,
std::map<std::string, unsigned int>& hashStations);
size_t getNrOfStations(std::vector<std::string>& filenames,
std::map<std::string, size_t>& hashStations);
int create1DFile(const std::vector< std::vector<MeteoData> >& data);
int writeHeader(std::ofstream &file, const std::vector< std::vector<MeteoData> >& stations, const std::string parameter_name);
void open2DFile(const std::vector< std::vector<MeteoData> >& stations,
const std::string& fileprefix, const std::string& label, const double& year,
std::ofstream& file);
int write2DmeteoFile(const std::vector< std::vector<MeteoData> >& data, const unsigned int& parindex,
int write2DmeteoFile(const std::vector< std::vector<MeteoData> >& data, const unsigned int& parindex,
const std::string& filename, const std::string& label);
void write2DMeteo(const std::vector< std::vector<MeteoData> >& data);
......
......@@ -146,8 +146,8 @@ bool IOManager::read_filtered_cache(const Date& start_date, const Date& end_date
{
if ((start_date >= fcache_start) && (end_date <= fcache_end)){
//it's already in the filtered_cache, so just copy the requested slice
for (unsigned int ii=0; ii<filtered_cache.size(); ii++){
unsigned int startpos = IOUtils::seek(start_date, filtered_cache[ii], false);
for (size_t ii=0; ii<filtered_cache.size(); ii++){
size_t startpos = IOUtils::seek(start_date, filtered_cache[ii], false);
if (startpos == IOUtils::npos){
if (filtered_cache[ii].size() > 0){
if (filtered_cache[ii][0].date <= end_date){
......@@ -158,8 +158,8 @@ bool IOManager::read_filtered_cache(const Date& start_date, const Date& end_date
if (startpos != IOUtils::npos){
vec_meteo.push_back(vector<MeteoData>());
size_t index = vec_meteo.size()-1;
for (unsigned int jj=startpos; jj<filtered_cache[ii].size(); jj++){
const size_t index = vec_meteo.size()-1;
for (size_t jj=startpos; jj<filtered_cache[ii].size(); jj++){
const MeteoData& md = filtered_cache[ii][jj];
if (md.date <= end_date){
vec_meteo[index].push_back(md);
......@@ -197,8 +197,8 @@ size_t IOManager::getMeteoData(const Date& i_date, METEO_TIMESERIE& vecMeteo)
//1. Check whether user wants raw data or processed data
if (processing_level == IOManager::raw){
rawio.readMeteoData(i_date-Duration(0.001, 0.), i_date+Duration(0.001, 0.), vec_cache);
for (unsigned int ii=0; ii<vec_cache.size(); ii++){
unsigned int index = IOUtils::seek(i_date, vec_cache[ii], true);
for (size_t ii=0; ii<vec_cache.size(); ii++){
const size_t index = IOUtils::seek(i_date, vec_cache[ii], true);
if (index != IOUtils::npos)
vecMeteo.push_back(vec_cache[ii][index]); //Insert station into vecMeteo
}
......@@ -218,15 +218,15 @@ size_t IOManager::getMeteoData(const Date& i_date, METEO_TIMESERIE& vecMeteo)
getMeteoData(i_date-proc_properties.time_before, i_date+proc_properties.time_after, vec_cache);
//vec_cache is either filtered or unfiltered, in any case it is wise to buffer it
for (unsigned int ii=0; ii<vec_cache.size(); ii++){//resampling for every station
for (size_t ii=0; ii<vec_cache.size(); ii++){//resampling for every station
if ((IOManager::resampled & processing_level) == IOManager::resampled){
//cout << "Resampling data for station " << ii << " (" << vec_cache[ii].size() << " elements)" << endl;
unsigned int position = meteoprocessor.resample(i_date, vec_cache[ii]);
const size_t position = meteoprocessor.resample(i_date, vec_cache[ii]);
if (position != IOUtils::npos)
vecMeteo.push_back(vec_cache[ii][position]);
} else { //only filtering activated
unsigned int index = IOUtils::seek(i_date, vec_cache[ii], true);
const size_t index = IOUtils::seek(i_date, vec_cache[ii], true);
if (index != IOUtils::npos)
vecMeteo.push_back(vec_cache[ii][index]); //Insert station into vecMeteo
}
......
......@@ -167,8 +167,8 @@ class MeteoData {
std::map<std::string, double> extraparameters; ///<All non-standard meteo parameters will end up in this map
std::map<std::string, double*> mapParameterByName; ///<Associate name and meteo parameter
std::map<unsigned int, double*> meteoparam; ///<Associate an unsigned int with every meteo parameter
std::map<unsigned int, std::string> meteoparamname; ///<Associate a name with every meteo parameter
std::map<size_t, double*> meteoparam; ///<Associate an unsigned int with every meteo parameter
std::map<size_t, std::string> meteoparamname; ///<Associate a name with every meteo parameter
size_t nrOfAllParameters;
bool resampled; ///<set this to true if MeteoData is result of resampling
......
......@@ -27,7 +27,7 @@ namespace mio {
//Quake3 fast 1/x² approximation
// For Magic Derivation see: Chris Lomont http://www.lomont.org/Math/Papers/2003/InvSqrt.pdf
// Credited to Greg Walsh.
// 32 Bit float magic number
// 32 Bit float magic number
#define SQRT_MAGIC_D 0x5f3759df
#define SQRT_MAGIC_F 0x5f375a86
......@@ -35,7 +35,7 @@ namespace mio {
//on a large scale interpolation test on TA, max relative error is 1e-6
inline float invSqrt(const float x) {
const float xhalf = 0.5f*x;
union {
// get bits for floating value
float x;
......@@ -48,7 +48,7 @@ inline float invSqrt(const float x) {
inline double invSqrt(const double x) {
const double xhalf = 0.5f*x;
union {
// get bits for floating value
float x;
......@@ -77,7 +77,7 @@ double Interpol2D::getReferenceAltitude(const DEMObject& dem)
if(dem.min_altitude!=IOUtils::nodata && dem.max_altitude!=IOUtils::nodata) {
//we use the median elevation as the reference elevation for reprojections
ref_altitude = 0.5 * (dem.min_altitude+dem.max_altitude);
ref_altitude = 0.5 * (dem.min_altitude+dem.max_altitude);
} else {
//since there is nothing else that we can do, we use an arbitrary elevation
ref_altitude = 1500.;
......@@ -148,7 +148,7 @@ void Interpol2D::getNeighbors(const double& x, const double& y,
std::vector< std::pair<double, unsigned int> >& list)
{
if(list.size()>0) list.clear();
for(unsigned int i=0; i<vecStations.size(); i++) {
const Coords& position = vecStations[i].position;
const double DX = x-position.getEasting();
......@@ -246,7 +246,7 @@ int Interpol2D::BiLinRegression(const std::vector<double>& in_X, const std::vect
//Now, the core interpolation functions: they project a given parameter to a reference altitude, given a constant lapse rate
//example: Ta projected to 1500m with a rate of -0.0065K/m
/**
* @brief Projects a given parameter to another elevation:
* @brief Projects a given parameter to another elevation:
* This implementation keeps the value constant as a function of the elevation.
* This interface has to follow the interface of *LapseRateProjectPtr
* @param value original value
......@@ -261,7 +261,7 @@ double Interpol2D::ConstProject(const double& value, const double&, const double
}
/**
* @brief Projects a given parameter to another elevation:
* @brief Projects a given parameter to another elevation:
* This implementation assumes a linear dependency of the value as a function of the elevation.
* This interface has to follow the interface of *LapseRateProjectPtr
* @param value original value
......@@ -280,7 +280,7 @@ double Interpol2D::LinProject(const double& value, const double& altitude, const
}
/**
* @brief Projects a given parameter to another elevation:
* @brief Projects a given parameter to another elevation:
* This implementation assumes a 2 segments linear dependency of the value as a function of the elevation.
* It uses Interpol2D::bilin_inflection as the inflection point altitude. This interface has to follow the interface of *LapseRateProjectPtr
* @param value original value
......@@ -302,7 +302,7 @@ double Interpol2D::BiLinProject(const double& value, const double& altitude, con
}
/**
* @brief Projects a given parameter to another elevation:
* @brief Projects a given parameter to another elevation:
* This implementation assumes that the value increases by a given fraction as a function of the elevation.
* This interface has to follow the interface of *LapseRateProjectPtr
* @param value original value
......@@ -322,7 +322,7 @@ double Interpol2D::FracProject(const double& value, const double& altitude, cons
//Filling Functions
/**
* @brief Grid filling function:
* @brief Grid filling function:
* This implementation builds a standard air pressure as a function of the elevation
* @param dem array of elevations (dem)
* @param grid 2D array to fill
......@@ -367,7 +367,7 @@ void Interpol2D::constantGrid2DFill(const double& value, const DEMObject& dem, G
}
/**
* @brief Grid filling function:
* @brief Grid filling function:
* This implementation fills a flat grid with a constant value and then reprojects it to the terrain's elevation.
* for example, the air temperature measured at one point at 1500m would be given as value, the 1500m as altitude and the dem would allow to reproject this temperature on the full DEM using the detrending function provided as pointer (with its previously calculated coefficients).
* @param value value to put in the grid
......@@ -377,7 +377,7 @@ void Interpol2D::constantGrid2DFill(const double& value, const DEMObject& dem, G
* @param funcptr detrending function pointer (that uses the detrending coefficients)
* @param grid 2D array to fill
*/
void Interpol2D::constantLapseGrid2DFill(const double& value, const double& altitude,
void Interpol2D::constantLapseGrid2DFill(const double& value, const double& altitude,
const DEMObject& dem, const std::vector<double>& vecCoefficients,
const LapseRateProjectPtr& funcptr, Grid2DObject& grid)
{
......@@ -401,11 +401,11 @@ double Interpol2D::IDWCore(const double& x, const double& y, const std::vector<d
const std::vector<StationData>& vecStations_in)
{
//The value at any given cell is the sum of the weighted contribution from each source
const unsigned int n_stations=vecStations_in.size();
const size_t n_stations=vecStations_in.size();
double parameter=0., norm=0.;
const double scale = 1.e6;
for (unsigned int i=0; i<n_stations; i++) {
for (size_t i=0; i<n_stations; i++) {
/*const double weight=1./(HorizontalDistance(x, y, vecStations_in[i].position.getEasting(),
vecStations_in[i].position.getNorthing()) + 1e-6);*/
const Coords& position = vecStations_in[i].position;
......@@ -420,7 +420,7 @@ double Interpol2D::IDWCore(const double& x, const double& y, const std::vector<d
}
/**
* @brief Grid filling function:
* @brief Grid filling function:
* This implementation fills a flat grid using Inverse Distance Weighting and then reproject it to the terrain's elevation.
* for example, the air temperatures measured at several stations would be given as values, the stations altitude and positions
* as positions and projected to a flat grid. Afterward, the grid would be reprojected to the correct elevation as given
......@@ -438,13 +438,13 @@ void Interpol2D::LapseIDW(const std::vector<double>& vecData_in, const std::vect
Grid2DObject& grid)
{ //multiple source stations: lapse rate projection, IDW Krieging, re-projection
const double ref_altitude = getReferenceAltitude(dem);
const unsigned int n_stations=vecStations_in.size();
const size_t n_stations=vecStations_in.size();
grid.set(dem.ncols, dem.nrows, dem.cellsize, dem.llcorner);
std::vector<double> vecTref(vecStations_in.size(), 0.0); // init to 0.0
for (unsigned int i=0; i<n_stations; i++) {
vecTref[i] = funcptr(vecData_in[i], vecStations_in[i].position.getAltitude(),
for (size_t i=0; i<n_stations; i++) {
vecTref[i] = funcptr(vecData_in[i], vecStations_in[i].position.getAltitude(),
ref_altitude, vecCoefficients);
}
......@@ -467,9 +467,9 @@ void Interpol2D::LapseIDW(const std::vector<double>& vecData_in, const std::vect
}
/**
* @brief Grid filling function:
* @brief Grid filling function:
* Similar to Interpol2D::LapseIDW but using a limited number of stations for each cell. We also assume a two segments regression for altitude detrending with
* a fixed 1200m above sea level inflection point.
* a fixed 1200m above sea level inflection point.
* @param vecData_in input values to use for the IDW