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  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
Reply | Threaded
Open this post in threaded view
|

Re: regression-kriging and co-kriging

edzer
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

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

pEpkey.asc (2K) Download Attachment