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