21#ifndef NRLIB_REGULARSURFACE_HPP
22#define NRLIB_REGULARSURFACE_HPP
30#include "../grid/grid2d.hpp"
31#include "../iotools/stringtools.hpp"
43 RegularSurface(
double x0,
double y0,
double lx,
double ly,
size_t nx,
size_t ny,
44 const A&
value = A());
63 A
GetZ(
double x,
double y)
const;
73 double y_min,
double y_max)
const;
78 for (i=this->
begin();i<this->
end();i++)
84 for (i=this->
begin();i<this->
end();i++)
85 if (*i != missing_val_)
91 for (i=this->
begin();i<this->
end();i++)
92 if (*i != missing_val_)
98 for (i=this->
begin();i<this->
end();i++)
99 if (*i != missing_val_)
112 A
MinNode(
size_t &i,
size_t &j)
const;
114 A
MaxNode(
size_t &i,
size_t &j)
const;
117 A
Avg(
int& n_nodes)
const;
134 inline size_t FindI(
double x)
const;
138 inline size_t FindJ(
double y)
const;
143 inline void FindIndex(
double x,
double y,
size_t& i,
size_t& j)
const;
146 void FindContIndex(
double x,
double y,
double& i,
double& j)
const;
160 {
return x_min_ + i*dx_; }
163 {
return y_min_ + j*dy_; }
165 void GetXY(
size_t i,
size_t j,
double &
x,
double & y)
const {
178 this->
GetIJ(index, i, j);
182 void GetNode(
size_t i,
size_t j,
double&
x,
double& y,
double&
z)
const {
189 double GetXMax()
const {
return x_min_ + lx_; }
190 double GetYMax()
const {
return y_min_ + ly_; }
191 double GetDX()
const {
return dx_; }
192 double GetDY()
const {
return dy_; }
197 double lx,
double ly);
200 void Resize(
size_t ni,
size_t nj,
const A& val = A());
204 return val == missing_val_;
214 return (*
this)(i, j) == missing_val_;
219 (*this)(i, j) = missing_val_;
243 double CellRelX(
double x)
const;
244 double CellRelY(
double y)
const;
260 inline static A GetBilinearZ(A corners[4],
double u,
double v);
262 inline A GetBilinearZMissingCorners(A corners[4],
double u,
double v)
const;
276 missing_val_(static_cast<A>(-999.0))
282 double lx,
double ly,
283 size_t nx,
size_t ny,
290 missing_val_(static_cast<A>(-999.0))
292 dx_ = (nx > 1) ? lx / (nx-1) : 1.0;
293 dy_ = (ny > 1) ? ly / (ny-1) : 1.0;
299 double lx,
double ly,
306 missing_val_(static_cast<A>(-999.0))
308 size_t ni = grid.
GetNI();
309 size_t nj = grid.
GetNJ();
310 dx_ = (ni > 1) ? lx / (ni-1) : 1.0;
311 dy_ = (nj > 1) ? ly / (nj-1) : 1.0;
320 missing_val_(static_cast<A>(-999.0))
331 A
b = corners[1] - a;
332 A c = corners[2] - a;
333 A d = corners[3] - a -
b - c;
335 return static_cast<A
>(a +
b*u + c*v + d*u*v);
340A RegularSurface<A>::GetBilinearZMissingCorners(A corners[4],
double u,
double v)
const
344 weights[0] = (1.0-u) * (1.0-v);
345 weights[1] = u * (1.0 - v);
346 weights[2] = (1.0-u) * v;
350 double sum_weights = 0.0;
352 for (
size_t i = 0; i < 4; ++i) {
353 if (!IsMissing(corners[i])) {
354 sum += weights[i] * corners[i];
355 sum_weights += weights[i];
359 if (sum_weights == 0.0) {
360 return GetMissingValue();
362 return static_cast<A
>(sum / sum_weights);
370 GetCorners(
x, y, corner);
374 for (
int pix = 0; pix < 4; pix++)
376 if (IsMissing(corner[pix]))
380 if (n_missing == 4) {
381 return(missing_val_);
384 double x1 = CellRelX(
x);
385 double y1 = CellRelY(y);
387 if (n_missing == 0) {
388 return GetBilinearZ(corner, x1, y1);
391 return GetBilinearZMissingCorners(corner, x1, y1);
401 GetCorners(
x, y, corner);
403 bool missing =
false;
406 while (pix < 4 && !missing) {
407 if (IsMissing(corner[pix]))
413 return(missing_val_);
415 double x1 = CellRelX(
x);
416 double y1 = CellRelY(y);
417 return GetBilinearZ(corner, x1, y1);
425 if (
x < GetXMin() ||
x > GetXMax() || y < GetYMin() || y > GetYMax())
434 for (
size_t i = 0; i < this->GetN(); i++) {
435 if ((*
this)(i) != missing_val_) {
442 (*
this)(i) = missing_val_;
452 for (
size_t i = 0; i < this->GetN(); i++) {
453 if ((*
this)(i) != missing_val_) {
460 (*
this)(i) = missing_val_;
469 for (
size_t i = 0; i < this->GetN(); i++) {
470 if ((*
this)(i) != missing_val_) {
477 (*
this)(i) = missing_val_;
486 for (
size_t i = 0; i < this->GetN(); i++) {
487 if ((*
this)(i) != missing_val_) {
494 (*
this)(i) = missing_val_;
504 A minVal = (*this)(0);
505 typename std::vector<A>::const_iterator i;
506 for (i = this->
begin(); i < this->end(); i++) {
507 if ((IsMissing(minVal) || (*i) < minVal) && IsMissing(*i) ==
false)
516 A minVal = (*this)(0);
517 typename std::vector<A>::const_iterator it;
518 size_t min_index = 0;
519 for (it = this->
begin(); it < this->end(); it++) {
520 if ((IsMissing(minVal) || (*it) < minVal) && IsMissing(*it) ==
false) {
522 min_index =
static_cast<size_t> (it - this->
begin());
525 this->GetIJ(min_index, i, j);
532 A maxVal = (*this)(0);
533 typename std::vector<A>::const_iterator i;
534 for (i = this->
begin(); i < this->end(); i++) {
535 if ((IsMissing(maxVal) || (*i) > maxVal) && IsMissing(*i) ==
false)
544 A maxVal = (*this)(0);
545 typename std::vector<A>::const_iterator it;
546 size_t max_index = 0;
547 for (it = this->
begin(); it < this->end(); it++) {
548 if ((IsMissing(maxVal) || (*it) > maxVal) && IsMissing(*it) ==
false) {
550 max_index =
static_cast<size_t> (it - this->
begin());
553 this->GetIJ(max_index, i, j);
568 A sum =
static_cast<A
>(0);
570 typename std::vector<A>::const_iterator i;
571 for (i = this->
begin(); i < this->end(); i++) {
572 if (!IsMissing(*i)) {
590 double y_min,
double y_max)
const
592 if (x_min < GetXMin() || x_max > GetXMax() ||
593 y_min < GetYMin() || y_max > GetYMax()) {
604 assert (
x >= GetXMin() &&
x <= GetXMax());
605 return static_cast<size_t>( (
x-GetXMin()) / GetDX() );
612 assert (y >= GetYMin() && y <= GetYMax());
613 return static_cast<size_t>( (y-GetYMin()) / GetDY() );
620 assert(
x >= GetXMin() &&
x <= GetXMax() &&
621 y >= GetYMin() && y <= GetYMax());
623 i =
static_cast<size_t>( (
x-GetXMin()) / GetDX() );
624 j =
static_cast<size_t>( (y-GetYMin()) / GetDY() );
631 assert(
x >= GetXMin() &&
x <= GetXMax() &&
632 y >= GetYMin() && y <= GetYMax());
634 i = (
x-GetXMin()) / GetDX() ;
635 j = (y-GetYMin()) / GetDY() ;
642 i =
static_cast<int>(std::floor( (
x-GetXMin()) / GetDX() ));
643 j =
static_cast<int>(std::floor( (y-GetYMin()) / GetDY() ));
651 int ni =
static_cast<int>(this->GetNI());
652 int nj =
static_cast<int>(this->GetNJ());
653 FindGeneralIndex(
x + 0.5*dx_, y + 0.5*dy_, tmp_i, tmp_j);
656 else if (tmp_i >= ni)
661 else if (tmp_j >= nj)
664 i =
static_cast<size_t>(tmp_i);
665 j =
static_cast<size_t>(tmp_j);
675 int ni =
static_cast<int>(this->GetNI());
676 int nj =
static_cast<int>(this->GetNJ());
677 FindGeneralIndex(
x, y, i , j);
679 if (i+1 >= 0 && i < ni && j+1 >= 0 && j < nj) {
680 if (i >= 0 && j >= 0)
681 corners[0] = (*this)(i, j);
683 corners[0] = missing_val_;
684 if (i+1 < ni && j >= 0)
685 corners[1] = (*this)(i+1, j);
687 corners[1] = missing_val_;
688 if (i >= 0 && j+1 < nj)
689 corners[2] = (*this)(i, j+1);
691 corners[2] = missing_val_;
692 if (i+1 < ni && j+1 < nj)
693 corners[3] = (*this)(i+1, j+1);
695 corners[3] = missing_val_;
698 corners[0] = missing_val_;
699 corners[1] = missing_val_;
700 corners[2] = missing_val_;
701 corners[3] = missing_val_;
713 if (!IsMissingAt(i-1,j-1)) {
714 sum += (*this)(i-1,j-1);
718 if (j < this->GetNJ() - 1) {
719 if (!IsMissingAt(i-1,j+1)) {
720 sum += (*this)(i-1,j+1);
724 if (!IsMissingAt(i-1,j)) {
725 sum += (*this)(i-1,j);
729 if (i < this->GetNI() - 1) {
731 if (!IsMissingAt(i+1,j-1)) {
732 sum += (*this)(i+1,j-1);
736 if (j < this->GetNJ() - 1) {
737 if (!IsMissingAt(i+1,j+1)) {
738 sum += (*this)(i+1,j+1);
742 if (!IsMissingAt(i+1,j)) {
743 sum += (*this)(i+1,j);
749 if (!IsMissingAt(i,j-1)) {
750 sum += (*this)(i,j-1);
754 if (j < this->GetNJ() - 1) {
755 if (!IsMissingAt(i,j+1)) {
756 sum += (*this)(i,j+1);
761 double avg = missing_val_;
762 bool non_missing =
false;
769 return (non_missing);
775 double lx,
double ly)
783 dx_ = (ni > 1) ? lx / (ni - 1) : 1.0;
784 dy_ = (nj > 1) ? ly / (nj - 1) : 1.0;
792 dx_ = (nx > 1) ? lx_ / (nx-1) : 1.0;
793 dy_ = (ny > 1) ? ly_ / (ny-1) : 1.0;
804 throw FileFormatError(
"Failed to determine file format for surface file: " + filename);
822 +
" as non-rotated grid is currently not supported.");
827 +
". Rotated grids are not supported.");
844 throw FileFormatError(
"Writing of surface to file " + filename +
" on format "
854 return std::modf((
x-GetXMin())/GetDX(), &i);
859double RegularSurface<A>::CellRelY(
double y)
const
862 return std::modf((y-GetYMin())/GetDY(), &j);
870 std::swap(x_min_, other.x_min_);
871 std::swap(y_min_, other.y_min_);
872 std::swap(lx_, other.lx_);
873 std::swap(ly_, other.ly_);
874 std::swap(dx_, other.dx_);
875 std::swap(dy_, other.dy_);
876 std::swap(name_, other.name_);
877 std::swap(missing_val_, other.missing_val_);
const cJSON *const b
Definition: cJSON.h:251
const char *const name
Definition: cJSON.h:258
int index
Definition: cJSON.h:168
int count
Definition: cJSON.h:212
char const int const cJSON_bool format
Definition: cJSON.h:161
const char *const string
Definition: cJSON.h:170
Definition: grid2d.hpp:32
virtual void Resize(size_t ni, size_t nj, const A &val=A())
Definition: grid2d.hpp:101
std::vector< A >::iterator iterator
Definition: grid2d.hpp:34
void Swap(Grid2D< A > &other)
Definition: grid2d.hpp:174
iterator end()
Definition: grid2d.hpp:56
void GetIJ(size_t index, size_t &i, size_t &j) const
Definition: grid2d.hpp:153
size_t GetNI() const
Definition: grid2d.hpp:61
iterator begin()
Definition: grid2d.hpp:55
size_t GetNJ() const
Definition: grid2d.hpp:62
Definition: regularsurface.hpp:39
A GetZInside(double x, double y) const
Definition: regularsurface.hpp:398
void GetXY(size_t i, size_t j, double &x, double &y) const
Definition: regularsurface.hpp:165
void FindContIndex(double x, double y, double &i, double &j) const
Definition: regularsurface.hpp:629
bool IsMissing(A val) const
Check if grid value is missing.
Definition: regularsurface.hpp:203
size_t FindJ(double y) const
Definition: regularsurface.hpp:610
bool MultiplyNonConform(const Surface< A > *s2)
Definition: regularsurface.hpp:467
bool IsInsideSurface(double x, double y) const
Checks if point is inside definition area for surface.
Definition: regularsurface.hpp:423
void Subtract(A c)
Definition: regularsurface.hpp:89
double GetDX() const
Definition: regularsurface.hpp:191
void SetDimensions(double x_min, double y_min, double lx, double ly)
Definition: regularsurface.hpp:774
double GetXMin() const
Definition: regularsurface.hpp:187
bool DivideNonConform(const Surface< A > *s2)
Definition: regularsurface.hpp:484
double GetYMin() const
Definition: regularsurface.hpp:188
Surface< A > * Clone() const
Generate a copy of the underlying object.
Definition: regularsurface.hpp:57
bool CreateNeighbourAvg(size_t i, size_t j)
Definition: regularsurface.hpp:707
double GetDY() const
Definition: regularsurface.hpp:192
size_t FindI(double x) const
Definition: regularsurface.hpp:602
void Add(A c)
Definition: regularsurface.hpp:82
bool EnclosesRectangle(double x_min, double x_max, double y_min, double y_max) const
Definition: regularsurface.hpp:589
A GetMissingValue() const
Definition: regularsurface.hpp:222
void SetMissingValue(A missing_val)
Definition: regularsurface.hpp:223
A Avg() const
Definition: regularsurface.hpp:558
void ReadFromFile(const std::string &filename, SurfaceFileFormat format=SURF_UNKNOWN)
Read surface file on given format.
Definition: regularsurface.hpp:798
A GetZ(double x, double y) const
Definition: regularsurface.hpp:367
void SetMissing(size_t i, size_t j)
Set missing.
Definition: regularsurface.hpp:218
double GetXMax() const
Definition: regularsurface.hpp:189
void Resize(size_t ni, size_t nj, const A &val=A())
Resize grid. Overrides Grid2D's resize.
Definition: regularsurface.hpp:789
void WriteToFile(const std::string &filename, SurfaceFileFormat format=SURF_STORM_BINARY) const
Definition: regularsurface.hpp:833
void GetCorners(double x, double y, A corners[4]) const
Definition: regularsurface.hpp:670
A Max() const
Definition: regularsurface.hpp:530
double GetYMax() const
Definition: regularsurface.hpp:190
bool AddNonConform(const Surface< A > *s2)
Definition: regularsurface.hpp:432
void FindNearestNodeIndex(double x, double y, size_t &i, size_t &j) const
Definition: regularsurface.hpp:648
double GetLengthY() const
Definition: regularsurface.hpp:194
double GetY(size_t j) const
Definition: regularsurface.hpp:162
RegularSurface()
Definition: regularsurface.hpp:268
void SetName(const std::string &name)
Definition: regularsurface.hpp:226
double GetLengthX() const
Definition: regularsurface.hpp:193
void FindGeneralIndex(double x, double y, int &i, int &j) const
Definition: regularsurface.hpp:640
void Swap(RegularSurface< A > &other)
Definition: regularsurface.hpp:867
void GetNode(size_t index, double &x, double &y, double &z) const
Definition: regularsurface.hpp:176
bool SubtractNonConform(const Surface< A > *s2)
Definition: regularsurface.hpp:450
A MinNode(size_t &i, size_t &j) const
Definition: regularsurface.hpp:514
void FindIndex(double x, double y, size_t &i, size_t &j) const
Definition: regularsurface.hpp:618
A Min() const
Definition: regularsurface.hpp:502
void Assign(A c)
Sets all values on the surface to a constant value.
Definition: regularsurface.hpp:76
A MaxNode(size_t &i, size_t &j) const
Definition: regularsurface.hpp:542
const std::string & GetName() const
Definition: regularsurface.hpp:225
double GetX(size_t i) const
Definition: regularsurface.hpp:159
bool IsMissingAt(size_t i, size_t j) const
Check if grid value is missing.
Definition: regularsurface.hpp:213
void Multiply(A c)
Definition: regularsurface.hpp:96
void GetNode(size_t i, size_t j, double &x, double &y, double &z) const
Definition: regularsurface.hpp:182
void GetXY(size_t index, double &x, double &y) const
Definition: regularsurface.hpp:170
Definition: surface.hpp:29
virtual bool IsMissing(A) const
Definition: surface.hpp:41
virtual A GetZ(double x, double y) const =0
Definition: exception.hpp:31
void WriteIrapClassicAsciiSurf(const RegularSurface< A > &surf, double angle, const std::string &filename)
Definition: surfaceio.hpp:650
void ReadIrapClassicAsciiSurf(const std::string &filename, RegularSurface< A > &surface, double &angle)
Definition: surfaceio.hpp:180
void ReadSgriSurf(const std::string &filename, RegularSurface< A > &surface, double &angle)
Definition: surfaceio.hpp:262
SurfaceFileFormat FindSurfaceFileType(const std::string &filename)
Find type of file.
void WriteStormBinarySurf(const RegularSurface< A > &surf, const std::string &filename)
Definition: surfaceio.hpp:700
SurfaceFileFormat
Definition: surfaceio.hpp:38
@ SURF_STORM_BINARY
Definition: surfaceio.hpp:41
@ SURF_SGRI
Definition: surfaceio.hpp:42
@ SURF_UNKNOWN
Definition: surfaceio.hpp:39
@ SURF_IRAP_CLASSIC_ASCII
Definition: surfaceio.hpp:40
std::string ToString(const T obj, int precision=-99999)
Definition: stringtools.hpp:177
void ReadStormBinarySurf(const std::string &filename, RegularSurface< A > &surface)
Definition: surfaceio.hpp:119
x y * z
Definition: exprtk.hpp:9663
T value(details::expression_node< T > *n)
Definition: exprtk.hpp:12955
x y t t *t x y t t t x y t t t x *y t *t t x *y t *t t x y t t t x y t t t x(y+z)