MeteoIODoc 2.11.0
Array3D.h
Go to the documentation of this file.
1// SPDX-License-Identifier: LGPL-3.0-or-later
2/***********************************************************************************/
3/* Copyright 2009 WSL Institute for Snow and Avalanche Research SLF-DAVOS */
4/***********************************************************************************/
5/* This file is part of MeteoIO.
6 MeteoIO is free software: you can redistribute it and/or modify
7 it under the terms of the GNU Lesser General Public License as published by
8 the Free Software Foundation, either version 3 of the License, or
9 (at your option) any later version.
10
11 MeteoIO is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 GNU Lesser General Public License for more details.
15
16 You should have received a copy of the GNU Lesser General Public License
17 along with MeteoIO. If not, see <http://www.gnu.org/licenses/>.
18*/
19#ifndef ARRAY3D_H
20#define ARRAY3D_H
21
22#include <meteoio/IOUtils.h>
24
25#include <vector>
26#include <limits>
27#include <iostream>
28#include <numeric>
29#include <algorithm>
30
31//forward declaration
32namespace mio { template <class T> class Array3D; template <class T> class Array3DProxy2; }
34
35namespace mio {
43template <class T> class Array3DProxy {
44 public:
45 friend class Array3D<T>;
46 Array3DProxy2<T> operator[](const size_t& i_any) {
47 return Array3DProxy2<T>(array3D, anx, i_any);
48 }
49
50 private:
51 Array3DProxy(Array3D<T>& i_array3D, const size_t& i_anx) : array3D(i_array3D), anx(i_anx){}
52 Array3D<T>& array3D;
53 const size_t anx;
54};
55
63template <class T> class Array3DProxy2 {
64 public:
65 friend class Array3DProxy<T>;
66 T& operator[](const size_t& i_anz) {
67 return array3D(anx, any, i_anz);
68 }
69
70 private:
71 Array3DProxy2(Array3D<T>& i_array3D, const size_t& i_anx,
72 const size_t& i_any) : array3D(i_array3D), anx(i_anx), any(i_any){}
73 Array3D<T>& array3D;
74 const size_t anx;
75 const size_t any;
76};
77
87template<class T> class Array3D {
88 public:
90
103 Array3D(const Array3D<T>& i_array3D,
104 const size_t& i_nx, const size_t& i_ny, const size_t& i_nz,
105 const size_t& i_ncols, const size_t& i_nrows, const size_t& i_ndepth);
106
113 Array3D(const size_t& anx, const size_t& any, const size_t& anz);
114
122 Array3D(const size_t& anx, const size_t& any, const size_t& anz, const T& init);
123
136 void subset(const Array3D<T>& i_array3D,
137 const size_t& i_nx, const size_t& i_ny, const size_t& i_nz,
138 const size_t& i_ncols, const size_t& i_nrows, const size_t& i_ndepth);
139
152 void fill(const Array3D<T>& i_array3D,
153 const size_t& i_nx, const size_t& i_ny, const size_t& i_nz,
154 const size_t& i_ncols, const size_t& i_nrows, const size_t& i_ndepth);
155
156 void fill(const Array3D<T>& i_array3D, const size_t& i_nx, const size_t& i_ny, const size_t& i_nz);
157
165 void insert(const Array2D<T>& layer, const size_t& depth);
166
172 void setKeepNodata(const bool i_keep_nodata);
173
178 bool getKeepNodata() const;
179
180 void resize(const size_t& anx, const size_t& any, const size_t& anz);
181 void resize(const size_t& anx, const size_t& any, const size_t& anz, const T& init);
182 void size(size_t& anx, size_t& any, size_t& anz) const;
183 size_t size() const;
184 size_t getNx() const;
185 size_t getNy() const;
186 size_t getNz() const;
187 void clear();
188 bool empty() const;
189
194 T getMin() const;
199 T getMax() const;
204 T getMean() const;
211 size_t getCount() const;
216 const Array3D<T> getAbs() const;
217 void abs();
218
219 const std::string toString() const;
220 template<class P> friend std::ostream& operator<<(std::ostream& os, const Array3D<P>& array);
221 template<class P> friend std::istream& operator>>(std::istream& is, Array3D<P>& array);
222
223 bool checkEpsilonEquality(const Array3D<double>& rhs, const double& epsilon) const;
224 static bool checkEpsilonEquality(const Array3D<double>& rhs1, const Array3D<double>& rhs2, const double& epsilon);
225
226 T& operator ()(const size_t& i);
227 const T operator ()(const size_t& i) const;
228 T& operator ()(const size_t& x, const size_t& y, const size_t& z);
229 const T operator ()(const size_t& x, const size_t& y, const size_t& z) const;
230 Array3DProxy<T> operator[](const size_t& i);
231
232 Array3D<T>& operator =(const T& value);
233
234 Array3D<T>& operator+=(const T& rhs);
235 const Array3D<T> operator+(const T& rhs) const;
237 const Array3D<T> operator+(const Array3D<T>& rhs) const;
238
239 Array3D<T>& operator-=(const T& rhs);
240 const Array3D<T> operator-(const T& rhs) const;
242 const Array3D<T> operator-(const Array3D<T>& rhs) const;
243
244 Array3D<T>& operator*=(const T& rhs);
245 const Array3D<T> operator*(const T& rhs) const;
247 const Array3D<T> operator*(const Array3D<T>& rhs) const;
248
249 Array3D<T>& operator/=(const T& rhs);
250 const Array3D<T> operator/(const T& rhs) const;
252 const Array3D<T> operator/(const Array3D<T>& rhs) const;
253
254 bool operator==(const Array3D<T>&) const;
255 bool operator!=(const Array3D<T>&) const;
256
257 protected:
258 std::vector<T> vecData;
259 size_t nx;
260 size_t ny;
261 size_t nz;
262 size_t nxny; //nx times ny
264};
265
266template<class T> inline T& Array3D<T>::operator()(const size_t& i) {
267#ifndef NOSAFECHECKS
268 return vecData.at(i);
269#else
270 return vecData[i];
271#endif
272}
273
274template<class T> inline const T Array3D<T>::operator()(const size_t& i) const {
275#ifndef NOSAFECHECKS
276 return vecData.at(i);
277#else
278 return vecData[i];
279#endif
280}
281
282template<class T> inline T& Array3D<T>::operator()(const size_t& x, const size_t& y, const size_t& z) {
283#ifndef NOSAFECHECKS
284 if ((x >= nx) || (y >= ny) || (z >= nz)) {
285 std::ostringstream ss;
286 ss << "Trying to access array(" << x << "," << y << "," << z << ")";
287 ss << " while array is (" << nx << "," << ny << "," << nz << ")";
288 throw IndexOutOfBoundsException(ss.str(), AT);
289 }
290#endif
291
292 //ROW-MAJOR alignment of the vector: fully C-compatible memory layout
293 return vecData[x + y*nx + z*nxny];
294}
295
296template<class T> inline const T Array3D<T>::operator()(const size_t& x, const size_t& y, const size_t& z) const {
297#ifndef NOSAFECHECKS
298 if ((x >= nx) || (y >= ny) || (z >= nz)) {
299 std::ostringstream ss;
300 ss << "Trying to access array(" << x << "," << y << "," << z << ")";
301 ss << " while array is (" << nx << "," << ny << "," << nz << ")";
302 throw IndexOutOfBoundsException(ss.str(), AT);
303 }
304#endif
305 return vecData[x + y*nx + z*nxny];
306}
307
308template<class T> Array3DProxy<T> Array3D<T>::operator[](const size_t& i) {
309 return Array3DProxy<T>(*this, i);
310}
311
312template<class T> Array3D<T>::Array3D() : vecData(), nx(0), ny(0), nz(0), nxny(0), keep_nodata(true)
313{
314}
315
316template<class T> Array3D<T>::Array3D(const Array3D<T>& i_array3D,
317 const size_t& i_nx, const size_t& i_ny, const size_t& i_nz,
318 const size_t& i_ncols, const size_t& i_nrows, const size_t& i_ndepth)
319 : vecData(i_ncols*i_nrows*i_ndepth), nx(i_ncols), ny(i_nrows), nz(i_ndepth), nxny(i_ncols*i_nrows), keep_nodata(true)
320{
321 subset(i_array3D, i_nx, i_ny, i_nz, i_ncols, i_nrows, i_ndepth);
322}
323
324template<class T> void Array3D<T>::subset(const Array3D<T>& i_array3D,
325 const size_t& i_nx, const size_t& i_ny, const size_t& i_nz,
326 const size_t& i_ncols, const size_t& i_nrows, const size_t& i_ndepth)
327{
328
329 if (((i_nx+i_ncols) > i_array3D.nx) || ((i_ny+i_nrows) > i_array3D.ny) || ((i_nz+i_ndepth) > i_array3D.nz)) {
330 std::ostringstream ss;
331 ss << "Trying to cut an array of size (" << i_array3D.nx << "," << i_array3D.ny << "," << i_array3D.nz << ") ";
332 ss << "to size (" << i_ncols << "," << i_nrows << "," << i_ndepth << ") ";
333 ss << "starting at (" << i_nx << "," << i_ny << "," << i_nz << ")";
334 throw IndexOutOfBoundsException(ss.str(), AT);
335 }
336
337 if ((i_ncols == 0) || (i_nrows == 0) || (i_ndepth == 0)) //the space has to make sense
338 throw IndexOutOfBoundsException("Trying to cut an array into a null sized array!", AT);
339
340 resize(i_ncols, i_nrows, i_ndepth); //create new Array3D object
341 //Copy by value subspace
342 for (size_t ii=0; ii<nz; ii++) {
343 for (size_t jj=0; jj<ny; jj++) {
344 for (size_t kk=0; kk<nx; kk++) {
345 //Running through the vector in order of memory alignment
346 operator()(kk,jj,ii) = i_array3D(i_nx+kk, i_ny+jj, i_nz+ii);
347 }
348 }
349 }
350}
351
352template<class T> void Array3D<T>::fill(const Array3D<T>& i_array3D, const size_t& i_nx, const size_t& i_ny, const size_t& i_nz)
353{
354 size_t i_ncols, i_nrows, i_ndepth;
355 i_array3D.size(i_ncols, i_nrows, i_ndepth);
356 fill(i_array3D, i_nx, i_ny, i_nz, i_ncols, i_nrows, i_ndepth);
357}
358
359template<class T> void Array3D<T>::fill(const Array3D<T>& i_array3D,
360 const size_t& i_nx, const size_t& i_ny, const size_t& i_nz,
361 const size_t& i_ncols, const size_t& i_nrows, const size_t& i_ndepth)
362{
363
364 if (((i_nx+i_ncols) > i_array3D.nx) || ((i_ny+i_nrows) > i_array3D.ny) || ((i_nz+i_ndepth) > i_array3D.nz)) {
365 std::ostringstream ss;
366 ss << "Filling an array of size (" << nx << "," << ny << "," << nz << ") ";
367 ss << "with an array of size (" << i_ncols << "," << i_nrows << "," << i_ndepth << ") ";
368 ss << "starting at (" << i_nx << "," << i_ny << "," << i_nz << ")";
369 throw IndexOutOfBoundsException(ss.str(), AT);
370 }
371
372 if ((i_ncols == 0) || (i_nrows == 0) || (i_ndepth == 0)) //the space has to make sense
373 throw IndexOutOfBoundsException("Filling an array with a null sized array!", AT);
374
375 //Copy by value subspace
376 for (size_t ii=i_nz; ii<(i_nz+i_ndepth); ii++) {
377 const size_t iz = ii-i_nz;
378 for (size_t jj=i_ny; jj<(i_ny+i_nrows); jj++) {
379 const size_t iy = jj-i_ny;
380 for (size_t kk=i_nx; kk<(i_nx+i_ncols); kk++) {
381 const size_t ix = kk-i_nx;
382 operator()(kk,jj,ii) = i_array3D(ix, iy, iz);
383 }
384 }
385 }
386}
387
388template<class T> void Array3D<T>::insert(const Array2D<T>& layer, const size_t& depth)
389{
390 if (depth>nz) {
391 std::ostringstream ss;
392 ss << "Trying to insert layer at depth " << depth << " ";
393 ss << "in an array of size (" << nx << "," << ny << "," << nz << ") ";
394 throw IndexOutOfBoundsException(ss.str(), AT);
395 }
396 if (layer.getNx()!=nx || layer.getNy()!=ny) {
397 std::ostringstream ss;
398 ss << "Trying to insert layer of size (" << layer.getNx() << "," << layer.getNy() << ") ";
399 ss << "in an array of size (" << nx << "," << ny << "," << nz << ") ";
400 throw IndexOutOfBoundsException(ss.str(), AT);
401 }
402 //copy plane in the correct position
403 for (size_t jj=0; jj<ny; jj++) {
404 for (size_t ii=0; ii<nx; ii++) {
405 operator()(ii,jj,depth) = layer(ii, jj);
406 }
407 }
408}
409
410template<class T> Array3D<T>::Array3D(const size_t& anx, const size_t& any, const size_t& anz)
411 : vecData(anx*any*anz), nx(anx), ny(any), nz(anz), nxny(anx*any), keep_nodata(true)
412{
413 //resize(anx, any, anz);
414}
415
416template<class T> Array3D<T>::Array3D(const size_t& anx, const size_t& any, const size_t& anz, const T& init)
417 : vecData(anx*any*anz, init), nx(anx), ny(any), nz(anz), nxny(anx*any), keep_nodata(true)
418{
419 //resize(anx, any, anz, init);
420}
421
422template<class T> void Array3D<T>::setKeepNodata(const bool i_keep_nodata) {
423 keep_nodata = i_keep_nodata;
424}
425
426template<class T> bool Array3D<T>::getKeepNodata() const {
427 return keep_nodata;
428}
429
430template<class T> void Array3D<T>::resize(const size_t& anx, const size_t& any, const size_t& anz) {
431 clear(); //we won't be able to "rescue" old values, so we reset the whole vector
432 vecData.resize(anx*any*anz);
433 nx = anx;
434 ny = any;
435 nz = anz;
436 nxny = nx*ny;
437}
438
439template<class T> void Array3D<T>::resize(const size_t& anx, const size_t& any, const size_t& anz, const T& init) {
440 clear(); //we won't be able to "rescue" old values, so we reset the whole vector
441 vecData.resize(anx*any*anz, init);
442 nx = anx;
443 ny = any;
444 nz = anz;
445 nxny = nx*ny;
446}
447
448template<class T> void Array3D<T>::size(size_t& anx, size_t& any, size_t& anz) const {
449 anx=nx;
450 any=ny;
451 anz=nz;
452}
453
454template<class T> size_t Array3D<T>::size() const {
455 return nx*ny*nz;
456}
457
458template<class T> size_t Array3D<T>::getNx() const {
459 return nx;
460}
461
462template<class T> size_t Array3D<T>::getNy() const {
463 return ny;
464}
465
466template<class T> size_t Array3D<T>::getNz() const {
467 return nz;
468}
469
470template<class T> void Array3D<T>::clear() {
471 vecData.clear();
472 nx = ny = nz = nxny = 0;
473}
474
475template<class T> bool Array3D<T>::empty() const {
476 return (nx==0 && ny==0 && nz==0);
477}
478
479template<class T> const std::string Array3D<T>::toString() const {
480 std::ostringstream os;
481 os << "<array3d>\n";
482 for (size_t kk=0; kk<nz; kk++) {
483 os << "depth[" << kk << "]\n";
484 for (size_t jj=0; jj<ny; jj++) {
485 for (size_t ii=0; ii<nx; ii++) {
486 os << operator()(ii,jj,kk) << " ";
487 }
488 os << "\n";
489 }
490 }
491 os << "</array3d>\n";
492 return os.str();
493}
494
495template<class P> std::ostream& operator<<(std::ostream& os, const Array3D<P>& array) {
496 os.write(reinterpret_cast<const char*>(&array.keep_nodata), sizeof(array.keep_nodata));
497 os.write(reinterpret_cast<const char*>(&array.nx), sizeof(array.nx));
498 os.write(reinterpret_cast<const char*>(&array.ny), sizeof(array.ny));
499 os.write(reinterpret_cast<const char*>(&array.nz), sizeof(array.nz));
500 os.write(reinterpret_cast<const char*>(&array.vecData[0]), static_cast<std::streamsize>(array.nx*array.ny*array.nz*sizeof(P)));
501 return os;
502}
503
504template<class P> std::istream& operator>>(std::istream& is, Array3D<P>& array) {
505 is.read(reinterpret_cast<char*>(&array.keep_nodata), sizeof(array.keep_nodata));
506 is.read(reinterpret_cast<char*>(&array.nx), sizeof(array.nx));
507 is.read(reinterpret_cast<char*>(&array.ny), sizeof(array.ny));
508 is.read(reinterpret_cast<char*>(&array.nz), sizeof(array.nz));
509 array.nxny = array.nx * array.ny;
510 array.vecData.resize(array.nx*array.ny*array.nz);
511 is.read(reinterpret_cast<char*>(&array.vecData[0]), static_cast<std::streamsize>(array.nx*array.ny*array.nz*sizeof(P))); //30 times faster than assign() or copy()
512 return is;
513}
514
515template<class T> T Array3D<T>::getMin() const {
516
517 T min = std::numeric_limits<T>::max();
518 const size_t nxyz = ny*nx*nz;
519
520 if (keep_nodata==false) {
521 min = *min_element(vecData.begin(), vecData.end());
522 if (min!=std::numeric_limits<T>::max()) return min;
523 else return (T)IOUtils::nodata;
524 } else {
525 for (size_t jj=0; jj<nxyz; jj++) {
526 const T val = vecData[jj];
527 if (val!=IOUtils::nodata && val<min) min=val;
528 }
529 if (min!=std::numeric_limits<T>::max()) return min;
530 else return (T)IOUtils::nodata;
531 }
532}
533
534template<class T> T Array3D<T>::getMax() const {
535
536 T max = -std::numeric_limits<T>::max();
537 const size_t nxyz = ny*nx*nz;
538
539 if (keep_nodata==false) {
540 max = *max_element(vecData.begin(), vecData.end());
541 if (max!=-std::numeric_limits<T>::max()) return max;
542 else return (T)IOUtils::nodata;
543 } else {
544 for (size_t jj=0; jj<nxyz; jj++) {
545 const T val = vecData[jj];
546 if (val!=IOUtils::nodata && val>max) max=val;
547 }
548 if (max!=-std::numeric_limits<T>::max()) return max;
549 else return (T)IOUtils::nodata;
550 }
551}
552
553template<class T> T Array3D<T>::getMean() const {
554
555 T mean = 0;
556 const size_t nxyz = nx*ny*nz;
557
558 if (keep_nodata==false) {
559 if (nxyz>0) return std::accumulate(vecData.begin(), vecData.end(), 0.) / (T)(nxyz);
560 else return (T)IOUtils::nodata;
561 } else {
562 size_t count = 0;
563 for (size_t jj=0; jj<nxyz; jj++) {
564 const T val = vecData[jj];
565 if (val!=IOUtils::nodata) {
566 mean += val;
567 count++;
568 }
569 }
570 if (count>0) return mean/(T)(count);
571 else return (T)IOUtils::nodata;
572 }
573}
574
575template<class T> size_t Array3D<T>::getCount() const
576{
577 const size_t nxyz = nx*ny*nz;
578
579 if (keep_nodata==false) {
580 return (size_t)nxyz;
581 } else {
582 size_t count = 0;
583 for (size_t ii=0; ii<nxyz; ii++) {
584 if (vecData[ii]!=IOUtils::nodata) count++;
585 }
586 return count;
587 }
588}
589
590template<class T> void Array3D<T>::abs() {
591 if (std::numeric_limits<T>::is_signed) {
592 const size_t nxyz = nx*ny*nz;
593 if (keep_nodata==false) {
594 for (size_t jj=0; jj<nxyz; jj++) {
595 T& val = vecData[jj];
596 if (val<0) val=-val;
597 }
598 } else {
599 for (size_t jj=0; jj<nxyz; jj++) {
600 T& val = vecData[jj];
601 if (val<0 && val!=IOUtils::nodata) val=-val;
602 }
603 }
604 }
605}
606
607template<class T> const Array3D<T> Array3D<T>::getAbs() const {
608 Array3D<T> result(*this); //make a copy
609 result.abs(); //already implemented
610
611 return result;
612}
613
614//arithmetic operators
615template<class T> bool Array3D<T>::checkEpsilonEquality(const Array3D<double>& rhs, const double& epsilon) const {
616 if (nx!=rhs.nx || ny!=rhs.ny || nz!=rhs.nz) return false;
617
618 const size_t nxyz = nx*ny*nz;
619 for (size_t jj=0; jj<nxyz; jj++)
620 if (IOUtils::checkEpsilonEquality(vecData[jj], rhs.vecData[jj], epsilon)==false) return false;
621
622 return true;
623}
624
625template<class T> bool Array3D<T>::checkEpsilonEquality(const Array3D<double>& rhs1, const Array3D<double>& rhs2, const double& epsilon) {
626 return rhs1.checkEpsilonEquality(rhs2, epsilon);
627}
628
629template<class T> Array3D<T>& Array3D<T>::operator=(const T& value) {
630 std::fill(vecData.begin(), vecData.end(), value);
631 return *this;
632}
633
634template<class T> Array3D<T>& Array3D<T>::operator+=(const Array3D<T>& rhs)
635{
636 //They have to have equal size
637 if ((rhs.nx != nx) || (rhs.ny != ny) || (rhs.nz != nz)) {
638 std::ostringstream ss;
639 ss << "Trying to add two Array3D objects with different dimensions: ";
640 ss << "(" << nx << "," << ny << "," << nz << ") + (" << rhs.nx << "," << rhs.ny << "," << rhs.nz << ")";
641 throw IOException(ss.str(), AT);
642 }
643 //Add to every single member of the Array3D<T>
644 const size_t nxyz = nx*ny*nz;
645 if (keep_nodata==false) {
646 for (size_t jj=0; jj<nxyz; jj++) {
647 vecData[jj] += rhs(jj);
648 }
649 } else {
650 for (size_t jj=0; jj<nxyz; jj++) {
651 if (vecData[jj]==IOUtils::nodata || rhs(jj)==IOUtils::nodata)
652 vecData[jj] = IOUtils::nodata;
653 else
654 vecData[jj] += rhs(jj);
655 }
656 }
657
658 return *this;
659}
660
661template<class T> const Array3D<T> Array3D<T>::operator+(const Array3D<T>& rhs) const
662{
663 Array3D<T> result(*this); //make a copy
664 result += rhs; //already implemented
665
666 return result;
667}
668
669template<class T> Array3D<T>& Array3D<T>::operator+=(const T& rhs)
670{
671 if (rhs==0.) return *this;
672
673 //Add to every single member of the Array3D<T>
674 const size_t nxyz = nx*ny*nz;
675 if (keep_nodata==false) {
676 for (size_t jj=0; jj<nxyz; jj++) {
677 vecData[jj] += rhs;
678 }
679 } else {
680 for (size_t jj=0; jj<nxyz; jj++) {
681 if (vecData[jj]!=IOUtils::nodata)
682 vecData[jj] += rhs;
683 }
684 }
685
686 return *this;
687}
688
689template<class T> const Array3D<T> Array3D<T>::operator+(const T& rhs) const
690{
691 Array3D<T> result(*this);
692 result += rhs; //already implemented
693
694 return result;
695}
696
697template<class T> Array3D<T>& Array3D<T>::operator-=(const Array3D<T>& rhs)
698{
699 //They have to have equal size
700 if ((rhs.nx != nx) || (rhs.ny != ny) || (rhs.nz != nz)) {
701 std::ostringstream ss;
702 ss << "Trying to substract two Array3D objects with different dimensions: ";
703 ss << "(" << nx << "," << ny << "," << nz << ") - (" << rhs.nx << "," << rhs.ny << "," << rhs.nz << ")";
704 throw IOException(ss.str(), AT);
705 }
706 //Substract to every single member of the Array3D<T>
707 const size_t nxyz = nx*ny*nz;
708 if (keep_nodata==false) {
709 for (size_t jj=0; jj<nxyz; jj++) {
710 vecData[jj] -= rhs(jj);
711 }
712 } else {
713 for (size_t jj=0; jj<nxyz; jj++) {
714 if (vecData[jj]==IOUtils::nodata || rhs(jj)==IOUtils::nodata)
715 vecData[jj] = IOUtils::nodata;
716 else
717 vecData[jj] -= rhs(jj);
718 }
719 }
720
721 return *this;
722}
723
724template<class T> const Array3D<T> Array3D<T>::operator-(const Array3D<T>& rhs) const
725{
726 Array3D<T> result(*this); //make a copy
727 result -= rhs; //already implemented
728
729 return result;
730}
731
732template<class T> Array3D<T>& Array3D<T>::operator-=(const T& rhs)
733{
734 *this += -rhs;
735 return *this;
736}
737
738template<class T> const Array3D<T> Array3D<T>::operator-(const T& rhs) const
739{
740 Array3D<T> result(*this);
741 result += -rhs; //already implemented
742
743 return result;
744}
745
746template<class T> Array3D<T>& Array3D<T>::operator*=(const Array3D<T>& rhs)
747{
748 //They have to have equal size
749 if ((rhs.nx != nx) || (rhs.ny != ny) || (rhs.nz != nz)) {
750 std::ostringstream ss;
751 ss << "Trying to multiply two Array3D objects with different dimensions: ";
752 ss << "(" << nx << "," << ny << "," << nz << ") * (" << rhs.nx << "," << rhs.ny << "," << rhs.nz << ")";
753 throw IOException(ss.str(), AT);
754 }
755 //Multiply every single member of the Array3D<T>
756 const size_t nxyz = nx*ny*nz;
757 if (keep_nodata==false) {
758 for (size_t jj=0; jj<nxyz; jj++) {
759 vecData[jj] *= rhs(jj);
760 }
761 } else {
762 for (size_t jj=0; jj<nxyz; jj++) {
763 if (vecData[jj]==IOUtils::nodata || rhs(jj)==IOUtils::nodata)
764 vecData[jj] = IOUtils::nodata;
765 else
766 vecData[jj] *= rhs(jj);
767 }
768 }
769
770 return *this;
771}
772
773template<class T> const Array3D<T> Array3D<T>::operator*(const Array3D<T>& rhs) const
774{
775 Array3D<T> result(*this); //make a copy
776 result *= rhs; //already implemented
777
778 return result;
779}
780
781template<class T> Array3D<T>& Array3D<T>::operator*=(const T& rhs)
782{
783 if (rhs==1.) return *this;
784
785 //Multiply every single member of the Array3D<T>
786 const size_t nxyz = nx*ny*nz;
787 if (keep_nodata==false) {
788 for (size_t jj=0; jj<nxyz; jj++) {
789 vecData[jj] *= rhs;
790 }
791 } else {
792 for (size_t jj=0; jj<nxyz; jj++) {
793 if (vecData[jj]!=IOUtils::nodata)
794 vecData[jj] *= rhs;
795 }
796 }
797
798 return *this;
799}
800
801template<class T> const Array3D<T> Array3D<T>::operator*(const T& rhs) const
802{
803 Array3D<T> result(*this);
804 result *= rhs; //already implemented
805
806 return result;
807}
808
809template<class T> Array3D<T>& Array3D<T>::operator/=(const Array3D<T>& rhs)
810{
811 //They have to have equal size
812 if ((rhs.nx != nx) || (rhs.ny != ny) || (rhs.nz != nz)) {
813 std::ostringstream ss;
814 ss << "Trying to divide two Array3D objects with different dimensions: ";
815 ss << "(" << nx << "," << ny << "," << nz << ") / (" << rhs.nx << "," << rhs.ny << "," << rhs.nz << ")";
816 throw IOException(ss.str(), AT);
817 }
818 //Divide every single member of the Array3D<T>
819 const size_t nxyz = nx*ny*nz;
820 if (keep_nodata==false) {
821 for (size_t jj=0; jj<nxyz; jj++) {
822 vecData[jj] /= rhs(jj);
823 }
824 } else {
825 for (size_t jj=0; jj<nxyz; jj++) {
826 if (vecData[jj]==IOUtils::nodata || rhs(jj)==IOUtils::nodata)
827 vecData[jj] = IOUtils::nodata;
828 else
829 vecData[jj] /= rhs(jj);
830 }
831 }
832
833 return *this;
834}
835
836template<class T> const Array3D<T> Array3D<T>::operator/(const Array3D<T>& rhs) const
837{
838 Array3D<T> result(*this); //make a copy
839 result /= rhs; //already implemented
840
841 return result;
842}
843
844template<class T> Array3D<T>& Array3D<T>::operator/=(const T& rhs)
845{
846 *this *= (1./rhs);
847 return *this;
848}
849
850template<class T> const Array3D<T> Array3D<T>::operator/(const T& rhs) const
851{
852 Array3D<T> result(*this);
853 result *= (1./rhs); //already implemented
854
855 return result;
856}
857
858template<class T> bool Array3D<T>::operator==(const Array3D<T>& in) const {
859 const size_t in_nx=in.getNx(), in_ny=in.getNy(), in_nz=in.getNz();
860
861 if (nx!=in_nx || ny!=in_ny || nz!=in_nz)
862 return false;
863
864 const size_t nxyz = nx*ny*nz;
865 for (size_t jj=0; jj<nxyz; jj++)
866 if ( !IOUtils::checkEpsilonEquality( vecData[jj] , in.vecData[jj], 1e-6) ) return false;
867
868 return true;
869}
870
871template<class T> bool Array3D<T>::operator!=(const Array3D<T>& in) const {
872 return !(*this==in);
873}
874
875} //end namespace mio
876
877#endif
#define AT
Definition: IOExceptions.h:28
The template class Array2D is a 2D Array (Matrix) able to hold any type of object as datatype....
Definition: Array2D.h:65
size_t getNx() const
Definition: Array2D.h:384
size_t getNy() const
Definition: Array2D.h:388
The template class Array3D is a 3D Array (Tensor) able to hold any type of object as datatype....
Definition: Array3D.h:87
bool operator!=(const Array3D< T > &) const
Operator that tests for inequality.
Definition: Array3D.h:871
static bool checkEpsilonEquality(const Array3D< double > &rhs1, const Array3D< double > &rhs2, const double &epsilon)
Definition: Array3D.h:625
const Array3D< T > operator+(const T &rhs) const
Definition: Array3D.h:689
bool checkEpsilonEquality(const Array3D< double > &rhs, const double &epsilon) const
Definition: Array3D.h:615
void clear()
Definition: Array3D.h:470
size_t nz
Definition: Array3D.h:261
const Array3D< T > operator/(const Array3D< T > &rhs) const
Definition: Array3D.h:836
void fill(const Array3D< T > &i_array3D, const size_t &i_nx, const size_t &i_ny, const size_t &i_nz, const size_t &i_ncols, const size_t &i_nrows, const size_t &i_ndepth)
A method that can be used to insert a subplane into an existing Array3D object that is passed as i_ar...
Definition: Array3D.h:359
T getMean() const
returns the mean value contained in the grid
Definition: Array3D.h:553
void insert(const Array2D< T > &layer, const size_t &depth)
insert a 2D layer in an Array3D object The (nx ,ny) dimensions of the 2D and the 3D arrays must match...
Definition: Array3D.h:388
Array3D< T > & operator-=(const Array3D< T > &rhs)
Definition: Array3D.h:697
size_t getNy() const
Definition: Array3D.h:462
Array3D< T > & operator=(const T &value)
Definition: Array3D.h:629
T getMax() const
returns the maximum value contained in the grid
Definition: Array3D.h:534
const Array3D< T > getAbs() const
returns the grid of the absolute value of values contained in the grid
Definition: Array3D.h:607
bool keep_nodata
Definition: Array3D.h:263
void subset(const Array3D< T > &i_array3D, const size_t &i_nx, const size_t &i_ny, const size_t &i_nz, const size_t &i_ncols, const size_t &i_nrows, const size_t &i_ndepth)
Definition: Array3D.h:324
const Array3D< T > operator-(const T &rhs) const
Definition: Array3D.h:738
Array3D(const Array3D< T > &i_array3D, const size_t &i_nx, const size_t &i_ny, const size_t &i_nz, const size_t &i_ncols, const size_t &i_nrows, const size_t &i_ndepth)
Definition: Array3D.h:316
void size(size_t &anx, size_t &any, size_t &anz) const
Definition: Array3D.h:448
friend std::istream & operator>>(std::istream &is, Array3D< P > &array)
Definition: Array3D.h:504
size_t nxny
Definition: Array3D.h:262
size_t getCount() const
returns the number of points contained in the grid. If setNodataHandling(IOUtils::RAW_NODATA),...
Definition: Array3D.h:575
void abs()
Definition: Array3D.h:590
size_t getNx() const
Definition: Array3D.h:458
Array3D< T > & operator+=(const Array3D< T > &rhs)
Definition: Array3D.h:634
bool empty() const
Definition: Array3D.h:475
Array3DProxy< T > operator[](const size_t &i)
Definition: Array3D.h:308
size_t size() const
Definition: Array3D.h:454
std::vector< T > vecData
The actual objects are stored in a one-dimensional vector.
Definition: Array3D.h:258
size_t getNz() const
Definition: Array3D.h:466
Array3D< T > & operator/=(const T &rhs)
Definition: Array3D.h:844
void resize(const size_t &anx, const size_t &any, const size_t &anz)
Definition: Array3D.h:430
void resize(const size_t &anx, const size_t &any, const size_t &anz, const T &init)
Definition: Array3D.h:439
Array3D< T > & operator-=(const T &rhs)
Definition: Array3D.h:732
size_t nx
Definition: Array3D.h:259
bool getKeepNodata() const
get how to process nodata values (ie: as nodata or as normal numbers)
Definition: Array3D.h:426
const Array3D< T > operator*(const T &rhs) const
Definition: Array3D.h:801
T & operator()(const size_t &i)
Definition: Array3D.h:266
Array3D(const size_t &anx, const size_t &any, const size_t &anz, const T &init)
Definition: Array3D.h:416
const Array3D< T > operator-(const Array3D< T > &rhs) const
Definition: Array3D.h:724
Array3D< T > & operator*=(const Array3D< T > &rhs)
Definition: Array3D.h:746
const Array3D< T > operator*(const Array3D< T > &rhs) const
Definition: Array3D.h:773
void setKeepNodata(const bool i_keep_nodata)
set how to process nodata values (ie: as nodata or as normal numbers)
Definition: Array3D.h:422
const Array3D< T > operator/(const T &rhs) const
Definition: Array3D.h:850
T getMin() const
returns the minimum value contained in the grid
Definition: Array3D.h:515
const std::string toString() const
Definition: Array3D.h:479
Array3D()
Definition: Array3D.h:312
friend std::ostream & operator<<(std::ostream &os, const Array3D< P > &array)
Definition: Array3D.h:495
Array3D(const size_t &anx, const size_t &any, const size_t &anz)
Definition: Array3D.h:410
Array3D< T > & operator*=(const T &rhs)
Definition: Array3D.h:781
void fill(const Array3D< T > &i_array3D, const size_t &i_nx, const size_t &i_ny, const size_t &i_nz)
Definition: Array3D.h:352
bool operator==(const Array3D< T > &) const
Operator that tests for equality.
Definition: Array3D.h:858
const Array3D< T > operator+(const Array3D< T > &rhs) const
Definition: Array3D.h:661
size_t ny
Definition: Array3D.h:260
Array3D< T > & operator+=(const T &rhs)
Definition: Array3D.h:669
Array3D< T > & operator/=(const Array3D< T > &rhs)
Definition: Array3D.h:809
The basic exception class adjusted for the needs of SLF software.
Definition: IOExceptions.h:40
thrown when an index is out of bounds
Definition: IOExceptions.h:106
static const double e
Definition: Meteoconst.h:68
const double nodata
This is the internal nodata value.
Definition: IOUtils.h:75
size_t count(const std::string &input, const std::string &search)
count how many times a substring appears in a string
Definition: IOUtils.cc:258
bool checkEpsilonEquality(const double &val1, const double &val2, const double &epsilon)
Check whether two values are equal regarding a certain epsilon environment (within certain radius of ...
Definition: IOUtils.h:121
Definition: Config.cc:31
std::istream & operator>>(std::istream &is, Config &cfg)
Definition: Config.cc:468
std::ostream & operator<<(std::ostream &os, const Config &cfg)
Definition: Config.cc:433