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

Last change on this file since 16560 was 16560, checked in by Mathieu Morlighem, 11 years ago

merged trunk-jpl and trunk for revision 16554
File size: 16.6 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 "./Variogram.h"
26#include "../../toolkits/toolkits.h"
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*/
79 _printf0_("Generating quadtree with a maximum box size " << minlength << " (depth=" << maxdepth << ")... ");
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 }
109 _printf0_("done\n");
110 _printf0_("Initial number of observations: " << n << "\n");
111 _printf0_(" Final number of observations: " << this->quadtree->NbObs << "\n");
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;
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 /*Find all observations that are in radius*/
143 this->quadtree->RangeSearch(&indices,&nobs,x_interp,y_interp,radius);
144 for (i=0;i<nobs;i++){
145 observation=dynamic_cast<Observation*>(this->GetObjectByOffset(indices[i]));
146 h2 = (observation->x-x_interp)*(observation->x-x_interp) + (observation->y-y_interp)*(observation->y-y_interp);
147 if(i==0){
148 hmin2 = h2;
149 index = indices[i];
150 }
151 else{
152 if(h2<hmin2){
153 hmin2 = h2;
154 index = indices[i];
155 }
156 }
157 }
158
159 /*Assign output pointer*/
160 if(index>=0){
161 observation=dynamic_cast<Observation*>(this->GetObjectByOffset(index));
162 *px=observation->x;
163 *py=observation->y;
164 *pobs=observation->value;
165 }
166 else{
167
168 *px=UNDEF;
169 *py=UNDEF;
170 *pobs=UNDEF;
171 }
172 xDelete<int>(indices);
173
174}/*}}}*/
175/*FUNCTION Observations::Distances{{{*/
176void Observations::Distances(IssmPDouble* distances,IssmPDouble *x,IssmPDouble *y,int n,IssmPDouble radius){
177
178 IssmPDouble xi,yi,obs;
179
180 for(int i=0;i<n;i++){
181 this->ClosestObservation(&xi,&yi,&obs,x[i],y[i],radius);
182 if(xi==UNDEF && yi==UNDEF)
183 distances[i]=UNDEF;
184 else
185 distances[i]=sqrt( (x[i]-xi)*(x[i]-xi) + (y[i]-yi)*(y[i]-yi) );
186 }
187}/*}}}*/
188/*FUNCTION Observations::ObservationList(IssmPDouble **px,IssmPDouble **py,IssmPDouble **pobs,int* pnobs,IssmPDouble x_interp,IssmPDouble y_interp,IssmPDouble radius,int maxdata){{{*/
189void Observations::ObservationList(IssmPDouble **px,IssmPDouble **py,IssmPDouble **pobs,int* pnobs,IssmPDouble x_interp,IssmPDouble y_interp,IssmPDouble radius,int maxdata){
190
191 /*Output and Intermediaries*/
192 bool stop;
193 int nobs,tempnobs,i,j,k,n,counter;
194 IssmPDouble h2,radius2;
195 int *indices = NULL;
196 int *tempindices = NULL;
197 IssmPDouble *dists = NULL;
198 IssmPDouble *x = NULL;
199 IssmPDouble *y = NULL;
200 IssmPDouble *obs = NULL;
201 Observation *observation = NULL;
202
203 /*If radius is not provided or is 0, return all observations*/
204 if(radius==0) radius=this->quadtree->root->length;
205
206 /*Compute radius square*/
207 radius2 = radius*radius;
208
209 /*Find all observations that are in radius*/
210 this->quadtree->RangeSearch(&tempindices,&tempnobs,x_interp,y_interp,radius);
211 if(tempnobs){
212 indices = xNew<int>(tempnobs);
213 dists = xNew<IssmPDouble>(tempnobs);
214 }
215 nobs = 0;
216 for (i=0;i<tempnobs;i++){
217 observation=dynamic_cast<Observation*>(this->GetObjectByOffset(tempindices[i]));
218 h2 = (observation->x-x_interp)*(observation->x-x_interp) + (observation->y-y_interp)*(observation->y-y_interp);
219
220 if(nobs==maxdata && h2>radius2) continue;
221 if(nobs<=maxdata){
222 indices[nobs] = tempindices[i];
223 dists[nobs] = h2;
224 nobs++;
225 }
226 if(nobs==1) continue;
227
228 /*Sort all dists up to now*/
229 n=nobs-1;
230 stop = false;
231 for(k=0;k<n-1;k++){
232 if(h2<dists[k]){
233 counter=1;
234 for(int jj=k;jj<n;jj++){
235 j = n-counter;
236 dists[j+1] = dists[j];
237 indices[j+1] = indices[j];
238 counter++;
239 }
240 dists[k] = h2;
241 indices[k] = tempindices[i];
242 stop = true;
243 break;
244 }
245 if(stop) break;
246 }
247 }
248 xDelete<IssmPDouble>(dists);
249 xDelete<int>(tempindices);
250
251 if(nobs){
252 /*Allocate vectors*/
253 x = xNew<IssmPDouble>(nobs);
254 y = xNew<IssmPDouble>(nobs);
255 obs = xNew<IssmPDouble>(nobs);
256
257 /*Loop over all observations and fill in x, y and obs*/
258 for (i=0;i<nobs;i++){
259 observation=dynamic_cast<Observation*>(this->GetObjectByOffset(indices[i]));
260 observation->WriteXYObs(&x[i],&y[i],&obs[i]);
261 }
262 }
263
264 /*Assign output pointer*/
265 xDelete<int>(indices);
266 *px=x;
267 *py=y;
268 *pobs=obs;
269 *pnobs=nobs;
270}/*}}}*/
271/*FUNCTION Observations::ObservationList(IssmPDouble **px,IssmPDouble **py,IssmPDouble **pobs,int* pnobs){{{*/
272void Observations::ObservationList(IssmPDouble **px,IssmPDouble **py,IssmPDouble **pobs,int* pnobs){
273
274 /*Output and Intermediaries*/
275 int nobs;
276 IssmPDouble *x = NULL;
277 IssmPDouble *y = NULL;
278 IssmPDouble *obs = NULL;
279 Observation *observation = NULL;
280
281 nobs = this->Size();
282
283 if(nobs){
284 x = xNew<IssmPDouble>(nobs);
285 y = xNew<IssmPDouble>(nobs);
286 obs = xNew<IssmPDouble>(nobs);
287 for(int i=0;i<this->Size();i++){
288 observation=dynamic_cast<Observation*>(this->GetObjectByOffset(i));
289 observation->WriteXYObs(&x[i],&y[i],&obs[i]);
290 }
291 }
292
293 /*Assign output pointer*/
294 *px=x;
295 *py=y;
296 *pobs=obs;
297 *pnobs=nobs;
298}/*}}}*/
299/*FUNCTION Observations::InterpolationIDW{{{*/
300void Observations::InterpolationIDW(IssmPDouble *pprediction,IssmPDouble x_interp,IssmPDouble y_interp,IssmPDouble radius,int mindata,int maxdata,IssmPDouble power){
301
302 /*Intermediaries*/
303 int i,n_obs;
304 IssmPDouble prediction;
305 IssmPDouble numerator,denominator,h,weight;
306 IssmPDouble *x = NULL;
307 IssmPDouble *y = NULL;
308 IssmPDouble *obs = NULL;
309
310 /*Some checks*/
311 _assert_(maxdata>0);
312 _assert_(pprediction);
313 _assert_(power>0);
314
315 /*If radius is not provided or is 0, return all observations*/
316 if(radius==0) radius=this->quadtree->root->length;
317
318 /*Get list of observations for current point*/
319 this->ObservationList(&x,&y,&obs,&n_obs,x_interp,y_interp,radius,maxdata);
320
321 /*If we have less observations than mindata, return UNDEF*/
322 if(n_obs<mindata){
323 prediction = UNDEF;
324 }
325 else{
326 numerator = 0.;
327 denominator = 0.;
328 for(i=0;i<n_obs;i++){
329 h = sqrt( (x[i]-x_interp)*(x[i]-x_interp) + (y[i]-y_interp)*(y[i]-y_interp));
330 if (h<0.0000001){
331 numerator = obs[i];
332 denominator = 1.;
333 break;
334 }
335 weight = 1./pow(h,power);
336 numerator += weight*obs[i];
337 denominator += weight;
338 }
339 prediction = numerator/denominator;
340 }
341
342 /*clean-up*/
343 *pprediction = prediction;
344 xDelete<IssmPDouble>(x);
345 xDelete<IssmPDouble>(y);
346 xDelete<IssmPDouble>(obs);
347}/*}}}*/
348/*FUNCTION Observations::InterpolationKriging{{{*/
349void Observations::InterpolationKriging(IssmPDouble *pprediction,IssmPDouble *perror,IssmPDouble x_interp,IssmPDouble y_interp,IssmPDouble radius,int mindata,int maxdata,Variogram* variogram){
350
351 /*Intermediaries*/
352 int i,j,n_obs;
353 IssmPDouble prediction,error;
354 IssmPDouble numerator,denominator,ratio;
355 IssmPDouble *x = NULL;
356 IssmPDouble *y = NULL;
357 IssmPDouble *obs = NULL;
358 IssmPDouble *Gamma = NULL;
359 IssmPDouble *GinvG0 = NULL;
360 IssmPDouble *Ginv1 = NULL;
361 IssmPDouble *GinvZ = NULL;
362 IssmPDouble *gamma0 = NULL;
363 IssmPDouble *ones = NULL;
364
365 /*Some checks*/
366 _assert_(mindata>0 && maxdata>0);
367 _assert_(pprediction && perror);
368
369 /*If radius is not provided or is 0, return all observations*/
370 if(radius==0) radius=this->quadtree->root->length;
371
372 /*Get list of observations for current point*/
373 this->ObservationList(&x,&y,&obs,&n_obs,x_interp,y_interp,radius,maxdata);
374
375 /*If we have less observations than mindata, return UNDEF*/
376 if(n_obs<mindata){
377 *pprediction = -999.0;
378 *perror = -999.0;
379 return;
380 }
381
382 /*Allocate intermediary matrix and vectors*/
383 Gamma = xNew<IssmPDouble>(n_obs*n_obs);
384 gamma0 = xNew<IssmPDouble>(n_obs);
385 ones = xNew<IssmPDouble>(n_obs);
386
387 /*First: Create semivariogram matrix for observations*/
388 for(i=0;i<n_obs;i++){
389 for(j=0;j<=i;j++){
390 //Gamma[i*n_obs+j] = variogram->SemiVariogram(x[i]-x[j],y[i]-y[j]);
391 Gamma[i*n_obs+j] = variogram->Covariance(x[i]-x[j],y[i]-y[j]);
392 Gamma[j*n_obs+i] = Gamma[i*n_obs+j];
393 }
394 }
395 for(i=0;i<n_obs;i++) ones[i]=1;
396
397 /*Get semivariogram vector associated to this location*/
398 //for(i=0;i<n_obs;i++) gamma0[i] = variogram->SemiVariogram(x[i]-x_interp,y[i]-y_interp);
399 for(i=0;i<n_obs;i++) gamma0[i] = variogram->Covariance(x[i]-x_interp,y[i]-y_interp);
400
401 /*Solve the three linear systems*/
402#if _HAVE_GSL_
403 DenseGslSolve(&GinvG0,Gamma,gamma0,n_obs); // Gamma^-1 gamma0
404 DenseGslSolve(&Ginv1, Gamma,ones,n_obs); // Gamma^-1 ones
405 DenseGslSolve(&GinvZ, Gamma,obs,n_obs); // Gamma^-1 Z
406#else
407 _error_("GSL is required");
408#endif
409
410 /*Prepare predictor*/
411 numerator=-1.; denominator=0.;
412 for(i=0;i<n_obs;i++) numerator +=GinvG0[i];
413 for(i=0;i<n_obs;i++) denominator+=Ginv1[i];
414 ratio=numerator/denominator;
415
416 prediction = 0.;
417 error = - numerator*numerator/denominator;
418 for(i=0;i<n_obs;i++) prediction += (gamma0[i]-ratio)*GinvZ[i];
419 for(i=0;i<n_obs;i++) error += gamma0[i]*GinvG0[i];
420
421 /*clean-up*/
422 *pprediction = prediction;
423 *perror = error;
424 xDelete<IssmPDouble>(x);
425 xDelete<IssmPDouble>(y);
426 xDelete<IssmPDouble>(obs);
427 xDelete<IssmPDouble>(Gamma);
428 xDelete<IssmPDouble>(gamma0);
429 xDelete<IssmPDouble>(ones);
430 xDelete<IssmPDouble>(GinvG0);
431 xDelete<IssmPDouble>(Ginv1);
432 xDelete<IssmPDouble>(GinvZ);
433
434}/*}}}*/
435/*FUNCTION Observations::InterpolationNearestNeighbor{{{*/
436void Observations::InterpolationNearestNeighbor(IssmPDouble *pprediction,IssmPDouble x_interp,IssmPDouble y_interp,IssmPDouble radius){
437
438 /*Intermediaries*/
439 IssmPDouble x,y,obs;
440
441 /*Get clostest observation*/
442 this->ClosestObservation(&x,&y,&obs,x_interp,y_interp,radius);
443
444 /*Assign output pointer*/
445 *pprediction = obs;
446}/*}}}*/
447/*FUNCTION Observations::InterpolationV4{{{*/
448void Observations::InterpolationV4(IssmPDouble *pprediction,IssmPDouble x_interp,IssmPDouble y_interp,IssmPDouble radius,int mindata,int maxdata){
449 /* Reference: David T. Sandwell, Biharmonic spline interpolation of GEOS-3
450 * and SEASAT altimeter data, Geophysical Research Letters, 2, 139-142,
451 * 1987. Describes interpolation using value or gradient of value in any
452 * dimension.*/
453
454 /*Intermediaries*/
455 int i,j,n_obs;
456 IssmPDouble prediction,h;
457 IssmPDouble *x = NULL;
458 IssmPDouble *y = NULL;
459 IssmPDouble *obs = NULL;
460 IssmPDouble *Green = NULL;
461 IssmPDouble *weights = NULL;
462 IssmPDouble *g = NULL;
463
464 /*Some checks*/
465 _assert_(maxdata>0);
466 _assert_(pprediction);
467
468 /*If radius is not provided or is 0, return all observations*/
469 if(radius==0) radius=this->quadtree->root->length;
470
471 /*Get list of observations for current point*/
472 this->ObservationList(&x,&y,&obs,&n_obs,x_interp,y_interp,radius,maxdata);
473
474 /*If we have less observations than mindata, return UNDEF*/
475 if(n_obs<mindata || n_obs<2){
476 prediction = UNDEF;
477 }
478 else{
479
480 /*Allocate intermediary matrix and vectors*/
481 Green = xNew<IssmPDouble>(n_obs*n_obs);
482 g = xNew<IssmPDouble>(n_obs);
483
484 /*First: distance vector*/
485 for(i=0;i<n_obs;i++){
486 h = sqrt( (x[i]-x_interp)*(x[i]-x_interp) + (y[i]-y_interp)*(y[i]-y_interp) );
487 if(h>0){
488 g[i] = h*h*(log(h)-1.);
489 }
490 else{
491 g[i] = 0.;
492 }
493 }
494
495 /*Build Green function matrix*/
496 for(i=0;i<n_obs;i++){
497 for(j=0;j<=i;j++){
498 h = sqrt( (x[i]-x[j])*(x[i]-x[j]) + (y[i]-y[j])*(y[i]-y[j]) );
499 if(h>0){
500 Green[j*n_obs+i] = h*h*(log(h)-1.);
501 }
502 else{
503 Green[j*n_obs+i] = 0.;
504 }
505 Green[i*n_obs+j] = Green[j*n_obs+i];
506 }
507 }
508
509 /*Compute weights*/
510#if _HAVE_GSL_
511 DenseGslSolve(&weights,Green,obs,n_obs); // Green^-1 obs
512#else
513 _error_("GSL is required");
514#endif
515
516 /*Interpolate*/
517 prediction = 0;
518 for(i=0;i<n_obs;i++) prediction += weights[i]*g[i];
519
520 }
521
522 /*clean-up*/
523 *pprediction = prediction;
524 xDelete<IssmPDouble>(x);
525 xDelete<IssmPDouble>(y);
526 xDelete<IssmPDouble>(obs);
527 xDelete<IssmPDouble>(Green);
528 xDelete<IssmPDouble>(g);
529 xDelete<IssmPDouble>(weights);
530}/*}}}*/
531/*FUNCTION Observations::QuadtreeColoring{{{*/
532void Observations::QuadtreeColoring(IssmPDouble* A,IssmPDouble *x,IssmPDouble *y,int n){
533
534 int xi,yi,level;
535
536 for(int i=0;i<n;i++){
537 this->quadtree->IntergerCoordinates(&xi,&yi,x[i],y[i]);
538 this->quadtree->QuadtreeDepth(&level,xi,yi);
539 A[i]=(IssmPDouble)level;
540 }
541
542}/*}}}*/
543/*FUNCTION Observations::Variomap{{{*/
544void Observations::Variomap(IssmPDouble* gamma,IssmPDouble *x,int n){
545
546 /*Output and Intermediaries*/
547 int i,j,k;
548 IssmPDouble distance;
549 Observation *observation1 = NULL;
550 Observation *observation2 = NULL;
551
552 IssmPDouble *counter = xNew<IssmPDouble>(n);
553 for(j=0;j<n;j++) counter[j] = 0.0;
554 for(j=0;j<n;j++) gamma[j] = 0.0;
555
556 for(i=0;i<this->Size();i++){
557 observation1=dynamic_cast<Observation*>(this->GetObjectByOffset(i));
558
559 for(j=i+1;j<this->Size();j++){
560 observation2=dynamic_cast<Observation*>(this->GetObjectByOffset(j));
561
562 distance=sqrt(pow(observation1->x - observation2->x,2) + pow(observation1->y - observation2->y,2));
563 if(distance>x[n-1]) continue;
564
565 int index = int(distance/(x[1]-x[0]));
566 if(index>n-1) index = n-1;
567 if(index<0) index = 0;
568
569 gamma[index] += 1./2.*pow(observation1->value - observation2->value,2);
570 counter[index] += 1.;
571 }
572 }
573
574 /*Normalize semivariogram*/
575 gamma[0]=0.;
576 for(k=0;k<n;k++){
577 if(counter[k]) gamma[k] = gamma[k]/counter[k];
578 }
579
580 /*Assign output pointer*/
581 xDelete<IssmPDouble>(counter);
582}/*}}}*/
Note: See TracBrowser for help on using the repository browser.