WSL/SLF GitLab Repository

IOUtils.cc 23.1 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
/***********************************************************************************/
/*  Copyright 2009 WSL Institute for Snow and Avalanche Research    SLF-DAVOS      */
/***********************************************************************************/
/* This file is part of MeteoIO.
    MeteoIO is free software: you can redistribute it and/or modify
    it under the terms of the GNU Lesser General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    MeteoIO is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU Lesser General Public License for more details.

    You should have received a copy of the GNU Lesser General Public License
    along with MeteoIO.  If not, see <http://www.gnu.org/licenses/>.
*/
18
#include <cmath>
19
#include <cstring>
20
#include <ctype.h>
21
#include <algorithm>
22

23
#include <meteoio/IOUtils.h>
24
#include <meteoio/MathOptim.h>
25
#include <meteoio/Config.h>    // to avoid forward declaration hell
26
#include <meteoio/dataClasses/MeteoData.h> // to avoid forward declaration hell
27
#include <meteoio/dataClasses/CoordsAlgorithms.h>
28
29
30

namespace mio {

31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
#ifdef _MSC_VER
//This is C99, Microsoft should move on and suppport it, it is almost 15 years old!!
double round(const double& x) {
	//middle value point test
	if (ceil(x+0.5) == floor(x+0.5)) {
		const int a = (int)ceil(x);
		if (a%2 == 0) {
			return ceil(x);
		} else {
			return floor(x);
		}
	} else {
		return floor(x+0.5);
	}
}
#endif

48
std::string getLibVersion() {
49
	std::ostringstream ss;
50
51
52
53
	ss << _VERSION << " compiled on " << __DATE__ << " " << __TIME__;
	return ss.str();
}

54
55
56
namespace IOUtils {

double bearing_to_angle(const double& bearing) {
57
	return (fmod(360.-bearing+90., 360.)*Cst::to_rad);
58
59
}

60
double angle_to_bearing(const double& angle) {
61
	return (fmod( 90.-angle*Cst::to_deg+360. , 360. ));
62
63
}

64
double bearing(std::string bearing_str)
65
{
66
67
68
	trim(bearing_str);
	toUpper(bearing_str);

69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
	if (bearing_str=="N") return 0.;
	if (bearing_str=="NNE") return 22.5;
	if (bearing_str=="NE") return 45.;
	if (bearing_str=="ENE") return 67.5;
	if (bearing_str=="E") return 90.;
	if (bearing_str=="ESE") return 112.5;
	if (bearing_str=="SE") return 135.;
	if (bearing_str=="SSE") return 157.5;
	if (bearing_str=="S") return 180.;
	if (bearing_str=="SSW") return 202.5;
	if (bearing_str=="SW") return 225.;
	if (bearing_str=="WSW") return 247.5;
	if (bearing_str=="W") return 270.;
	if (bearing_str=="WNW") return 292.5;
	if (bearing_str=="NW") return 315.;
	if (bearing_str=="NNW") return 337.5;
85
86
87
88
89

	//no match
	return nodata;
}

90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
std::string bearing(double bearing)
{
	if (bearing==nodata) return "n/a";

	bearing = fmod( fmod(bearing, 360.)+360., 360.); //from -infty to +infty back to [0, 360]

	if (bearing<=11.25 || bearing>348.75) return "N";
	if (bearing<=33.75) return "NNE";
	if (bearing<=56.25) return "NE";
	if (bearing<=78.75) return "ENE";
	if (bearing<=101.25) return "E";
	if (bearing<=123.75) return "ESE";
	if (bearing<=146.25) return "SE";
	if (bearing<=168.75) return "SSE";
	if (bearing<=191.25) return "S";
	if (bearing<=213.75) return "SSW";
	if (bearing<=236.25) return "SW";
	if (bearing<=258.75) return "WSW";
	if (bearing<=281.25) return "W";
	if (bearing<=303.75) return "WNW";
	if (bearing<=326.25) return "NW";
	if (bearing<=348.75) return "NNW";

	//should not be reached
	return "";
}

117
void stripComments(std::string& str)
118
{
119
	const size_t found = str.find_first_of("#;");
120
121
122
123
124
	if (found != std::string::npos){
		str.erase(found); //rest of line disregarded
	}
}

125
void trim(std::string& str)
126
{
127
128
129
	const std::string whitespaces(" \t\f\v\n\r");
	const size_t startpos = str.find_first_not_of(whitespaces); // Find the first character position after excluding leading blank spaces
	const size_t endpos = str.find_last_not_of(whitespaces); // Find the first character position from reverse af
130
131

	// if all spaces or empty return an empty string
132
	if(startpos!=std::string::npos && endpos!=std::string::npos) {
133
134
		str.erase(endpos+1); //right trim
		str.erase(0, startpos); //left trim
135
	} else {
136
		str.clear();
137
138
139
	}
}

140
void toUpper(std::string& str) {
141
	std::transform(str.begin(), str.end(), str.begin(), ::toupper);
142
143
}

144
void toLower(std::string& str) {
145
	std::transform(str.begin(), str.end(), str.begin(), ::tolower);
146
147
}

148
149
150
151
152
std::string strToUpper(std::string str) {
	//based on http://cpp-next.com/archive/2009/08/want-speed-pass-by-value/
	//it is better to let the compiler copy (or not!) the parameters
	std::transform(str.begin(), str.end(), str.begin(), ::toupper);
	return str;
153
154
}

155
156
157
std::string strToLower(std::string str) {
	std::transform(str.begin(), str.end(), str.begin(), ::tolower);
	return str;
158
159
}

160
bool isNumeric(std::string str, const unsigned int& nBase)
161
{
162
163
	trim(str); //delete trailing and leading whitespaces and tabs
	std::istringstream iss(str);
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179

	if( nBase == 10 ) {
		double tmp;
		iss >> tmp;
	} else if( nBase == 8 || nBase == 16 ) {
		int tmp;
		iss >> ( ( nBase == 8 ) ? std::oct : std::hex ) >> tmp;
	} else
		return false;

	if( !iss ) //the conversion failed
		return false;

	return ( iss.rdbuf()->in_avail() == 0 ); //true if nothing was left after conversion
}

180
bool readKeyValuePair(const std::string& in_line, const std::string& delimiter, std::string &key, std::string &value, const bool& setToUpperCase)
181
{
182
	const size_t pos = ((delimiter==" ") || (delimiter=="\t"))? in_line.find_first_of(" \t", 0) : in_line.find(delimiter); //first occurence of delimiter
183
184

	if(pos != std::string::npos) { //ignore in_lines that are empty or without '='
185
186
		key = in_line.substr(0, pos);
		value = in_line.substr(pos + 1);
187

188
189
		trim(key);
		trim(value);
190

191
		if (key.empty() || value.empty()) {
192
193
194
			return false;
		}

195
		if (setToUpperCase)
196
			toUpper(key);
197
	} else {
198
199
		key="";
		value="";
200
201
202
203
204
205
		return false;
	}

	return true;
}

206
std::string getLogName() {
207
208
209
210
211
212
213
214
	char *tmp;

	if((tmp=getenv("USERNAME"))==NULL) { //Windows & Unix
		if((tmp=getenv("LOGNAME"))==NULL) { //Unix
			tmp=getenv("USER"); //Windows & Unix
		}
	}

215
216
	if(tmp==NULL) return std::string("N/A");
	return std::string(tmp);
217
218
}

219
void readKeyValueHeader(std::map<std::string,std::string>& headermap,
220
221
                        std::istream& fin, const size_t& linecount,
                        const std::string& delimiter, const bool& keep_case)
222
{
223
	size_t linenr = 0;
224
	std::string line;
225
226

	//make a test for end of line encoding:
227
	const char eol = getEoln(fin);
228

229
	for (size_t ii=0; ii< linecount; ii++){
230
		if (std::getline(fin, line, eol)) {
231
			std::string key, value;
232
			linenr++;
233
234
			const bool result = readKeyValuePair(line, delimiter, key, value);
			if(result) {
235
236
				if(!keep_case) headermap[ strToLower(key) ] = value;
				else headermap[key] = value;
237
			} else { //  means if ((key == "") || (value==""))
238
				std::ostringstream out;
239
240
241
242
243
244
245
246
247
				out << "Invalid key value pair in line: " << linenr << " of header";
				throw IOException(out.str(), AT);
			}
		} else {
			throw InvalidFormatException("Premature EOF while reading Header", AT);
		}
	}
}

248
size_t readLineToVec(const std::string& line_in, std::vector<double>& vec_data)
249
250
251
252
253
254
255
256
257
258
259
{
	vec_data.clear();
	std::istringstream iss(line_in); //construct inputstream with the string line as input
	iss.setf(std::ios::fixed);
	iss.precision(std::numeric_limits<double>::digits10);

	double tmp;
	while (!iss.eof()) {
		iss >> std::skipws >> tmp;

		if (iss.fail()) {
260
			std::ostringstream ss;
261
262
263
264
265
266
267
268
269
			ss << "Can not read column " << vec_data.size()+1 << " in data line \"" << line_in << "\"";
			throw InvalidFormatException(ss.str(), AT);
		}
		vec_data.push_back(tmp);
	}

	return vec_data.size();
}

270
size_t readLineToVec(const std::string& line_in, std::vector<std::string>& vecString)
271
272
273
274
275
276
277
278
{
	vecString.clear();
	std::istringstream iss(line_in); //construct inputstream with the string line as input

	std::string tmp_string;
	while (!iss.eof()) {
		iss >> std::skipws >> tmp_string;

279
		if (!tmp_string.empty()) {
280
			vecString.push_back(tmp_string);
281
			tmp_string.clear();
282
283
284
285
286
287
		}
	}

	return vecString.size();
}

288
size_t readLineToVec(const std::string& line_in, std::vector<std::string>& vecString, const char& delim)
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
{
	vecString.clear();
	std::string tmp_string;
	std::istringstream iss(line_in);

	while (getline(iss, tmp_string, delim)){
		vecString.push_back(tmp_string);
	}

	return vecString.size();
}

// generic template function convertString must be defined in the header

const char ALPHANUM[] = "abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ0123456789";
304
const char NUM[] = "0123456789";
305

306
template<> bool convertString<std::string>(std::string& t, const std::string& str, std::ios_base& (*f)(std::ios_base&))
307
308
{
	(void)f;
309
310
	t = str;
	trim(t); //delete trailing and leading whitespaces and tabs
311
312
313
	return true;
}

314
template<> bool convertString<bool>(bool& t, const std::string& str, std::ios_base& (*f)(std::ios_base&))
315
{
316
	std::string s(str);
317
318
	trim(s); //delete trailing and leading whitespaces and tabs

319
	if (toupper(s[0])=='T' || toupper(s[0])=='Y') {
320
		t = true;
321
	} else if (toupper(s[0])=='F' || toupper(s[0])=='N') {
322
323
324
325
326
327
328
329
330
331
332
		t = false;
	} else {
		std::istringstream iss(s);
		int i;
		iss >> f >> i; //Convert first part of stream with the formatter (e.g. std::dec, std::oct)
		if (iss.fail()) {//Conversion failed
			return false;
		}
		t = (i != 0);
	}

333
	const std::string::size_type pos = s.find_first_not_of(ALPHANUM);
334
335
336
	if (pos != std::string::npos) {
		std::string tmp = s.substr(pos);
		trim(tmp);
337
		if (!tmp.empty() && tmp[0] != '#' && tmp[0] != ';') {//if line holds more than one value it's invalid
338
339
340
341
342
343
344
			return false;
		}
	}

	return true;
}

345
346
347
348
349
template<> bool convertString<double>(double& t, const std::string& str, std::ios_base& (*f)(std::ios_base&))
{
	if (f == std::dec) {
		//First check if string is empty
		const char* start = str.c_str();
350
		while(*start && isspace(*start)) start++;
351
352
353
354
355
356
357
358
359
360
361
362
		if (*start == '\0' || *start == '#' || *start == ';') { // line empty or comment
			t = static_cast<double> (nodata);
			return true;
		}

		//string is not empty
		char* end;
		t = strtod(str.c_str(), &end); //would strip leading whitespaces, but already done

		if (*end == '\0') { //conversion successful
			return true;
		} else { // conversion might have worked, let's check what is left
363
			while((*end != '\0') && isspace(*end)) end++;
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398

			if (*end == '\0' || *end == '#' || *end == ';') { // we allow the number to be followed by a comment
				return true;
			}

			return false; // Invalid string to convert to double
		}
	}

	std::string s(str);
	trim(s); //delete trailing and leading whitespaces and tabs
	if (s.empty()) {
		t = static_cast<double> (nodata);
		return true;
	}

	std::istringstream iss(s);
	iss.setf(std::ios::fixed);
	iss.precision(std::numeric_limits<double>::digits10); //try to read values with maximum precision
	iss >> f >> t; //Convert first part of stream with the formatter (e.g. std::dec, std::oct)

	if (iss.fail()) {
		//Conversion failed
		return false;
	}
	std::string tmp;
	getline(iss,  tmp); //get rest of line, if any
	trim(tmp);
	if (!tmp.empty() && tmp[0] != '#' && tmp[0] != ';') {
		//if line holds more than one value it's invalid
		return false;
	}
	return true;
}

399
template<> bool convertString<unsigned int>(unsigned int& t, const std::string& str, std::ios_base& (*f)(std::ios_base&))
400
{
401
	std::string s(str);
402
	trim(s); //delete trailing and leading whitespaces and tabs
403
	if (s.empty()) {
404
405
406
407
408
409
410
411
412
413
414
		t = unodata;
		return true;
	} else {
		std::istringstream iss(s);
		iss.setf(std::ios::fixed);
		iss.precision(std::numeric_limits<double>::digits10); //try to read values with maximum precision
		iss >> f >> t; //Convert first part of stream with the formatter (e.g. std::dec, std::oct)
		if (iss.fail()) {
			//Conversion failed
			return false;
		}
415
		std::string tmp;
416
417
		getline(iss,  tmp); //get rest of line, if any
		trim(tmp);
418
		if (!tmp.empty() && tmp[0] != '#' && tmp[0] != ';') {
419
420
421
422
423
424
425
			//if line holds more than one value it's invalid
			return false;
		}
		return true;
	}
}

426
bool convertString(Date& t, const std::string& str, const double& time_zone, std::ios_base& (*f)(std::ios_base&))
427
{
428
	std::string s(str);
429
430
431
	trim(s); //delete trailing and leading whitespaces and tabs

	(void)f;
432
	int year;
433
434
	unsigned int month, day, hour, minute;
	double second;
435
	char rest[32] = "";
436

437
438
439
440
441
442
	const char *c_str = s.c_str();
	if (sscanf(c_str, "%d-%u-%u %u:%u:%lg%31s", &year, &month, &day, &hour, &minute, &second, rest) >= 6) {
		std::string timezone_iso(rest);
		stripComments(timezone_iso);
		const double tz = (timezone_iso.empty())? time_zone : Date::parseTimeZone(timezone_iso);
		if(tz==nodata) return false;
443
		t.setDate(year, month, day, hour, minute, second, tz);
444
445
446
447
448
449
450
		return true;

	} else if (sscanf(c_str, "%d-%u-%uT%u:%u:%lg%31s", &year, &month, &day, &hour, &minute, &second, rest) >= 6) { //ISO
		std::string timezone_iso(rest);
		stripComments(timezone_iso);
		const double tz = (timezone_iso.empty())? time_zone : Date::parseTimeZone(timezone_iso);
		if(tz==nodata) return false;
451
		t.setDate(year, month, day, hour, minute, second, tz);
452
453
454
455
456
457
458
		return true;

	} else if (sscanf(c_str, "%d-%u-%u %u:%u%31s", &year, &month, &day, &hour, &minute, rest) >= 5) {
		std::string timezone_iso(rest);
		stripComments(timezone_iso);
		const double tz = (timezone_iso.empty())? time_zone : Date::parseTimeZone(timezone_iso);
		if(tz==nodata) return false;
Mathias Bavay's avatar
Mathias Bavay committed
459
		t.setDate(year, month, day, hour, minute, static_cast<unsigned>(0), tz);
460
461
462
463
464
465
466
		return true;

	} else if (sscanf(c_str, "%d-%u-%uT%u:%u%31s", &year, &month, &day, &hour, &minute, rest) >= 5) {
		std::string timezone_iso(rest);
		stripComments(timezone_iso);
		const double tz = (timezone_iso.empty())? time_zone : Date::parseTimeZone(timezone_iso);
		if(tz==nodata) return false;
Mathias Bavay's avatar
Mathias Bavay committed
467
		t.setDate(year, month, day, hour, minute, static_cast<unsigned>(0), tz);
468
469
470
471
472
473
474
		return true;

	} else if (sscanf(c_str, "%d-%u-%u%31s", &year, &month, &day, rest) >= 3) {
		std::string timezone_iso(rest);
		stripComments(timezone_iso);
		const double tz = (timezone_iso.empty())? time_zone : Date::parseTimeZone(timezone_iso);
		if(tz==nodata) return false;
Mathias Bavay's avatar
Mathias Bavay committed
475
		t.setDate(year, month, day, static_cast<unsigned>(0), static_cast<unsigned>(0), static_cast<unsigned>(0), tz);
476
477
478
479
480
481
482
		return true;

	} else if (sscanf(c_str, "%u:%u%31s", &hour, &minute, rest) >= 2) {
		std::string timezone_iso(rest);
		stripComments(timezone_iso);
		const double tz = (timezone_iso.empty())? time_zone : Date::parseTimeZone(timezone_iso);
		if(tz==nodata) return false;
Mathias Bavay's avatar
Mathias Bavay committed
483
		t.setDate( (static_cast<double>(hour))/24. + (static_cast<double>(minute))/24./60. , tz);
484
485
		return true;

486
	} else {
487
		//try to read purely numerical date, potentially surrounded by other chars
488
489
490
491
492
		//and potentially containing an ISO time zone string
		const size_t in_len = s.length();

		//extract date/time
		const size_t date_beg = s.find_first_of(NUM);
493
		if (date_beg==npos || date_beg==in_len) return false;
494
		size_t date_end = s.find_first_not_of(NUM, date_beg+1);
495
		if (date_end==npos) date_end = in_len;
496
497
498
499
		const std::string date = s.substr(date_beg, date_end-date_beg);

		//parse date/time
		const size_t date_len = date.length();
500
501
502
503
504
505
		if (date_len<10 || date_len>14) return false;
		if (convertString(year,date.substr(0,4))==false) return false;
		if (convertString(month,date.substr(4,2))==false) return false;
		if (convertString(day,date.substr(6,2))==false) return false;
		if (convertString(hour,date.substr(8,2))==false) return false;
		if (date_len==10)
506
507
			minute=0;
		else {
508
			if (date_len>=12) {
509
				if( convertString(minute,date.substr(10,2))==false ) return false;
510
511
			} else
				return false;
512
			if (date_len==12)
513
514
				second=0;
			else {
515
516
				if (date_len==14) {
					if (convertString(second,date.substr(12,2))==false) return false;
517
518
519
520
				} else
					return false;
			}
		}
521

522
523
524
		//extract potential ISO time zone string
		double tz = time_zone;
		const size_t tz_beg = s.find_first_of("+-", date_end);
525
		if (tz_beg!=npos && tz_beg!=in_len) {
526
			size_t tz_end = s.find_first_not_of("0123456789:", date_end+1);
527
			if (tz_end==npos) tz_end = in_len;
528
529
530
			const std::string timezone_iso = s.substr(tz_beg, tz_end-tz_beg);
			if(!timezone_iso.empty()) tz = Date::parseTimeZone(timezone_iso);
		}
531

532
		t.setDate( year, month, day, hour, minute, second, tz );
533
534
535
536
537
	}

	return true;
}

538
template<> bool convertString<Coords>(Coords& t, const std::string& str, std::ios_base& (*f)(std::ios_base&))
539
{
540
	std::string s(str);
541
542
543
544
545
	trim(s); //delete trailing and leading whitespaces and tabs

	(void)f;
	double lat, lon;
	try {
546
		CoordsAlgorithms::parseLatLon(s, lat, lon);
547
	} catch(const IOException&) {
548
549
550
551
552
553
554
555
		return false;
	}
	t.setLatLon(lat, lon, nodata);

	return true;
}


556
void getProjectionParameters(const Config& cfg, std::string& coordin, std::string& coordinparam,
557
                                      std::string& coordout, std::string& coordoutparam) {
558
	cfg.getValue("COORDSYS", "Input", coordin);
559
560
561
	cfg.getValue("COORDPARAM", "Input", coordinparam, IOUtils::nothrow);
	cfg.getValue("COORDSYS", "Output", coordout, IOUtils::nothrow);
	cfg.getValue("COORDPARAM", "Output", coordoutparam, IOUtils::nothrow);
562
563
}

564
void getTimeZoneParameters(const Config& cfg, double& tz_in, double& tz_out) {
565
566
	cfg.getValue("TIME_ZONE", "Input", tz_in, IOUtils::nothrow);
	cfg.getValue("TIME_ZONE", "Output", tz_out, IOUtils::nothrow);
567
}
Cyril Perot's avatar
Cyril Perot committed
568

569
570
//returns index of element, if element does not exist it returns closest index after soughtdate
//the element needs to be an exact hit or embedded between two measurments
571
size_t seek(const Date& soughtdate, const std::vector<MeteoData>& vecM, const bool& exactmatch)
572
{
573
574
	if (vecM.empty() || soughtdate > vecM.back().date || soughtdate < vecM.front().date) {
		//the sought date is not contained in the vector, return npos
575
		return npos;
576
577
	}

578
	const size_t max_idx = vecM.size()-1; //obviously, the index must be <= max_idx
579

580
581
582
583
584
585
	//since usually the sampling rate is quite constant, try to guess where our point
	//should be and provide a much smaller search interval around it
	const double start_date = vecM.front().date.getJulian(true);
	const double end_date = vecM.back().date.getJulian(true);
	const double curr_date = soughtdate.getJulian(true);
	const double raw_pos = (curr_date-start_date) / (end_date-start_date) * static_cast<double>(max_idx); //always >=0
Mathias Bavay's avatar
Mathias Bavay committed
586
587
	const size_t start_idx = static_cast<size_t>( floor(raw_pos*.9) );
	const size_t end_idx = std::min( static_cast<size_t>( ceil(raw_pos*1.1) ), max_idx);
588

589
590
591
	//first and last index of the search interval, either using our initial guess or the full vector
	size_t first = (curr_date >= vecM[start_idx].date.getJulian(true))? start_idx : 0;
	size_t last = (curr_date <= vecM[end_idx].date.getJulian(true))? end_idx : max_idx;
592

593
594
	//if we reach this point: the date is spanned by the buffer and there are at least two elements
	if (exactmatch){
595
596
		//perform binary search
		while (first <= last) {
597
			const size_t mid = (first + last) / 2;  // compute mid point
Cyril Perot's avatar
Cyril Perot committed
598
			if (soughtdate > vecM[mid].date)
599
				first = mid + 1;                   // repeat search in top half
Cyril Perot's avatar
Cyril Perot committed
600
			else if (soughtdate < vecM[mid].date)
601
				last = mid - 1;                    // repeat search in bottom half
Cyril Perot's avatar
Cyril Perot committed
602
603
			else
				return mid;                        // found it. return position
604
605
		}
	} else {
606
607
		//perform binary search
		while (first <= last) {
608
			const size_t mid = (first + last) / 2;  // compute mid point
Cyril Perot's avatar
Cyril Perot committed
609

610
			if (mid < max_idx) {
611
612
				if ((soughtdate > vecM[mid].date) && (soughtdate < vecM[mid+1].date))
					return mid+1;
613
			}
614

Cyril Perot's avatar
Cyril Perot committed
615
			if (soughtdate > vecM[mid].date)
616
				first = mid + 1;                   // repeat search in top half
Cyril Perot's avatar
Cyril Perot committed
617
			else if (soughtdate < vecM[mid].date)
618
619
				last = mid - 1;                    // repeat search in bottom half
			else
Cyril Perot's avatar
Cyril Perot committed
620
				return mid;                        // found it. return position
621
622
623
		}
	}

624
	return npos;
625
626
}

627
void getArraySliceParams(const size_t& dimx, const size_t& nbworkers, const size_t &wk, size_t& startx, size_t& nx)
628
629
{
	if(nbworkers>dimx) {
630
		std::ostringstream ss;
631
632
633
634
		ss << "Can not split " << dimx << " columns in " << nbworkers << " bands!";
		throw InvalidArgumentException(ss.str(), AT);
	}

635
636
	const size_t nx_avg = dimx / nbworkers;
	const size_t remainder = dimx % nbworkers;
637
638
639
640
641
642
643
644
645
646

	if(wk<=remainder) { //distribute remainder, 1 extra column per worker, on first workers
		nx = nx_avg+1;
		startx = (wk-1)*nx;
	} else { //all remainder has been distributed, we now attribute a normal number of columns
		nx = nx_avg;
		startx = remainder*(nx+1) + (wk-1-remainder)*nx;
	}
}

647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
double unitsPrefix(const char& prefix)
{
	if (prefix == 'f') {
		return 1e-15;
	} else if (prefix == 'p') {
		return 1e-12;
	} else if (prefix == 'n') {
		return 1e-9;
	} else if (prefix == 'u') {
		return 1e-6;
	} else if (prefix == 'm') {
		return 1e-3;
	} else if (prefix == 'c') {
		return 1e-2;
	} else if (prefix == 'd') {
		return 1e-1;
	} else if (prefix == 'h') {
		return 1e2;
	} else if (prefix == 'k') {
		return 1e3;
	} else if (prefix == 'M') {
		return 1e6;
	} else if (prefix == 'G') {
		return 1e9;
	} else if (prefix == 'T') {
		return 1e12;
	} else if (prefix == 'P') {
		return 1e15;
	}

	const std::string prefix_str( 1, prefix );
	throw IOException("Invalid unit prefix '"+prefix_str+"'", AT);
}

681
682
//currently, only the most simple ase of units are handled. Composite units
//such as 'W/m2 <-> mW/cm2' are NOT handled.
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
double unitsConversion(const double& val, std::string unitIn, std::string unitOut)
{
	if (val==IOUtils::nodata)
		return IOUtils::nodata;
	if (unitIn.empty() || unitOut.empty())
		throw InvalidArgumentException("Units can not be empty!", AT);

	if (unitIn=="degK" || unitIn=="°K" || unitIn=="Kelvin")
		unitIn = "K";
	if (unitOut=="degK" || unitOut=="°K" || unitOut=="Kelvin")
		unitOut = "K";
	if (unitIn=="degC" || unitIn=="Celsius")
		unitIn = "°C";
	if (unitOut=="degC" || unitOut=="Celsius")
		unitOut = "°C";
	if (unitIn=="degF" || unitIn=="Fahrenheit")
		unitIn = "°F";
	if (unitOut=="degF" || unitOut=="Fahrenheit")
		unitOut = "°F";

	if (unitIn=="°C" && unitOut=="K") {
		return (val+Cst::t_water_triple_pt);
	} else if (unitIn=="K" && unitOut=="°C") {
		return (val-Cst::t_water_triple_pt);
	} else if (unitIn=="K" && unitOut=="°F") {
		return ((val-Cst::t_water_triple_pt)*1.8+32.);
	} else if (unitIn=="°F" && unitOut=="K") {
		return ((val-32.)/1.8+Cst::t_water_triple_pt);
	}  else if (unitIn=="°F" && unitOut=="°C") {
		return ((val-32.)/1.8);
	}  else if (unitIn=="°C" && unitOut=="°F") {
		return (val*1.8+32.);
	} else {
716
717
718
		//extract the unit prefix
		const double inPrefix_factor = (isalpha(unitIn[0]) && isalpha(unitIn[1]))? unitsPrefix( unitIn[0] ) : 1;
		const double outPrefix_factor = (isalpha(unitOut[0]) && isalpha(unitOut[1]))? unitsPrefix( unitOut[0] ) : 1;
719

720
721
722
723
724
725
726
727
728
729
730
		//extract the unit exponent
		const char in_last_char = unitIn[ unitIn.size()-1 ];
		const char out_last_char = unitOut[ unitOut.size()-1 ];
		const unsigned char inExponent = (isdigit(in_last_char))? static_cast<unsigned char>( in_last_char-'0' ) : static_cast<unsigned char>( 1 );
		const unsigned char outExponent = (isdigit(out_last_char))? static_cast<unsigned char>( out_last_char-'0' ) : static_cast<unsigned char>( 1 );

		//compute the input and output units factor
		const double inFactor = (inExponent==1)? inPrefix_factor : Optim::fastPow(inPrefix_factor, inExponent);
		const double outFactor = (outExponent==1)? outPrefix_factor : Optim::fastPow(outPrefix_factor, outExponent);

		const double ratio = inFactor / outFactor;
731
732
		return val*ratio;
	}
733
	//throw ConversionFailedException("Unable to perform unit conversion.", AT);
734
735
}

736
} //namespace IOUtils
737
} //namespace