MeteoIODoc  MeteoIODoc-2.6.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Array4D.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 ARRAY4D_H
19 #define ARRAY4D_H
20 
21 #include <meteoio/IOUtils.h>
22 #include <meteoio/IOExceptions.h>
23 
24 #include <vector>
25 #include <limits>
26 #include <iostream>
27 
28 namespace mio {
29 
36 template<class T> class Array4D {
37  public:
38  Array4D();
39 
54  Array4D(const Array4D<T>& i_array4D,
55  const size_t& i_nw, const size_t& i_nx, const size_t& i_ny, const size_t& i_nz,
56  const size_t& i_sizeW, const size_t& i_sizeX, const size_t& i_sizeY, const size_t& i_sizeZ);
57 
58  Array4D(const size_t& anw, const size_t& anx, const size_t& any, const size_t& anz);
59 
68  Array4D(const size_t& i_sizeW, const size_t& i_sizeX, const size_t& i_sizeY, const size_t& i_sizeZ, const T& init);
69 
84  void subset(const Array4D<T>& i_array4D,
85  const size_t& i_nw, const size_t& i_nx, const size_t& i_ny, const size_t& i_nz,
86  const size_t& i_sizeW, const size_t& i_sizeX, const size_t& i_sizeY, const size_t& i_sizeZ);
87 
102  void fill(const Array4D<T>& i_array4D,
103  const size_t& i_nw, const size_t& i_nx, const size_t& i_ny, const size_t& i_nz,
104  const size_t& i_sizeW, const size_t& i_sizeX, const size_t& i_sizeY, const size_t& i_sizeZ);
105 
106  void fill(const Array4D<T>& i_array4D, const size_t& i_nw, const size_t& i_nx, const size_t& i_ny, const size_t& i_nz);
107 
113  void setKeepNodata(const bool i_keep_nodata);
114 
119  bool getKeepNodata();
120 
121  void resize(const size_t& anw, const size_t& anx, const size_t& any, const size_t& anz);
122  void resize(const size_t& anw, const size_t& anx, const size_t& any, const size_t& anz, const T& init);
123  void size(size_t& anw, size_t& anx, size_t& any, size_t& anz) const;
124  size_t size() const;
125  size_t getNw() const;
126  size_t getNx() const;
127  size_t getNy() const;
128  size_t getNz() const;
129  void clear();
130  bool empty() const;
131 
136  T getMin() const;
141  T getMax() const;
146  T getMean() const;
153  size_t getCount() const;
158  const Array4D<T> getAbs() const;
159  void abs();
160 
161  const std::string toString() const;
162  template<class P> friend std::iostream& operator<<(std::iostream& os, const Array4D<P>& array);
163  template<class P> friend std::iostream& operator>>(std::iostream& is, Array4D<P>& array);
164 
165  bool checkEpsilonEquality(const Array4D<double>& rhs, const double& epsilon) const;
166  static bool checkEpsilonEquality(const Array4D<double>& rhs1, const Array4D<double>& rhs2, const double& epsilon);
167 
168  T& operator ()(const size_t& i);
169  const T operator ()(const size_t& i) const;
170  T& operator ()(const size_t& w, const size_t& x, const size_t& y, const size_t& z);
171  const T operator ()(const size_t& w, const size_t& x, const size_t& y, const size_t& z) const;
172 
174  Array4D<T>& operator =(const T& value);
175 
176  Array4D<T>& operator+=(const T& rhs);
177  const Array4D<T> operator+(const T& rhs) const;
178  Array4D<T>& operator+=(const Array4D<T>& rhs);
179  const Array4D<T> operator+(const Array4D<T>& rhs) const;
180 
181  Array4D<T>& operator-=(const T& rhs);
182  const Array4D<T> operator-(const T& rhs) const;
183  Array4D<T>& operator-=(const Array4D<T>& rhs);
184  const Array4D<T> operator-(const Array4D<T>& rhs) const;
185 
186  Array4D<T>& operator*=(const T& rhs);
187  const Array4D<T> operator*(const T& rhs) const;
188  Array4D<T>& operator*=(const Array4D<T>& rhs);
189  const Array4D<T> operator*(const Array4D<T>& rhs) const;
190 
191  Array4D<T>& operator/=(const T& rhs);
192  const Array4D<T> operator/(const T& rhs) const;
193  Array4D<T>& operator/=(const Array4D<T>& rhs);
194  const Array4D<T> operator/(const Array4D<T>& rhs) const;
195 
196  bool operator==(const Array4D<T>&) const;
197  bool operator!=(const Array4D<T>&) const;
198 
199  protected:
200  std::vector<T> vecData;
201  size_t nw;
202  size_t nx;
203  size_t ny;
204  size_t nz;
205  size_t nwnx; //xw time nx
206  size_t nwnxny; //xw time nx times ny
208 };
209 
210 template<class T> inline T& Array4D<T>::operator()(const size_t& i) {
211 #ifndef NOSAFECHECKS
212  return vecData.at(i);
213 #else
214  return vecData[i];
215 #endif
216 }
217 
218 template<class T> inline const T Array4D<T>::operator()(const size_t& i) const {
219 #ifndef NOSAFECHECKS
220  return vecData.at(i);
221 #else
222  return vecData[i];
223 #endif
224 }
225 
226 template<class T> inline T& Array4D<T>::operator()(const size_t& w, const size_t& x, const size_t& y, const size_t& z) {
227 #ifndef NOSAFECHECKS
228  if ((w >= nw) || (x >= nx) || (y >= ny) || (z >= nz)) {
229  std::stringstream ss;
230  ss << "Trying to access array(" << w << "," << x << "," << y << "," << z << ")";
231  ss << " while array is ("<< nw << "," << nx << "," << ny << "," << nz << ")";
232  throw IndexOutOfBoundsException(ss.str(), AT);
233  }
234 #endif
235 
236  //ROW-MAJOR alignment of the vector: fully C-compatible memory layout
237  return vecData[w + x*nw + y*nwnx + z*nwnxny];
238 }
239 
240 template<class T> inline const T Array4D<T>::operator()(const size_t& w, const size_t& x, const size_t& y, const size_t& z) const {
241 #ifndef NOSAFECHECKS
242  if ((w >= nw) || (x >= nx) || (y >= ny) || (z >= nz)) {
243  std::stringstream ss;
244  ss << "Trying to access array("<< w << "," << x << "," << y << "," << z << ")";
245  ss << " while array is ("<< nw << "," << nx << "," << ny << "," << nz << ")";
246  throw IndexOutOfBoundsException(ss.str(), AT);
247  }
248 #endif
249  return vecData[w + nw*x + y*nwnx + z*nwnxny];
250 }
251 
252 template<class T> Array4D<T>::Array4D() : vecData(), nw(0), nx(0), ny(0), nz(0), nwnx(0), nwnxny(0), keep_nodata(true)
253 {
254 }
255 
256 template<class T> Array4D<T>::Array4D(const Array4D<T>& i_array4D,
257  const size_t& i_nw,const size_t& i_nx, const size_t& i_ny, const size_t& i_nz,
258  const size_t& i_sizeW, const size_t& i_sizeX, const size_t& i_sizeY, const size_t& i_sizeZ)
259  : vecData(i_sizeW*i_sizeX*i_sizeY*i_sizeZ), nw(i_sizeW), nx(i_sizeX), ny(i_sizeY), nz(i_sizeZ), nwnx(i_sizeW*i_sizeX), nwnxny(i_sizeW*i_sizeX*i_sizeY), keep_nodata(true)
260 {
261  subset(i_array4D, i_nw, i_nx, i_ny, i_nz, i_sizeW, i_sizeX, i_sizeY, i_sizeZ);
262 }
263 
264 template<class T> void Array4D<T>::subset(const Array4D<T>& i_array4D,
265  const size_t& i_nw, const size_t& i_nx, const size_t& i_ny, const size_t& i_nz,
266  const size_t& i_sizeW, const size_t& i_sizeX, const size_t& i_sizeY, const size_t& i_sizeZ)
267 {
268 
269  if (((i_nw+i_sizeW) > i_array4D.nw) || ((i_nx+i_sizeX) > i_array4D.nx) || ((i_ny+i_sizeY) > i_array4D.ny) || ((i_nz+i_sizeZ) > i_array4D.nz)) {
270  std::stringstream ss;
271  ss << "Trying to cut an array of size ("<< i_array4D.nw << "," << i_array4D.nx << "," << i_array4D.ny << "," << i_array4D.nz << ") ";
272  ss << "to size (" << i_sizeW << "," << i_sizeX << "," << i_sizeY << "," << i_sizeZ << ") ";
273  ss << "starting at ("<< i_nw << "," << i_nx << "," << i_ny << "," << i_nz << ")";
274  throw IndexOutOfBoundsException(ss.str(), AT);
275  }
276 
277  if ((i_sizeW == 0) || (i_sizeX == 0) || (i_sizeY == 0) || (i_sizeZ == 0)) //the space has to make sense
278  throw IndexOutOfBoundsException("Trying to cut an array into a null sized array!", AT);
279 
280  resize(i_sizeW, i_sizeX, i_sizeY, i_sizeZ); //create new Array4D object
281  //Copy by value subspace
282  for (size_t ii=0; ii<nz; ii++) {
283  for (size_t jj=0; jj<ny; jj++) {
284  for (size_t kk=0; kk<nx; kk++) {
285  for (size_t ll=0; ll<nw; ll++)
286  //Running through the vector in order of memory alignment
287  operator()(ll,kk,jj,ii) = i_array4D(i_nw+ll, i_nx+kk, i_ny+jj, i_nz+ii);
288  }
289  }
290  }
291 }
292 
293 template<class T> void Array4D<T>::fill(const Array4D<T>& i_array4D, const size_t& i_nw, const size_t& i_nx, const size_t& i_ny, const size_t& i_nz)
294 {
295  size_t i_sizeW, i_sizeX, i_sizeY, i_sizeZ;
296  i_array4D.size(i_sizeW, i_sizeX, i_sizeY, i_sizeZ);
297  fill(i_array4D, i_nw, i_nx, i_ny, i_nz, i_sizeW, i_sizeX, i_sizeY, i_sizeZ);
298 }
299 
300 template<class T> void Array4D<T>::fill(const Array4D<T>& i_array4D,
301  const size_t& i_nw, const size_t& i_nx, const size_t& i_ny, const size_t& i_nz,
302  const size_t& i_sizeW, const size_t& i_sizeX, const size_t& i_sizeY, const size_t& i_sizeZ)
303 {
304 
305  if (((i_nw+i_sizeW) > i_array4D.nw) || ((i_nx+i_sizeX) > i_array4D.nx) || ((i_ny+i_sizeY) > i_array4D.ny) || ((i_nz+i_sizeZ) > i_array4D.nz)) {
306  std::stringstream ss;
307  ss << "Filling an array of size (" << nw << "," << nx << "," << ny << "," << nz << ") ";
308  ss << "with an array of size (" << i_sizeW << "," << i_sizeX << "," << i_sizeY << "," << i_sizeZ << ") ";
309  ss << "starting at (" << i_nw << "," << i_nx << "," << i_ny << "," << i_nz << ")";
310  throw IndexOutOfBoundsException(ss.str(), AT);
311  }
312 
313  if ((i_sizeW == 0) || (i_sizeX == 0) || (i_sizeY == 0) || (i_sizeZ == 0)) //the space has to make sense
314  throw IndexOutOfBoundsException("Filling an array with a null sized array!", AT);
315 
316  //Copy by value subspace
317  for (size_t ii=i_nz; ii<(i_nz+i_sizeZ); ii++) {
318  for (size_t jj=i_ny; jj<(i_ny+i_sizeY); jj++) {
319  for (size_t kk=i_nx; kk<(i_nx+i_sizeX); kk++) {
320  for (size_t ll=i_nw; ll<(i_nw+i_sizeW); ll++) {
321  const size_t iw = ll-i_nw;
322  const size_t ix = kk-i_nx;
323  const size_t iy = jj-i_ny;
324  const size_t iz = ii-i_nz;
325  operator()(ll,kk,jj,ii) = i_array4D(iw,ix, iy, iz);
326  }
327  }
328  }
329  }
330 }
331 
332 template<class T> Array4D<T>::Array4D(const size_t& anw, const size_t& anx, const size_t& any, const size_t& anz)
333  : vecData(anw*anx*any*anz), nw(anw), nx(anx), ny(any), nz(anz), nwnx(anw*anx), nwnxny(anw*anx*any), keep_nodata(true)
334 {
335  //resize(anw, anx, any, anz);
336 }
337 
338 template<class T> Array4D<T>::Array4D(const size_t& anw, const size_t& anx, const size_t& any, const size_t& anz, const T& init)
339  : vecData(anw*anx*any*anz, init), nw(anw), nx(anx), ny(any), nz(anz), nwnx(anw*anx), nwnxny(anw*anx*any), keep_nodata(true)
340 {
341  //resize(anw, anx, any, anz, init);
342 }
343 
344 template<class T> void Array4D<T>::setKeepNodata(const bool i_keep_nodata) {
345  keep_nodata = i_keep_nodata;
346 }
347 
348 template<class T> bool Array4D<T>::getKeepNodata() {
349  return keep_nodata;
350 }
351 
352 template<class T> void Array4D<T>::resize(const size_t& anw, const size_t& anx, const size_t& any, const size_t& anz) {
353  clear(); //we won't be able to "rescue" old values, so we reset the whole vector
354  vecData.resize(anw*anx*any*anz);
355  nw = anw;
356  nx = anx;
357  ny = any;
358  nz = anz;
359  nwnx = nw*nx;
360  nwnxny = nw*nx*ny;
361 }
362 
363 template<class T> void Array4D<T>::resize(const size_t& anw, const size_t& anx, const size_t& any, const size_t& anz, const T& init) {
364  clear(); //we won't be able to "rescue" old values, so we reset the whole vector
365  vecData.resize(anw*anx*any*anz, init);
366  nw = anw;
367  nx = anx;
368  ny = any;
369  nz = anz;
370  nwnx = nw*nx;
371  nwnxny = nw*nx*ny;
372 }
373 
374 template<class T> void Array4D<T>::size(size_t& anw, size_t& anx, size_t& any, size_t& anz) const {
375  anw=nw;
376  anx=nx;
377  any=ny;
378  anz=nz;
379 }
380 
381 template<class T> size_t Array4D<T>::size() const {
382  return nw*nx*ny*nz;
383 }
384 
385 template<class T> size_t Array4D<T>::getNw() const {
386  return nw;
387 }
388 
389 template<class T> size_t Array4D<T>::getNx() const {
390  return nx;
391 }
392 
393 template<class T> size_t Array4D<T>::getNy() const {
394  return ny;
395 }
396 
397 template<class T> size_t Array4D<T>::getNz() const {
398  return nz;
399 }
400 
401 template<class T> void Array4D<T>::clear() {
402  vecData.clear();
403  nw = nx = ny = nz = nwnx = nwnxny = 0;
404 }
405 
406 template<class T> bool Array4D<T>::empty() const {
407  return (nw==0 && nx==0 && ny==0 && nz==0);
408 }
409 
410 template<class T> const std::string Array4D<T>::toString() const {
411  std::stringstream os;
412  os << "<array4d>\n";
413  for (size_t ll=0; ll<nw; ll++) {
414  os << "dim4[" << ll << "]\n";
415  for (size_t kk=0; kk<nz; kk++) {
416  os << "depth[" << kk << "]\n";
417  for (size_t ii=0; ii<nx; ii++) {
418  for (size_t jj=0; jj<ny; jj++) {
419  os << operator()(ii,jj,kk,ll) << " ";
420  }
421  os << "\n";
422  }
423  }
424  }
425  os << "</array4d>\n";
426  return os.str();
427 }
428 
429 template<class P> std::iostream& operator<<(std::iostream& os, const Array4D<P>& array) {
430  os.write(reinterpret_cast<const char*>(&array.keep_nodata), sizeof(array.keep_nodata));
431  os.write(reinterpret_cast<const char*>(&array.nx), sizeof(array.nx));
432  os.write(reinterpret_cast<const char*>(&array.ny), sizeof(array.ny));
433  os.write(reinterpret_cast<const char*>(&array.nz), sizeof(array.nz));
434  os.write(reinterpret_cast<const char*>(&array.nw), sizeof(array.nw));
435  os.write(reinterpret_cast<const char*>(&array.vecData[0]), static_cast<std::streamsize>(array.nx*array.ny*array.nz*array.nw*sizeof(P)));
436  return os;
437 }
438 
439 template<class P> std::iostream& operator>>(std::iostream& is, Array4D<P>& array) {
440  is.read(reinterpret_cast<char*>(&array.keep_nodata), sizeof(array.keep_nodata));
441  is.read(reinterpret_cast<char*>(&array.nx), sizeof(array.nx));
442  is.read(reinterpret_cast<char*>(&array.ny), sizeof(array.ny));
443  is.read(reinterpret_cast<char*>(&array.nz), sizeof(array.nz));
444  is.read(reinterpret_cast<char*>(&array.nw), sizeof(array.nw));
445  array.vecData.resize(array.nx*array.ny*array.nz*array.nw);
446  is.read(reinterpret_cast<char*>(&array.vecData[0]), static_cast<std::streamsize>(array.nx*array.ny*array.nz*array.nw*sizeof(P))); //30 times faster than assign() or copy()
447  return is;
448 }
449 
450 template<class T> T Array4D<T>::getMin() const {
451 
452  T min = std::numeric_limits<T>::max();
453  const size_t nwyz = nwnxny*nz;
454 
455  if (keep_nodata==false) {
456  for (size_t jj=0; jj<nwyz; jj++) {
457  const T val = vecData[jj];
458  if (val<min) min=val;
459  }
460  return min;
461  } else {
462  for (size_t jj=0; jj<nwyz; jj++) {
463  const T val = vecData[jj];
464  if (val!=IOUtils::nodata && val<min) min=val;
465  }
466  if (min!=std::numeric_limits<T>::max()) return min;
467  else return (T)IOUtils::nodata;
468  }
469 }
470 
471 template<class T> T Array4D<T>::getMax() const {
472 
473  T max = -std::numeric_limits<T>::max();
474  const size_t nwyz = nwnxny*nz;
475 
476  if (keep_nodata==false) {
477  for (size_t jj=0; jj<nwyz; jj++) {
478  const T val = vecData[jj];
479  if (val>max) max=val;
480  }
481  return max;
482  } else {
483  for (size_t jj=0; jj<nwyz; jj++) {
484  const T val = vecData[jj];
485  if (val!=IOUtils::nodata && val>max) max=val;
486  }
487  if (max!=-std::numeric_limits<T>::max()) return max;
488  else return (T)IOUtils::nodata;
489  }
490 }
491 
492 template<class T> T Array4D<T>::getMean() const {
493 
494  T mean = 0;
495  const size_t nwyz = nwnxny*nz;
496 
497  if (keep_nodata==false) {
498  for (size_t jj=0; jj<nwyz; jj++) {
499  const T val = vecData[jj];
500  mean += val;
501  }
502  if (nwyz>0) return mean/(T)(nwyz);
503  else return (T)0;
504  } else {
505  size_t count = 0;
506  for (size_t jj=0; jj<nwyz; jj++) {
507  const T val = vecData[jj];
508  if (val!=IOUtils::nodata) {
509  mean += val;
510  count++;
511  }
512  }
513  if (count>0) return mean/(T)(count);
514  else return (T)IOUtils::nodata;
515  }
516 }
517 
518 template<class T> size_t Array4D<T>::getCount() const
519 {
520  const size_t nwyz = nwnxny*nz;
521 
522  if (keep_nodata==false) {
523  return (size_t)nwyz;
524  } else {
525  size_t count = 0;
526  for (size_t ii=0; ii<nwyz; ii++) {
527  if (vecData[ii]!=IOUtils::nodata) count++;
528  }
529  return count;
530  }
531 }
532 
533 template<class T> void Array4D<T>::abs() {
534  if (std::numeric_limits<T>::is_signed) {
535  const size_t nwyz = nwnxny*nz;
536  if (keep_nodata==false) {
537  for (size_t jj=0; jj<nwyz; jj++) {
538  T& val = vecData[jj];
539  if (val<0) val=-val;
540  }
541  } else {
542  for (size_t jj=0; jj<nwyz; jj++) {
543  T& val = vecData[jj];
544  if (val<0 && val!=IOUtils::nodata) val=-val;
545  }
546  }
547  }
548 }
549 
550 
551 template<class T> const Array4D<T> Array4D<T>::getAbs() const {
552  Array4D<T> result(*this); //make a copy
553  result.abs(); //already implemented
554 
555  return result;
556 }
557 
558 //arithmetic operators
559 template<class T> bool Array4D<T>::checkEpsilonEquality(const Array4D<double>& rhs, const double& epsilon) const {
560  if (nw!=rhs.nw || nx!=rhs.nx || ny!=rhs.ny || nz!=rhs.nz) return false;
561 
562  const size_t nwyz = nwnxny*nz;
563  for (size_t jj=0; jj<nwyz; jj++)
564  if (IOUtils::checkEpsilonEquality(vecData[jj], rhs.vecData[jj], epsilon)==false) return false;
565 
566  return true;
567 }
568 
569 template<class T> bool Array4D<T>::checkEpsilonEquality(const Array4D<double>& rhs1, const Array4D<double>& rhs2, const double& epsilon) {
570  return rhs1.checkEpsilonEquality(rhs2, epsilon);
571 }
572 
573 template<class T> Array4D<T>& Array4D<T>::operator=(const Array4D<T>& source) {
574  if (this != &source) {
575  keep_nodata = source.keep_nodata;
576  nw = source.nw;
577  nx = source.nx;
578  ny = source.ny;
579  nz = source.nz;
580  nwnx = source.nwnx;
581  nwnxny = source.nwnxny;
582  vecData = source.vecData;
583  }
584  return *this;
585 }
586 
587 template<class T> Array4D<T>& Array4D<T>::operator=(const T& value) {
588  std::fill(vecData.begin(), vecData.end(), value);
589  return *this;
590 }
591 
592 template<class T> Array4D<T>& Array4D<T>::operator+=(const Array4D<T>& rhs)
593 {
594  //They have to have equal size
595  if ((rhs.nw != nw) || (rhs.nx != nx) || (rhs.ny != ny) || (rhs.nz != nz)) {
596  std::stringstream ss;
597  ss << "Trying to add two Array4D objects with different dimensions: ";
598  ss << "(" << nw << "," << nx << "," << ny << "," << nz << ") + (" << rhs.nw << "," << rhs.nx << "," << rhs.ny << "," << rhs.nz << ")";
599  throw IOException(ss.str(), AT);
600  }
601  //Add to every single member of the Array4D<T>
602  const size_t nwyz = nwnxny*nz;
603  if (keep_nodata==false) {
604  for (size_t jj=0; jj<nwyz; jj++) {
605  vecData[jj] += rhs(jj);
606  }
607  } else {
608  for (size_t jj=0; jj<nwyz; jj++) {
609  if (vecData[jj]==IOUtils::nodata || rhs(jj)==IOUtils::nodata)
610  vecData[jj] = IOUtils::nodata;
611  else
612  vecData[jj] += rhs(jj);
613  }
614  }
615 
616  return *this;
617 }
618 
619 template<class T> const Array4D<T> Array4D<T>::operator+(const Array4D<T>& rhs) const
620 {
621  Array4D<T> result(*this); //make a copy
622  result += rhs; //already implemented
623 
624  return result;
625 }
626 
627 template<class T> Array4D<T>& Array4D<T>::operator+=(const T& rhs)
628 {
629  if (rhs==0.) return *this;
630 
631  //Add to every single member of the Array4D<T>
632  const size_t nwyz = nwnxny*nz;
633  if (keep_nodata==false) {
634  for (size_t jj=0; jj<nwyz; jj++) {
635  vecData[jj] += rhs;
636  }
637  } else {
638  for (size_t jj=0; jj<nwyz; jj++) {
639  if (vecData[jj]!=IOUtils::nodata)
640  vecData[jj] += rhs;
641  }
642  }
643 
644  return *this;
645 }
646 
647 template<class T> const Array4D<T> Array4D<T>::operator+(const T& rhs) const
648 {
649  Array4D<T> result(*this);
650  result += rhs; //already implemented
651 
652  return result;
653 }
654 
655 template<class T> Array4D<T>& Array4D<T>::operator-=(const Array4D<T>& rhs)
656 {
657  //They have to have equal size
658  if ((rhs.nw != nw) || (rhs.nx != nx) || (rhs.ny != ny) || (rhs.nz != nz)) {
659  std::stringstream ss;
660  ss << "Trying to substract two Array4D objects with different dimensions: ";
661  ss << "(" << nw << "," << nx << "," << ny << "," << nz << ") - (" << rhs.nw << "," << rhs.nx << "," << rhs.ny << "," << rhs.nz << ")";
662  throw IOException(ss.str(), AT);
663  }
664  //Substract to every single member of the Array4D<T>
665  const size_t nwyz = nwnxny*nz;
666  if (keep_nodata==false) {
667  for (size_t jj=0; jj<nwyz; jj++) {
668  vecData[jj] -= rhs(jj);
669  }
670  } else {
671  for (size_t jj=0; jj<nwyz; jj++) {
672  if (vecData[jj]==IOUtils::nodata || rhs(jj)==IOUtils::nodata)
673  vecData[jj] = IOUtils::nodata;
674  else
675  vecData[jj] -= rhs(jj);
676  }
677  }
678 
679  return *this;
680 }
681 
682 template<class T> const Array4D<T> Array4D<T>::operator-(const Array4D<T>& rhs) const
683 {
684  Array4D<T> result(*this); //make a copy
685  result -= rhs; //already implemented
686 
687  return result;
688 }
689 
690 template<class T> Array4D<T>& Array4D<T>::operator-=(const T& rhs)
691 {
692  *this += -rhs;
693  return *this;
694 }
695 
696 template<class T> const Array4D<T> Array4D<T>::operator-(const T& rhs) const
697 {
698  Array4D<T> result(*this);
699  result += -rhs; //already implemented
700 
701  return result;
702 }
703 
704 template<class T> Array4D<T>& Array4D<T>::operator*=(const Array4D<T>& rhs)
705 {
706  //They have to have equal size
707  if ((rhs.nw != nw) || (rhs.nx != nx) || (rhs.ny != ny) || (rhs.nz != nz)) {
708  std::stringstream ss;
709  ss << "Trying to multiply two Array4D objects with different dimensions: ";
710  ss << "("<< nw << "," << nx << "," << ny << "," << nz << ") * (" << rhs.nw << "," << rhs.nx << "," << rhs.ny << "," << rhs.nz << ")";
711  throw IOException(ss.str(), AT);
712  }
713  //Multiply every single member of the Array4D<T>
714  const size_t nwxyz = nwnxny*nz;
715  if (keep_nodata==false) {
716  for (size_t jj=0; jj<nwxyz; jj++) {
717  vecData[jj] *= rhs(jj);
718  }
719  } else {
720  for (size_t jj=0; jj<nwxyz; jj++) {
721  if (vecData[jj]==IOUtils::nodata || rhs(jj)==IOUtils::nodata)
722  vecData[jj] = IOUtils::nodata;
723  else
724  vecData[jj] *= rhs(jj);
725  }
726  }
727 
728  return *this;
729 }
730 
731 template<class T> const Array4D<T> Array4D<T>::operator*(const Array4D<T>& rhs) const
732 {
733  Array4D<T> result(*this); //make a copy
734  result *= rhs; //already implemented
735 
736  return result;
737 }
738 
739 template<class T> Array4D<T>& Array4D<T>::operator*=(const T& rhs)
740 {
741  if (rhs==1.) return *this;
742 
743  //Multiply every single member of the Array4D<T>
744  const size_t nwxyz = nwnxny*nz;
745  if (keep_nodata==false) {
746  for (size_t jj=0; jj<nwxyz; jj++) {
747  vecData[jj] *= rhs;
748  }
749  } else {
750  for (size_t jj=0; jj<nwxyz; jj++) {
751  if (vecData[jj]!=IOUtils::nodata)
752  vecData[jj] *= rhs;
753  }
754  }
755 
756  return *this;
757 }
758 
759 template<class T> const Array4D<T> Array4D<T>::operator*(const T& rhs) const
760 {
761  Array4D<T> result(*this);
762  result *= rhs; //already implemented
763 
764  return result;
765 }
766 
767 template<class T> Array4D<T>& Array4D<T>::operator/=(const Array4D<T>& rhs)
768 {
769  //They have to have equal size
770  if ((rhs.nw != nw) || (rhs.nx != nx) || (rhs.ny != ny) || (rhs.nz != nz)) {
771  std::stringstream ss;
772  ss << "Trying to divide two Array4D objects with different dimensions: ";
773  ss << "(" << nw << "," << nx << "," << ny << "," << nz << ") / (" << rhs.nw << "," << rhs.nx << "," << rhs.ny << "," << rhs.nz << ")";
774  throw IOException(ss.str(), AT);
775  }
776  //Divide every single member of the Array4D<T>
777  const size_t nwxyz = nwnxny*nz;
778  if (keep_nodata==false) {
779  for (size_t jj=0; jj<nwxyz; jj++) {
780  vecData[jj] /= rhs(jj);
781  }
782  } else {
783  for (size_t jj=0; jj<nwxyz; jj++) {
784  if (vecData[jj]==IOUtils::nodata || rhs(jj)==IOUtils::nodata)
785  vecData[jj] = IOUtils::nodata;
786  else
787  vecData[jj] /= rhs(jj);
788  }
789  }
790 
791  return *this;
792 }
793 
794 template<class T> const Array4D<T> Array4D<T>::operator/(const Array4D<T>& rhs) const
795 {
796  Array4D<T> result(*this); //make a copy
797  result /= rhs; //already implemented
798 
799  return result;
800 }
801 
802 template<class T> Array4D<T>& Array4D<T>::operator/=(const T& rhs)
803 {
804  *this *= (1./rhs);
805  return *this;
806 }
807 
808 template<class T> const Array4D<T> Array4D<T>::operator/(const T& rhs) const
809 {
810  Array4D<T> result(*this);
811  result *= (1./rhs); //already implemented
812 
813  return result;
814 }
815 
816 template<class T> bool Array4D<T>::operator==(const Array4D<T>& in) const {
817  const size_t in_nx=in.getNx(), in_ny=in.getNy(), in_nz=in.getNz(), in_nw=in.getNw();
818 
819  if (nx!=in_nx || ny!=in_ny || nz!=in_nz || nw!=in_nw)
820  return false;
821 
822  const size_t nwxyz = nx*ny*nz*nw;
823  for (size_t jj=0; jj<nwxyz; jj++)
824  if ( !IOUtils::checkEpsilonEquality( vecData[jj] , in.vecData[jj], 1e-6) ) return false;
825 
826  return true;
827 }
828 
829 template<class T> bool Array4D<T>::operator!=(const Array4D<T>& in) const {
830  return !(*this==in);
831 }
832 
833 } //end namespace mio
834 
835 #endif
void resize(const size_t &anw, const size_t &anx, const size_t &any, const size_t &anz)
Definition: Array4D.h:352
bool operator==(const Array4D< T > &) const
Operator that tests for equality.
Definition: Array4D.h:816
bool keep_nodata
Definition: Array4D.h:207
bool getKeepNodata()
get how to process nodata values (ie: as nodata or as normal numbers)
Definition: Array4D.h:348
void abs()
Definition: Array4D.h:533
void size(size_t &anw, size_t &anx, size_t &any, size_t &anz) const
Definition: Array4D.h:374
size_t getNw() const
Definition: Array4D.h:385
bool checkEpsilonEquality(const Array4D< double > &rhs, const double &epsilon) const
Definition: Array4D.h:559
size_t size() const
Definition: Array4D.h:381
bool operator!=(const Array4D< T > &) const
Operator that tests for inequality.
Definition: Array4D.h:829
Array4D< T > & operator*=(const T &rhs)
Definition: Array4D.h:739
size_t nx
Definition: Array4D.h:202
T & operator()(const size_t &i)
Definition: Array4D.h:210
Array4D< T > & operator-=(const T &rhs)
Definition: Array4D.h:690
size_t nw
Definition: Array4D.h:201
size_t getNz() const
Definition: Array4D.h:397
Array4D< T > & operator=(const Array4D< T > &)
Definition: Array4D.h:573
const Array4D< T > getAbs() const
returns the grid of the absolute value of values contained in the grid
Definition: Array4D.h:551
const Array4D< T > operator*(const T &rhs) const
Definition: Array4D.h:759
const Array4D< T > operator+(const T &rhs) const
Definition: Array4D.h:647
Array4D()
Definition: Array4D.h:252
std::vector< T > vecData
The actual objects are stored in a one-dimensional vector.
Definition: Array4D.h:200
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
size_t nwnxny
Definition: Array4D.h:206
The template class Array4D is a 4D Array (Tensor) able to hold any type of object as datatype...
Definition: Array4D.h:36
size_t getNx() const
Definition: Array4D.h:389
const Array4D< T > operator/(const T &rhs) const
Definition: Array4D.h:808
void fill(const Array4D< T > &i_array4D, const size_t &i_nw, const size_t &i_nx, const size_t &i_ny, const size_t &i_nz, const size_t &i_sizeW, const size_t &i_sizeX, const size_t &i_sizeY, const size_t &i_sizeZ)
A method that can be used to insert a subplane into an existing Array4D object that is passed as i_ar...
Definition: Array4D.h:300
T getMin() const
returns the minimum value contained in the grid
Definition: Array4D.h:450
bool empty() const
Definition: Array4D.h:406
size_t nz
Definition: Array4D.h:204
const std::string toString() const
Definition: Array4D.h:410
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: Array4D.h:518
#define AT
Definition: IOExceptions.h:29
friend std::iostream & operator>>(std::iostream &is, Array4D< P > &array)
Definition: Array4D.h:439
const double nodata
This is the internal nodata value.
Definition: IOUtils.h:60
size_t getNy() const
Definition: Array4D.h:393
void clear()
Definition: Array4D.h:401
T getMax() const
returns the maximum value contained in the grid
Definition: Array4D.h:471
size_t nwnx
Definition: Array4D.h:205
Array4D< T > & operator+=(const T &rhs)
Definition: Array4D.h:627
size_t ny
Definition: Array4D.h:203
thrown when an index is out of bounds
Definition: IOExceptions.h:108
const double e
Definition: Meteoconst.h:66
Array4D< T > & operator/=(const T &rhs)
Definition: Array4D.h:802
void setKeepNodata(const bool i_keep_nodata)
set how to process nodata values (ie: as nodata or as normal numbers)
Definition: Array4D.h:344
void subset(const Array4D< T > &i_array4D, const size_t &i_nw, const size_t &i_nx, const size_t &i_ny, const size_t &i_nz, const size_t &i_sizeW, const size_t &i_sizeX, const size_t &i_sizeY, const size_t &i_sizeZ)
Definition: Array4D.h:264
T getMean() const
returns the mean value contained in the grid
Definition: Array4D.h:492
The basic exception class adjusted for the needs of SLF software.
Definition: IOExceptions.h:41
std::iostream & operator>>(std::iostream &is, Config &cfg)
Definition: Config.cc:141
value::array array
Definition: picojson.h:194
const Array4D< T > operator-(const T &rhs) const
Definition: Array4D.h:696