Télécharger krig_01.dgibi

Retour à la liste

** Test du krigeage en dimension 2
 
* Sur la base sur la reference suivante :
* Edward H. Isaaks, R. Mohan Srivasta
* Applied Geostatistics, Oxford University Press, 1989 - 561 pages
*
* L'exemple traite ici est situe au chapitre 12, pages 290-296
* Consultable sur : https://www.geokniga.org/bookfiles/geokniga-anintroductiontoappliedgeostatistics.pdf
 
 
** Options generales
OPTI 'ECHO' 0 ;
itrac = FAUX ;
SAUT 2 'LIGNE' ;
 
 
 
************************************************************************
*        Donnees communes : points de mesures et points cibles         *
************************************************************************
 
** Valeurs des points de mesure
lxm   = PROG  61.  63.  64.  68.  71.  73.  75. ;
lym   = PROG 139. 140. 129. 128. 140. 141. 128. ;
lfm   = PROG 477. 696. 227. 646. 606. 791. 783. ;
 
** Variogramme (modele gaussien avec pepite, palier et portee)
c0    = 0. ;
c1    = 10. ;
a     = 10. ;
lh    = PROG 0. 'PAS' 0.01 20. ;
lg    = c0 + (c1 * (1. - (EXP (-3. * lh / a)))) ;
REMP lg 1 0. ;
gam   = EVOL 'VERT' 'MANU' lh lg ;
SI itrac ;
  DESS gam 'TITR' 'Variogramme (modele gaussien)' ;
FINSI ;
 
** Valeurs des points cibles (ou l'on va kriger)
xmin  = 60. ;
xmax  = 80. ;
ymin  = 120. ;
ymax  = 145. ;
dx    = 0.25 ;
lxc0  = PROG xmin 'PAS' dx xmax ;
lyc0  = PROG ymin 'PAS' dx ymax ;
 
** Valeurs au point de reference (pour tester le calcul)
xref  = 65. ;
yref  = 137. ;
fref  = 592.7 ;
vref  = 8.96 ;
 
 
 
************************************************************************
*      Cas 1 : krigeage avec mesures LISTREEL et cibles LISTREEL       *
*                            et variogramme                            *
************************************************************************
 
** Liste des coordonnes des cibles
lxc   = PROG ;
lyc   = PROG ;
REPE by (DIME lyc0) ;
  REPE bx (DIME lxc0) ;
    lxc   = lxc ET (EXTR lxc0 &bx) ;
    lyc   = lyc ET (EXTR lyc0 &by) ;
  FIN bx ;
FIN by ;
ncib  = DIME lxc ;
 
** Krigeage au points cibles (avec le variogramme)
t1                   = TABL ;
t1 . 'MESURES'       = TABL ;
t1 . 'MESURES' . 'X' = lxm ;
t1 . 'MESURES' . 'Y' = lym ;
t1 . 'MESURES' . 'FONC' = lfm ;
t1 . 'CIBLES'        = TABL ;
t1 . 'CIBLES'  . 'X' = lxc ;
t1 . 'CIBLES'  . 'Y' = lyc ;
t1 . 'COORDONNEES'   = MOTS 'X' 'Y' ;
t1 . 'COMPOSANTE'    = MOT 'FONC' ;
t1 . 'VARIOGRAMME'   = gam ;
TEMP 'ZERO' ;
est1 var1 = KRIG t1;
tps1  = TEMP 'HORL' ;
MESS 'CAS #1' ;
MESS '======' ;
MESS ;
MESS 'Krigeage sur' ' ' ncib ' points realise en' ' ' tps1 ' ms' ;
MESS ;
 
** Comparaison a la solution de reference
ixref = POSI xref 'DANS' lxc0 ;
iyref = POSI yref 'DANS' lyc0 ;
iref  = ixref + ((iyref - 1) * (DIME lxc0)) ;
xcal1 = EXTR lxc iref ;
ycal1 = EXTR lyc iref ;
fcal1 = EXTR est1 iref ;
vcal1 = EXTR var1 iref ;
errf1 = ((FLOT (@ARR fcal1 1)) - fref) / fref ;
errv1 = ((FLOT (@ARR vcal1 2)) - vref) / vref ;
MESS '+--------------------------+--------------------------+' ;
MESS '  Valeurs de reference     | Valeurs calculees         ' ;
MESS '+--------------------------+--------------------------+' ;
MESS '   x            y          |  x            y           ' ;
MESS ' ' xref ' ' yref         ' |' xcal1 ' ' ycal1        ' ' ;
MESS ' --------------------------+-------------------------- ' ;
MESS '  Fonction     Variance    | Fonction     Variance     ' ;
MESS ' ' fref ' ' vref         ' |' fcal1 ' ' vcal1        ' ' ;
MESS ' --------------------------+-------------------------- ' ;
MESS '                           | Erreur       Erreur       ' ;
MESS '                           |' errf1 ' ' errv1        ' ' ;
MESS '+--------------------------+--------------------------+' ;
SAUT 2 'LIGNE' ;
 
 
 
************************************************************************
*      Cas 2 : krigeage avec mesures LISTREEL et cibles CHPOINT        *
*                           et covariogramme                           *
************************************************************************
 
** Covariogramme
lc    = c1 * (EXP (-3. * lh / a)) ;
REMP lc 1 (c0 + c1) ;
cov   = EVOL 'VERT' 'MANU' lh lc ;
 
** Maillage de points cibles, ou l'on veut interpoler
OPTI 'DIME' 2 'ELEM' 'QUA4' 'DENS' dx ;
lig1  = DROI (xmin ymin) (xmax ymin) ;
mail  = lig1 TRAN (0. (ymax - ymin)) ;
ncib  = NBNO mail ;
con   = CONT mail ;
x y   = COOR mail ;
cib   = (NOMC 'X' x) ET (NOMC 'Y' y) ;
 
** Krigeage au points cibles (avec le covariogramme)
t2                   = TABL ;
t2 . 'MESURES'       = TABL ;
t2 . 'MESURES' . 'X' = lxm ;
t2 . 'MESURES' . 'Y' = lym ;
t2 . 'MESURES' . 'FONC' = lfm ;
t2 . 'CIBLES'        = cib ;
t2 . 'COORDONNEES'   = MOTS 'X' 'Y' ;
t2 . 'COMPOSANTE'    = MOT 'FONC' ;
t2 . 'COVARIOGRAMME' = cov ;
TEMP 'ZERO' ;
est2 var2 = KRIG t2 ;
tps2  = TEMP 'HORL' ;
MESS 'CAS #2' ;
MESS '======' ;
MESS ;
MESS 'Krigeage sur' ' ' ncib ' points realise en' ' ' tps2 ' ms' ;
MESS ;
 
** Comparaison a la solution de reference
pref  = mail POIN 'PROCHE' (xref yref) ;
xcal2 ycal2 = COOR pref ;
fcal2 = EXTR est2 'FONC' pref ;
vcal2 = EXTR var2 'FONC' pref ;
errf2 = ((FLOT (@ARR fcal2 1)) - fref) / fref ;
errv2 = ((FLOT (@ARR vcal2 2)) - vref) / vref ;
MESS '+--------------------------+--------------------------+' ;
MESS '  Valeurs de reference     | Valeurs calculees         ' ;
MESS '+--------------------------+--------------------------+' ;
MESS '   x            y          |  x            y           ' ;
MESS ' ' xref ' ' yref         ' |' xcal2 ' ' ycal2        ' ' ;
MESS ' --------------------------+-------------------------- ' ;
MESS '  Fonction     Variance    | Fonction     Variance     ' ;
MESS ' ' fref ' ' vref         ' |' fcal2 ' ' vcal2        ' ' ;
MESS ' --------------------------+-------------------------- ' ;
MESS '                           | Erreur       Erreur       ' ;
MESS '                           |' errf2 ' ' errv2        ' ' ;
MESS '+--------------------------+--------------------------+' ;
SAUT 2 'LIGNE' ;
 
 
 
************************************************************************
*       Cas 3 : krigeage avec mesures CHPOINT et cibles LISTREEL       *
*                            et variogramme                            *
************************************************************************
 
** Maillage et champs pour visualiser ces points de mesure
ptm   = POIN lxm lym ;
xm ym = COOR ptm ;
mes   = (NOMC 'X' xm) ET (NOMC 'Y' ym) ET
        (MANU 'CHPO' ptm 1 'FONC' lfm 'NATURE' 'DIFFUS') ;
 
** Krigeage au points cibles (avec le variogramme)
t3                   = TABL ;
t3 . 'MESURES'       = mes ;
t3 . 'CIBLES'        = TABL ;
t3 . 'CIBLES'  . 'X' = lxc ;
t3 . 'CIBLES'  . 'Y' = lyc ;
t3 . 'COORDONNEES'   = MOTS 'X' 'Y' ;
t3 . 'COMPOSANTE'    = MOT 'FONC' ;
t3 . 'VARIOGRAMME'   = gam ;
TEMP 'ZERO' ;
est3 var3 = KRIG t3 ;
tps3  = TEMP 'HORL' ;
MESS 'CAS #3' ;
MESS '======' ;
MESS ;
MESS 'Krigeage sur' ' ' ncib ' points realise en' ' ' tps3 ' ms' ;
MESS ;
 
** Comparaison a la solution de reference
xcal3 = xcal1 ;
ycal3 = xcal1 ;
fcal3 = EXTR est3 iref ;
vcal3 = EXTR var3 iref ;
errf3 = ((FLOT (@ARR fcal3 1)) - fref) / fref ;
errv3 = ((FLOT (@ARR vcal3 2)) - vref) / vref ;
MESS '+--------------------------+--------------------------+' ;
MESS '  Valeurs de reference     | Valeurs calculees         ' ;
MESS '+--------------------------+--------------------------+' ;
MESS '   x            y          |  x            y           ' ;
MESS ' ' xref ' ' yref         ' |' xcal3 ' ' ycal3        ' ' ;
MESS ' --------------------------+-------------------------- ' ;
MESS '  Fonction     Variance    | Fonction     Variance     ' ;
MESS ' ' fref ' ' vref         ' |' fcal3 ' ' vcal3        ' ' ;
MESS ' --------------------------+-------------------------- ' ;
MESS '                           | Erreur       Erreur       ' ;
MESS '                           |' errf3 ' ' errv3        ' ' ;
MESS '+--------------------------+--------------------------+' ;
SAUT 2 'LIGNE' ;
 
 
 
************************************************************************
*       Cas 4 : krigeage avec mesures CHPOINT et cibles CHPOINT        *
*                           et covariogramme                           *
************************************************************************
 
** Krigeage au points cibles (avec le covariogramme)
t4                   = TABL ;
t4 . 'MESURES'       = mes ;
t4 . 'CIBLES'        = cib ;
t4 . 'COORDONNEES'   = MOTS 'X' 'Y' ;
t4 . 'COMPOSANTE'    = MOT 'FONC' ;
t4 . 'COVARIOGRAMME' = cov ;
TEMP 'ZERO' ;
est4 var4 = KRIG t4 ;
tps4  = TEMP 'HORL' ;
MESS 'CAS #4' ;
MESS '======' ;
MESS ;
MESS 'Krigeage sur' ' ' ncib ' points realise en' ' ' tps4 ' ms' ;
MESS ;
 
** Comparaison a la solution de reference
xcal4 = xcal2 ;
ycal4 = ycal2 ;
fcal4 = EXTR est4 'FONC' pref ;
vcal4 = EXTR var4 'FONC' pref ;
errf4 = ((FLOT (@ARR fcal4 1)) - fref) / fref ;
errv4 = ((FLOT (@ARR vcal4 2)) - vref) / vref ;
MESS '+--------------------------+--------------------------+' ;
MESS '  Valeurs de reference     | Valeurs calculees         ' ;
MESS '+--------------------------+--------------------------+' ;
MESS '   x            y          |  x            y           ' ;
MESS ' ' xref ' ' yref         ' |' xcal4 ' ' ycal4        ' ' ;
MESS ' --------------------------+-------------------------- ' ;
MESS '  Fonction     Variance    | Fonction     Variance     ' ;
MESS ' ' fref ' ' vref         ' |' fcal4 ' ' vcal4        ' ' ;
MESS ' --------------------------+-------------------------- ' ;
MESS '                           | Erreur       Erreur       ' ;
MESS '                           |' errf4 ' ' errv4        ' ' ;
MESS '+--------------------------+--------------------------+' ;
SAUT 2 'LIGNE' ;
 
** Trace des champs
SI itrac ;
  TRAC est4 mail con 'TITR' 'Interpolation de la fonction par Krigeage' ;
  TRAC var4 mail con 'TITR' 'Variance de l''estimation par Krigeage' ;
FINSI ;
 
 
** Sortie en erreur si ecart trop grand
errmax = MAXI 'ABS' errf1 errf2 errf3 errf4
                    errv1 errv2 errv3 errv4 ;
MESS 'Erreur max. =' errmax ;
SI (errmax > 1.E-10) ;
  ERRE '--> Mauvaise interpolation :-(' ;
FINSI ;
 
 
FIN ;
 
 
 

© Cast3M 2003 - Tous droits réservés.
Mentions légales