regression-kriging and co-kriging

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

regression-kriging and co-kriging

Emanuele Barca
Dear Edzer,

maybe I found the solution. I found this in the predict function help:
"When a non-stationary (i.e., non-constant) mean is used, both for
simulation and prediction purposes the variogram model defined should be
that of the residual process, and not that of the raw observations"
Since my data were, actually, non-stationary, I applied the universal
co-kriging instead usual co-kriging.
now the maps of regression-kring and co-kriging are actually similar s
expected.
did I understand correctly the quoted sentence?

regards

emanuele barca
------------------------------

>
> Message: 2
> Date: Sat, 10 Aug 2019 10:41:38 +0200
> From: Edzer Pebesma <[hidden email]>
> To: [hidden email]
> Subject: Re: [R-sig-Geo] regression-kriging and co-kriging
> Message-ID: <[hidden email]>
> Content-Type: text/plain; charset="utf-8"
>
> Hard to tell from your script. Maybe give a reproducible example?
>
> On 8/6/19 1:07 PM, Emanuele Barca wrote:
>> Dear  r-sig-geo friends,
>>
>> I produced two maps garnered in the following way:
>>
>> # for regression-kriging
>> Piezo.map <-autoKrige(LivStat ~  Z, input_data = mydata.sp, new_data
>> = covariates,  model = "Ste")
>>
>> Piezork.pred <- Piezo.map$krige_output$var1.pred
>> Piezork.coords <- Piezo.map$krige_output@coords
>> Piezork.out <- as.data.frame(cbind(Piezork.coords, Piezork.pred))
>> colnames(Piezork.out)[1:2] <- c("X", "Y")
>> coordinates(Piezork.out) = ~ X + Y
>> gridded(Piezork.out) <- TRUE
>>
>> spplot(Piezork.out, main = list(label = "R-k Hydraulic head", cex =
>> 1.5))
>>
>> #for co-kriging
>> g <- gstat(id = "Piezo", formula = LivStat ~ 1, data = mydata.sp, set
>> = list(nocheck = 1))
>> g <- gstat(g, id = "Z", formula = Z ~ 1, data = mydata.sp, set =
>> list(nocheck = 1))
>>
>> v.g <- variogram(g)
>>
>> #g <- gstat(g, id = "Piezo", model = vgm(150, "Mat", 1350, 0.0, kappa
>> = 1.9), fill.all = T)#
>> g <- gstat(g, id = "Piezo", model = vgm(0.7, "Ste", 1300, 18, kappa =
>> 1.9), fill.all = T)#
>> g.fit <- fit.lmc(v.g, g, fit.lmc = TRUE, correct.diagonal = 1.01) #
>> fit multivariable variogram model , fit.lmc = TRUE, correct.diagonal =
>> 1.01
>> g.fit
>> plot(v.g, model = g.fit, main = "Fitted Variogram Models - Raw Data")#
>> #gridded(covariates) <- TRUE
>> g.cok <- predict(g.fit, newdata = covariates)#grid
>>
>> g.cok.pred <- g.cok@data$Piezo.pred
>> aaaa <- na.omit(g.cok.pred)
>> g.cok.coords <- g.cok@coords
>> g.cok.out <- as.data.frame(cbind(g.cok.coords, g.cok.pred))
>> colnames(g.cok.out)[1:2] <- c("X", "Y")
>> coordinates(g.cok.out) = ~ X + Y
>> gridded(g.cok.out) <- TRUE
>> spplot(g.cok.out, main = list(label = "Hydraulic head with
>> Co-kriging", cex = 1.5))
>>
>> ###########################################################################################################################
>>
>> I am unable to understand why the first map appears as a raster and
>> the second not, notwithstanding the fact that they are both computed
>> on the same "covariates" DEM???
>>
>> where is the mistake???
>>
>> regards
>>
>> emanuele
>>
>> ________________________________________________________
>> Emanuele Barca                               Researcher
>> Water Research Institute                       (IRSA-CNR)
>> Via De Blasio, 5                       70123 Bari (Italy)
>> Phone +39 080 5820535               Fax  +39 080 5313365
>> Mobile +39 340 3420689
>> _________________________________________________________
>>
>>
>>
>> ---
>> Questa e-mail è stata controllata per individuare virus con Avast
>> antivirus.
>> https://www.avast.com/antivirus
>>
>> [[alternative HTML version deleted]]
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> [hidden email]
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>
> --
> Edzer Pebesma
> Institute for Geoinformatics
> Heisenbergstrasse 2, 48151 Muenster, Germany
> Phone: +49 251 8333081
>
> -------------- next part --------------
> A non-text attachment was scrubbed...
> Name: pEpkey.asc
> Type: application/pgp-keys
> Size: 2472 bytes
> Desc: not available
> URL:
> <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20190810/b5cb9d77/attachment-0001.bin>
>
>
> ------------------------------
>
> Subject: Digest Footer
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
>
> ------------------------------
>
> End of R-sig-Geo Digest, Vol 192, Issue 7
> *****************************************

_______________________________________________
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: regression-kriging and co-kriging

edzer


On 8/12/19 8:21 PM, Emanuele Barca wrote:

> Dear Edzer,
>
> maybe I found the solution. I found this in the predict function help:
> "When a non-stationary (i.e., non-constant) mean is used, both for
> simulation and prediction purposes the variogram model defined should be
> that of the residual process, and not that of the raw observations"
> Since my data were, actually, non-stationary, I applied the universal
> co-kriging instead usual co-kriging.
> now the maps of regression-kring and co-kriging are actually similar s
> expected.
> did I understand correctly the quoted sentence?
I think so, but hard to be sure given the information you provide.

>
> regards
>
> emanuele barca
> ------------------------------
>>
>> Message: 2
>> Date: Sat, 10 Aug 2019 10:41:38 +0200
>> From: Edzer Pebesma <[hidden email]>
>> To: [hidden email]
>> Subject: Re: [R-sig-Geo] regression-kriging and co-kriging
>> Message-ID: <[hidden email]>
>> Content-Type: text/plain; charset="utf-8"
>>
>> Hard to tell from your script. Maybe give a reproducible example?
>>
>> On 8/6/19 1:07 PM, Emanuele Barca wrote:
>>> Dear  r-sig-geo friends,
>>>
>>> I produced two maps garnered in the following way:
>>>
>>> # for regression-kriging
>>> Piezo.map <-autoKrige(LivStat ~  Z, input_data = mydata.sp, new_data
>>> = covariates,  model = "Ste")
>>>
>>> Piezork.pred <- Piezo.map$krige_output$var1.pred
>>> Piezork.coords <- Piezo.map$krige_output@coords
>>> Piezork.out <- as.data.frame(cbind(Piezork.coords, Piezork.pred))
>>> colnames(Piezork.out)[1:2] <- c("X", "Y")
>>> coordinates(Piezork.out) = ~ X + Y
>>> gridded(Piezork.out) <- TRUE
>>>
>>> spplot(Piezork.out, main = list(label = "R-k Hydraulic head", cex =
>>> 1.5))
>>>
>>> #for co-kriging
>>> g <- gstat(id = "Piezo", formula = LivStat ~ 1, data = mydata.sp, set
>>> = list(nocheck = 1))
>>> g <- gstat(g, id = "Z", formula = Z ~ 1, data = mydata.sp, set =
>>> list(nocheck = 1))
>>>
>>> v.g <- variogram(g)
>>>
>>> #g <- gstat(g, id = "Piezo", model = vgm(150, "Mat", 1350, 0.0, kappa
>>> = 1.9), fill.all = T)#
>>> g <- gstat(g, id = "Piezo", model = vgm(0.7, "Ste", 1300, 18, kappa =
>>> 1.9), fill.all = T)#
>>> g.fit <- fit.lmc(v.g, g, fit.lmc = TRUE, correct.diagonal = 1.01) #
>>> fit multivariable variogram model , fit.lmc = TRUE, correct.diagonal
>>> = 1.01
>>> g.fit
>>> plot(v.g, model = g.fit, main = "Fitted Variogram Models - Raw Data")#
>>> #gridded(covariates) <- TRUE
>>> g.cok <- predict(g.fit, newdata = covariates)#grid
>>>
>>> g.cok.pred <- g.cok@data$Piezo.pred
>>> aaaa <- na.omit(g.cok.pred)
>>> g.cok.coords <- g.cok@coords
>>> g.cok.out <- as.data.frame(cbind(g.cok.coords, g.cok.pred))
>>> colnames(g.cok.out)[1:2] <- c("X", "Y")
>>> coordinates(g.cok.out) = ~ X + Y
>>> gridded(g.cok.out) <- TRUE
>>> spplot(g.cok.out, main = list(label = "Hydraulic head with
>>> Co-kriging", cex = 1.5))
>>>
>>> ###########################################################################################################################
>>>
>>>
>>> I am unable to understand why the first map appears as a raster and
>>> the second not, notwithstanding the fact that they are both computed
>>> on the same "covariates" DEM???
>>>
>>> where is the mistake???
>>>
>>> regards
>>>
>>> emanuele
>>>
>>> ________________________________________________________
>>> Emanuele Barca                               Researcher
>>> Water Research Institute                       (IRSA-CNR)
>>> Via De Blasio, 5                       70123 Bari (Italy)
>>> Phone +39 080 5820535               Fax  +39 080 5313365
>>> Mobile +39 340 3420689
>>> _________________________________________________________
>>>
>>>
>>>
>>> ---
>>> Questa e-mail è stata controllata per individuare virus con Avast
>>> antivirus.
>>> https://www.avast.com/antivirus
>>>
>>>     [[alternative HTML version deleted]]
>>>
>>> _______________________________________________
>>> R-sig-Geo mailing list
>>> [hidden email]
>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>
>>
>> --
>> Edzer Pebesma
>> Institute for Geoinformatics
>> Heisenbergstrasse 2, 48151 Muenster, Germany
>> Phone: +49 251 8333081
>>
>> -------------- next part --------------
>> A non-text attachment was scrubbed...
>> Name: pEpkey.asc
>> Type: application/pgp-keys
>> Size: 2472 bytes
>> Desc: not available
>> URL:
>> <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20190810/b5cb9d77/attachment-0001.bin>
>>
>>
>>
>> ------------------------------
>>
>> Subject: Digest Footer
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> [hidden email]
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>>
>> ------------------------------
>>
>> End of R-sig-Geo Digest, Vol 192, Issue 7
>> *****************************************
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
--
Edzer Pebesma
Institute for Geoinformatics
Heisenbergstrasse 2, 48151 Muenster, Germany
Phone: +49 251 8333081

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

pEpkey.asc (2K) Download Attachment