WSL/SLF GitLab Repository

Array1D.h 16.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
19
#ifndef ARRAY1D_H
#define ARRAY1D_H
20
21
22
23

#include <vector>
#include <limits>
#include <iostream>
24
#include <iterator>
25

26
27
28
#include <meteoio/IOUtils.h>
#include <meteoio/IOExceptions.h>

29
30
31
namespace mio {

/**
32
33
 * @class Array1D
 * @brief The template class Array1D is a 1D array (vector) able to hold any type of object as datatype.
34
 * If the compilation flag NOSAFECHECKS is used, bounds check is turned off (leading to increased performances).
35
 *
36
 * @ingroup data_str
37
38
39
40
 * @author Thomas Egger
 * @date   2009-05-02
 */

41
template<class T> class Array1D {
42
	public:
43
		Array1D(const size_t& asize=0);
44

45
46
47
48
49
		/**
		* A constructor that creates an array filled with constant values
		* @param asize size of the new array
		* @param init initial value to fill the array with
		*/
50
		Array1D(const size_t& asize, const T& init);
51

52
53
		/**
		* @brief set how to process nodata values (ie: as nodata or as normal numbers)
54
55
		* @param i_keep_nodata true means that NODATA is interpreted as NODATA, false means that it is a normal number
		By default, arrays keep nodata.
56
		*/
57
58
59
60
61
62
63
		void setKeepNodata(const bool i_keep_nodata);

		/**
		* @brief get how to process nodata values (ie: as nodata or as normal numbers)
		* @return true means that NODATA is interpreted as NODATA, false means that it is a normal number
		*/
		bool getKeepNodata();
64

65
66
		void size(size_t& nx) const;
		size_t getNx() const;
67

68
69
		void resize(const size_t& asize);
		void resize(const size_t& asize, const T& init);
70
		void clear();
71
		bool empty() const;
72

73
		void insertAt(const int& index, T e);
74
		void removeAt(const size_t& index);
75

76
77
78
79
		/**
		* @brief returns the minimum value contained in the grid
		* @return minimum value
		*/
80
		T getMin() const;
81
82
83
84
		/**
		* @brief returns the maximum value contained in the grid
		* @return maximum value
		*/
85
		T getMax() const;
86
87
88
89
		/**
		* @brief returns the mean value contained in the grid
		* @return mean value
		*/
90
		T getMean() const;
91
		/**
92
93
94
		* @brief returns the number of points contained in the grid.
		* If setNodataHandling(IOUtils::RAW_NODATA), then the number of points is the size of the grid.
		* If setNodataHandling(IOUtils::PARSE_NODATA), then it is the number of non-nodata values in the grid
95
		* @return count
96
		*/
97
		size_t getCount() const;
98
99
100
101
		/**
		* @brief returns the grid of the absolute value of values contained in the grid
		* @return grid of abs(grid)
		*/
102
		const Array1D<T> getAbs() const;
103
		void abs();
104

105

106
107
108
		const std::string toString() const;
		template<class P> friend std::iostream& operator<<(std::iostream& os, const Array1D<P>& array);
		template<class P> friend std::iostream& operator>>(std::iostream& is, Array1D<P>& array);
109

110
111
		bool checkEpsilonEquality(const Array1D<double>& rhs, const double& epsilon) const;
		static bool checkEpsilonEquality(const Array1D<double>& rhs1, const Array1D<double>& rhs2, const double& epsilon);
112

113
114
115
116
		T& operator [](const size_t& index);
		const T operator [](const size_t& index) const;
		T& operator ()(const size_t& index);
		const T operator ()(const size_t& index) const;
117

118
119
		Array1D<T>& operator =(const Array1D<T>&);
		Array1D<T>& operator =(const T& value);
120

121
		Array1D<T>& operator+=(const T& rhs);
122
		const Array1D<T> operator+(const T& rhs) const;
123
		Array1D<T>& operator+=(const Array1D<T>& rhs);
124
		const Array1D<T> operator+(const Array1D<T>& rhs) const;
125

126
		Array1D<T>& operator-=(const T& rhs);
127
		const Array1D<T> operator-(const T& rhs) const;
128
		Array1D<T>& operator-=(const Array1D<T>& rhs);
129
		const Array1D<T> operator-(const Array1D<T>& rhs) const;
130

131
		Array1D<T>& operator*=(const T& rhs);
132
		const Array1D<T> operator*(const T& rhs) const;
133
		Array1D<T>& operator*=(const Array1D<T>& rhs);
134
		const Array1D<T> operator*(const Array1D<T>& rhs) const;
135

136
		Array1D<T>& operator/=(const T& rhs);
137
		const Array1D<T> operator/(const T& rhs) const;
138
		Array1D<T>& operator/=(const Array1D<T>& rhs);
139
		const Array1D<T> operator/(const Array1D<T>& rhs) const;
140

141
142
143
		bool operator==(const Array1D<T>&) const; ///<Operator that tests for equality
		bool operator!=(const Array1D<T>&) const; ///<Operator that tests for inequality

144
145
	protected:
		std::vector<T> vecData; ///<the actual data structure, that holds the objects of type T
146
		size_t nx; ///<this is introduced to omit the costly vecData.size()
147
		bool keep_nodata;
148
149
};

150
template<class T> Array1D<T>::Array1D(const size_t& asize)
151
152
153
                  : vecData(asize), nx(asize), keep_nodata(true)
{
	//resize(asize);
154
155
}

156
template<class T> Array1D<T>::Array1D(const size_t& asize, const T& init)
157
158
159
                 : vecData(asize, init), nx(asize), keep_nodata(true)
{
	//resize(asize, init);
160
161
}

162
template<class T> void Array1D<T>::setKeepNodata(const bool i_keep_nodata) {
163
	keep_nodata = i_keep_nodata;
164
165
}

166
template<class T> bool Array1D<T>::getKeepNodata() {
167
168
169
	return keep_nodata;
}

170
template<class T> void Array1D<T>::size(size_t& o_nx) const {
171
172
	o_nx = nx;
}
173

174
template<class T> size_t Array1D<T>::getNx() const {
175
176
177
	return nx;
}

178
template<class T> void Array1D<T>::resize(const size_t& asize) {
179
	vecData.clear();
180
181
	vecData.resize(asize);
	nx = asize;
182
183
}

184
template<class T> void Array1D<T>::resize(const size_t& asize, const T& init) {
185
	vecData.clear();
186
187
	vecData.resize(asize, init);
	nx = asize;
188
189
}

190
template<class T> inline T& Array1D<T>::operator()(const size_t& index) {
191
192
#ifndef NOSAFECHECKS
	if (index >= nx) {
193
194
		std::stringstream ss;
		ss << "Trying to access array(" << index << ")";
195
		ss << " while array is (" << nx << ")";
196
		throw IndexOutOfBoundsException(ss.str(), AT);
197
198
199
200
201
	}
#endif
	return vecData[index];
}

202
template<class T> inline const T Array1D<T>::operator()(const size_t& index) const {
203
204
#ifndef NOSAFECHECKS
	if (index >= nx) {
205
206
		std::stringstream ss;
		ss << "Trying to access array(" << index << ")";
207
		ss << " while array is (" << nx << ")";
208
		throw IndexOutOfBoundsException(ss.str(), AT);
209
210
211
212
213
	}
#endif
	return vecData[index];
}

214
template<class T> inline T& Array1D<T>::operator [](const size_t& index) {
215
#ifndef NOSAFECHECKS
216
217
	return vecData.at(index);
#else
218
	return vecData[index];
219
#endif
220
221
}

222
template<class T> inline const T Array1D<T>::operator [](const size_t& index) const {
223
#ifndef NOSAFECHECKS
224
225
	return vecData.at(index);
#else
226
	return vecData[index];
227
#endif
228
229
}

230
template<class T> void Array1D<T>::clear() {
231
232
233
234
	vecData.clear();
	nx = 0;
}

235
template<class T> bool Array1D<T>::empty() const {
236
237
238
	return (nx==0);
}

239
240
template<class T> const std::string Array1D<T>::toString() const {
	std::stringstream os;
241
	os << "<array1d>\n";
242
	for(size_t ii=0; ii<nx; ii++) {
243
		os << vecData[ii] << " ";
244
245
	}
	os << "\n</array1d>\n";
246
247
248
249
250
251
	return os.str();
}

template<class P> std::iostream& operator<<(std::iostream& os, const Array1D<P>& array) {
	os.write(reinterpret_cast<const char*>(&array.keep_nodata), sizeof(array.keep_nodata));
	os.write(reinterpret_cast<const char*>(&array.nx), sizeof(array.nx));
252
	os.write(reinterpret_cast<const char*>(&array.vecData[0]), static_cast<std::streamsize>(array.nx*sizeof(P)));
253
254
255
	return os;
}

256
257
258
259
template<class P> std::iostream& operator>>(std::iostream& is, Array1D<P>& array) {
	is.read(reinterpret_cast<char*>(&array.keep_nodata), sizeof(array.keep_nodata));
	is.read(reinterpret_cast<char*>(&array.nx), sizeof(array.nx));
	array.vecData.resize(array.nx);
260
	is.read(reinterpret_cast<char*>(&array.vecData[0]), static_cast<std::streamsize>(array.nx*sizeof(P))); //30 times faster than assign() or copy()
261
	return is;
262
263
}

264
template<class T> void Array1D<T>::insertAt(const int& index, T e) {
265
266
267
268
269
270
271
	if (index < 0) {
		vecData.push_back(e);
                nx++;
	} else if ((index >= 0) && (index < (int)vecData.size())) {
		vecData.insert(vecData.begin() + index, e);
		nx++;
	} else {
272
		std::stringstream ss;
273
		ss << "Inserting an element at (" << index << ") in an array of size [" << nx << "]";
274
		throw IndexOutOfBoundsException(ss.str(), AT);
275
276
277
	}
}

278
template<class T> void Array1D<T>::removeAt(const size_t& index) {
279
280
281
282
283
284
	if (index < vecData.size()) {
		vecData.erase(vecData.begin()+index);
		nx--;
	}
}

285
template<class T> T Array1D<T>::getMin() const {
286
287
288

	T min = std::numeric_limits<T>::max();

289
	if(keep_nodata==false) {
290
		for (size_t ii=0; ii<nx; ii++) {
291
292
293
294
			const T val = vecData[ii];
			if(val<min) min=val;
		}
		return min;
295
	} else {
296
		for (size_t ii=0; ii<nx; ii++) {
297
298
299
300
301
302
303
304
			const T val = vecData[ii];
			if(val!=IOUtils::nodata && val<min) min=val;
		}
		if(min!=std::numeric_limits<T>::max()) return min;
		else return (T)IOUtils::nodata;
	}
}

305
template<class T> T Array1D<T>::getMax() const {
306
307
308

	T max = -std::numeric_limits<T>::max();

309
	if(keep_nodata==false) {
310
		for (size_t ii=0; ii<nx; ii++) {
311
312
313
314
			const T val = vecData[ii];
			if(val>max) max=val;
		}
		return max;
315
	} else {
316
		for (size_t ii=0; ii<nx; ii++) {
317
318
319
320
321
322
323
324
			const T val = vecData[ii];
			if(val!=IOUtils::nodata && val>max) max=val;
		}
		if(max!=-std::numeric_limits<T>::max()) return max;
		else return (T)IOUtils::nodata;
	}
}

325
template<class T> T Array1D<T>::getMean() const {
326
327
328

	T mean = 0;

329
	if(keep_nodata==false) {
330
		for (size_t ii=0; ii<nx; ii++) {
331
332
333
			const T val = vecData[ii];
			mean += val;
		}
334
		const size_t count = nx;
335
336
		if(count>0) return mean/(T)(count);
		else return (T)0;
337
	} else {
338
339
		size_t count = 0;
		for (size_t ii=0; ii<nx; ii++) {
340
341
342
343
344
345
346
347
348
349
350
			const T val = vecData[ii];
			if(val!=IOUtils::nodata) {
				mean += val;
				count++;
			}
		}
		if(count>0) return mean/(T)(count);
		else return (T)IOUtils::nodata;
	}
}

351
template<class T> size_t Array1D<T>::getCount() const
352
{
353
	if(keep_nodata==false) {
354
		return (size_t)nx;
355
	} else {
356
		size_t count = 0;
357
		for (size_t ii=0; ii<nx; ii++) {
358
			if(vecData[ii]!=IOUtils::nodata) count++;
359
		}
360
		return count;
361
362
363
	}
}

364
template<class T> void Array1D<T>::abs() {
365
	if(std::numeric_limits<T>::is_signed) {
366
		if(keep_nodata==false) {
367
			for (size_t ii=0; ii<nx; ii++) {
368
				T& val = vecData[ii];
369
370
				if(val<0) val=-val;
			}
371
		} else {
372
			for (size_t ii=0; ii<nx; ii++) {
373
				T& val = vecData[ii];
374
375
376
377
378
379
				if(val<0 && val!=IOUtils::nodata) val=-val;
			}
		}
	}
}

380
template<class T> const Array1D<T> Array1D<T>::getAbs() const {
381
	Array1D<T> result(*this); //make a copy
382
	result.abs(); //already implemented
383
384
385
386

	return result;
}

387
//arithmetic operators
388
template<class T> bool Array1D<T>::checkEpsilonEquality(const Array1D<double>& rhs, const double& epsilon) const {
389
390
	if(nx!=rhs.nx) return false;

391
	for (size_t jj=0; jj<nx; jj++)
392
393
394
395
396
		if(IOUtils::checkEpsilonEquality(vecData[jj], rhs.vecData[jj], epsilon)==false) return false;

	return true;
}

397
template<class T> bool Array1D<T>::checkEpsilonEquality(const Array1D<double>& rhs1, const Array1D<double>& rhs2, const double& epsilon) {
398
399
400
	return rhs1.checkEpsilonEquality(rhs2, epsilon);
}

401
template<class T> Array1D<T>& Array1D<T>::operator=(const Array1D<T>& source) {
402
403
404
	if(this != &source) {
		vecData = source.vecData;
		nx = source.nx;
405
		keep_nodata = source.keep_nodata;
406
407
408
409
	}
	return *this;
}

410
411
template<class T> Array1D<T>& Array1D<T>::operator=(const T& value) {
	//reset every single member of the Array1D<T>
412
413
414
415
	std::fill(vecData.begin(), vecData.end(), value);
	return *this;
}

416
template<class T> Array1D<T>& Array1D<T>::operator+=(const Array1D<T>& rhs)
417
418
{
	//They have to have equal size
419
420
	if (rhs.nx != nx) {
		std::stringstream ss;
421
		ss << "Trying to add two Array1D objects with different dimensions: ";
422
423
424
		ss << "(" << nx << ") + (" << rhs.nx << ")";
		throw IOException(ss.str(), AT);
	}
425

426
	//Add to every single member of the Array1D<T>
427
	if(keep_nodata==false) {
428
		for (size_t ii=0; ii<nx; ii++) {
429
			vecData[ii] += rhs(ii);
430
431
		}
	} else {
432
		for (size_t ii=0; ii<nx; ii++) {
433
434
			if(vecData[ii]==IOUtils::nodata || rhs(ii)==IOUtils::nodata)
				vecData[ii] = IOUtils::nodata;
435
			else
436
				vecData[ii] += rhs(ii);
437
		}
438
	}
439
440
441
442

	return *this;
}

443
template<class T> const Array1D<T> Array1D<T>::operator+(const Array1D<T>& rhs) const
444
{
445
	Array1D<T> result(*this); //make a copy
446
447
448
449
450
	result += rhs; //already implemented

	return result;
}

451
template<class T> Array1D<T>& Array1D<T>::operator+=(const T& rhs)
452
{
453
	//Add to every single member of the Array1D<T>
454
	if(keep_nodata==false) {
455
		for (size_t ii=0; ii<nx; ii++) {
456
			vecData[ii] += rhs;
457
458
		}
	} else {
459
		for (size_t ii=0; ii<nx; ii++) {
460
461
			if(vecData[ii]!=IOUtils::nodata)
				vecData[ii] += rhs;
462
		}
463
464
465
466
467
	}

	return *this;
}

468
template<class T> const Array1D<T> Array1D<T>::operator+(const T& rhs) const
469
{
470
	Array1D<T> result(*this);
471
472
473
474
475
	result += rhs; //already implemented

	return result;
}

476
template<class T> Array1D<T>& Array1D<T>::operator-=(const Array1D<T>& rhs)
477
478
{
	//They have to have equal size
479
480
	if (rhs.nx != nx) {
		std::stringstream ss;
481
		ss << "Trying to substract two Array1D objects with different dimensions: ";
482
483
484
		ss << "(" << nx << ") - (" << rhs.nx << ")";
		throw IOException(ss.str(), AT);
	}
485

486
	//Substract to every single member of the Array1D<T>
487
	if(keep_nodata==false) {
488
		for (size_t ii=0; ii<nx; ii++) {
489
			vecData[ii] -= rhs(ii);
490
491
		}
	} else {
492
		for (size_t ii=0; ii<nx; ii++) {
493
494
			if(vecData[ii]==IOUtils::nodata || rhs(ii)==IOUtils::nodata)
				vecData[ii] = IOUtils::nodata;
495
			else
496
				vecData[ii] -= rhs(ii);
497
		}
498
499
500
501
502
	}

	return *this;
}

503
template<class T> const Array1D<T> Array1D<T>::operator-(const Array1D<T>& rhs) const
504
{
505
	Array1D<T> result(*this); //make a copy
506
507
508
509
510
	result -= rhs; //already implemented

	return result;
}

511
template<class T> Array1D<T>& Array1D<T>::operator-=(const T& rhs)
512
{
513
	*this += -rhs;
514
515
516
	return *this;
}

517
template<class T> const Array1D<T> Array1D<T>::operator-(const T& rhs) const
518
{
519
	Array1D<T> result(*this);
520
	result += -rhs; //already implemented
521
522
523
524

	return result;
}

525
template<class T> Array1D<T>& Array1D<T>::operator*=(const Array1D<T>& rhs)
526
527
{
	//They have to have equal size
528
529
	if (rhs.nx != nx){
		std::stringstream ss;
530
		ss << "Trying to multiply two Array1D objects with different dimensions: ";
531
532
533
		ss << "(" << nx << ") * (" << rhs.nx << ")";
		throw IOException(ss.str(), AT);
	}
534
	//Multiply every single member of the Array1D<T>
535
	if(keep_nodata==false) {
536
		for (size_t ii=0; ii<nx; ii++) {
537
			vecData[ii] *= rhs(ii);
538
539
		}
	} else {
540
		for (size_t ii=0; ii<nx; ii++) {
541
542
			if(vecData[ii]==IOUtils::nodata || rhs(ii)==IOUtils::nodata)
				vecData[ii] = IOUtils::nodata;
543
			else
544
				vecData[ii] *= rhs(ii);
545
		}
546
547
548
549
550
	}

	return *this;
}

551
template<class T> const Array1D<T> Array1D<T>::operator*(const Array1D<T>& rhs) const
552
{
553
	Array1D<T> result(*this); //make a copy
554
555
556
557
558
	result *= rhs; //already implemented

	return result;
}

559
template<class T> Array1D<T>& Array1D<T>::operator*=(const T& rhs)
560
{
561
	//Multiply every single member of the Array1D<T>
562
	if(keep_nodata==false) {
563
		for (size_t ii=0; ii<nx; ii++) {
564
			vecData[ii] *= rhs;
565
566
		}
	} else {
567
		for (size_t ii=0; ii<nx; ii++) {
568
569
			if(vecData[ii]!=IOUtils::nodata)
				vecData[ii] *= rhs;
570
		}
571
572
573
574
575
	}

	return *this;
}

576
template<class T> const Array1D<T> Array1D<T>::operator*(const T& rhs) const
577
{
578
	Array1D<T> result(*this);
579
580
581
582
583
	result *= rhs; //already implemented

	return result;
}

584
template<class T> Array1D<T>& Array1D<T>::operator/=(const Array1D<T>& rhs)
585
586
{
	//They have to have equal size
587
588
	if (rhs.nx != nx){
		std::stringstream ss;
589
		ss << "Trying to divide two Array1D objects with different dimensions: ";
590
591
592
		ss << "(" << nx << ") / (" << rhs.nx << ")";
		throw IOException(ss.str(), AT);
	}
593
	//Divide every single member of the Array1D<T>
594
	if(keep_nodata==false) {
595
		for (size_t ii=0; ii<nx; ii++) {
596
			vecData[ii] /= rhs(ii);
597
598
		}
	} else {
599
		for (size_t ii=0; ii<nx; ii++) {
600
601
			if(vecData[ii]==IOUtils::nodata || rhs(ii)==IOUtils::nodata)
				vecData[ii] = IOUtils::nodata;
602
			else
603
				vecData[ii] /= rhs(ii);
604
		}
605
606
607
608
609
	}

	return *this;
}

610
template<class T> const Array1D<T> Array1D<T>::operator/(const Array1D<T>& rhs) const
611
{
612
	Array1D<T> result(*this); //make a copy
613
614
615
616
617
	result /= rhs; //already implemented

	return result;
}

618
template<class T> Array1D<T>& Array1D<T>::operator/=(const T& rhs)
619
{
620
	*this *= (1./rhs);
621
622
623
	return *this;
}

624
template<class T> const Array1D<T> Array1D<T>::operator/(const T& rhs) const
625
{
626
	Array1D<T> result(*this);
627
	result *= 1./rhs; //already implemented
628
629
630
631

	return result;
}

632
template<class T> bool Array1D<T>::operator==(const Array1D<T>& in) const {
633
	const size_t in_nx = in.getNx();
634
635
636
637

	if(nx!=in_nx)
		return false;

638
	for(size_t jj=0; jj<nx; jj++)
639
640
641
642
643
644
645
646
647
		if( !IOUtils::checkEpsilonEquality( vecData[jj] , in.vecData[jj], 1e-6) ) return false;

	return true;
}

template<class T> bool Array1D<T>::operator!=(const Array1D<T>& in) const {
	return !(*this==in);
}

648
649
650
} //end namespace mio

#endif