Ice Sheet System Model  4.18
Code documentation
Public Member Functions | Private Attributes
Observations Class Reference

Declaration of Observations class. More...

#include <Observations.h>

Inheritance diagram for Observations:
DataSet

Public Member Functions

 Observations ()
 
 Observations (IssmDouble *observations_list, IssmDouble *x, IssmDouble *y, int n, Options *options)
 
 ~Observations ()
 
void InitCovertree (IssmDouble *observations_list, IssmDouble *x, IssmDouble *y, int n, Options *options)
 
void InitQuadtree (IssmDouble *observations_list, IssmDouble *x, IssmDouble *y, int n, Options *options)
 
void ClosestObservation (IssmDouble *px, IssmDouble *py, IssmDouble *pobs, IssmDouble x_interp, IssmDouble y_interp, IssmDouble radius)
 
void ClosestObservationCovertree (IssmDouble *px, IssmDouble *py, IssmDouble *pobs, IssmDouble x_interp, IssmDouble y_interp, IssmDouble radius)
 
void ClosestObservationQuadtree (IssmDouble *px, IssmDouble *py, IssmDouble *pobs, IssmDouble x_interp, IssmDouble y_interp, IssmDouble radius)
 
void Distances (IssmPDouble *distances, IssmPDouble *x, IssmPDouble *y, int n, IssmPDouble radius)
 
void InterpolationIDW (IssmDouble *pprediction, IssmDouble x_interp, IssmDouble y_interp, IssmDouble radius, int mindata, int maxdata, IssmDouble power)
 
void InterpolationV4 (IssmDouble *pprediction, IssmDouble x_interp, IssmDouble y_interp, IssmDouble radius, int mindata, int maxdata)
 
void InterpolationKriging (IssmDouble *pprediction, IssmDouble *perror, IssmDouble x_interp, IssmDouble y_interp, IssmDouble radius, int mindata, int maxdata, Variogram *variogram)
 
void InterpolationNearestNeighbor (IssmDouble *pprediction, IssmDouble x_interp, IssmDouble y_interp, IssmDouble radius)
 
void ObservationList (IssmDouble **px, IssmDouble **py, IssmDouble **pobs, int *pnobs)
 
void ObservationList (IssmDouble **px, IssmDouble **py, IssmDouble **pobs, int *pnobs, IssmDouble x_interp, IssmDouble y_interp, IssmDouble radius, int maxdata)
 
void ObservationListCovertree (IssmDouble **px, IssmDouble **py, IssmDouble **pobs, int *pnobs, IssmDouble x_interp, IssmDouble y_interp, IssmDouble radius, int maxdata)
 
void ObservationListQuadtree (IssmDouble **px, IssmDouble **py, IssmDouble **pobs, int *pnobs, IssmDouble x_interp, IssmDouble y_interp, IssmDouble radius, int maxdata)
 
void QuadtreeColoring (IssmDouble *A, IssmDouble *x, IssmDouble *y, int n)
 
void Variomap (IssmDouble *gamma, IssmDouble *x, int n)
 
- Public Member Functions inherited from DataSet
 DataSet ()
 
 DataSet (int enum_type)
 
 ~DataSet ()
 
void Marshall (char **pmarshalled_data, int *pmarshalled_data_size, int marshall_direction)
 
int GetEnum ()
 
int GetEnum (int offset)
 
void Echo ()
 
void DeepEcho ()
 
int AddObject (Object *object)
 
int DeleteObject (int id)
 
int Size ()
 
void clear ()
 
Object * GetObjectByOffset (int offset)
 
Object * GetObjectById (int *poffset, int eid)
 
void Presort ()
 
void Sort ()
 
DataSet * Copy (void)
 
int DeleteObject (Object *object)
 

Private Attributes

int treetype
 
Quadtree * quadtree
 
Covertree * covertree
 

Additional Inherited Members

- Data Fields inherited from DataSet
std::vector< Object * > objects
 
int enum_type
 
int sorted
 
int presorted
 
int numsorted
 
int * sorted_ids
 
int * id_offsets
 

Detailed Description

Declaration of Observations class.

Declaration of Observations class. Observations are vector lists (Containers) of Observation objects.

Definition at line 16 of file Observations.h.

Constructor & Destructor Documentation

◆ Observations() [1/2]

Observations::Observations ( )

Definition at line 33 of file Observations.cpp.

33  {/*{{{*/
34  this->treetype = 0;
35  this->quadtree = NULL;
36  this->covertree = NULL;
37  return;
38 }

◆ Observations() [2/2]

Observations::Observations ( IssmDouble *  observations_list,
IssmDouble *  x,
IssmDouble *  y,
int  n,
Options *  options 
)

Definition at line 40 of file Observations.cpp.

40  {/*{{{*/
41 
42  /*Check that there are observations*/
43  if(n<=0) _error_("No observation found");
44 
45  /*Get tree type (FIXME)*/
46  IssmDouble dtree = 0.;
47  options->Get(&dtree,"treetype",1.);
48  this->treetype = reCast<int>(dtree);
49  switch(this->treetype){
50  case 1:
51  this->covertree = NULL;
52  this->InitQuadtree(observations_list,x,y,n,options);
53  break;
54  case 2:
55  this->quadtree = NULL;
56  this->InitCovertree(observations_list,x,y,n,options);
57  break;
58  default:
59  _error_("Tree type "<<this->treetype<<" not supported yet (1: quadtree, 2: covertree)");
60  }
61 }

◆ ~Observations()

Observations::~Observations ( )

Definition at line 63 of file Observations.cpp.

63  {/*{{{*/
64  switch(this->treetype){
65  case 1:
66  delete this->quadtree;
67  break;
68  case 2:
69  delete this->covertree;
70  break;
71  default:
72  _printf_("Tree type "<<this->treetype<<" not supported yet (1: quadtree, 2: covertree)");
73  }
74  return;
75 }

Member Function Documentation

◆ InitCovertree()

void Observations::InitCovertree ( IssmDouble *  observations_list,
IssmDouble *  x,
IssmDouble *  y,
int  n,
Options *  options 
)

Definition at line 79 of file Observations.cpp.

79  {/*{{{*/
80 
81  /*Intermediaries*/
82  IssmPDouble minspacing,mintrimming,maxtrimming;
83 
84  /*Checks*/
85  _assert_(n);
86 
87  /*Get trimming limits*/
88  options->Get(&mintrimming,"mintrimming",-1.e+21);
89  options->Get(&maxtrimming,"maxtrimming",+1.e+21);
90  options->Get(&minspacing,"minspacing",0.01);
91  if(minspacing<=0) _error_("minspacing must > 0");
92 
93  /*Get maximum distance between 2 points
94  * maxDist should be the maximum distance that any two points
95  * can have between each other. IE p.distance(q) < maxDist for all
96  * p,q that you will ever try to insert. The cover tree may be invalid
97  * if an inaccurate maxDist is given.*/
98  IssmPDouble xmin = x[0];
99  IssmPDouble xmax = x[0];
100  IssmPDouble ymin = y[0];
101  IssmPDouble ymax = y[0];
102  for(int i=1;i<n;i++){
103  if(x[i]<xmin) xmin=x[i];
104  if(x[i]>xmax) xmax=x[i];
105  if(y[i]<ymin) ymin=y[i];
106  if(y[i]>ymax) ymax=y[i];
107  }
108  IssmPDouble maxDist = sqrt(pow(xmax-xmin,2)+pow(ymax-ymin,2));
109  IssmPDouble base = 2.;
110  int maxdepth = ceilf(log(maxDist)/log(base));
111 
112  _printf0_("Generating covertree with a maximum depth " << maxdepth <<"... ");
113  this->covertree=new Covertree(maxdepth);
114 
115  for(int i=0;i<n;i++){
116 
117  /*First check limits*/
118  if(observations_list[i]>maxtrimming) continue;
119  if(observations_list[i]<mintrimming) continue;
120 
121  /*Second, check that this observation is not too close from another one*/
122  Observation newobs = Observation(x[i],y[i],observations_list[i]);
123  if(i>0 && this->covertree->getRoot()){
124  /*Get closest obs and see if it is too close*/
125  std::vector<Observation> kNN=(this->covertree->kNearestNeighbors(newobs,1));
126  Observation oldobs = (*kNN.begin());
127  if(oldobs.distance(newobs)<minspacing) continue;
128  }
129 
130  this->covertree->insert(newobs);
131  }
132  _printf0_("done\n");
133 }

◆ InitQuadtree()

void Observations::InitQuadtree ( IssmDouble *  observations_list,
IssmDouble *  x,
IssmDouble *  y,
int  n,
Options *  options 
)

Definition at line 135 of file Observations.cpp.

135  {/*{{{*/
136 
137  /*Intermediaries*/
138  int i,maxdepth,level,counter,index;
139  int xi,yi;
140  IssmPDouble xmin,xmax,ymin,ymax;
141  IssmPDouble offset,minlength,minspacing,mintrimming,maxtrimming;
142  Observation *observation = NULL;
143 
144  /*Checks*/
145  _assert_(n);
146 
147  /*Get extrema*/
148  xmin=x[0]; ymin=y[0];
149  xmax=x[0]; ymax=y[0];
150  for(i=1;i<n;i++){
151  xmin=min(xmin,x[i]); ymin=min(ymin,y[i]);
152  xmax=max(xmax,x[i]); ymax=max(ymax,y[i]);
153  }
154  offset=0.05*(xmax-xmin); xmin-=offset; xmax+=offset;
155  offset=0.05*(ymax-ymin); ymin-=offset; ymax+=offset;
156 
157  /*Get trimming limits*/
158  options->Get(&mintrimming,"mintrimming",-1.e+21);
159  options->Get(&maxtrimming,"maxtrimming",+1.e+21);
160  options->Get(&minspacing,"minspacing",0.01);
161  if(minspacing<=0) _error_("minspacing must > 0");
162 
163  /*Get Minimum box size*/
164  if(options->GetOption("boxlength")){
165  options->Get(&minlength,"boxlength");
166  if(minlength<=0)_error_("boxlength should be a positive number");
167  maxdepth=reCast<int,IssmPDouble>(log(max(xmax-xmin,ymax-ymin)/minlength +1)/log(2.0));
168  }
169  else{
170  maxdepth = 30;
171  minlength=max(xmax-xmin,ymax-ymin)/IssmPDouble((1L<<maxdepth)-1);
172  }
173 
174  /*Initialize Quadtree*/
175  _printf0_("Generating quadtree with a maximum box size " << minlength << " (depth=" << maxdepth << ")... ");
176  this->quadtree = new Quadtree(xmin,xmax,ymin,ymax,maxdepth);
177 
178  /*Add observations one by one*/
179  counter = 0;
180  for(i=0;i<n;i++){
181 
182  /*First check limits*/
183  if(observations_list[i]>maxtrimming) continue;
184  if(observations_list[i]<mintrimming) continue;
185 
186  /*Second, check that this observation is not too close from another one*/
187  this->quadtree->ClosestObs(&index,x[i],y[i]);
188  if(index>=0){
189  observation=xDynamicCast<Observation*>(this->GetObjectByOffset(index));
190  if(pow(observation->x-x[i],2)+pow(observation->y-y[i],2) < minspacing) continue;
191  }
192 
193  this->quadtree->IntergerCoordinates(&xi,&yi,x[i],y[i]);
194  this->quadtree->QuadtreeDepth2(&level,xi,yi);
195  if((int)level <= maxdepth){
196  observation = new Observation(x[i],y[i],xi,yi,counter++,observations_list[i]);
197  this->quadtree->Add(observation);
198  this->AddObject(observation);
199  }
200  else{
201  /*We need to average with the current observations*/
202  this->quadtree->AddAndAverage(x[i],y[i],observations_list[i]);
203  }
204  }
205  _printf0_("done\n");
206  _printf0_("Initial number of observations: " << n << "\n");
207  _printf0_(" Final number of observations: " << this->quadtree->NbObs << "\n");
208 }

◆ ClosestObservation()

void Observations::ClosestObservation ( IssmDouble *  px,
IssmDouble *  py,
IssmDouble *  pobs,
IssmDouble  x_interp,
IssmDouble  y_interp,
IssmDouble  radius 
)

Definition at line 212 of file Observations.cpp.

212  {/*{{{*/
213 
214  switch(this->treetype){
215  case 1:
216  this->ClosestObservationQuadtree(px,py,pobs,x_interp,y_interp,radius);
217  break;
218  case 2:
219  this->ClosestObservationCovertree(px,py,pobs,x_interp,y_interp,radius);
220  break;
221  default:
222  _error_("Tree type "<<this->treetype<<" not supported yet (1: quadtree, 2: covertree)");
223  }
224 
225 }/*}}}*/

◆ ClosestObservationCovertree()

void Observations::ClosestObservationCovertree ( IssmDouble *  px,
IssmDouble *  py,
IssmDouble *  pobs,
IssmDouble  x_interp,
IssmDouble  y_interp,
IssmDouble  radius 
)

Definition at line 226 of file Observations.cpp.

226  {/*{{{*/
227 
228  IssmPDouble hmin = UNDEF;
229 
230  if(this->covertree->getRoot()){
231  /*Get closest obs and see if it is too close*/
232  Observation newobs = Observation(x_interp,y_interp,0.);
233  std::vector<Observation> kNN=(this->covertree->kNearestNeighbors(newobs,1));
234  Observation observation = (*kNN.begin());
235  hmin = observation.distance(newobs);
236  if(hmin<=radius){
237  *px = observation.x;
238  *py = observation.y;
239  *pobs = observation.value;
240  return;
241  }
242  }
243 
244  *px = UNDEF;
245  *py = UNDEF;
246  *pobs = UNDEF;
247 }/*}}}*/

◆ ClosestObservationQuadtree()

void Observations::ClosestObservationQuadtree ( IssmDouble *  px,
IssmDouble *  py,
IssmDouble *  pobs,
IssmDouble  x_interp,
IssmDouble  y_interp,
IssmDouble  radius 
)

Definition at line 248 of file Observations.cpp.

248  {/*{{{*/
249 
250  /*Output and Intermediaries*/
251  int nobs,i,index;
252  IssmPDouble hmin,h2,hmin2;
253  int *indices = NULL;
254  Observation *observation = NULL;
255 
256  /*If radius is not provided or is 0, return all observations*/
257  if(radius==0) radius=this->quadtree->root->length;
258 
259  /*For CPPcheck*/
260  hmin = 2*radius;
261 
262  /*First, find closest point in Quadtree (fast but might not be the true closest obs)*/
263  this->quadtree->ClosestObs(&index,x_interp,y_interp);
264  if(index>=0){
265  observation=xDynamicCast<Observation*>(this->GetObjectByOffset(index));
266  hmin = sqrt((observation->x-x_interp)*(observation->x-x_interp) + (observation->y-y_interp)*(observation->y-y_interp));
267  if(hmin<radius) radius=hmin;
268  }
269 
270  /*Find all observations that are in radius*/
271  this->quadtree->RangeSearch(&indices,&nobs,x_interp,y_interp,radius);
272  for (i=0;i<nobs;i++){
273  observation=xDynamicCast<Observation*>(this->GetObjectByOffset(indices[i]));
274  h2 = (observation->x-x_interp)*(observation->x-x_interp) + (observation->y-y_interp)*(observation->y-y_interp);
275  if(i==0){
276  hmin2 = h2;
277  index = indices[i];
278  }
279  else{
280  if(h2<hmin2){
281  hmin2 = h2;
282  index = indices[i];
283  }
284  }
285  }
286 
287  /*Assign output pointer*/
288  if(nobs || hmin==radius){
289  observation=xDynamicCast<Observation*>(this->GetObjectByOffset(index));
290  *px = observation->x;
291  *py = observation->y;
292  *pobs = observation->value;
293  }
294  else{
295  *px = UNDEF;
296  *py = UNDEF;
297  *pobs = UNDEF;
298  }
299  xDelete<int>(indices);
300 
301 }/*}}}*/

◆ Distances()

void Observations::Distances ( IssmPDouble *  distances,
IssmPDouble *  x,
IssmPDouble *  y,
int  n,
IssmPDouble  radius 
)

Definition at line 302 of file Observations.cpp.

302  {/*{{{*/
303 
304  IssmPDouble xi,yi,obs;
305 
306  for(int i=0;i<n;i++){
307  this->ClosestObservation(&xi,&yi,&obs,x[i],y[i],radius);
308  if(xi==UNDEF && yi==UNDEF){
309  distances[i]=UNDEF;
310  }
311  else{
312  distances[i]=sqrt( (x[i]-xi)*(x[i]-xi) + (y[i]-yi)*(y[i]-yi) );
313  }
314  }
315 }/*}}}*/

◆ InterpolationIDW()

void Observations::InterpolationIDW ( IssmDouble *  pprediction,
IssmDouble  x_interp,
IssmDouble  y_interp,
IssmDouble  radius,
int  mindata,
int  maxdata,
IssmDouble  power 
)

Definition at line 481 of file Observations.cpp.

481  {/*{{{*/
482 
483  /*Intermediaries*/
484  int i,n_obs;
485  IssmPDouble prediction;
486  IssmPDouble numerator,denominator,h,weight;
487  IssmPDouble *x = NULL;
488  IssmPDouble *y = NULL;
489  IssmPDouble *obs = NULL;
490 
491  /*Some checks*/
492  _assert_(maxdata>0);
493  _assert_(pprediction);
494  _assert_(power>0);
495 
496  /*Get list of observations for current point*/
497  this->ObservationList(&x,&y,&obs,&n_obs,x_interp,y_interp,radius,maxdata);
498 
499  /*If we have less observations than mindata, return UNDEF*/
500  if(n_obs<mindata){
501  prediction = UNDEF;
502  }
503  else{
504  numerator = 0.;
505  denominator = 0.;
506  for(i=0;i<n_obs;i++){
507  h = sqrt( (x[i]-x_interp)*(x[i]-x_interp) + (y[i]-y_interp)*(y[i]-y_interp));
508  if (h<0.0000001){
509  numerator = obs[i];
510  denominator = 1.;
511  break;
512  }
513  weight = 1./pow(h,power);
514  numerator += weight*obs[i];
515  denominator += weight;
516  }
517  prediction = numerator/denominator;
518  }
519 
520  /*clean-up*/
521  *pprediction = prediction;
522  xDelete<IssmPDouble>(x);
523  xDelete<IssmPDouble>(y);
524  xDelete<IssmPDouble>(obs);
525 }/*}}}*/

◆ InterpolationV4()

void Observations::InterpolationV4 ( IssmDouble *  pprediction,
IssmDouble  x_interp,
IssmDouble  y_interp,
IssmDouble  radius,
int  mindata,
int  maxdata 
)

Definition at line 610 of file Observations.cpp.

610  {/*{{{*/
611  /* Reference: David T. Sandwell, Biharmonic spline interpolation of GEOS-3
612  * and SEASAT altimeter data, Geophysical Research Letters, 2, 139-142,
613  * 1987. Describes interpolation using value or gradient of value in any
614  * dimension.*/
615 
616  /*Intermediaries*/
617  int i,j,n_obs;
618  IssmPDouble prediction,h;
619  IssmPDouble *x = NULL;
620  IssmPDouble *y = NULL;
621  IssmPDouble *obs = NULL;
622  IssmPDouble *Green = NULL;
623  IssmPDouble *weights = NULL;
624  IssmPDouble *g = NULL;
625 
626  /*Some checks*/
627  _assert_(maxdata>0);
628  _assert_(pprediction);
629 
630  /*Get list of observations for current point*/
631  this->ObservationList(&x,&y,&obs,&n_obs,x_interp,y_interp,radius,maxdata);
632 
633  /*If we have less observations than mindata, return UNDEF*/
634  if(n_obs<mindata || n_obs<2){
635  prediction = UNDEF;
636  }
637  else{
638 
639  /*Allocate intermediary matrix and vectors*/
640  Green = xNew<IssmPDouble>(n_obs*n_obs);
641  g = xNew<IssmPDouble>(n_obs);
642 
643  /*First: distance vector*/
644  for(i=0;i<n_obs;i++){
645  h = sqrt( (x[i]-x_interp)*(x[i]-x_interp) + (y[i]-y_interp)*(y[i]-y_interp) );
646  if(h>0){
647  g[i] = h*h*(log(h)-1.);
648  }
649  else{
650  g[i] = 0.;
651  }
652  }
653 
654  /*Build Green function matrix*/
655  for(i=0;i<n_obs;i++){
656  for(j=0;j<=i;j++){
657  h = sqrt( (x[i]-x[j])*(x[i]-x[j]) + (y[i]-y[j])*(y[i]-y[j]) );
658  if(h>0){
659  Green[j*n_obs+i] = h*h*(log(h)-1.);
660  }
661  else{
662  Green[j*n_obs+i] = 0.;
663  }
664  Green[i*n_obs+j] = Green[j*n_obs+i];
665  }
666  /*Zero diagonal (should be done already, but just in case)*/
667  Green[i*n_obs+i] = 0.;
668  }
669 
670  /*Compute weights*/
671 #if _HAVE_GSL_
672  DenseGslSolve(&weights,Green,obs,n_obs); // Green^-1 obs
673 #else
674  _error_("GSL is required");
675 #endif
676 
677  /*Interpolate*/
678  prediction = 0;
679  for(i=0;i<n_obs;i++) prediction += weights[i]*g[i];
680 
681  }
682 
683  /*clean-up*/
684  *pprediction = prediction;
685  xDelete<IssmPDouble>(x);
686  xDelete<IssmPDouble>(y);
687  xDelete<IssmPDouble>(obs);
688  xDelete<IssmPDouble>(Green);
689  xDelete<IssmPDouble>(g);
690  xDelete<IssmPDouble>(weights);
691 }/*}}}*/

◆ InterpolationKriging()

void Observations::InterpolationKriging ( IssmDouble *  pprediction,
IssmDouble *  perror,
IssmDouble  x_interp,
IssmDouble  y_interp,
IssmDouble  radius,
int  mindata,
int  maxdata,
Variogram *  variogram 
)

Definition at line 526 of file Observations.cpp.

526  {/*{{{*/
527 
528  /*Intermediaries*/
529  int i,j,n_obs;
530  IssmPDouble prediction,error;
531  IssmPDouble *x = NULL;
532  IssmPDouble *y = NULL;
533  IssmPDouble *obs = NULL;
534  IssmPDouble *Lambda = NULL;
535 
536  /*Some checks*/
537  _assert_(mindata>0 && maxdata>0);
538  _assert_(pprediction && perror);
539 
540  /*Get list of observations for current point*/
541  this->ObservationList(&x,&y,&obs,&n_obs,x_interp,y_interp,radius,maxdata);
542 
543  /*If we have less observations than mindata, return UNDEF*/
544  if(n_obs<mindata){
545  *pprediction = -999.0;
546  *perror = -999.0;
547  return;
548  }
549 
550  /*Allocate intermediary matrix and vectors*/
551  IssmPDouble* A = xNew<IssmPDouble>((n_obs+1)*(n_obs+1));
552  IssmPDouble* B = xNew<IssmPDouble>(n_obs+1);
553 
554  IssmDouble unbias = variogram->Covariance(0.,0.);
555  /*First: Create semivariogram matrix for observations*/
556  for(i=0;i<n_obs;i++){
557  //printf("%g %g ==> %g\n",x[i],y[i],sqrt(pow(x[i]-x_interp,2)+pow(y[i]-y_interp,2)));
558  for(j=0;j<=i;j++){
559  A[i*(n_obs+1)+j] = variogram->Covariance(x[i]-x[j],y[i]-y[j]);
560  A[j*(n_obs+1)+i] = A[i*(n_obs+1)+j];
561  }
562  A[i*(n_obs+1)+n_obs] = unbias;
563  //A[i*(n_obs+1)+n_obs] = 1.;
564  }
565  for(i=0;i<n_obs;i++) A[n_obs*(n_obs+1)+i]=unbias;
566  //for(i=0;i<n_obs;i++) A[n_obs*(n_obs+1)+i]=1.;
567  A[n_obs*(n_obs+1)+n_obs] = 0.;
568 
569  /*Get semivariogram vector associated to this location*/
570  for(i=0;i<n_obs;i++) B[i] = variogram->Covariance(x[i]-x_interp,y[i]-y_interp);
571  B[n_obs] = unbias;
572  //B[n_obs] = 1.;
573 
574  /*Solve the three linear systems*/
575 #if _HAVE_GSL_
576  DenseGslSolve(&Lambda,A,B,n_obs+1); // Gamma^-1 Z
577 #else
578  _error_("GSL is required");
579 #endif
580 
581  /*Compute predictor*/
582  prediction = 0.;
583  for(i=0;i<n_obs;i++) prediction += Lambda[i]*obs[i];
584 
585  /*Compute error (GSLIB p15 eq II.14)*/
586  error = variogram->Covariance(0.,0.)*(1. - Lambda[n_obs]);;
587  for(i=0;i<n_obs;i++) error += -Lambda[i]*B[i];
588 
589  /*clean-up*/
590  *pprediction = prediction;
591  *perror = error;
592  xDelete<IssmPDouble>(x);
593  xDelete<IssmPDouble>(y);
594  xDelete<IssmPDouble>(obs);
595  xDelete<IssmPDouble>(A);
596  xDelete<IssmPDouble>(B);
597  xDelete<IssmPDouble>(Lambda);
598 }/*}}}*/

◆ InterpolationNearestNeighbor()

void Observations::InterpolationNearestNeighbor ( IssmDouble *  pprediction,
IssmDouble  x_interp,
IssmDouble  y_interp,
IssmDouble  radius 
)

Definition at line 599 of file Observations.cpp.

599  {/*{{{*/
600 
601  /*Intermediaries*/
602  IssmPDouble x,y,obs;
603 
604  /*Get clostest observation*/
605  this->ClosestObservation(&x,&y,&obs,x_interp,y_interp,radius);
606 
607  /*Assign output pointer*/
608  *pprediction = obs;
609 }/*}}}*/

◆ ObservationList() [1/2]

void Observations::ObservationList ( IssmDouble **  px,
IssmDouble **  py,
IssmDouble **  pobs,
int *  pnobs 
)

Definition at line 316 of file Observations.cpp.

316  {/*{{{*/
317 
318  /*Output and Intermediaries*/
319  int nobs;
320  IssmPDouble *x = NULL;
321  IssmPDouble *y = NULL;
322  IssmPDouble *obs = NULL;
323  Observation *observation = NULL;
324 
325  nobs = this->Size();
326 
327  if(nobs){
328  x = xNew<IssmPDouble>(nobs);
329  y = xNew<IssmPDouble>(nobs);
330  obs = xNew<IssmPDouble>(nobs);
331  for(int i=0;i<this->Size();i++){
332  observation=xDynamicCast<Observation*>(this->GetObjectByOffset(i));
333  observation->WriteXYObs(&x[i],&y[i],&obs[i]);
334  }
335  }
336 
337  /*Assign output pointer*/
338  *px=x;
339  *py=y;
340  *pobs=obs;
341  *pnobs=nobs;
342 }/*}}}*/

◆ ObservationList() [2/2]

void Observations::ObservationList ( IssmDouble **  px,
IssmDouble **  py,
IssmDouble **  pobs,
int *  pnobs,
IssmDouble  x_interp,
IssmDouble  y_interp,
IssmDouble  radius,
int  maxdata 
)

Definition at line 343 of file Observations.cpp.

343  {/*{{{*/
344 
345  switch(this->treetype){
346  case 1:
347  this->ObservationListQuadtree(px,py,pobs,pnobs,x_interp,y_interp,radius,maxdata);
348  break;
349  case 2:
350  this->ObservationListCovertree(px,py,pobs,pnobs,x_interp,y_interp,radius,maxdata);
351  break;
352  default:
353  _error_("Tree type "<<this->treetype<<" not supported yet (1: quadtree, 2: covertree)");
354  }
355 }/*}}}*/

◆ ObservationListCovertree()

void Observations::ObservationListCovertree ( IssmDouble **  px,
IssmDouble **  py,
IssmDouble **  pobs,
int *  pnobs,
IssmDouble  x_interp,
IssmDouble  y_interp,
IssmDouble  radius,
int  maxdata 
)

Definition at line 356 of file Observations.cpp.

356  {/*{{{*/
357 
358  double *x = NULL;
359  double *y = NULL;
360  double *obs = NULL;
361  Observation observation=Observation(x_interp,y_interp,0.);
362  std::vector<Observation> kNN;
363 
364  kNN=(this->covertree->kNearestNeighbors(observation, maxdata));
365  //cout << "kNN's size: " << kNN.size() << " (maxdata = " <<maxdata<<")"<<endl;
366 
367  //kNN is sort from closest to farthest neighbor
368  //searches for the first neighbor that is out of radius
369  //deletes and resizes the kNN vector
370  vector<Observation>::iterator it;
371  if(radius>0.){
372  for (it = kNN.begin(); it != kNN.end(); ++it) {
373  //(*it).print();
374  //cout << "\n" << (*it).distance(observation) << endl;
375  if ((*it).distance(observation) > radius) {
376  break;
377  }
378  }
379  kNN.erase(it, kNN.end());
380  }
381 
382  /*Allocate vectors*/
383  x = new double[kNN.size()];
384  y = new double[kNN.size()];
385  obs = new double[kNN.size()];
386 
387  /*Loop over all observations and fill in x, y and obs*/
388  int i = 0;
389  for(it = kNN.begin(); it != kNN.end(); ++it) {
390  (*it).WriteXYObs((*it), &x[i], &y[i], &obs[i]);
391  i++;
392  }
393 
394  *px=x;
395  *py=y;
396  *pobs=obs;
397  *pnobs = kNN.size();
398 }/*}}}*/

◆ ObservationListQuadtree()

void Observations::ObservationListQuadtree ( IssmDouble **  px,
IssmDouble **  py,
IssmDouble **  pobs,
int *  pnobs,
IssmDouble  x_interp,
IssmDouble  y_interp,
IssmDouble  radius,
int  maxdata 
)

Definition at line 399 of file Observations.cpp.

399  {/*{{{*/
400 
401  /*Output and Intermediaries*/
402  bool stop;
403  int nobs,tempnobs,i,j,k,n,counter;
404  IssmPDouble h2,radius2;
405  int *indices = NULL;
406  int *tempindices = NULL;
407  IssmPDouble *dists = NULL;
408  IssmPDouble *x = NULL;
409  IssmPDouble *y = NULL;
410  IssmPDouble *obs = NULL;
411  Observation *observation = NULL;
412 
413  /*If radius is not provided or is 0, return all observations*/
414  if(radius==0.) radius=this->quadtree->root->length*2.;
415 
416  /*Compute radius square*/
417  radius2 = radius*radius;
418 
419  /*Find all observations that are in radius*/
420  this->quadtree->RangeSearch(&tempindices,&tempnobs,x_interp,y_interp,radius);
421  if(tempnobs){
422  indices = xNew<int>(tempnobs);
423  dists = xNew<IssmPDouble>(tempnobs);
424  }
425  nobs = 0;
426  for(i=0;i<tempnobs;i++){
427  observation=xDynamicCast<Observation*>(this->GetObjectByOffset(tempindices[i]));
428  h2 = (observation->x-x_interp)*(observation->x-x_interp) + (observation->y-y_interp)*(observation->y-y_interp);
429 
430  if(nobs==maxdata && h2>radius2) continue;
431  if(nobs<maxdata){
432  indices[nobs] = tempindices[i];
433  dists[nobs] = h2;
434  nobs++;
435  }
436  if(nobs==1) continue;
437 
438  /*Sort all dists up to now*/
439  n=nobs-1;
440  stop = false;
441  for(k=0;k<n-1;k++){
442  if(h2<dists[k]){
443  counter=1;
444  for(int jj=k;jj<n;jj++){
445  j = n-counter;
446  dists[j+1] = dists[j];
447  indices[j+1] = indices[j];
448  counter++;
449  }
450  dists[k] = h2;
451  indices[k] = tempindices[i];
452  stop = true;
453  break;
454  }
455  if(stop) break;
456  }
457  }
458  xDelete<IssmPDouble>(dists);
459  xDelete<int>(tempindices);
460 
461  if(nobs){
462  /*Allocate vectors*/
463  x = xNew<IssmPDouble>(nobs);
464  y = xNew<IssmPDouble>(nobs);
465  obs = xNew<IssmPDouble>(nobs);
466 
467  /*Loop over all observations and fill in x, y and obs*/
468  for(i=0;i<nobs;i++){
469  observation=xDynamicCast<Observation*>(this->GetObjectByOffset(indices[i]));
470  observation->WriteXYObs(&x[i],&y[i],&obs[i]);
471  }
472  }
473 
474  /*Assign output pointer*/
475  xDelete<int>(indices);
476  *px=x;
477  *py=y;
478  *pobs=obs;
479  *pnobs=nobs;
480 }/*}}}*/

◆ QuadtreeColoring()

void Observations::QuadtreeColoring ( IssmDouble *  A,
IssmDouble *  x,
IssmDouble *  y,
int  n 
)

Definition at line 692 of file Observations.cpp.

692  {/*{{{*/
693 
694  if(this->treetype!=1) _error_("Tree type is not quadtree");
695  int xi,yi,level;
696 
697  for(int i=0;i<n;i++){
698  this->quadtree->IntergerCoordinates(&xi,&yi,x[i],y[i]);
699  this->quadtree->QuadtreeDepth(&level,xi,yi);
700  A[i]=(IssmPDouble)level;
701  }
702 
703 }/*}}}*/

◆ Variomap()

void Observations::Variomap ( IssmDouble *  gamma,
IssmDouble *  x,
int  n 
)

Definition at line 704 of file Observations.cpp.

704  {/*{{{*/
705 
706  /*Output and Intermediaries*/
707  int i,j,k;
708  IssmPDouble distance;
709  Observation *observation1 = NULL;
710  Observation *observation2 = NULL;
711 
712  IssmPDouble *counter = xNew<IssmPDouble>(n);
713  for(j=0;j<n;j++) counter[j] = 0.0;
714  for(j=0;j<n;j++) gamma[j] = 0.0;
715 
716  for(i=0;i<this->Size();i++){
717  observation1=xDynamicCast<Observation*>(this->GetObjectByOffset(i));
718 
719  for(j=i+1;j<this->Size();j++){
720  observation2=xDynamicCast<Observation*>(this->GetObjectByOffset(j));
721 
722  distance=sqrt(pow(observation1->x - observation2->x,2) + pow(observation1->y - observation2->y,2));
723  if(distance>x[n-1]) continue;
724 
725  int index = int(distance/(x[1]-x[0]));
726  if(index>n-1) index = n-1;
727  if(index<0) index = 0;
728 
729  gamma[index] += 1./2.*pow(observation1->value - observation2->value,2);
730  counter[index] += 1.;
731  }
732  }
733 
734  /*Normalize semivariogram*/
735  gamma[0]=0.;
736  for(k=0;k<n;k++){
737  if(counter[k]) gamma[k] = gamma[k]/counter[k];
738  }
739 
740  /*Assign output pointer*/
741  xDelete<IssmPDouble>(counter);
742 }/*}}}*/

Field Documentation

◆ treetype

int Observations::treetype
private

Definition at line 19 of file Observations.h.

◆ quadtree

Quadtree* Observations::quadtree
private

Definition at line 20 of file Observations.h.

◆ covertree

Covertree* Observations::covertree
private

Definition at line 21 of file Observations.h.

The documentation for this class was generated from the following files:
DataSet::Size
int Size()
Definition: DataSet.cpp:399
Observation::distance
double distance(const Observation &ob) const
Definition: Observation.cpp:78
Covertree::insert
void insert(const Observation &newObservation)
Definition: Covertree.cpp:49
_assert_
#define _assert_(ignore)
Definition: exceptions.h:37
IssmDouble
double IssmDouble
Definition: types.h:37
Observation
Definition: Observation.h:10
_printf0_
#define _printf0_(StreamArgs)
Definition: Print.h:29
Observation::x
double x
Definition: Observation.h:13
DataSet::AddObject
int AddObject(Object *object)
Definition: DataSet.cpp:252
_printf_
#define _printf_(StreamArgs)
Definition: Print.h:22
Observations::covertree
Covertree * covertree
Definition: Observations.h:21
Quadtree::QuadtreeDepth
void QuadtreeDepth(int *A, int xi, int yi)
Definition: Quadtree.cpp:410
Quadtree::root
QuadtreeBox * root
Definition: Quadtree.h:48
Covertree
Definition: Covertree.h:8
Observations::treetype
int treetype
Definition: Observations.h:19
Variogram::Covariance
virtual double Covariance(double deltax, double deltay)=0
Covertree::getRoot
CoverTreeNode * getRoot() const
Definition: Covertree.cpp:367
Observations::ClosestObservation
void ClosestObservation(IssmDouble *px, IssmDouble *py, IssmDouble *pobs, IssmDouble x_interp, IssmDouble y_interp, IssmDouble radius)
Definition: Observations.cpp:212
Observations::ClosestObservationQuadtree
void ClosestObservationQuadtree(IssmDouble *px, IssmDouble *py, IssmDouble *pobs, IssmDouble x_interp, IssmDouble y_interp, IssmDouble radius)
Definition: Observations.cpp:248
Observation::WriteXYObs
void WriteXYObs(const Observation &ob, double *px, double *py, double *pobs)
Definition: Observation.cpp:90
Observations::ObservationListCovertree
void ObservationListCovertree(IssmDouble **px, IssmDouble **py, IssmDouble **pobs, int *pnobs, IssmDouble x_interp, IssmDouble y_interp, IssmDouble radius, int maxdata)
Definition: Observations.cpp:356
Observations::InitCovertree
void InitCovertree(IssmDouble *observations_list, IssmDouble *x, IssmDouble *y, int n, Options *options)
Definition: Observations.cpp:79
UNDEF
#define UNDEF
Definition: constants.h:8
Observations::quadtree
Quadtree * quadtree
Definition: Observations.h:20
Quadtree::AddAndAverage
void AddAndAverage(double x, double y, double value)
Definition: Quadtree.cpp:197
Observations::ClosestObservationCovertree
void ClosestObservationCovertree(IssmDouble *px, IssmDouble *py, IssmDouble *pobs, IssmDouble x_interp, IssmDouble y_interp, IssmDouble radius)
Definition: Observations.cpp:226
Observations::InitQuadtree
void InitQuadtree(IssmDouble *observations_list, IssmDouble *x, IssmDouble *y, int n, Options *options)
Definition: Observations.cpp:135
Observation::y
double y
Definition: Observation.h:13
Quadtree::QuadtreeBox::length
double length
Definition: Quadtree.h:21
Options::Get
void Get(OptionType *pvalue, const char *name)
Definition: Options.h:21
Observations::ObservationList
void ObservationList(IssmDouble **px, IssmDouble **py, IssmDouble **pobs, int *pnobs)
Definition: Observations.cpp:316
Options::GetOption
Option * GetOption(const char *name)
Definition: Options.cpp:67
Quadtree::ClosestObs
void ClosestObs(int *pindex, double x, double y)
Definition: Quadtree.cpp:245
Quadtree::RangeSearch
void RangeSearch(int **pindices, int *pnobs, double x, double y, double range)
Definition: Quadtree.cpp:393
_error_
#define _error_(StreamArgs)
Definition: exceptions.h:49
Quadtree::Add
void Add(Observation *observation)
Definition: Quadtree.cpp:115
DataSet::GetObjectByOffset
Object * GetObjectByOffset(int offset)
Definition: DataSet.cpp:334
min
IssmDouble min(IssmDouble a, IssmDouble b)
Definition: extrema.cpp:14
Quadtree::IntergerCoordinates
void IntergerCoordinates(int *xi, int *yi, double x, double y)
Definition: Quadtree.cpp:303
Observation::value
double value
Definition: Observation.h:17
Quadtree::NbObs
int NbObs
Definition: Quadtree.h:50
Quadtree::QuadtreeDepth2
void QuadtreeDepth2(int *A, int xi, int yi)
Definition: Quadtree.cpp:437
Covertree::kNearestNeighbors
std::vector< Observation > kNearestNeighbors(const Observation &p, const unsigned int &k) const
Definition: Covertree.cpp:155
DenseGslSolve
void DenseGslSolve(IssmPDouble **pX, IssmPDouble *A, IssmPDouble *B, int n)
Definition: DenseGslSolve.cpp:25
Quadtree
Definition: Quadtree.h:7
max
IssmDouble max(IssmDouble a, IssmDouble b)
Definition: extrema.cpp:24
IssmPDouble
IssmDouble IssmPDouble
Definition: types.h:38
Observations::ObservationListQuadtree
void ObservationListQuadtree(IssmDouble **px, IssmDouble **py, IssmDouble **pobs, int *pnobs, IssmDouble x_interp, IssmDouble y_interp, IssmDouble radius, int maxdata)
Definition: Observations.cpp:399
Generated on Thu Jul 2 2020 08:09:29 for Ice Sheet System Model by   doxygen 1.8.19