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

Last change on this file since 24686 was 24686, checked in by Mathieu Morlighem, 5 years ago

merged trunk-jpl and trunk for revision 24684
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 _printf_("Tree type "<<this->treetype<<" not supported yet (1: quadtree, 2: covertree)");
73 }
74 return;
75}
76/*}}}*/
77
78/*Initialize data structures*/
79void Observations::InitCovertree(IssmPDouble* observations_list,IssmPDouble* x,IssmPDouble* y,int n,Options* options){/*{{{*/
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}
134/*}}}*/
135void Observations::InitQuadtree(IssmPDouble* observations_list,IssmPDouble* x,IssmPDouble* y,int n,Options* options){/*{{{*/
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}
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::ClosestObservationCovertree(IssmPDouble *px,IssmPDouble *py,IssmPDouble *pobs,IssmPDouble x_interp,IssmPDouble y_interp,IssmPDouble radius){/*{{{*/
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}/*}}}*/
248void Observations::ClosestObservationQuadtree(IssmPDouble *px,IssmPDouble *py,IssmPDouble *pobs,IssmPDouble x_interp,IssmPDouble y_interp,IssmPDouble radius){/*{{{*/
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}/*}}}*/
302void Observations::Distances(IssmPDouble* distances,IssmPDouble *x,IssmPDouble *y,int n,IssmPDouble radius){/*{{{*/
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}/*}}}*/
316void Observations::ObservationList(IssmPDouble **px,IssmPDouble **py,IssmPDouble **pobs,int* pnobs){/*{{{*/
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}/*}}}*/
343void Observations::ObservationList(IssmPDouble **px,IssmPDouble **py,IssmPDouble **pobs,int* pnobs,IssmPDouble x_interp,IssmPDouble y_interp,IssmPDouble radius,int maxdata){/*{{{*/
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}/*}}}*/
356void Observations::ObservationListCovertree(double **px,double **py,double **pobs,int* pnobs,double x_interp,double y_interp,double radius,int maxdata){/*{{{*/
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}/*}}}*/
399void Observations::ObservationListQuadtree(IssmPDouble **px,IssmPDouble **py,IssmPDouble **pobs,int* pnobs,IssmPDouble x_interp,IssmPDouble y_interp,IssmPDouble radius,int maxdata){/*{{{*/
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}/*}}}*/
481void Observations::InterpolationIDW(IssmPDouble *pprediction,IssmPDouble x_interp,IssmPDouble y_interp,IssmPDouble radius,int mindata,int maxdata,IssmPDouble power){/*{{{*/
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}/*}}}*/
526void Observations::InterpolationKriging(IssmPDouble *pprediction,IssmPDouble *perror,IssmPDouble x_interp,IssmPDouble y_interp,IssmPDouble radius,int mindata,int maxdata,Variogram* variogram){/*{{{*/
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}/*}}}*/
599void Observations::InterpolationNearestNeighbor(IssmPDouble *pprediction,IssmPDouble x_interp,IssmPDouble y_interp,IssmPDouble radius){/*{{{*/
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}/*}}}*/
610void Observations::InterpolationV4(IssmPDouble *pprediction,IssmPDouble x_interp,IssmPDouble y_interp,IssmPDouble radius,int mindata,int maxdata){/*{{{*/
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}/*}}}*/
692void Observations::QuadtreeColoring(IssmPDouble* A,IssmPDouble *x,IssmPDouble *y,int n){/*{{{*/
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}/*}}}*/
704void Observations::Variomap(IssmPDouble* gamma,IssmPDouble *x,int n){/*{{{*/
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}/*}}}*/
Note: See TracBrowser for help on using the repository browser.