Changeset 12419

Timestamp:

06/15/12 11:41:06 (13 years ago)

Author:

Mathieu Morlighem

Message:

Moved Kriging to method of Observations

Location:

issm/trunk-jpl/src/c/Container

Files:

• Elements.h (modified) (1 diff)
• Observations.cpp (modified) (4 diffs)
• Observations.h (modified) (2 diffs)

issm/trunk-jpl/src/c/Container/Elements.h

@@ -15 +15 @@
 class DataSet;
 class Inputs;
-
 
 class Elements: public DataSet{

issm/trunk-jpl/src/c/Container/Observations.cpp

@@ -20 +20 @@
 #include "../shared/shared.h"
 #include "../include/include.h"
+#include "../modules/modules.h"
 #include "../EnumDefinitions/EnumDefinitions.h"
 #include "../io/io.h"
@@ -115 +116 @@
 
 /*Methods*/
-/*FUNCTION Observations::ObservationList{{{*/
+/*FUNCTION Observations::ClosestObservation{{{*/
+void Observations::ClosestObservation(double *px,double *py,double *pobs,double x_interp,double y_interp,double radius){
+
+        /*Output and Intermediaries*/
+        bool         stop;
+        int          nobs,i,index;
+        double       h2,hmin2,radius2;
+        int         *indices      = NULL;
+        Observation *observation  = NULL;
+
+        /*If radius is not provided or is 0, return all observations*/
+        if(radius==0) radius=this->quadtree->root->length;
+
+        /*Compute radius square*/
+        radius2 = radius*radius;
+
+        /*Find all observations that are in radius*/
+        this->quadtree->RangeSearch(&indices,&nobs,x_interp,y_interp,radius);
+        for (i=0;i<nobs;i++){
+                observation=(Observation*)this->GetObjectByOffset(indices[i]);
+                h2 = (observation->x-x_interp)*(observation->x-x_interp) + (observation->y-y_interp)*(observation->y-y_interp);
+
+                if(i==0){
+                        hmin2 = h2;
+                        index = i;
+                }
+                else{
+                        if(h2<hmin2){
+                                hmin2 = h2;
+                                index = i;
+                        }
+                }
+        }  
+
+        /*Assign output pointer*/
+        if(!nobs){
+                *px=UNDEF;
+                *py=UNDEF;
+                *pobs=UNDEF;
+        }
+        else{
+                observation=(Observation*)this->GetObjectByOffset(indices[index]);
+                *px=observation->x;
+                *py=observation->y;
+                *pobs=observation->value;
+        }
+        xfree((void**)&indices);
+
+}/*}}}*/
+/*FUNCTION Observations::ObservationList(double **px,double **py,double **pobs,int* pnobs,double x_interp,double y_interp,double radius,int maxdata){{{*/
 void Observations::ObservationList(double **px,double **py,double **pobs,int* pnobs,double x_interp,double y_interp,double radius,int maxdata){
 
@@ -130 +180 @@
         Observation *observation  = NULL;
 
+        /*If radius is not provided or is 0, return all observations*/
+        if(radius==0) radius=this->quadtree->root->length;
+
         /*Compute radius square*/
-        if(radius==0) radius=this->quadtree->root->length;
         radius2 = radius*radius;
 
@@ -196 +248 @@
         *pnobs=nobs;
 }/*}}}*/
+/*FUNCTION Observations::ObservationList(double **px,double **py,double **pobs,int* pnobs){{{*/
+void Observations::ObservationList(double **px,double **py,double **pobs,int* pnobs){
+
+        /*Output and Intermediaries*/
+        int          nobs;
+        double      *x            = NULL;
+        double      *y            = NULL;
+        double      *obs          = NULL;
+        Observation *observation  = NULL;
+
+        nobs = this->Size();
+
+        if(nobs){
+                x   = (double*)xmalloc(nobs*sizeof(double));
+                y   = (double*)xmalloc(nobs*sizeof(double));
+                obs = (double*)xmalloc(nobs*sizeof(double));
+                for(int i=0;i<this->Size();i++){
+                        observation=(Observation*)this->GetObjectByOffset(i);
+                        observation->WriteXYObs(&x[i],&y[i],&obs[i]);
+                }
+        }
+
+        /*Assign output pointer*/
+        *px=x;
+        *py=y;
+        *pobs=obs;
+        *pnobs=nobs;
+}/*}}}*/
+/*FUNCTION Observations::InterpolationKriging{{{*/
+void Observations::InterpolationKriging(double *pprediction,double *perror,double x_interp,double y_interp,double radius,int mindata,int maxdata,Variogram* variogram){
+
+        /*Intermediaries*/
+        int           i,j,n_obs;
+        double        prediction,error;
+        double        numerator,denominator,ratio;
+        double       *x            = NULL;
+        double       *y            = NULL;
+        double       *obs          = NULL;
+        double       *Gamma        = NULL;
+        double       *GinvG0       = NULL;
+        double       *Ginv1        = NULL;
+        double       *GinvZ        = NULL;
+        double       *gamma0       = NULL;
+        double       *ones         = NULL;
+
+        /*Some checks*/
+        _assert_(mindata>0 && maxdata>0);
+        _assert_(pprediction && perror);
+
+        /*If radius is not provided or is 0, return all observations*/
+        if(radius==0) radius=this->quadtree->root->length;
+
+        /*Get list of observations for current point*/
+        this->ObservationList(&x,&y,&obs,&n_obs,x_interp,y_interp,radius,maxdata);
+
+        /*If we have less observations than mindata, return UNDEF*/
+        if(n_obs<mindata){
+                *pprediction = -999.0; 
+                *perror      = -999.0; 
+                return;
+        }
+
+        /*Allocate intermediary matrix and vectors*/
+        Gamma  = (double*)xmalloc(n_obs*n_obs*sizeof(double));
+        gamma0 = (double*)xmalloc(n_obs*sizeof(double));
+        ones   = (double*)xmalloc(n_obs*sizeof(double));
+
+        /*First: Create semivariogram matrix for observations*/
+        for(i=0;i<n_obs;i++){
+                for(j=0;j<=i;j++){
+                        //Gamma[i*n_obs+j] = variogram->SemiVariogram(x[i]-x[j],y[i]-y[j]);
+                        Gamma[i*n_obs+j] = variogram->Covariance(x[i]-x[j],y[i]-y[j]);
+                        Gamma[j*n_obs+i] = Gamma[i*n_obs+j];
+                }
+        }
+        for(i=0;i<n_obs;i++) ones[i]=1;
+
+        /*Get semivariogram vector associated to this location*/
+        //for(i=0;i<n_obs;i++) gamma0[i] = variogram->SemiVariogram(x[i]-x_interp,y[i]-y_interp);
+        for(i=0;i<n_obs;i++) gamma0[i] = variogram->Covariance(x[i]-x_interp,y[i]-y_interp);
+
+        /*Solve the three linear systems*/
+#if _HAVE_GSL_
+        SolverxGsl(&GinvG0,Gamma,gamma0,n_obs); // Gamma^-1 gamma0
+        SolverxGsl(&Ginv1, Gamma,ones,n_obs);   // Gamma^-1 ones
+        SolverxGsl(&GinvZ, Gamma,obs,n_obs);    // Gamma^-1 Z
+#else
+        _error_("GSL is required");
+#endif
+
+        /*Prepare predictor*/
+        numerator=-1.; denominator=0.;
+        for(i=0;i<n_obs;i++) numerator  +=GinvG0[i];
+        for(i=0;i<n_obs;i++) denominator+=Ginv1[i];
+        ratio=numerator/denominator;
+
+        prediction = 0.;
+        error      = - numerator*numerator/denominator;
+        for(i=0;i<n_obs;i++) prediction += (gamma0[i]-ratio)*GinvZ[i];
+        for(i=0;i<n_obs;i++) error += gamma0[i]*GinvG0[i];
+
+        /*clean-up*/
+        *pprediction = prediction;
+        *perror = error;
+        xfree((void**)&x);
+        xfree((void**)&y);
+        xfree((void**)&obs);
+        xfree((void**)&Gamma);
+        xfree((void**)&gamma0);
+        xfree((void**)&ones);
+        xfree((void**)&GinvG0);
+        xfree((void**)&Ginv1);
+        xfree((void**)&GinvZ);
+
+}/*}}}*/
 /*FUNCTION Observations::QuadtreeColoring{{{*/
 void Observations::QuadtreeColoring(double* A,double *x,double *y,int n){

issm/trunk-jpl/src/c/Container/Observations.h

@@ -8 +8 @@
 class Obsevration;
 class Quadtree;
+class Variogram;
 class Options;
 
@@ -23 +24 @@
 
                 /*Methods*/
+                void ClosestObservation(double *px,double *py,double *pobs,double x_interp,double y_interp,double radius);
+                void InterpolationKriging(double *pprediction,double *perror,double x_interp,double y_interp,double radius,int mindata,int maxdata,Variogram* variogram);
+                void ObservationList(double **px,double **py,double **pobs,int* pnobs);
                 void ObservationList(double **px,double **py,double **pobs,int* pnobs,double x_interp,double y_interp,double radius,int maxdata);
                 void QuadtreeColoring(double* A,double *x,double *y,int n);