Changeset 12298


Ignore:
Timestamp:
05/29/12 11:25:53 (13 years ago)
Author:
Mathieu Morlighem
Message:

Added many options to Kriging now equivalent to gslib

Location:
issm/trunk-jpl/src
Files:
4 edited

Legend:

Unmodified
Added
Removed

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

    r12292 r12298  
    3636
    3737        /*Intermediaries*/
    38         int          i,maxdepth,level,counter;
     38        int          i,j,maxdepth,level,counter,index;
    3939        int          xi,yi;
    4040        double       xmin,xmax,ymin,ymax;
    41         double       offset,minlength;
     41        double       offset,minlength,minspacing,mintrimming,maxtrimming;
     42        int         *indices     = NULL;
    4243        Observation *observation = NULL;
    4344
    … …  
    5253        offset=0.05*(ymax-ymin); ymin-=offset; ymax+=offset;
    5354
     55        /*Get trimming limits*/
     56        options->Get(&mintrimming,"mintrimming",-1.e+21);
     57        options->Get(&maxtrimming,"maxtrimming",+1.e+21);
     58        options->Get(&minspacing,"minspacing",0.01);
     59        if(minspacing<=0) _error_("minspacing must > 0");
     60
    5461        /*Get Minimum box size*/
    5562        if(options->GetOption("boxlength")){
    … …  
    7077        counter = 0;
    7178        for(i=0;i<n;i++){
     79
     80                /*First check limits*/
     81                if(observations_list[i]>maxtrimming) continue;
     82                if(observations_list[i]<mintrimming) continue;
     83
     84                /*First check that this observation is not too close from another one*/
     85                this->quadtree->ClosestObs(&index,x[i],y[i]);
     86                if(index>=0){
     87                        observation=(Observation*)this->GetObjectByOffset(index);
     88                        if(pow(observation->x-x[i],2)+pow(observation->y-y[i],2) < minspacing) continue;
     89                }
     90
    7291                this->quadtree->IntergerCoordinates(&xi,&yi,x[i],y[i]);
    7392                this->quadtree->QuadtreeDepth2(&level,xi,yi);
    … …  
    83102        }
    84103        printf("done\n");
     104        printf("Initial number of observations: %i\n",n);
     105        printf("  Final number of observations: %i\n",this->quadtree->NbObs);
    85106}
    86107/*}}}*/
    … …  
    98119        /*Output and Intermediaries*/
    99120        bool         stop;
    100         int          nobs,i,j,k,n,counter;
     121        int          nobs,tempnobs,i,j,k,n,counter;
    101122        double       h2,radius2;
    102123        int         *indices      = NULL;
     124        int         *tempindices  = NULL;
    103125        double      *dists        = NULL;
    104126        double      *x            = NULL;
    … …  
    112134
    113135        /*Find all observations that are in radius*/
    114         indices = (int*)xmalloc(this->Size()*sizeof(int));
    115         dists   = (double*)xmalloc(this->Size()*sizeof(double));
    116         nobs    = 0;
    117 
    118         for (i=0;i<this->Size();i++){
    119                 observation=(Observation*)this->GetObjectByOffset(i);
     136        this->quadtree->RangeSearch(&tempindices,&tempnobs,x_interp,y_interp,radius);
     137        if(tempnobs){
     138                indices = (int*)xmalloc(tempnobs*sizeof(int));
     139                dists   = (double*)xmalloc(tempnobs*sizeof(double));
     140        }
     141        nobs = 0;
     142        for (i=0;i<tempnobs;i++){
     143                observation=(Observation*)this->GetObjectByOffset(tempindices[i]);
    120144                h2 = (observation->x-x_interp)*(observation->x-x_interp) + (observation->y-y_interp)*(observation->y-y_interp);
    121145
    122146                if(nobs==maxdata && h2>radius2) continue;
    123147                if(nobs<=maxdata){
    124                         indices[nobs]   = i;
     148                        indices[nobs]   = tempindices[i];
    125149                        dists[nobs]     = h2;
    126150                        nobs++;
    … …  
    141165                                }
    142166                                dists[k]   = h2;
    143                                 indices[k] = i;
     167                                indices[k] = tempindices[i];
    144168                                stop = true;
    145169                                break;
    … …  
    149173        } 
    150174        xfree((void**)&dists);
     175        xfree((void**)&tempindices);
    151176
    152177        if(nobs){

  • issm/trunk-jpl/src/c/objects/Kriging/Quadtree.cpp

    r12273 r12298  
    250250        }
    251251}/*}}}*/
     252/*FUNCTION Quadtree::ClosestObs{{{1*/
     253void Quadtree::ClosestObs(int *pindex,double x,double y){
     254
     255        QuadtreeBox **pbox = NULL;
     256        QuadtreeBox  *box  = NULL;
     257        int           xi,yi;
     258        int           level,levelbin;
     259        int           index = -1;
     260        double        length,length2;
     261
     262        /*Get integer coodinates*/
     263        this->IntergerCoordinates(&xi,&yi,x,y);
     264
     265        /*Initialize levels*/
     266        level    = 0;
     267        levelbin = (1L<<this->MaxDepth);// = 2^30
     268
     269        /*Get inital box (the largest)*/
     270        pbox=&root;
     271
     272        /*Find the smallest box where this point is located*/
     273        while((box=*pbox) && (box->nbitems<0)){
     274
     275                levelbin>>=1; level+=1;
     276
     277                pbox = &box->box[IJ(xi,yi,levelbin)];
     278        }
     279
     280        /*Add obervation in this box (should be full)*/
     281        if(box && box->nbitems>0){
     282                index  = box->obs[0]->index;
     283                length = pow(box->obs[0]->x - x,2.) + pow(box->obs[0]->y - y,2.);
     284                for(int i=1;i<box->nbitems;i++){
     285                        length2 = pow(box->obs[i]->x - x,2.) + pow(box->obs[i]->y - y,2.);
     286                        if(length2<length){
     287                                index  = box->obs[i]->index;
     288                                length = length2;
     289                        }
     290                }
     291        }
     292
     293        *pindex=index;
     294}/*}}}*/
    252295/*FUNCTION Quadtree::Echo{{{1*/
    253296void  Quadtree::Echo(void){
    … …  
    453496
    454497        /*Allocate indices (maximum by default*/
    455         indices = (int*)xmalloc(this->NbObs*sizeof(int));
     498        if(this->NbObs) indices = (int*)xmalloc(this->NbObs*sizeof(int));
    456499        nobs = 0;
    457500
    458         this->root->RangeSearch(indices,&nobs,x,y,range);
     501        if(this->root) this->root->RangeSearch(indices,&nobs,x,y,range);
    459502
    460503        /*Clean-up and return*/

  • issm/trunk-jpl/src/c/objects/Kriging/Quadtree.h

    r12273 r12298  
    3434
    3535                                /*Methods*/
    36                                 int IsWithinRange(double x,double y,double range);
    37                                 void RangeSearch(int* indices,int *pnobs,double x,double y,double range);
    38                                 void WriteObservations(int* indices,int *pnobs);
     36                                int          IsWithinRange(double x,double y,double range);
     37                                void         RangeSearch(int *indices,int *pnobs,double x,double y,double range);
     38                                void         WriteObservations(int *indices,int *pnobs);
    3939
    4040                };
    … …  
    5454                void         Add(Observation *observation);
    5555                void         AddAndAverage(double x,double y,double value);
     56                void         ClosestObs(int *pindex,double x,double y);
    5657                void         DeepEcho(void);
    5758                void         Echo(void);

  • issm/trunk-jpl/src/modules/Kriging/Kriging.cpp

    r12251 r12298  
    5555void KrigingUsage(void){
    5656        _printf_(true,"\n");
    57         _printf_(true,"   usage: predictions=%s(x,y,observations,x_interp,y_interp);\n",__FUNCT__);
     57        _printf_(true,"   usage: predictions=%s(x,y,observations,x_interp,y_interp,'options');\n",__FUNCT__);
     58        _printf_(true,"   available options:\n");
     59        _printf_(true,"      -'model': Available variogram models 'gaussian' (default),'spherical','power','exponential'\n");
     60        _printf_(true,"         -'nugget': nugget effect (default 0.2)\n");
     61        _printf_(true,"         -'range':  for gaussian, spherical and exponential models (default sqrt(3))\n");
     62        _printf_(true,"         -'sill':   for gaussian, spherical and exponential models (default 1)\n");
     63        _printf_(true,"         -'slope':  for power model (default 1)\n");
     64        _printf_(true,"         -'power':  for power model (default 1)\n");
     65        _printf_(true,"      -'searchradius': search radius for each prediction (default is observations span)\n");
     66        _printf_(true,"      -'boxlength':    minimum length of quadtree boxes (useful to decrease the number of observations)\n");
     67        _printf_(true,"      -'maxdata':      minimum number of observations for a prediction (default is 50)\n");
     68        _printf_(true,"      -'mindata':      maximum number of observations for a prediction (default is 1)\n");
     69        _printf_(true,"      -'maxtrimming':  maximum trimming value (default is -1.e+21)\n");
     70        _printf_(true,"      -'mintrimming':  minimum trimming value (default is +1.e+21)\n");
     71        _printf_(true,"      -'minspacing':   minimum distance between observation (default is 0.01)\n");
    5872        _printf_(true,"\n");
    5973}
Note: See TracChangeset for help on using the changeset viewer.