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