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

Last change on this file since 15396 was 15396, checked in by Mathieu Morlighem, 12 years ago

merged trunk-jpl and trunk for revision 15394
File size: 16.7 KB
RevLine 
[14996]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
[15012]18#include "../Options/Options.h"
[14996]19#include "./Observations.h"
20#include "./Observation.h"
[15067]21#include "../../datastructures/datastructures.h"
[14996]22#include "../../shared/shared.h"
23
24#include "./Quadtree.h"
25#include "./Variogram.h"
[14997]26#include "../../toolkits/toolkits.h"
[14996]27
28using namespace std;
29/*}}}*/
30
31/*Object constructors and destructor*/
32/*FUNCTION Observations::Observations(){{{*/
33Observations::Observations(){
34 this->quadtree = NULL;
35 return;
36}
37/*}}}*/
38/*FUNCTION Observations::Observations(IssmPDouble* observations_list,IssmPDouble* x,IssmPDouble* y,int n,Options* options){{{*/
39Observations::Observations(IssmPDouble* observations_list,IssmPDouble* x,IssmPDouble* y,int n,Options* options){
40
41 /*Intermediaries*/
42 int i,maxdepth,level,counter,index;
43 int xi,yi;
44 IssmPDouble xmin,xmax,ymin,ymax;
45 IssmPDouble offset,minlength,minspacing,mintrimming,maxtrimming;
46 Observation *observation = NULL;
47
48 /*Check that observations is not empty*/
49 if(n==0) _error_("No observation found");
50
51 /*Get extrema*/
52 xmin=x[0]; ymin=y[0];
53 xmax=x[0]; ymax=y[0];
54 for(i=1;i<n;i++){
55 xmin=min(xmin,x[i]); ymin=min(ymin,y[i]);
56 xmax=max(xmax,x[i]); ymax=max(ymax,y[i]);
57 }
58 offset=0.05*(xmax-xmin); xmin-=offset; xmax+=offset;
59 offset=0.05*(ymax-ymin); ymin-=offset; ymax+=offset;
60
61 /*Get trimming limits*/
62 options->Get(&mintrimming,"mintrimming",-1.e+21);
63 options->Get(&maxtrimming,"maxtrimming",+1.e+21);
64 options->Get(&minspacing,"minspacing",0.01);
65 if(minspacing<=0) _error_("minspacing must > 0");
66
67 /*Get Minimum box size*/
68 if(options->GetOption("boxlength")){
69 options->Get(&minlength,"boxlength");
70 if(minlength<=0)_error_("boxlength should be a positive number");
71 maxdepth=reCast<int,IssmPDouble>(log(max(xmax-xmin,ymax-ymin)/minlength +1)/log(2.0));
72 }
73 else{
74 maxdepth = 30;
75 minlength=max(xmax-xmin,ymax-ymin)/IssmPDouble((1L<<maxdepth)-1);
76 }
77
78 /*Initialize Quadtree*/
[15100]79 _printf0_("Generating quadtree with a maximum box size " << minlength << " (depth=" << maxdepth << ")... ");
[14996]80 this->quadtree = new Quadtree(xmin,xmax,ymin,ymax,maxdepth);
81
82 /*Add observations one by one*/
83 counter = 0;
84 for(i=0;i<n;i++){
85
86 /*First check limits*/
87 if(observations_list[i]>maxtrimming) continue;
88 if(observations_list[i]<mintrimming) continue;
89
90 /*First check that this observation is not too close from another one*/
91 this->quadtree->ClosestObs(&index,x[i],y[i]);
92 if(index>=0){
93 observation=dynamic_cast<Observation*>(this->GetObjectByOffset(index));
94 if(pow(observation->x-x[i],2)+pow(observation->y-y[i],2) < minspacing) continue;
95 }
96
97 this->quadtree->IntergerCoordinates(&xi,&yi,x[i],y[i]);
98 this->quadtree->QuadtreeDepth2(&level,xi,yi);
99 if((int)level <= maxdepth){
100 observation = new Observation(x[i],y[i],xi,yi,counter++,observations_list[i]);
101 this->quadtree->Add(observation);
102 this->AddObject(observation);
103 }
104 else{
105 /*We need to average with the current observations*/
106 this->quadtree->AddAndAverage(x[i],y[i],observations_list[i]);
107 }
108 }
[15104]109 _printf0_("done\n");
[15100]110 _printf0_("Initial number of observations: " << n << "\n");
111 _printf0_(" Final number of observations: " << this->quadtree->NbObs << "\n");
[14996]112}
113/*}}}*/
114/*FUNCTION Observations::~Observations(){{{*/
115Observations::~Observations(){
116 delete quadtree;
117 return;
118}
119/*}}}*/
120
121/*Methods*/
122/*FUNCTION Observations::ClosestObservation{{{*/
123void Observations::ClosestObservation(IssmPDouble *px,IssmPDouble *py,IssmPDouble *pobs,IssmPDouble x_interp,IssmPDouble y_interp,IssmPDouble radius){
124
125 /*Output and Intermediaries*/
126 int nobs,i,index;
127 IssmPDouble hmin,h2,hmin2,radius2;
128 int *indices = NULL;
129 Observation *observation = NULL;
130
131 /*If radius is not provided or is 0, return all observations*/
132 if(radius==0) radius=this->quadtree->root->length;
133
134 /*First, find closest point in Quadtree (fast but might not be the true closest obs)*/
135 this->quadtree->ClosestObs(&index,x_interp,y_interp);
136 if(index>=0){
137 observation=dynamic_cast<Observation*>(this->GetObjectByOffset(index));
138 hmin = sqrt((observation->x-x_interp)*(observation->x-x_interp) + (observation->y-y_interp)*(observation->y-y_interp));
139 if(hmin<radius) radius=hmin;
140 }
141
142 /*Compute radius square*/
143 radius2 = radius*radius;
144
145 /*Find all observations that are in radius*/
146 this->quadtree->RangeSearch(&indices,&nobs,x_interp,y_interp,radius);
147 for (i=0;i<nobs;i++){
148 observation=dynamic_cast<Observation*>(this->GetObjectByOffset(indices[i]));
149 h2 = (observation->x-x_interp)*(observation->x-x_interp) + (observation->y-y_interp)*(observation->y-y_interp);
150 if(i==0){
151 hmin2 = h2;
152 index = indices[i];
153 }
154 else{
155 if(h2<hmin2){
156 hmin2 = h2;
157 index = indices[i];
158 }
159 }
160 }
161
162 /*Assign output pointer*/
163 if(index>=0){
164 observation=dynamic_cast<Observation*>(this->GetObjectByOffset(index));
165 *px=observation->x;
166 *py=observation->y;
167 *pobs=observation->value;
168 }
169 else{
170
171 *px=UNDEF;
172 *py=UNDEF;
173 *pobs=UNDEF;
174 }
175 xDelete<int>(indices);
176
177}/*}}}*/
178/*FUNCTION Observations::Distances{{{*/
179void Observations::Distances(IssmPDouble* distances,IssmPDouble *x,IssmPDouble *y,int n,IssmPDouble radius){
180
181 IssmPDouble xi,yi,obs;
182
183 for(int i=0;i<n;i++){
184 this->ClosestObservation(&xi,&yi,&obs,x[i],y[i],radius);
185 if(xi==UNDEF && yi==UNDEF)
186 distances[i]=UNDEF;
187 else
188 distances[i]=sqrt( (x[i]-xi)*(x[i]-xi) + (y[i]-yi)*(y[i]-yi) );
189 }
190}/*}}}*/
191/*FUNCTION Observations::ObservationList(IssmPDouble **px,IssmPDouble **py,IssmPDouble **pobs,int* pnobs,IssmPDouble x_interp,IssmPDouble y_interp,IssmPDouble radius,int maxdata){{{*/
192void Observations::ObservationList(IssmPDouble **px,IssmPDouble **py,IssmPDouble **pobs,int* pnobs,IssmPDouble x_interp,IssmPDouble y_interp,IssmPDouble radius,int maxdata){
193
194 /*Output and Intermediaries*/
195 bool stop;
196 int nobs,tempnobs,i,j,k,n,counter;
197 IssmPDouble h2,radius2;
198 int *indices = NULL;
199 int *tempindices = NULL;
200 IssmPDouble *dists = NULL;
201 IssmPDouble *x = NULL;
202 IssmPDouble *y = NULL;
203 IssmPDouble *obs = NULL;
204 Observation *observation = NULL;
205
206 /*If radius is not provided or is 0, return all observations*/
207 if(radius==0) radius=this->quadtree->root->length;
208
209 /*Compute radius square*/
210 radius2 = radius*radius;
211
212 /*Find all observations that are in radius*/
213 this->quadtree->RangeSearch(&tempindices,&tempnobs,x_interp,y_interp,radius);
214 if(tempnobs){
215 indices = xNew<int>(tempnobs);
216 dists = xNew<IssmPDouble>(tempnobs);
217 }
218 nobs = 0;
219 for (i=0;i<tempnobs;i++){
220 observation=dynamic_cast<Observation*>(this->GetObjectByOffset(tempindices[i]));
221 h2 = (observation->x-x_interp)*(observation->x-x_interp) + (observation->y-y_interp)*(observation->y-y_interp);
222
223 if(nobs==maxdata && h2>radius2) continue;
224 if(nobs<=maxdata){
225 indices[nobs] = tempindices[i];
226 dists[nobs] = h2;
227 nobs++;
228 }
229 if(nobs==1) continue;
230
231 /*Sort all dists up to now*/
232 n=nobs-1;
233 stop = false;
234 for(k=0;k<n-1;k++){
235 if(h2<dists[k]){
236 counter=1;
237 for(int jj=k;jj<n;jj++){
238 j = n-counter;
239 dists[j+1] = dists[j];
240 indices[j+1] = indices[j];
241 counter++;
242 }
243 dists[k] = h2;
244 indices[k] = tempindices[i];
245 stop = true;
246 break;
247 }
248 if(stop) break;
249 }
250 }
251 xDelete<IssmPDouble>(dists);
252 xDelete<int>(tempindices);
253
254 if(nobs){
255 /*Allocate vectors*/
256 x = xNew<IssmPDouble>(nobs);
257 y = xNew<IssmPDouble>(nobs);
258 obs = xNew<IssmPDouble>(nobs);
259
260 /*Loop over all observations and fill in x, y and obs*/
261 for (i=0;i<nobs;i++){
262 observation=dynamic_cast<Observation*>(this->GetObjectByOffset(indices[i]));
263 observation->WriteXYObs(&x[i],&y[i],&obs[i]);
264 }
265 }
266
267 /*Assign output pointer*/
268 xDelete<int>(indices);
269 *px=x;
270 *py=y;
271 *pobs=obs;
272 *pnobs=nobs;
273}/*}}}*/
274/*FUNCTION Observations::ObservationList(IssmPDouble **px,IssmPDouble **py,IssmPDouble **pobs,int* pnobs){{{*/
275void Observations::ObservationList(IssmPDouble **px,IssmPDouble **py,IssmPDouble **pobs,int* pnobs){
276
277 /*Output and Intermediaries*/
278 int nobs;
279 IssmPDouble *x = NULL;
280 IssmPDouble *y = NULL;
281 IssmPDouble *obs = NULL;
282 Observation *observation = NULL;
283
284 nobs = this->Size();
285
286 if(nobs){
287 x = xNew<IssmPDouble>(nobs);
288 y = xNew<IssmPDouble>(nobs);
289 obs = xNew<IssmPDouble>(nobs);
290 for(int i=0;i<this->Size();i++){
291 observation=dynamic_cast<Observation*>(this->GetObjectByOffset(i));
292 observation->WriteXYObs(&x[i],&y[i],&obs[i]);
293 }
294 }
295
296 /*Assign output pointer*/
297 *px=x;
298 *py=y;
299 *pobs=obs;
300 *pnobs=nobs;
301}/*}}}*/
302/*FUNCTION Observations::InterpolationIDW{{{*/
303void Observations::InterpolationIDW(IssmPDouble *pprediction,IssmPDouble x_interp,IssmPDouble y_interp,IssmPDouble radius,int mindata,int maxdata,IssmPDouble power){
304
305 /*Intermediaries*/
306 int i,n_obs;
307 IssmPDouble prediction;
308 IssmPDouble numerator,denominator,h,weight;
309 IssmPDouble *x = NULL;
310 IssmPDouble *y = NULL;
311 IssmPDouble *obs = NULL;
312
313 /*Some checks*/
314 _assert_(maxdata>0);
315 _assert_(pprediction);
316 _assert_(power>0);
317
318 /*If radius is not provided or is 0, return all observations*/
319 if(radius==0) radius=this->quadtree->root->length;
320
321 /*Get list of observations for current point*/
322 this->ObservationList(&x,&y,&obs,&n_obs,x_interp,y_interp,radius,maxdata);
323
324 /*If we have less observations than mindata, return UNDEF*/
325 if(n_obs<mindata){
326 prediction = UNDEF;
327 }
328 else{
329 numerator = 0.;
330 denominator = 0.;
331 for(i=0;i<n_obs;i++){
332 h = sqrt( (x[i]-x_interp)*(x[i]-x_interp) + (y[i]-y_interp)*(y[i]-y_interp));
333 if (h<0.0000001){
334 numerator = obs[i];
335 denominator = 1.;
336 break;
337 }
338 weight = 1./pow(h,power);
339 numerator += weight*obs[i];
340 denominator += weight;
341 }
342 prediction = numerator/denominator;
343 }
344
345 /*clean-up*/
346 *pprediction = prediction;
347 xDelete<IssmPDouble>(x);
348 xDelete<IssmPDouble>(y);
349 xDelete<IssmPDouble>(obs);
350}/*}}}*/
351/*FUNCTION Observations::InterpolationKriging{{{*/
352void Observations::InterpolationKriging(IssmPDouble *pprediction,IssmPDouble *perror,IssmPDouble x_interp,IssmPDouble y_interp,IssmPDouble radius,int mindata,int maxdata,Variogram* variogram){
353
354 /*Intermediaries*/
355 int i,j,n_obs;
356 IssmPDouble prediction,error;
357 IssmPDouble numerator,denominator,ratio;
358 IssmPDouble *x = NULL;
359 IssmPDouble *y = NULL;
360 IssmPDouble *obs = NULL;
361 IssmPDouble *Gamma = NULL;
362 IssmPDouble *GinvG0 = NULL;
363 IssmPDouble *Ginv1 = NULL;
364 IssmPDouble *GinvZ = NULL;
365 IssmPDouble *gamma0 = NULL;
366 IssmPDouble *ones = NULL;
367
368 /*Some checks*/
369 _assert_(mindata>0 && maxdata>0);
370 _assert_(pprediction && perror);
371
372 /*If radius is not provided or is 0, return all observations*/
373 if(radius==0) radius=this->quadtree->root->length;
374
375 /*Get list of observations for current point*/
376 this->ObservationList(&x,&y,&obs,&n_obs,x_interp,y_interp,radius,maxdata);
377
378 /*If we have less observations than mindata, return UNDEF*/
379 if(n_obs<mindata){
380 *pprediction = -999.0;
381 *perror = -999.0;
382 return;
383 }
384
385 /*Allocate intermediary matrix and vectors*/
386 Gamma = xNew<IssmPDouble>(n_obs*n_obs);
387 gamma0 = xNew<IssmPDouble>(n_obs);
388 ones = xNew<IssmPDouble>(n_obs);
389
390 /*First: Create semivariogram matrix for observations*/
391 for(i=0;i<n_obs;i++){
392 for(j=0;j<=i;j++){
393 //Gamma[i*n_obs+j] = variogram->SemiVariogram(x[i]-x[j],y[i]-y[j]);
394 Gamma[i*n_obs+j] = variogram->Covariance(x[i]-x[j],y[i]-y[j]);
395 Gamma[j*n_obs+i] = Gamma[i*n_obs+j];
396 }
397 }
398 for(i=0;i<n_obs;i++) ones[i]=1;
399
400 /*Get semivariogram vector associated to this location*/
401 //for(i=0;i<n_obs;i++) gamma0[i] = variogram->SemiVariogram(x[i]-x_interp,y[i]-y_interp);
402 for(i=0;i<n_obs;i++) gamma0[i] = variogram->Covariance(x[i]-x_interp,y[i]-y_interp);
403
404 /*Solve the three linear systems*/
405#if _HAVE_GSL_
406 DenseGslSolve(&GinvG0,Gamma,gamma0,n_obs); // Gamma^-1 gamma0
407 DenseGslSolve(&Ginv1, Gamma,ones,n_obs); // Gamma^-1 ones
408 DenseGslSolve(&GinvZ, Gamma,obs,n_obs); // Gamma^-1 Z
409#else
410 _error_("GSL is required");
411#endif
412
413 /*Prepare predictor*/
414 numerator=-1.; denominator=0.;
415 for(i=0;i<n_obs;i++) numerator +=GinvG0[i];
416 for(i=0;i<n_obs;i++) denominator+=Ginv1[i];
417 ratio=numerator/denominator;
418
419 prediction = 0.;
420 error = - numerator*numerator/denominator;
421 for(i=0;i<n_obs;i++) prediction += (gamma0[i]-ratio)*GinvZ[i];
422 for(i=0;i<n_obs;i++) error += gamma0[i]*GinvG0[i];
423
424 /*clean-up*/
425 *pprediction = prediction;
426 *perror = error;
427 xDelete<IssmPDouble>(x);
428 xDelete<IssmPDouble>(y);
429 xDelete<IssmPDouble>(obs);
430 xDelete<IssmPDouble>(Gamma);
431 xDelete<IssmPDouble>(gamma0);
432 xDelete<IssmPDouble>(ones);
433 xDelete<IssmPDouble>(GinvG0);
434 xDelete<IssmPDouble>(Ginv1);
435 xDelete<IssmPDouble>(GinvZ);
436
437}/*}}}*/
438/*FUNCTION Observations::InterpolationNearestNeighbor{{{*/
439void Observations::InterpolationNearestNeighbor(IssmPDouble *pprediction,IssmPDouble x_interp,IssmPDouble y_interp,IssmPDouble radius){
440
441 /*Intermediaries*/
442 IssmPDouble x,y,obs;
443
444 /*Get clostest observation*/
445 this->ClosestObservation(&x,&y,&obs,x_interp,y_interp,radius);
446
447 /*Assign output pointer*/
448 *pprediction = obs;
449}/*}}}*/
450/*FUNCTION Observations::InterpolationV4{{{*/
451void Observations::InterpolationV4(IssmPDouble *pprediction,IssmPDouble x_interp,IssmPDouble y_interp,IssmPDouble radius,int mindata,int maxdata){
452 /* Reference: David T. Sandwell, Biharmonic spline interpolation of GEOS-3
453 * and SEASAT altimeter data, Geophysical Research Letters, 2, 139-142,
454 * 1987. Describes interpolation using value or gradient of value in any
455 * dimension.*/
456
457 /*Intermediaries*/
458 int i,j,n_obs;
459 IssmPDouble prediction,h;
460 IssmPDouble *x = NULL;
461 IssmPDouble *y = NULL;
462 IssmPDouble *obs = NULL;
463 IssmPDouble *Green = NULL;
464 IssmPDouble *weights = NULL;
465 IssmPDouble *g = NULL;
466
467 /*Some checks*/
468 _assert_(maxdata>0);
469 _assert_(pprediction);
470
471 /*If radius is not provided or is 0, return all observations*/
472 if(radius==0) radius=this->quadtree->root->length;
473
474 /*Get list of observations for current point*/
475 this->ObservationList(&x,&y,&obs,&n_obs,x_interp,y_interp,radius,maxdata);
476
477 /*If we have less observations than mindata, return UNDEF*/
478 if(n_obs<mindata || n_obs<2){
479 prediction = UNDEF;
480 }
481 else{
482
483 /*Allocate intermediary matrix and vectors*/
484 Green = xNew<IssmPDouble>(n_obs*n_obs);
485 g = xNew<IssmPDouble>(n_obs);
486
487 /*First: distance vector*/
488 for(i=0;i<n_obs;i++){
489 h = sqrt( (x[i]-x_interp)*(x[i]-x_interp) + (y[i]-y_interp)*(y[i]-y_interp) );
490 if(h>0){
491 g[i] = h*h*(log(h)-1.);
492 }
493 else{
494 g[i] = 0.;
495 }
496 }
497
498 /*Build Green function matrix*/
499 for(i=0;i<n_obs;i++){
500 for(j=0;j<=i;j++){
501 h = sqrt( (x[i]-x[j])*(x[i]-x[j]) + (y[i]-y[j])*(y[i]-y[j]) );
502 if(h>0){
503 Green[j*n_obs+i] = h*h*(log(h)-1.);
504 }
505 else{
506 Green[j*n_obs+i] = 0.;
507 }
508 Green[i*n_obs+j] = Green[j*n_obs+i];
509 }
510 }
511
512 /*Compute weights*/
513#if _HAVE_GSL_
514 DenseGslSolve(&weights,Green,obs,n_obs); // Green^-1 obs
515#else
516 _error_("GSL is required");
517#endif
518
519 /*Interpolate*/
520 prediction = 0;
521 for(i=0;i<n_obs;i++) prediction += weights[i]*g[i];
522
523 }
524
525 /*clean-up*/
526 *pprediction = prediction;
527 xDelete<IssmPDouble>(x);
528 xDelete<IssmPDouble>(y);
529 xDelete<IssmPDouble>(obs);
530 xDelete<IssmPDouble>(Green);
531 xDelete<IssmPDouble>(g);
532 xDelete<IssmPDouble>(weights);
533}/*}}}*/
534/*FUNCTION Observations::QuadtreeColoring{{{*/
535void Observations::QuadtreeColoring(IssmPDouble* A,IssmPDouble *x,IssmPDouble *y,int n){
536
537 int xi,yi,level;
538
539 for(int i=0;i<n;i++){
540 this->quadtree->IntergerCoordinates(&xi,&yi,x[i],y[i]);
541 this->quadtree->QuadtreeDepth(&level,xi,yi);
542 A[i]=(IssmPDouble)level;
543 }
544
545}/*}}}*/
546/*FUNCTION Observations::Variomap{{{*/
547void Observations::Variomap(IssmPDouble* gamma,IssmPDouble *x,int n){
548
549 /*Output and Intermediaries*/
550 int i,j,k;
551 IssmPDouble distance;
552 Observation *observation1 = NULL;
553 Observation *observation2 = NULL;
554
555 IssmPDouble *counter = xNew<IssmPDouble>(n);
556 for(j=0;j<n;j++) counter[j] = 0.0;
557 for(j=0;j<n;j++) gamma[j] = 0.0;
558
559 for(i=0;i<this->Size();i++){
560 observation1=dynamic_cast<Observation*>(this->GetObjectByOffset(i));
561
562 for(j=i+1;j<this->Size();j++){
563 observation2=dynamic_cast<Observation*>(this->GetObjectByOffset(j));
564
565 distance=sqrt(pow(observation1->x - observation2->x,2.) + pow(observation1->y - observation2->y,2.));
566 if(distance>x[n-1]) continue;
567
568 int index = int(distance/(x[1]-x[0]));
569 if(index>n-1) index = n-1;
570 if(index<0) index = 0;
571
572 gamma[index] += 1./2.*pow(observation1->value - observation2->value,2.);
573 counter[index] += 1.;
574 }
575 }
576
577 /*Normalize semivariogram*/
578 gamma[0]=0.;
579 for(k=0;k<n;k++){
580 if(counter[k]) gamma[k] = gamma[k]/counter[k];
581 }
582
583 /*Assign output pointer*/
584 xDelete<IssmPDouble>(counter);
585}/*}}}*/
Note: See TracBrowser for help on using the repository browser.