Problems with gstat gaussian variogram and cross validation when fix nugget at 0

classic Classic list List threaded Threaded
2 messages Options
Reply | Threaded
Open this post in threaded view
|

Problems with gstat gaussian variogram and cross validation when fix nugget at 0

Stefano Menichetti
Dear all, I have a question about fixing nugget at 0 in a Gaussian
variogram using *gstat *package

I have encountered a serious problem on a variogram of piezometric data
that give "strange" results if I want to fix nugget to 0.

This is the file and how to import:

/library(gstat)/
//file =>
https://drive.google.com/open?id=1AOilEhR8G8YANVuFIBUcmEg-6UG1Ajrr//

/lay1 =� read.delim(paste(path,*file*,sep=""),sep="\t",dec=",")//
//SpData = na.exclude(subset(lay1, select = c("X","Y","Observed"))) //
//colnames(SpData)=c("X","Y","Z")/

The initial variogram:

/coordinates(SpData) = ~X+Y//
//var = variogram(Z~1,SpData)/

show a slow initial rise and so suggest a *Gaussian *type.

After a first variogram with:

/m = vgm(200,"Gau",3000,0)//
//
//plot(var, model=m, //
//main = paste(m$model[2],"//
//Nug=",m$psill[1],", Sill=",round(m$psill[2]),",
Range=",round(m$range[2]),sep ="")//
//)/

I have obtained a fitted variogram like this,

/m = fit.variogram(var, m)/

with a reasonable cross-validation results:

/krg_cv = krige.cv(Z~1, SpData, model = m)//
//plot(krg_cv$observed,krg_cv$var1.pred, xlab="Z observed", ylab="Z
predicted")//
//abline(a=0,b=1)//
//paste("media residui =",mean(krg_cv$residual))//
//paste("RMS residui =",mean(krg_cv$residual^2)^0.5)/

Resulting map it is very smooth, despite the small, it seems, nugget and
so *I have tried to fix nugget at 0* (exact interpolation).

But cross validation results with

/m = vgm(216,"Gau",3640,*0*)/

are dramatically wrong and also very different from cv results with

/m = vgm(216,"Gau",3640,*1e-3*)/

What happens ? I wrong in something ? Why nugget paraneter at 0 or at
small value influnce so hard the results ?

Thank you very much in advance,


Stefano


Dott. Geol. Stefano Menichetti
=============================================================
ARPAT - Agenzia per la Protezione dell'Ambiente della Toscana
Settore Tecnico SIRA - Sistema Informativo Regionale Ambientale via Nicola Porpora 22 - 50144 Firenze
tel. 0553206333 cell. 3668217978 (3383550147) fax 0553206410
[hidden email]  http:\\www.arpat.toscana.it  http:\\sira.arpat.toscana.it


        [[alternative HTML version deleted]]


_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: Problems with gstat gaussian variogram and cross validation when fix nugget at 0

Jon.SKOIEN
Stefano,

Fixing the nugget at zero for the Gaussian varigoram can often cause numerical problems in kriging. If you have two observations that are close in space, their modelled semivariance will be so small that you get numerical issues with the covariance matrix. In some cases it will be singular, in your case it could be inverted, but the absolute values of the weights will be very high for some of the cross-validation locations. It seems you have two stations that are very close to each other, if you remove one of these, you can see how the results will be more similar (still not the same) to the model with a small nugget effect.

Best,
Jon

--
Jon Olav Skøien

European Commission
Joint Research Centre – JRC.E.1
Disaster Risk Management Unit
Building 26b 1/144 | Via E.Fermi 2749, I-21027 Ispra (VA) Italy, TP 267
[hidden email]
Tel:  +39 0332 789205
Disclaimer: Views expressed in this email are those of the individual and do not necessarily represent official views of the European Commission.


________________________________________
From: R-sig-Geo [[hidden email]] on behalf of Stefano Menichetti [[hidden email]]
Sent: 12 October 2018 16:36
To: [hidden email]
Subject: [R-sig-Geo] Problems with gstat gaussian variogram and cross validation when fix nugget at 0

Dear all, I have a question about fixing nugget at 0 in a Gaussian
variogram using *gstat *package

I have encountered a serious problem on a variogram of piezometric data
that give "strange" results if I want to fix nugget to 0.

This is the file and how to import:

/library(gstat)/
//file =>
https://drive.google.com/open?id=1AOilEhR8G8YANVuFIBUcmEg-6UG1Ajrr//

/lay1 =� read.delim(paste(path,*file*,sep=""),sep="\t",dec=",")//
//SpData = na.exclude(subset(lay1, select = c("X","Y","Observed"))) //
//colnames(SpData)=c("X","Y","Z")/

The initial variogram:

/coordinates(SpData) = ~X+Y//
//var = variogram(Z~1,SpData)/

show a slow initial rise and so suggest a *Gaussian *type.

After a first variogram with:

/m = vgm(200,"Gau",3000,0)//
//
//plot(var, model=m, //
//main = paste(m$model[2],"//
//Nug=",m$psill[1],", Sill=",round(m$psill[2]),",
Range=",round(m$range[2]),sep ="")//
//)/

I have obtained a fitted variogram like this,

/m = fit.variogram(var, m)/

with a reasonable cross-validation results:

/krg_cv = krige.cv(Z~1, SpData, model = m)//
//plot(krg_cv$observed,krg_cv$var1.pred, xlab="Z observed", ylab="Z
predicted")//
//abline(a=0,b=1)//
//paste("media residui =",mean(krg_cv$residual))//
//paste("RMS residui =",mean(krg_cv$residual^2)^0.5)/

Resulting map it is very smooth, despite the small, it seems, nugget and
so *I have tried to fix nugget at 0* (exact interpolation).

But cross validation results with

/m = vgm(216,"Gau",3640,*0*)/

are dramatically wrong and also very different from cv results with

/m = vgm(216,"Gau",3640,*1e-3*)/

What happens ? I wrong in something ? Why nugget paraneter at 0 or at
small value influnce so hard the results ?

Thank you very much in advance,


Stefano


Dott. Geol. Stefano Menichetti
=============================================================
ARPAT - Agenzia per la Protezione dell'Ambiente della Toscana
Settore Tecnico SIRA - Sistema Informativo Regionale Ambientale via Nicola Porpora 22 - 50144 Firenze
tel. 0553206333 cell. 3668217978 (3383550147) fax 0553206410
[hidden email]  http:\\www.arpat.toscana.it  http:\\sira.arpat.toscana.it


        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo