source: issm/trunk/src/c/classes/kriging/Observations.cpp@ 20500

Last change on this file since 20500 was 20500, checked in by Mathieu Morlighem, 9 years ago

merged trunk-jpl and trunk for revision 20497
File size: 21.0 KB
Line 
1/*
2 * \file Observations.cpp
3 * \brief: Implementation of Observations class, derived from DataSet class.
4 */
5
6/*Headers: {{{*/
7#ifdef HAVE_CONFIG_H
8 #include <config.h>
9#else
10#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
11#endif
12
13#include <vector>
14#include <functional>
15#include <algorithm>
16#include <iostream>
17
18#include "../Options/Options.h"
19#include "./Observations.h"
20#include "./Observation.h"
21#include "../../datastructures/datastructures.h"
22#include "../../shared/shared.h"
23
24#include "./Quadtree.h"
25#include "./Covertree.h"
26#include "./Variogram.h"
27#include "../../toolkits/toolkits.h"
28
29using namespace std;
30/*}}}*/
31
32/*Object constructors and destructor*/
33Observations::Observations(){/*{{{*/
34 this->treetype = 0;
35 this->quadtree = NULL;
36 this->covertree = NULL;
37 return;
38}
39/*}}}*/
40Observations::Observations(IssmPDouble* observations_list,IssmPDouble* x,IssmPDouble* y,int n,Options* options){/*{{{*/
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}
62/*}}}*/
63Observations::~Observations(){/*{{{*/
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 _error_("Tree type "<<this->treetype<<" not supported yet (1: quadtree, 2: covertree)");
73 }
74 return;
75}
76/*}}}*/
77
78/*Initialize data structures*/
79void Observations::InitQuadtree(IssmPDouble* observations_list,IssmPDouble* x,IssmPDouble* y,int n,Options* options){/*{{{*/
80
81 /*Intermediaries*/
82 int i,maxdepth,level,counter,index;
83 int xi,yi;
84 IssmPDouble xmin,xmax,ymin,ymax;
85 IssmPDouble offset,minlength,minspacing,mintrimming,maxtrimming;
86 Observation *observation = NULL;
87
88 /*Checks*/
89 _assert_(n);
90
91 /*Get extrema*/
92 xmin=x[0]; ymin=y[0];
93 xmax=x[0]; ymax=y[0];
94 for(i=1;i<n;i++){
95 xmin=min(xmin,x[i]); ymin=min(ymin,y[i]);
96 xmax=max(xmax,x[i]); ymax=max(ymax,y[i]);
97 }
98 offset=0.05*(xmax-xmin); xmin-=offset; xmax+=offset;
99 offset=0.05*(ymax-ymin); ymin-=offset; ymax+=offset;
100
101 /*Get trimming limits*/
102 options->Get(&mintrimming,"mintrimming",-1.e+21);
103 options->Get(&maxtrimming,"maxtrimming",+1.e+21);
104 options->Get(&minspacing,"minspacing",0.01);
105 if(minspacing<=0) _error_("minspacing must > 0");
106
107 /*Get Minimum box size*/
108 if(options->GetOption("boxlength")){
109 options->Get(&minlength,"boxlength");
110 if(minlength<=0)_error_("boxlength should be a positive number");
111 maxdepth=reCast<int,IssmPDouble>(log(max(xmax-xmin,ymax-ymin)/minlength +1)/log(2.0));
112 }
113 else{
114 maxdepth = 30;
115 minlength=max(xmax-xmin,ymax-ymin)/IssmPDouble((1L<<maxdepth)-1);
116 }
117
118 /*Initialize Quadtree*/
119 _printf0_("Generating quadtree with a maximum box size " << minlength << " (depth=" << maxdepth << ")... ");
120 this->quadtree = new Quadtree(xmin,xmax,ymin,ymax,maxdepth);
121
122 /*Add observations one by one*/
123 counter = 0;
124 for(i=0;i<n;i++){
125
126 /*First check limits*/
127 if(observations_list[i]>maxtrimming) continue;
128 if(observations_list[i]<mintrimming) continue;
129
130 /*Second, check that this observation is not too close from another one*/
131 this->quadtree->ClosestObs(&index,x[i],y[i]);
132 if(index>=0){
133 observation=xDynamicCast<Observation*>(this->GetObjectByOffset(index));
134 if(pow(observation->x-x[i],2)+pow(observation->y-y[i],2) < minspacing) continue;
135 }
136
137 this->quadtree->IntergerCoordinates(&xi,&yi,x[i],y[i]);
138 this->quadtree->QuadtreeDepth2(&level,xi,yi);
139 if((int)level <= maxdepth){
140 observation = new Observation(x[i],y[i],xi,yi,counter++,observations_list[i]);
141 this->quadtree->Add(observation);
142 this->AddObject(observation);
143 }
144 else{
145 /*We need to average with the current observations*/
146 this->quadtree->AddAndAverage(x[i],y[i],observations_list[i]);
147 }
148 }
149 _printf0_("done\n");
150 _printf0_("Initial number of observations: " << n << "\n");
151 _printf0_(" Final number of observations: " << this->quadtree->NbObs << "\n");
152}
153/*}}}*/
154void Observations::InitCovertree(IssmPDouble* observations_list,IssmPDouble* x,IssmPDouble* y,int n,Options* options){/*{{{*/
155
156 /*Intermediaries*/
157 IssmPDouble minspacing,mintrimming,maxtrimming;
158
159 /*Checks*/
160 _assert_(n);
161
162 /*Get trimming limits*/
163 options->Get(&mintrimming,"mintrimming",-1.e+21);
164 options->Get(&maxtrimming,"maxtrimming",+1.e+21);
165 options->Get(&minspacing,"minspacing",0.01);
166 if(minspacing<=0) _error_("minspacing must > 0");
167
168 /*Get maximum distance between 2 points
169 * maxDist should be the maximum distance that any two points
170 * can have between each other. IE p.distance(q) < maxDist for all
171 * p,q that you will ever try to insert. The cover tree may be invalid
172 * if an inaccurate maxDist is given.*/
173 IssmPDouble xmin = x[0];
174 IssmPDouble xmax = x[0];
175 IssmPDouble ymin = y[0];
176 IssmPDouble ymax = y[0];
177 for(int i=1;i<n;i++){
178 if(x[i]<xmin) xmin=x[i];
179 if(x[i]>xmax) xmax=x[i];
180 if(y[i]<ymin) ymin=y[i];
181 if(y[i]>ymax) ymax=y[i];
182 }
183 IssmPDouble maxDist = sqrt(pow(xmax-xmin,2)+pow(ymax-ymin,2));
184 IssmPDouble base = 2.;
185 int maxdepth = ceilf(log(maxDist)/log(base));
186
187 _printf0_("Generating covertree with a maximum depth " << maxdepth <<"... ");
188 this->covertree=new Covertree(maxdepth);
189
190 for(int i=0;i<n;i++){
191
192 /*First check limits*/
193 if(observations_list[i]>maxtrimming) continue;
194 if(observations_list[i]<mintrimming) continue;
195
196 /*Second, check that this observation is not too close from another one*/
197 Observation newobs = Observation(x[i],y[i],observations_list[i]);
198 if(i>0 && this->covertree->getRoot()){
199 /*Get closest obs and see if it is too close*/
200 std::vector<Observation> kNN=(this->covertree->kNearestNeighbors(newobs,1));
201 Observation oldobs = (*kNN.begin());
202 if(oldobs.distance(newobs)<minspacing) continue;
203 }
204
205 this->covertree->insert(newobs);
206 }
207 _printf0_("done\n");
208}
209/*}}}*/
210
211/*Methods*/
212void Observations::ClosestObservation(IssmPDouble *px,IssmPDouble *py,IssmPDouble *pobs,IssmPDouble x_interp,IssmPDouble y_interp,IssmPDouble radius){/*{{{*/
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}/*}}}*/
226void Observations::ClosestObservationQuadtree(IssmPDouble *px,IssmPDouble *py,IssmPDouble *pobs,IssmPDouble x_interp,IssmPDouble y_interp,IssmPDouble radius){/*{{{*/
227
228 /*Output and Intermediaries*/
229 int nobs,i,index;
230 IssmPDouble hmin,h2,hmin2;
231 int *indices = NULL;
232 Observation *observation = NULL;
233
234 /*If radius is not provided or is 0, return all observations*/
235 if(radius==0) radius=this->quadtree->root->length;
236
237 /*First, find closest point in Quadtree (fast but might not be the true closest obs)*/
238 this->quadtree->ClosestObs(&index,x_interp,y_interp);
239 if(index>=0){
240 observation=xDynamicCast<Observation*>(this->GetObjectByOffset(index));
241 hmin = sqrt((observation->x-x_interp)*(observation->x-x_interp) + (observation->y-y_interp)*(observation->y-y_interp));
242 if(hmin<radius) radius=hmin;
243 }
244
245 /*Find all observations that are in radius*/
246 this->quadtree->RangeSearch(&indices,&nobs,x_interp,y_interp,radius);
247 for (i=0;i<nobs;i++){
248 observation=xDynamicCast<Observation*>(this->GetObjectByOffset(indices[i]));
249 h2 = (observation->x-x_interp)*(observation->x-x_interp) + (observation->y-y_interp)*(observation->y-y_interp);
250 if(i==0){
251 hmin2 = h2;
252 index = indices[i];
253 }
254 else{
255 if(h2<hmin2){
256 hmin2 = h2;
257 index = indices[i];
258 }
259 }
260 }
261
262 /*Assign output pointer*/
263 if(nobs || hmin==radius){
264 observation=xDynamicCast<Observation*>(this->GetObjectByOffset(index));
265 *px = observation->x;
266 *py = observation->y;
267 *pobs = observation->value;
268 }
269 else{
270 *px = UNDEF;
271 *py = UNDEF;
272 *pobs = UNDEF;
273 }
274 xDelete<int>(indices);
275
276}/*}}}*/
277void Observations::ClosestObservationCovertree(IssmPDouble *px,IssmPDouble *py,IssmPDouble *pobs,IssmPDouble x_interp,IssmPDouble y_interp,IssmPDouble radius){/*{{{*/
278
279 IssmPDouble hmin = UNDEF;
280
281 if(this->covertree->getRoot()){
282 /*Get closest obs and see if it is too close*/
283 Observation newobs = Observation(x_interp,y_interp,0.);
284 std::vector<Observation> kNN=(this->covertree->kNearestNeighbors(newobs,1));
285 Observation observation = (*kNN.begin());
286 hmin = observation.distance(newobs);
287 if(hmin<=radius){
288 *px = observation.x;
289 *py = observation.y;
290 *pobs = observation.value;
291 return;
292 }
293 }
294
295 *px = UNDEF;
296 *py = UNDEF;
297 *pobs = UNDEF;
298}/*}}}*/
299void Observations::Distances(IssmPDouble* distances,IssmPDouble *x,IssmPDouble *y,int n,IssmPDouble radius){/*{{{*/
300
301 IssmPDouble xi,yi,obs;
302
303 for(int i=0;i<n;i++){
304 this->ClosestObservation(&xi,&yi,&obs,x[i],y[i],radius);
305 if(xi==UNDEF && yi==UNDEF){
306 distances[i]=UNDEF;
307 }
308 else{
309 distances[i]=sqrt( (x[i]-xi)*(x[i]-xi) + (y[i]-yi)*(y[i]-yi) );
310 }
311 }
312}/*}}}*/
313void Observations::ObservationList(IssmPDouble **px,IssmPDouble **py,IssmPDouble **pobs,int* pnobs){/*{{{*/
314
315 /*Output and Intermediaries*/
316 int nobs;
317 IssmPDouble *x = NULL;
318 IssmPDouble *y = NULL;
319 IssmPDouble *obs = NULL;
320 Observation *observation = NULL;
321
322 nobs = this->Size();
323
324 if(nobs){
325 x = xNew<IssmPDouble>(nobs);
326 y = xNew<IssmPDouble>(nobs);
327 obs = xNew<IssmPDouble>(nobs);
328 for(int i=0;i<this->Size();i++){
329 observation=xDynamicCast<Observation*>(this->GetObjectByOffset(i));
330 observation->WriteXYObs(&x[i],&y[i],&obs[i]);
331 }
332 }
333
334 /*Assign output pointer*/
335 *px=x;
336 *py=y;
337 *pobs=obs;
338 *pnobs=nobs;
339}/*}}}*/
340void Observations::ObservationList(IssmPDouble **px,IssmPDouble **py,IssmPDouble **pobs,int* pnobs,IssmPDouble x_interp,IssmPDouble y_interp,IssmPDouble radius,int maxdata){/*{{{*/
341
342 switch(this->treetype){
343 case 1:
344 this->ObservationListQuadtree(px,py,pobs,pnobs,x_interp,y_interp,radius,maxdata);
345 break;
346 case 2:
347 this->ObservationListCovertree(px,py,pobs,pnobs,x_interp,y_interp,radius,maxdata);
348 break;
349 default:
350 _error_("Tree type "<<this->treetype<<" not supported yet (1: quadtree, 2: covertree)");
351 }
352}/*}}}*/
353void Observations::ObservationListQuadtree(IssmPDouble **px,IssmPDouble **py,IssmPDouble **pobs,int* pnobs,IssmPDouble x_interp,IssmPDouble y_interp,IssmPDouble radius,int maxdata){/*{{{*/
354
355 /*Output and Intermediaries*/
356 bool stop;
357 int nobs,tempnobs,i,j,k,n,counter;
358 IssmPDouble h2,radius2;
359 int *indices = NULL;
360 int *tempindices = NULL;
361 IssmPDouble *dists = NULL;
362 IssmPDouble *x = NULL;
363 IssmPDouble *y = NULL;
364 IssmPDouble *obs = NULL;
365 Observation *observation = NULL;
366
367 /*If radius is not provided or is 0, return all observations*/
368 if(radius==0.) radius=this->quadtree->root->length*2.;
369
370 /*Compute radius square*/
371 radius2 = radius*radius;
372
373 /*Find all observations that are in radius*/
374 this->quadtree->RangeSearch(&tempindices,&tempnobs,x_interp,y_interp,radius);
375 if(tempnobs){
376 indices = xNew<int>(tempnobs);
377 dists = xNew<IssmPDouble>(tempnobs);
378 }
379 nobs = 0;
380 for(i=0;i<tempnobs;i++){
381 observation=xDynamicCast<Observation*>(this->GetObjectByOffset(tempindices[i]));
382 h2 = (observation->x-x_interp)*(observation->x-x_interp) + (observation->y-y_interp)*(observation->y-y_interp);
383
384 if(nobs==maxdata && h2>radius2) continue;
385 if(nobs<maxdata){
386 indices[nobs] = tempindices[i];
387 dists[nobs] = h2;
388 nobs++;
389 }
390 if(nobs==1) continue;
391
392 /*Sort all dists up to now*/
393 n=nobs-1;
394 stop = false;
395 for(k=0;k<n-1;k++){
396 if(h2<dists[k]){
397 counter=1;
398 for(int jj=k;jj<n;jj++){
399 j = n-counter;
400 dists[j+1] = dists[j];
401 indices[j+1] = indices[j];
402 counter++;
403 }
404 dists[k] = h2;
405 indices[k] = tempindices[i];
406 stop = true;
407 break;
408 }
409 if(stop) break;
410 }
411 }
412 xDelete<IssmPDouble>(dists);
413 xDelete<int>(tempindices);
414
415 if(nobs){
416 /*Allocate vectors*/
417 x = xNew<IssmPDouble>(nobs);
418 y = xNew<IssmPDouble>(nobs);
419 obs = xNew<IssmPDouble>(nobs);
420
421 /*Loop over all observations and fill in x, y and obs*/
422 for(i=0;i<nobs;i++){
423 observation=xDynamicCast<Observation*>(this->GetObjectByOffset(indices[i]));
424 observation->WriteXYObs(&x[i],&y[i],&obs[i]);
425 }
426 }
427
428 /*Assign output pointer*/
429 xDelete<int>(indices);
430 *px=x;
431 *py=y;
432 *pobs=obs;
433 *pnobs=nobs;
434}/*}}}*/
435void Observations::ObservationListCovertree(double **px,double **py,double **pobs,int* pnobs,double x_interp,double y_interp,double radius,int maxdata){/*{{{*/
436
437 double *x = NULL;
438 double *y = NULL;
439 double *obs = NULL;
440 Observation observation=Observation(x_interp,y_interp,0.);
441 std::vector<Observation> kNN;
442
443 kNN=(this->covertree->kNearestNeighbors(observation, maxdata));
444 //cout << "kNN's size: " << kNN.size() << " (maxdata = " <<maxdata<<")"<<endl;
445
446 //kNN is sort from closest to farthest neighbor
447 //searches for the first neighbor that is out of radius
448 //deletes and resizes the kNN vector
449 vector<Observation>::iterator it;
450 if(radius>0.){
451 for (it = kNN.begin(); it != kNN.end(); ++it) {
452 //(*it).print();
453 //cout << "\n" << (*it).distance(observation) << endl;
454 if ((*it).distance(observation) > radius) {
455 break;
456 }
457 }
458 kNN.erase(it, kNN.end());
459 }
460
461 /*Allocate vectors*/
462 x = new double[kNN.size()];
463 y = new double[kNN.size()];
464 obs = new double[kNN.size()];
465
466 /*Loop over all observations and fill in x, y and obs*/
467 int i = 0;
468 for(it = kNN.begin(); it != kNN.end(); ++it) {
469 (*it).WriteXYObs((*it), &x[i], &y[i], &obs[i]);
470 i++;
471 }
472
473 *px=x;
474 *py=y;
475 *pobs=obs;
476 *pnobs = kNN.size();
477}/*}}}*/
478void Observations::InterpolationIDW(IssmPDouble *pprediction,IssmPDouble x_interp,IssmPDouble y_interp,IssmPDouble radius,int mindata,int maxdata,IssmPDouble power){/*{{{*/
479
480 /*Intermediaries*/
481 int i,n_obs;
482 IssmPDouble prediction;
483 IssmPDouble numerator,denominator,h,weight;
484 IssmPDouble *x = NULL;
485 IssmPDouble *y = NULL;
486 IssmPDouble *obs = NULL;
487
488 /*Some checks*/
489 _assert_(maxdata>0);
490 _assert_(pprediction);
491 _assert_(power>0);
492
493 /*Get list of observations for current point*/
494 this->ObservationList(&x,&y,&obs,&n_obs,x_interp,y_interp,radius,maxdata);
495
496 /*If we have less observations than mindata, return UNDEF*/
497 if(n_obs<mindata){
498 prediction = UNDEF;
499 }
500 else{
501 numerator = 0.;
502 denominator = 0.;
503 for(i=0;i<n_obs;i++){
504 h = sqrt( (x[i]-x_interp)*(x[i]-x_interp) + (y[i]-y_interp)*(y[i]-y_interp));
505 if (h<0.0000001){
506 numerator = obs[i];
507 denominator = 1.;
508 break;
509 }
510 weight = 1./pow(h,power);
511 numerator += weight*obs[i];
512 denominator += weight;
513 }
514 prediction = numerator/denominator;
515 }
516
517 /*clean-up*/
518 *pprediction = prediction;
519 xDelete<IssmPDouble>(x);
520 xDelete<IssmPDouble>(y);
521 xDelete<IssmPDouble>(obs);
522}/*}}}*/
523void Observations::InterpolationKriging(IssmPDouble *pprediction,IssmPDouble *perror,IssmPDouble x_interp,IssmPDouble y_interp,IssmPDouble radius,int mindata,int maxdata,Variogram* variogram){/*{{{*/
524
525 /*Intermediaries*/
526 int i,j,n_obs;
527 IssmPDouble prediction,error;
528 IssmPDouble *x = NULL;
529 IssmPDouble *y = NULL;
530 IssmPDouble *obs = NULL;
531 IssmPDouble *Lambda = NULL;
532
533 /*Some checks*/
534 _assert_(mindata>0 && maxdata>0);
535 _assert_(pprediction && perror);
536
537 /*Get list of observations for current point*/
538 this->ObservationList(&x,&y,&obs,&n_obs,x_interp,y_interp,radius,maxdata);
539
540 /*If we have less observations than mindata, return UNDEF*/
541 if(n_obs<mindata){
542 *pprediction = -999.0;
543 *perror = -999.0;
544 return;
545 }
546
547 /*Allocate intermediary matrix and vectors*/
548 IssmPDouble* A = xNew<IssmPDouble>((n_obs+1)*(n_obs+1));
549 IssmPDouble* B = xNew<IssmPDouble>(n_obs+1);
550
551 IssmDouble unbias = variogram->Covariance(0.,0.);
552 /*First: Create semivariogram matrix for observations*/
553 for(i=0;i<n_obs;i++){
554 //printf("%g %g ==> %g\n",x[i],y[i],sqrt(pow(x[i]-x_interp,2)+pow(y[i]-y_interp,2)));
555 for(j=0;j<=i;j++){
556 A[i*(n_obs+1)+j] = variogram->Covariance(x[i]-x[j],y[i]-y[j]);
557 A[j*(n_obs+1)+i] = A[i*(n_obs+1)+j];
558 }
559 A[i*(n_obs+1)+n_obs] = unbias;
560 //A[i*(n_obs+1)+n_obs] = 1.;
561 }
562 for(i=0;i<n_obs;i++) A[n_obs*(n_obs+1)+i]=unbias;
563 //for(i=0;i<n_obs;i++) A[n_obs*(n_obs+1)+i]=1.;
564 A[n_obs*(n_obs+1)+n_obs] = 0.;
565
566 /*Get semivariogram vector associated to this location*/
567 for(i=0;i<n_obs;i++) B[i] = variogram->Covariance(x[i]-x_interp,y[i]-y_interp);
568 B[n_obs] = unbias;
569 //B[n_obs] = 1.;
570
571 /*Solve the three linear systems*/
572#if _HAVE_GSL_
573 DenseGslSolve(&Lambda,A,B,n_obs+1); // Gamma^-1 Z
574#else
575 _error_("GSL is required");
576#endif
577
578 /*Compute predictor*/
579 prediction = 0.;
580 for(i=0;i<n_obs;i++) prediction += Lambda[i]*obs[i];
581
582 /*Compute error (GSLIB p15 eq II.14)*/
583 error = variogram->Covariance(0.,0.)*(1. - Lambda[n_obs]);;
584 for(i=0;i<n_obs;i++) error += -Lambda[i]*B[i];
585
586 /*clean-up*/
587 *pprediction = prediction;
588 *perror = error;
589 xDelete<IssmPDouble>(x);
590 xDelete<IssmPDouble>(y);
591 xDelete<IssmPDouble>(obs);
592 xDelete<IssmPDouble>(A);
593 xDelete<IssmPDouble>(B);
594 xDelete<IssmPDouble>(Lambda);
595}/*}}}*/
596void Observations::InterpolationNearestNeighbor(IssmPDouble *pprediction,IssmPDouble x_interp,IssmPDouble y_interp,IssmPDouble radius){/*{{{*/
597
598 /*Intermediaries*/
599 IssmPDouble x,y,obs;
600
601 /*Get clostest observation*/
602 this->ClosestObservation(&x,&y,&obs,x_interp,y_interp,radius);
603
604 /*Assign output pointer*/
605 *pprediction = obs;
606}/*}}}*/
607void Observations::InterpolationV4(IssmPDouble *pprediction,IssmPDouble x_interp,IssmPDouble y_interp,IssmPDouble radius,int mindata,int maxdata){/*{{{*/
608 /* Reference: David T. Sandwell, Biharmonic spline interpolation of GEOS-3
609 * and SEASAT altimeter data, Geophysical Research Letters, 2, 139-142,
610 * 1987. Describes interpolation using value or gradient of value in any
611 * dimension.*/
612
613 /*Intermediaries*/
614 int i,j,n_obs;
615 IssmPDouble prediction,h;
616 IssmPDouble *x = NULL;
617 IssmPDouble *y = NULL;
618 IssmPDouble *obs = NULL;
619 IssmPDouble *Green = NULL;
620 IssmPDouble *weights = NULL;
621 IssmPDouble *g = NULL;
622
623 /*Some checks*/
624 _assert_(maxdata>0);
625 _assert_(pprediction);
626
627 /*Get list of observations for current point*/
628 this->ObservationList(&x,&y,&obs,&n_obs,x_interp,y_interp,radius,maxdata);
629
630 /*If we have less observations than mindata, return UNDEF*/
631 if(n_obs<mindata || n_obs<2){
632 prediction = UNDEF;
633 }
634 else{
635
636 /*Allocate intermediary matrix and vectors*/
637 Green = xNew<IssmPDouble>(n_obs*n_obs);
638 g = xNew<IssmPDouble>(n_obs);
639
640 /*First: distance vector*/
641 for(i=0;i<n_obs;i++){
642 h = sqrt( (x[i]-x_interp)*(x[i]-x_interp) + (y[i]-y_interp)*(y[i]-y_interp) );
643 if(h>0){
644 g[i] = h*h*(log(h)-1.);
645 }
646 else{
647 g[i] = 0.;
648 }
649 }
650
651 /*Build Green function matrix*/
652 for(i=0;i<n_obs;i++){
653 for(j=0;j<=i;j++){
654 h = sqrt( (x[i]-x[j])*(x[i]-x[j]) + (y[i]-y[j])*(y[i]-y[j]) );
655 if(h>0){
656 Green[j*n_obs+i] = h*h*(log(h)-1.);
657 }
658 else{
659 Green[j*n_obs+i] = 0.;
660 }
661 Green[i*n_obs+j] = Green[j*n_obs+i];
662 }
663 /*Zero diagonal (should be done already, but just in case)*/
664 Green[i*n_obs+i] = 0.;
665 }
666
667 /*Compute weights*/
668#if _HAVE_GSL_
669 DenseGslSolve(&weights,Green,obs,n_obs); // Green^-1 obs
670#else
671 _error_("GSL is required");
672#endif
673
674 /*Interpolate*/
675 prediction = 0;
676 for(i=0;i<n_obs;i++) prediction += weights[i]*g[i];
677
678 }
679
680 /*clean-up*/
681 *pprediction = prediction;
682 xDelete<IssmPDouble>(x);
683 xDelete<IssmPDouble>(y);
684 xDelete<IssmPDouble>(obs);
685 xDelete<IssmPDouble>(Green);
686 xDelete<IssmPDouble>(g);
687 xDelete<IssmPDouble>(weights);
688}/*}}}*/
689void Observations::QuadtreeColoring(IssmPDouble* A,IssmPDouble *x,IssmPDouble *y,int n){/*{{{*/
690
691 if(this->treetype!=1) _error_("Tree type is not quadtree");
692 int xi,yi,level;
693
694 for(int i=0;i<n;i++){
695 this->quadtree->IntergerCoordinates(&xi,&yi,x[i],y[i]);
696 this->quadtree->QuadtreeDepth(&level,xi,yi);
697 A[i]=(IssmPDouble)level;
698 }
699
700}/*}}}*/
701void Observations::Variomap(IssmPDouble* gamma,IssmPDouble *x,int n){/*{{{*/
702
703 /*Output and Intermediaries*/
704 int i,j,k;
705 IssmPDouble distance;
706 Observation *observation1 = NULL;
707 Observation *observation2 = NULL;
708
709 IssmPDouble *counter = xNew<IssmPDouble>(n);
710 for(j=0;j<n;j++) counter[j] = 0.0;
711 for(j=0;j<n;j++) gamma[j] = 0.0;
712
713 for(i=0;i<this->Size();i++){
714 observation1=xDynamicCast<Observation*>(this->GetObjectByOffset(i));
715
716 for(j=i+1;j<this->Size();j++){
717 observation2=xDynamicCast<Observation*>(this->GetObjectByOffset(j));
718
719 distance=sqrt(pow(observation1->x - observation2->x,2) + pow(observation1->y - observation2->y,2));
720 if(distance>x[n-1]) continue;
721
722 int index = int(distance/(x[1]-x[0]));
723 if(index>n-1) index = n-1;
724 if(index<0) index = 0;
725
726 gamma[index] += 1./2.*pow(observation1->value - observation2->value,2);
727 counter[index] += 1.;
728 }
729 }
730
731 /*Normalize semivariogram*/
732 gamma[0]=0.;
733 for(k=0;k<n;k++){
734 if(counter[k]) gamma[k] = gamma[k]/counter[k];
735 }
736
737 /*Assign output pointer*/
738 xDelete<IssmPDouble>(counter);
739}/*}}}*/
740
Note: See TracBrowser for help on using the repository browser.