WSL/SLF GitLab Repository

Array2D.h 22.3 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
/***********************************************************************************/
/*  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/>.
*/
#ifndef ARRAY2D_H
#define ARRAY2D_H

#include <meteoio/IOExceptions.h>
#include <meteoio/IOUtils.h>

#include <vector>
#include <limits>
#include <iostream>

namespace mio {

template <class T> class Array2D;

/**
 * @class Array2DProxy
 * @brief The template class Array2DProxy is a helper class for the template class Array2D
 *        with the purpose of adding the [][] operator to Array2D
 *
 * @author Thomas Egger
 */
template <class T> class Array2DProxy {
	public:
		friend class Array2D<T>;
		T& operator[](const unsigned int& j) {
			return array2D(anx, j);
		}

	private:
47
		Array2DProxy(Array2D<T>& i_array2D, const unsigned int& i_anx) : array2D(i_array2D), anx(i_anx){}
48
49
		Array2D<T>& array2D;
		const unsigned int anx;
50
};
51
52
53

/**
 * @class Array2D
54
 * @brief The template class Array2D is a 2D Array (Matrix) able to hold any type of object as datatype.
55
 * It relies on the Array2DProxy class to provide the [][] operator (slower than the (i,j) call).
56
 * If the compilation flag NOSAFECHECKS is used, bounds check is turned off (leading to increased performances).
57
 *
58
 * @ingroup data_str
59
60
61
62
63
 * @author Thomas Egger
 */
template<class T> class Array2D {
	public:
		Array2D();
64
65
66
67
68
69

		/**
		* A constructor that creates an array of a given size
		* @param anx number of columns of the new array
		* @param any number of rows of the new array
		*/
70
		Array2D(const unsigned int& anx, const unsigned int& any);
71
72
73
74
75
76
77

		/**
		* A constructor that creates an array filled with constant values
		* @param anx number of columns of the new array
		* @param any number of rows of the new array
		* @param init initial value to fill the array with
		*/
78
79
		Array2D(const unsigned int& anx, const unsigned int& any, const T& init);

80
81
		virtual ~Array2D();

82
83
		/**
		* A constructor that can be used to create an Array2D object that is contained in the
84
85
86
87
88
89
90
		* one passed as i_array2D argument. The resulting Array2D object is a by value copy of
		* a subplane of the plane spanned by the i_array2D
		* @param i_array2D array containing to extract the values from
		* @param i_nx lower left corner cell X index
		* @param i_ny lower left corner cell Y index
		* @param i_ncols number of columns of the new array
		* @param i_nrows number of rows of the new array
91
		*/
92
		Array2D(const Array2D<T>& i_array2D, const unsigned int& i_nx, const unsigned int& i_ny,
93
		        const unsigned int& i_ncols, const unsigned int& i_nrows);
94

95
		/**
96
		* @brief A method that can be used to cut out a subplane of an existing Array2D object
97
98
99
100
101
102
103
		* that is passed as i_array2D argument. The resulting Array2D object is a by value copy of
		* a subplane of the plane spanned by the i_array2D
		* @param i_array2D array containing to extract the values from
		* @param i_nx lower left corner cell X index
		* @param i_ny lower left corner cell Y index
		* @param i_ncols number of columns of the new array
		* @param i_nrows number of rows of the new array
104
		*/
105
106
		void subset(const Array2D<T>& i_array2D, const unsigned int& i_nx, const unsigned int& i_ny,
		            const unsigned int& i_ncols, const unsigned int& i_nrows);
107

108
109
110
111
112
113
114
115
116
117
118
119
120
		/**
		* @brief A method that can be used to insert a subplane into an existing Array2D object
		* that is passed as i_array2D argument. This is exactly the opposite of the subset method
		* an can be used to rebuild an array from subsets.
		* @param i_array2D array containing to subset values
		* @param i_nx lower left corner cell X index
		* @param i_ny lower left corner cell Y index
		* @param i_ncols number of columns of the new array
		* @param i_nrows number of rows of the new array
		*/
		void fill(const Array2D<T>& i_array2D, const unsigned int& i_nx, const unsigned int& i_ny,
		            const unsigned int& i_ncols, const unsigned int& i_nrows);

121
122
		void fill(const Array2D<T>& i_array2D, const unsigned int& i_nx, const unsigned int& i_ny);

123
124
		/**
		* @brief set how to process nodata values (ie: as nodata or as normal numbers)
125
126
		* @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.
127
		*/
128
129
130
131
132
133
134
		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();
135

136
137
138
		void resize(const unsigned int& nx, const unsigned int& ny);
		void resize(const unsigned int& nx, const unsigned int& ny, const T& init);
		void size(unsigned int& nx, unsigned int& ny) const;
139
140
		unsigned int getNx() const;
		unsigned int getNy() const;
141

142
143
144
		void clear();
		bool isEmpty() const;

145
146
147
148
		/**
		* @brief returns the minimum value contained in the grid
		* @return minimum value
		*/
149
		T getMin() const;
150
151
152
153
		/**
		* @brief returns the maximum value contained in the grid
		* @return maximum value
		*/
154
		T getMax() const;
155
156
157
158
		/**
		* @brief returns the mean value contained in the grid
		* @return mean value
		*/
159
		T getMean() const;
160
		/**
161
		* @brief returns the number of points contained in the grid.
162
163
		* 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
164
		* @return count
165
		*/
166
		size_t getCount() const;
167
168
169
170
		/**
		* @brief returns the grid of the absolute value of values contained in the grid
		* @return grid of abs(grid)
		*/
171
172
		const Array2D<T> getAbs() const;
		void abs();
173

174
175
		template<class P> friend std::ostream& operator<<(std::ostream& os, const Array2D<P>& array);

176
177
178
		bool checkEpsilonEquality(const Array2D<double>& rhs, const double& epsilon) const;
		static bool checkEpsilonEquality(const Array2D<double>& rhs1, const Array2D<double>& rhs2, const double& epsilon);

179
180
		T& operator ()(const unsigned int& x, const unsigned int& y);
		const T operator ()(const unsigned int& x, const unsigned int& y) const;
181
182
		T& operator ()(const unsigned int& i);
		const T operator ()(const unsigned int& i) const;
183
184
185
		Array2DProxy<T> operator[](const unsigned int& i);

		Array2D<T>& operator =(const Array2D<T>&);
186
		Array2D<T>& operator =(const T& value);
187

188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
		Array2D<T>& operator+=(const T& rhs);
		const Array2D<T> operator+(const T& rhs);
		Array2D<T>& operator+=(const Array2D<T>& rhs);
		const Array2D<T> operator+(const Array2D<T>& rhs);

		Array2D<T>& operator-=(const T& rhs);
		const Array2D<T> operator-(const T& rhs);
		Array2D<T>& operator-=(const Array2D<T>& rhs);
		const Array2D<T> operator-(const Array2D<T>& rhs);

		Array2D<T>& operator*=(const T& rhs);
		const Array2D<T> operator*(const T& rhs);
		Array2D<T>& operator*=(const Array2D<T>& rhs);
		const Array2D<T> operator*(const Array2D<T>& rhs);

		Array2D<T>& operator/=(const T& rhs);
		const Array2D<T> operator/(const T& rhs);
		Array2D<T>& operator/=(const Array2D<T>& rhs);
		const Array2D<T> operator/(const Array2D<T>& rhs);

	protected:
		std::vector<T> vecData;
		unsigned int nx;
		unsigned int ny;
212
		bool keep_nodata;
213
214
};

215
template<class T> inline T& Array2D<T>::operator()(const unsigned int& i) {
216
217
218
#ifndef NOSAFECHECKS
	return vecData.at(i);
#else
219
	return vecData[i];
220
#endif
221
222
}

223
template<class T> inline const T Array2D<T>::operator()(const unsigned int& i) const {
224
225
226
#ifndef NOSAFECHECKS
	return vecData.at(i);
#else
227
	return vecData[i];
228
#endif
229
}
230
template<class T> inline T& Array2D<T>::operator()(const unsigned int& x, const unsigned int& y) {
231
232
#ifndef NOSAFECHECKS
	if ((x >= nx) || (y >= ny)) {
233
234
		std::stringstream ss;
		ss << "Trying to access array(" << x << "," << y << ")";
235
		ss << " while array is (" << nx << "," << ny << ")";
236
		throw IndexOutOfBoundsException(ss.str(), AT);
237
238
	}
#endif
239
	//COLUMN-MAJOR alignment of the vector: fully C-compatible memory layout
240
	return vecData[x + y*nx];
241
242
}

243
template<class T> inline const T Array2D<T>::operator()(const unsigned int& x, const unsigned int& y) const {
244
245
#ifndef NOSAFECHECKS
	if ((x >= nx) || (y >= ny)) {
246
247
		std::stringstream ss;
		ss << "Trying to access array(" << x << "," << y << ")";
248
		ss << " while array is (" << nx << "," << ny << ")";
249
		throw IndexOutOfBoundsException(ss.str(), AT);
250
251
	}
#endif
252
	return vecData[x + y*nx];
253
254
255
}

template<class T> Array2DProxy<T> Array2D<T>::operator[](const unsigned int& i) {
256
	return Array2DProxy<T>(*this, i);
257
258
259
260
}

template<class T> Array2D<T>::Array2D() {
	nx = ny = 0;
261
	keep_nodata = true;
262
263
}

264
265
template<class T> Array2D<T>::~Array2D() { }

266
267
template<class T> Array2D<T>::Array2D(const Array2D<T>& i_array2D, const unsigned int& i_nx, const unsigned int& i_ny,
                                      const unsigned int& i_ncols, const unsigned int& i_nrows)
268
{
269
	keep_nodata = true;
270
	subset(i_array2D, i_nx, i_ny, i_ncols, i_nrows);
271
272
}

273
274
template<class T> void Array2D<T>::subset(const Array2D<T>& i_array2D, const unsigned int& i_nx, const unsigned int& i_ny,
                                          const unsigned int& i_ncols, const unsigned int& i_nrows)
275
{
276
277
278
279
280
281
	if (((i_nx+i_ncols) > i_array2D.nx) || ((i_ny+i_nrows) > i_array2D.ny)) {
		std::stringstream ss;
		ss << "Trying to cut an array of size (" << nx << "," << ny << ") ";
		ss << "to size (" << i_ncols << "," << i_nrows << ") starting at (" << i_nx << "," << i_ny << ")";
		throw IndexOutOfBoundsException(ss.str(), AT);
	}
282

283
	if ((i_ncols == 0) || (i_nrows == 0)) //the plane to copy has to make sense
284
		throw IndexOutOfBoundsException("Trying to cut an array into a null sized array!", AT);
285

286
	resize(i_ncols, i_nrows); //create new Array2D object
287
	//Copy by value subspace
288
289
	for (unsigned int jj=0; jj<ny; jj++) {
		for (unsigned int ii=0; ii<nx; ii++) {
290
			operator()(ii,jj) = i_array2D(i_nx+ii, i_ny+jj);
291
292
293
294
		}
	}
}

295
296
297
298
299
300
301
template<class T> void Array2D<T>::fill(const Array2D<T>& i_array2D, const unsigned int& i_nx, const unsigned int& i_ny)
{
	unsigned int i_ncols, i_nrows;
	i_array2D.size(i_ncols, i_nrows);
	fill(i_array2D, i_nx, i_ny, i_ncols, i_nrows);
}

302
303
304
template<class T> void Array2D<T>::fill(const Array2D<T>& i_array2D, const unsigned int& i_nx, const unsigned int& i_ny,
                                        const unsigned int& i_ncols, const unsigned int& i_nrows)
{
305
306
307
308
309
310
311
	if (((i_nx+i_ncols) > nx) || ((i_ny+i_nrows) > ny)) {
		std::stringstream ss;
		ss << "Filling an array of size (" << nx << "," << ny << ") ";
		ss << "with an array of size (" << i_ncols << "," << i_nrows << ") ";
		ss << "starting at (" << i_nx << "," << i_ny << ")";
		throw IndexOutOfBoundsException(ss.str(), AT);
	}
312
313

	if ((i_ncols == 0) || (i_nrows == 0)) //the plane to copy has to make sense
314
		throw IndexOutOfBoundsException("Filling an array with a null sized array!", AT);
315
316
317
318
319
320
321
322
323
324

	for(unsigned int jj=i_ny; jj<(i_ny+i_nrows); jj++) {
		for(unsigned int ii=i_nx; ii<(i_nx+i_ncols); ii++) {
			const unsigned int ix = ii-i_nx;
			const unsigned int iy = jj-i_ny;
			operator()(ii,jj) = i_array2D(ix, iy);
		}
	}
}

325
326
template<class T> Array2D<T>::Array2D(const unsigned int& anx, const unsigned int& any, const T& init) {
	nx = ny = 0;
327
	keep_nodata = true;
328
329
330
331
332
	resize(anx,any,init);
}

template<class T> Array2D<T>::Array2D(const unsigned int& anx, const unsigned int& any) {
	nx = ny = 0;
333
	keep_nodata = true;
334
335
336
	resize(anx,any);
}

337
338
339
340
341
342
template<class T> void Array2D<T>::setKeepNodata(const bool i_keep_nodata) {
	keep_nodata = i_keep_nodata;
}

template<class T> bool Array2D<T>::getKeepNodata() {
	return keep_nodata;
343
344
}

345
template<class T> void Array2D<T>::resize(const unsigned int& anx, const unsigned int& any) {
346
347
348
349
	clear(); //we won't be able to "rescue" old values, so we reset the whole vector
	vecData.resize(anx*any);
	nx = anx;
	ny = any;
350
351
352
}

template<class T> void Array2D<T>::resize(const unsigned int& anx, const unsigned int& any, const T& init) {
353
	clear(); //we won't be able to "rescue" old values, so we reset the whole vector
354
355
356
	vecData.resize(anx*any, init);
	nx = anx;
	ny = any;
357
358
}

359
template<class T> void Array2D<T>::size(unsigned int& anx, unsigned int& any) const {
360
361
362
363
	anx=nx;
	any=ny;
}

364
365
366
367
368
369
370
371
template<class T> unsigned int Array2D<T>::getNx() const {
	return nx;
}

template<class T> unsigned int Array2D<T>::getNy() const {
	return ny;
}

372
373
374
375
376
template<class T> void Array2D<T>::clear() {
	vecData.clear();
	nx=ny=0;
}

377
378
379
380
template<class T> bool Array2D<T>::isEmpty() const {
	return (nx==0 && ny==0);
}

381
382
template<class T> std::ostream& operator<<(std::ostream& os, const Array2D<T>& array) {
	os << "<array2d>\n";
383
	for(unsigned int jj=0; jj<array.ny; jj++) {
384
		unsigned int jnx = jj*array.nx;
385
		for (unsigned int ii=0; ii<array.nx; ii++) {
386
			os << array(ii+jnx) << " ";
387
388
389
390
391
392
393
		}
		os << "\n";
	}
	os << "</array2d>\n";
	return os;
}

394
template<class T> T Array2D<T>::getMin() const {
395
396
397

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

398
	const unsigned int nxy = ny*nx;
399
	if(keep_nodata==false) {
400
401
402
		for (unsigned int jj=0; jj<nxy; jj++) {
			const T val = operator()(jj);
			if(val<min) min=val;
403
404
		}
		return min;
405
	} else {
406
407
408
		for (unsigned int jj=0; jj<nxy; jj++) {
			const T val = operator()(jj);
			if(val!=IOUtils::nodata && val<min) min=val;
409
410
411
412
413
414
		}
		if(min!=std::numeric_limits<T>::max()) return min;
		else return (T)IOUtils::nodata;
	}
}

415
template<class T> T Array2D<T>::getMax() const {
416
417
418

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

419
	const unsigned int nxy = ny*nx;
420
	if(keep_nodata==false) {
421
422
423
		for (unsigned int jj=0; jj<nxy; jj++) {
			const T val = operator()(jj);
			if(val>max) max=val;
424
425
		}
		return max;
426
	} else {
427
428
429
		for (unsigned int jj=0; jj<nxy; jj++) {
			const T val = operator()(jj);
			if(val!=IOUtils::nodata && val>max) max=val;
430
431
432
433
434
435
		}
		if(max!=-std::numeric_limits<T>::max()) return max;
		else return (T)IOUtils::nodata;
	}
}

436
template<class T> T Array2D<T>::getMean() const {
437
438

	T mean = 0;
439
	const unsigned int nxy = nx*ny;
440

441
	if(keep_nodata==false) {
442
		for (unsigned int jj=0; jj<nxy; jj++) {
443
444
445
			const T val = operator()(jj);
			mean += val;
		}
446
		if(nxy>0) return mean/(T)(nxy);
447
		else return (T)0;
448
	} else {
449
		unsigned int count = 0;
450
451
452
453
454
		for (unsigned int jj=0; jj<nxy; jj++) {
			const T val = operator()(jj);
			if(val!=IOUtils::nodata) {
				mean += val;
				count++;
455
456
457
458
459
460
461
			}
		}
		if(count>0) return mean/(T)(count);
		else return (T)IOUtils::nodata;
	}
}

462
template<class T> size_t Array2D<T>::getCount() const
463
{
464
465
	const unsigned int nxy = nx*ny;

466
	if(keep_nodata==false) {
467
		return (size_t)nxy;
468
	} else {
469
		size_t count = 0;
470
		for (unsigned int ii=0; ii<nxy; ii++) {
471
			if(vecData[ii]!=IOUtils::nodata) count++;
472
		}
473
		return count;
474
475
476
	}
}

477
template<class T> void Array2D<T>::abs() {
478
479
	if(std::numeric_limits<T>::is_signed) {
		const unsigned int nxy = nx*ny;
480
		if(keep_nodata==false) {
481
482
483
484
			for (unsigned int ii=0; ii<nxy; ii++) {
				T& val = operator()(ii);
				if(val<0) val=-val;
			}
485
		} else {
486
487
488
489
490
491
492
493
			for (unsigned int ii=0; ii<nxy; ii++) {
				T& val = operator()(ii);
				if(val<0 && val!=IOUtils::nodata) val=-val;
			}
		}
	}
}

494
template<class T> const Array2D<T> Array2D<T>::getAbs() const {
495
	Array2D<T> result = *this; //make a copy
496
	result.abs(); //already implemented
497
498
499
500
501

	return result;
}


502
//arithmetic operators
503
504
505
506
507
508
509
510
511
512
513
514
515
516
template<class T> bool Array2D<T>::checkEpsilonEquality(const Array2D<double>& rhs, const double& epsilon) const {
	if(nx!=rhs.nx || ny!=rhs.ny) return false;

	const unsigned int nxy = nx*ny;
	for (unsigned int jj=0; jj<nxy; jj++)
		if(IOUtils::checkEpsilonEquality(vecData[jj], rhs.vecData[jj], epsilon)==false) return false;

	return true;
}

template<class T> bool Array2D<T>::checkEpsilonEquality(const Array2D<double>& rhs1, const Array2D<double>& rhs2, const double& epsilon) {
	return rhs1.checkEpsilonEquality(rhs2, epsilon);
}

517
518
template<class T> Array2D<T>& Array2D<T>::operator=(const Array2D<T>& source) {
	if(this != &source) {
519
		keep_nodata = source.keep_nodata;
520
521
522
523
524
525
526
		nx = source.nx;
		ny = source.ny;
		vecData = source.vecData;
	}
	return *this;
}

527
template<class T> Array2D<T>& Array2D<T>::operator=(const T& value) {
528
529
	std::fill(vecData.begin(), vecData.end(), value);
	//for(unsigned int i=0; i<vecData.size(); i++) vecData[i] = value;
530
531
532
	return *this;
}

533
534
535
template<class T> Array2D<T>& Array2D<T>::operator+=(const Array2D<T>& rhs)
{
	//They have to have equal size
536
537
538
539
540
541
	if ((rhs.nx != nx) || (rhs.ny != ny)) {
		std::stringstream ss;
		ss << "Trying to add two Array2D objects with different dimensions: ";
		ss << "(" << nx << "," << ny << ") + (" << rhs.nx << "," << rhs.ny << ")";
		throw IOException(ss.str(), AT);
	}
542

543
	const unsigned int nxy = nx*ny;
544
	//Add to every single member of the Array2D<T>
545
546
	if(keep_nodata==false) {
		for (unsigned int jj=0; jj<nxy; jj++)
547
			operator()(jj) += rhs(jj);
548
549
550
551
552
553
554
555
	} else {
		for (unsigned int jj=0; jj<nxy; jj++) {
			if(operator()(jj)==IOUtils::nodata || rhs(jj)==IOUtils::nodata)
				operator()(jj) = IOUtils::nodata;
			else
				operator()(jj) += rhs(jj);
		}
	}
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570

	return *this;
}

template<class T> const Array2D<T> Array2D<T>::operator+(const Array2D<T>& rhs)
{
	Array2D<T> result = *this; //make a copy
	result += rhs; //already implemented

	return result;
}

template<class T> Array2D<T>& Array2D<T>::operator+=(const T& rhs)
{
	//Add to every single member of the Array2D<T>
571
	const unsigned int nxy = nx*ny;
572
573
574
575
576
577
578
579
580
581

	if(keep_nodata==false) {
		for (unsigned int jj=0; jj<nxy; jj++)
			operator()(jj) += rhs;
	} else {
		for (unsigned int jj=0; jj<nxy; jj++) {
			if(operator()(jj)!=IOUtils::nodata)
				operator()(jj) += rhs;
		}
	}
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596

	return *this;
}

template<class T> const Array2D<T> Array2D<T>::operator+(const T& rhs)
{
	Array2D<T> result = *this;
	result += rhs; //already implemented

	return result;
}

template<class T> Array2D<T>& Array2D<T>::operator-=(const Array2D<T>& rhs)
{
	//They have to have equal size
597
598
599
600
601
602
	if ((rhs.nx != nx) || (rhs.ny != ny)){
		std::stringstream ss;
		ss << "Trying to substract two Array2D objects with different dimensions: ";
		ss << "(" << nx << "," << ny << ") - (" << rhs.nx << "," << rhs.ny << ")";
		throw IOException(ss.str(), AT);
	}
603
	//Substract to every single member of the Array2D<T>
604
	const unsigned int nxy = nx*ny;
605
606
607
608
609
610
611
612
613
614
615
616

	if(keep_nodata==false) {
		for (unsigned int jj=0; jj<nxy; jj++)
			operator()(jj) -= rhs(jj);
	} else {
		for (unsigned int jj=0; jj<nxy; jj++) {
			if(operator()(jj)==IOUtils::nodata || rhs(jj)==IOUtils::nodata)
				operator()(jj) = IOUtils::nodata;
			else
				operator()(jj) -= rhs(jj);
		}
	}
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631

	return *this;
}

template<class T> const Array2D<T> Array2D<T>::operator-(const Array2D<T>& rhs)
{
	Array2D<T> result = *this; //make a copy
	result -= rhs; //already implemented

	return result;
}

template<class T> Array2D<T>& Array2D<T>::operator-=(const T& rhs)
{
	//Substract to every single member of the Array2D<T>
632
	const unsigned int nxy = nx*ny;
633

634
635
636
637
638
639
640
641
642
	if(keep_nodata==false) {
		for (unsigned int jj=0; jj<nxy; jj++)
			operator()(jj) -= rhs;
	} else {
		for (unsigned int jj=0; jj<nxy; jj++) {
			if(operator()(jj)!=IOUtils::nodata)
				operator()(jj) -= rhs;
		}
	}
643
644
645
646
647
648
649
650
651
652
653
654
655
656
	return *this;
}

template<class T> const Array2D<T> Array2D<T>::operator-(const T& rhs)
{
	Array2D<T> result = *this;
	result -= rhs; //already implemented

	return result;
}

template<class T> Array2D<T>& Array2D<T>::operator*=(const Array2D<T>& rhs)
{
	//They have to have equal size
657
658
659
660
661
662
	if ((rhs.nx != nx) || (rhs.ny != ny)){
		std::stringstream ss;
		ss << "Trying to multiply two Array2D objects with different dimensions: ";
		ss << "(" << nx << "," << ny << ") * (" << rhs.nx << "," << rhs.ny << ")";
		throw IOException(ss.str(), AT);
	}
663
	//Add to every single member of the Array2D<T>
664
	const unsigned int nxy = nx*ny;
665
666
667
668
669
670
671
672
673
674
675
676

	if(keep_nodata==false) {
		for (unsigned int jj=0; jj<nxy; jj++)
			operator()(jj) *= rhs(jj);
	} else {
		for (unsigned int jj=0; jj<nxy; jj++) {
			if(operator()(jj)==IOUtils::nodata || rhs(jj)==IOUtils::nodata)
				operator()(jj) = IOUtils::nodata;
			else
				operator()(jj) *= rhs(jj);
		}
	}
677
678
679
680
681
682
683
684
685
686
687
688
689
690

	return *this;
}

template<class T> const Array2D<T> Array2D<T>::operator*(const Array2D<T>& rhs)
{
	Array2D<T> result = *this; //make a copy
	result *= rhs; //already implemented

	return result;
}

template<class T> Array2D<T>& Array2D<T>::operator*=(const T& rhs)
{
691
	//Multiply to every single member of the Array2D<T>
692
	const unsigned int nxy = nx*ny;
693
694
695

	if(keep_nodata==false) {
		for (unsigned int jj=0; jj<nxy; jj++)
696
			operator()(jj) *= rhs;
697
698
699
700
701
702
	} else {
		for (unsigned int jj=0; jj<nxy; jj++) {
			if(operator()(jj)!=IOUtils::nodata)
				operator()(jj) *= rhs;
		}
	}
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717

	return *this;
}

template<class T> const Array2D<T> Array2D<T>::operator*(const T& rhs)
{
	Array2D<T> result = *this;
	result *= rhs; //already implemented

	return result;
}

template<class T> Array2D<T>& Array2D<T>::operator/=(const Array2D<T>& rhs)
{
	//They have to have equal size
718
719
720
721
722
723
	if ((rhs.nx != nx) || (rhs.ny != ny)){
		std::stringstream ss;
		ss << "Trying to divide two Array2D objects with different dimensions: ";
		ss << "(" << nx << "," << ny << ") / (" << rhs.nx << "," << rhs.ny << ")";
		throw IOException(ss.str(), AT);
	}
724
	//Divide every single member of the Array2D<T>
725
	const unsigned int nxy = nx*ny;
726
727
728
729
730
731
732
733
734
735
736
737

	if(keep_nodata==false) {
		for (unsigned int jj=0; jj<nxy; jj++)
			operator()(jj) /= rhs(jj);
	} else {
		for (unsigned int jj=0; jj<nxy; jj++) {
			if(operator()(jj)==IOUtils::nodata || rhs(jj)==IOUtils::nodata)
				operator()(jj) = IOUtils::nodata;
			else
				operator()(jj) /= rhs(jj);
		}
	}
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752

	return *this;
}

template<class T> const Array2D<T> Array2D<T>::operator/(const Array2D<T>& rhs)
{
	Array2D<T> result = *this; //make a copy
	result /= rhs; //already implemented

	return result;
}

template<class T> Array2D<T>& Array2D<T>::operator/=(const T& rhs)
{
	//Divide every single member of the Array2D<T>
753
	const unsigned int nxy = nx*ny;
754
755
756
757
758
759
760
761
762
763

	if(keep_nodata==false) {
		for (unsigned int jj=0; jj<nxy; jj++)
			operator()(jj) /= rhs;
	} else {
		for (unsigned int jj=0; jj<nxy; jj++) {
			if(operator()(jj)!=IOUtils::nodata)
				operator()(jj) /= rhs;
		}
	}
764
765
766
767
768
769
770
771
772
773
774
775
776
777
	return *this;
}

template<class T> const Array2D<T> Array2D<T>::operator/(const T& rhs)
{
	Array2D<T> result = *this;
	result /= rhs; //already implemented

	return result;
}

} //end namespace mio

#endif