WSL/SLF GitLab Repository

IOUtils.cc 22.9 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
28
29

namespace mio {

30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
#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

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

53
54
55
namespace IOUtils {

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

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

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

68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
	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;
84
85
86
87
88

	//no match
	return nodata;
}

89
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
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 "";
}

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

124
void trim(std::string& str)
125
{
126
127
128
	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
129
130

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

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

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

147
148
149
150
151
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;
152
153
}

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

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

	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
}

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

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

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

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

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

	return true;
}

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

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

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

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

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

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

247
size_t readLineToVec(const std::string& line_in, std::vector<double>& vec_data)
248
249
250
251
252
253
254
255
256
257
258
{
	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()) {
259
			std::ostringstream ss;
260
261
262
263
264
265
266
267
268
			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();
}

269
size_t readLineToVec(const std::string& line_in, std::vector<std::string>& vecString)
270
271
272
273
274
275
276
277
{
	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;

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

	return vecString.size();
}

287
size_t readLineToVec(const std::string& line_in, std::vector<std::string>& vecString, const char& delim)
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
{
	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";
303
const char NUM[] = "0123456789";
304

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

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

318
	if (toupper(s[0])=='T' || toupper(s[0])=='Y') {
319
		t = true;
320
	} else if (toupper(s[0])=='F' || toupper(s[0])=='N') {
321
322
323
324
325
326
327
328
329
330
331
		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);
	}

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

	return true;
}

344
345
346
347
348
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();
349
		while(*start && isspace(*start)) start++;
350
351
352
353
354
355
356
357
358
359
360
361
		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
362
			while((*end != '\0') && isspace(*end)) end++;
363
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

			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;
}

398
template<> bool convertString<unsigned int>(unsigned int& t, const std::string& str, std::ios_base& (*f)(std::ios_base&))
399
{
400
	std::string s(str);
401
	trim(s); //delete trailing and leading whitespaces and tabs
402
	if (s.empty()) {
403
404
405
406
407
408
409
410
411
412
413
		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;
		}
414
		std::string tmp;
415
416
		getline(iss,  tmp); //get rest of line, if any
		trim(tmp);
417
		if (!tmp.empty() && tmp[0] != '#' && tmp[0] != ';') {
418
419
420
421
422
423
424
			//if line holds more than one value it's invalid
			return false;
		}
		return true;
	}
}

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

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

436
437
438
439
440
441
	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;
442
		t.setDate(year, month, day, hour, minute, second, tz);
443
444
445
446
447
448
449
		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;
450
		t.setDate(year, month, day, hour, minute, second, tz);
451
452
453
454
455
456
457
		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;
458
		t.setDate(year, month, day, hour, minute, (unsigned)0, tz);
459
460
461
462
463
464
465
		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;
466
		t.setDate(year, month, day, hour, minute, (unsigned)0, tz);
467
468
469
470
471
472
473
		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;
474
		t.setDate(year, month, day, (unsigned)0, (unsigned)0, (unsigned)0, tz);
475
476
477
478
479
480
481
482
483
484
		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;
		t.setDate( ((double)hour)/24. + ((double)minute)/24./60. , tz);
		return true;

485
	} else {
486
		//try to read purely numerical date, potentially surrounded by other chars
487
488
489
490
491
		//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);
492
		if (date_beg==npos || date_beg==in_len) return false;
493
		size_t date_end = s.find_first_not_of(NUM, date_beg+1);
494
		if (date_end==npos) date_end = in_len;
495
496
497
498
		const std::string date = s.substr(date_beg, date_end-date_beg);

		//parse date/time
		const size_t date_len = date.length();
499
500
501
502
503
504
		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)
505
506
			minute=0;
		else {
507
			if (date_len>=12) {
508
				if( convertString(minute,date.substr(10,2))==false ) return false;
509
510
			} else
				return false;
511
			if (date_len==12)
512
513
				second=0;
			else {
514
515
				if (date_len==14) {
					if (convertString(second,date.substr(12,2))==false) return false;
516
517
518
519
				} else
					return false;
			}
		}
520

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

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

	return true;
}

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

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

	return true;
}


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

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

568
569
//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
570
size_t seek(const Date& soughtdate, const std::vector<MeteoData>& vecM, const bool& exactmatch)
571
{
572
573
	if (vecM.empty() || soughtdate > vecM.back().date || soughtdate < vecM.front().date) {
		//the sought date is not contained in the vector, return npos
574
		return npos;
575
576
	}

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

579
580
581
582
583
584
	//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
585
	const size_t start_idx = (size_t)floor(raw_pos*.9);
586
	const size_t end_idx = std::min( (size_t)ceil(raw_pos*1.1), max_idx);
587

588
589
590
	//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;
591

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

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

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

623
	return npos;
624
625
}

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

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

	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;
	}
}

646
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
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);
}

680
681
//currently, only the most simple ase of units are handled. Composite units
//such as 'W/m2 <-> mW/cm2' are NOT handled.
682
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
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 {
715
716
717
		//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;
718

719
720
721
722
723
724
725
726
727
728
729
		//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;
730
731
732
733
734
		return val*ratio;
	}
	throw ConversionFailedException("Unable to perform unit conversion.", AT);
}

735
} //namespace IOUtils
736
} //namespace