Compactación del suelo
Olinda, R. A.
Noviembre, 2016
Universidad de San Carlos de Guatemala
Centro Universitario del Norte
CUNOR
Use R!
R Markdown
Sobre Programación básica en R, puede ser consultada Aqui
Comienza cargando las librerías y preparando los datos a utilizar.
library(geoR) #El análisis de los datos Geoestadísticos## --------------------------------------------------------------
## Analysis of Geostatistical Data
## For an Introduction to geoR go to http://www.leg.ufpr.br/geoR
## geoR version 1.7-5.2 (built on 2016-05-02) is now loaded
## --------------------------------------------------------------library(MASS) #Estadística descriptiva de
library(moments) #Estadística descriptiva
library(scatterplot3d) #Obtener gráficos de dispersión 3dPara realizar los análisis exploratorios a seguir, utilizaremos el paquete geoR e inicialmente necesitaremos crear un objeto del tipo geodata, como sigue.
Importacion del datos para el R
suelo= read.table('castro.txt', head=TRUE)
head(suelo)## Lat Lon X0.10 X10.15 X15.20 X20.25 X25.30 X30.35 X35.40
## 1 608039.6 7250137 2.06 2.55 1.83 1.54 1.12 0.84 0.80
## 2 608063.1 7250146 2.03 1.58 2.21 2.36 2.61 2.44 1.76
## 3 608086.6 7250154 1.67 1.38 1.61 1.79 2.13 1.46 1.36
## 4 608007.6 7250152 1.23 1.02 0.96 1.19 1.39 1.45 0.76
## 5 608031.1 7250161 0.57 0.75 1.22 1.69 2.12 2.16 2.34
## 6 608054.6 7250169 1.67 1.60 1.20 1.75 1.70 1.58 1.21
## Xcd Ycd
## 1 -49.93059 -24.85994
## 2 -49.93036 -24.85986
## 3 -49.93013 -24.85978
## 4 -49.93091 -24.85980
## 5 -49.93068 -24.85973
## 6 -49.93045 -24.85965Análise descritiva profundidade
geosuelo <- as.geodata(suelo, coords.col=c(1, 2), data.col=7)
class(geosuelo)## [1] "geodata"attach(geosuelo)Argumentos utilizados por la función as.geodata() anterior:
obj - Compactación del suelo (base de datos);
coords.col - 1:2 (columnas con la longitud y la latitud en UTM);
data.col - 7 (variable de interés).
Medidas de resumen son facilmente obtenidas, utilizando la función summary().
summary(geosuelo)## Number of data points: 358
##
## Coordinates summary
## Lat Lon
## min 607546.1 7250137
## max 608125.0 7250863
##
## Distance summary
## min max
## 24.98484 728.01244
##
## Data summary
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.490 1.362 1.695 1.786 2.090 3.990A seguir son mostrados dos gráficos que ayudan a visualizar los datos.
plot(geosuelo, low=TRUE)
scatterplot3d(geosuelo$coords[,1], geosuelo$coords[,2],
geosuelo$data,pch=1,xlab="Longitud",ylab="Latitud", zlab="compactación del suelo",color="red", main="Gráfico de Dispersión", zlim=c(0,5))
Generamos la grilla de trabajo se puede usar otra funcion expand.grid, en help pueden ver las opciones
bord = read.table('bord4.txt')
geosuelo$borders <- with(bord, cbind(V1,V2)) Otras formas de visualizar los datos, pueden ser obtenida de la siguinte manera.
par(mfrow=c(1,2))
points(geosuelo,xlab="Longitud",ylab="Latitud")
points(geosuelo,xlab="Longitud",ylab="Latitud",col=gray(seq(1,0.1,l=90)))
Estadística descriptiva de compactación del suelo
y1<-mean(suelo$X25.30)
y2<-min(suelo$X25.30)
y3<-max(suelo$X25.30)
y4<-sd(suelo$X25.30)
y5<-100*sd(suelo$X25.30)/mean(suelo$X25.30)
y6<-kurtosis(suelo$X25.30)
y7<-skewness(suelo$X25.30)
descriptiva<-data.frame(y1,y2,y3,y4,y5,y6,y7)
descriptiva## y1 y2 y3 y4 y5 y6 y7
## 1 1.785726 0.49 3.99 0.5901606 33.04877 4.175745 0.9699079Con el fin de hacer constante la variabilidad en los datos, se ha evaludado la posibilidad de transformarlos usando la familia de transformaciones potenciales de Box-Cox.
boxcox(geosuelo)
abline(v=-1,col="red")
La expresión de la variable transformada es:
\(datos^{trans}=\frac{datos^{0.5}-1}{0.5}\)
plot(geosuelo, low=T,lam=0) #Comprobación de los supuestos
Comprobando gráficamente la dependencia espacial con 100 simulaciones
par(mfrow=c(2,2))
var1 = variog(geosuelo, max.dist=350)## variog: computing omnidirectional variogramplot (var1)
var2 = variog(geosuelo, option="cloud")## variog: computing omnidirectional variogramplot (var2)
var3 = variog(geosuelo,uvec=seq(0,350,l=50))## variog: computing omnidirectional variogramplot(var3)
env.var = variog.mc.env(geosuelo, obj.v=var3, nsim=100) ## variog.env: generating 100 simulations by permutating data values
## variog.env: computing the empirical variogram for the 100 simulations
## variog.env: computing the envelopsplot(var3, env=env.var)
Si al menos un punto está fuera de “envelope” del simulación hay evidencia de que existe una dependencia espacial.
De acuerdo con lo que discutimos por la mañana, es necesario estudiar la superficie de tendencia.
plot(geosuelo,lowess=TRUE)#Superficie de tendencia constante.
plot(geosuelo,lowess=TRUE,lam=0) #Superficie de tendencia constante con los datos transformados
plot(geosuelo,lowess=TRUE,trend="1st") #Superficie de tendencia de 1er orden.
plot(geosuelo,lowess=TRUE,trend="1st",lam=0) #Superficie de tendencia de 1er orden con los datos transformados.
plot(geosuelo,lowess=TRUE,trend="2nd") #Superficie de tendencia de 2er orden.
plot(geosuelo,lowess=TRUE,trend=~coords[,2]) #Superficie de tendencia de 2er orden.
De acuerdo con los gráficos anteriores, con cual superficie de tendencia trabajaríamos?
Además de la superficie de tendencia, debemos entender la estructura de covarianza presente en los datos (esto si ella existe). A seguir son obtenidos los variogramas por medio de la función variog().
varc=variog(geosuelo, max.dist=350)## variog: computing omnidirectional variogramvarct=variog(geosuelo, lam=0, max.dist=350)## variog: computing omnidirectional variogramvar1st=variog(geosuelo, trend="1st", max.dist=350)## variog: computing omnidirectional variogramvar1stt=variog(geosuelo, lam=0,trend="1st", max.dist=350)## variog: computing omnidirectional variogramArgumentos utilizados por la función variog() anterior:
geodata - geosuelo (un objeto geodata);
trend - 1st (tipo de tendencia);
max.dist - 350 (distancia máxima).
Abajo se encuentran los comandos para plotar los variogramas.
par(mfrow=c(2,2))
plot(varc,main="Tendencia constante",xlab="Distancia en grados",ylab="Semivarianza")
plot(varct,main="Tendencia constante/transformación",xlab="Distancia en grados",ylab="Semivarianza")
plot(var1st,main="Tendencia de primer orden",xlab="Distancia en grados",ylab="Semivarianza")
plot(var1stt,main="Tendencia de primer orden/transformacion",xlab="Distancia en grados",ylab="Semivarianza")
Para concluir este tópico mostramos una función auxiliar que facilita la construcción de variogramas empíricos en direcciones diferentes que pueden ser utilizados para investigar anisotropía.
vario4c=variog4(geosuelo, max.dist=250)## variog: computing variogram for direction = 0 degrees (0 radians)
## tolerance angle = 22.5 degrees (0.393 radians)
## variog: computing variogram for direction = 45 degrees (0.785 radians)
## tolerance angle = 22.5 degrees (0.393 radians)
## variog: computing variogram for direction = 90 degrees (1.571 radians)
## tolerance angle = 22.5 degrees (0.393 radians)
## variog: computing variogram for direction = 135 degrees (2.356 radians)
## tolerance angle = 22.5 degrees (0.393 radians)
## variog: computing omnidirectional variogramvario4ct=variog4(geosuelo, lam=0, max.dist=250)## variog: computing variogram for direction = 0 degrees (0 radians)
## tolerance angle = 22.5 degrees (0.393 radians)
## variog: computing variogram for direction = 45 degrees (0.785 radians)
## tolerance angle = 22.5 degrees (0.393 radians)
## variog: computing variogram for direction = 90 degrees (1.571 radians)
## tolerance angle = 22.5 degrees (0.393 radians)
## variog: computing variogram for direction = 135 degrees (2.356 radians)
## tolerance angle = 22.5 degrees (0.393 radians)
## variog: computing omnidirectional variogramvario41st=variog4(geosuelo, trend="1st", max.dist=250)## variog: computing variogram for direction = 0 degrees (0 radians)
## tolerance angle = 22.5 degrees (0.393 radians)
## variog: computing variogram for direction = 45 degrees (0.785 radians)
## tolerance angle = 22.5 degrees (0.393 radians)
## variog: computing variogram for direction = 90 degrees (1.571 radians)
## tolerance angle = 22.5 degrees (0.393 radians)
## variog: computing variogram for direction = 135 degrees (2.356 radians)
## tolerance angle = 22.5 degrees (0.393 radians)
## variog: computing omnidirectional variogramvario41stt=variog4(geosuelo, lam=0,trend="1st", max.dist=250)## variog: computing variogram for direction = 0 degrees (0 radians)
## tolerance angle = 22.5 degrees (0.393 radians)
## variog: computing variogram for direction = 45 degrees (0.785 radians)
## tolerance angle = 22.5 degrees (0.393 radians)
## variog: computing variogram for direction = 90 degrees (1.571 radians)
## tolerance angle = 22.5 degrees (0.393 radians)
## variog: computing variogram for direction = 135 degrees (2.356 radians)
## tolerance angle = 22.5 degrees (0.393 radians)
## variog: computing omnidirectional variogrampar(mfrow=c(1,2))
plot(vario4c, omni=TRUE,legend=FALSE) #Incluyendo el variograma omnidireccional
legend("bottomright", c("omnid.", "0º","45º","90º","135º"), lty=2, col=c("black","red","green","blue"),xpd=TRUE,cex=0.8)
plot(vario4ct, omni=TRUE,legend=FALSE)
legend("bottomright", c("omnid.", "0º","45º","90º","135º"), lty=2, col=c("black","red","green","blue"),xpd=TRUE,cex=0.8)
# Gráfico con la firma, la fecha y la versión
mtext(side=3, at=par("usr")[3], adj=0.7,
cex=0.8, col="gray40", line=1,
text=paste("Ricardo Olinda --",
format(Sys.time(), "%d/%m/%Y %H:%M:%S --"),
R.version.string))
par(mfrow=c(1,2))
plot(vario41st, omni=TRUE,legend=FALSE)
legend("bottomright", c("omnid.", "0º","45º","90º","135º"), lty=2, col=c("black","red","green","blue"),xpd=TRUE,cex=0.8)
plot(vario41stt, omni=TRUE,legend=FALSE)
legend("bottomright", c("omnid.", "0º","45º","90º","135º"), lty=2, col=c("black","red","green","blue"),xpd=TRUE,cex=0.8)
# Gráfico con la firma, la fecha y la versión
mtext(side=3, at=par("usr")[3], adj=0.7,
cex=0.8, col="gray40", line=1,
text=paste("Ricardo Olinda --",
format(Sys.time(), "%d/%m/%Y %H:%M:%S --"),
R.version.string))
Variograma con diferentes gráficos en comparación con el efecto omnidireccional
plot(vario4c, omni=TRUE, same=FALSE)
Tendrás que estudiar adecuadamente la cuestión de la anisotropía en su datos.
Esto no le exime de la anisotropía efecto combinado + tendencia.
Sólo para probar las posibilidades, vamos a hacer el variograma hacia 45 grados.
va.45 <- variog(geosuelo, trend="1st",max.dist=250, dir=pi/4,estimator="classical")## variog: computing variogram for direction = 45 degrees (0.785 radians)
## tolerance angle = 22.5 degrees (0.393 radians)plot(va.45)
Variograma visual sólo para dar “disparos” iniciales para los parámetros semivariagrama
v1 = variog(geosuelo, max.dist=250,dir=pi/4,estimator="classical") ## variog: computing variogram for direction = 45 degrees (0.785 radians)
## tolerance angle = 22.5 degrees (0.393 radians)plot(v1)
## Ajuste (visual) del modelo para el variograma
#ef1 = eyefit(v1)
#summary(ef1)Estimación de los parámetros (máxima verosimilitud) los valores iniciales “ef1” sugerido visualmente al ver la función variograma del modelo ajustado por máxima verosimilitud (ML).
Las estimativas de máxima verosimilitud son obtenidas ajustando el modelo a los datos. Los modelos pueden ser comparados por medio de medidas de comparación de modelos, como, por ejemplo, el AIC y BIC.
ml1 <- likfit(geosuelo, ini=c(0.30,80.13), nug=0.04,
cov.model= "matern",kappa=0.5)## ---------------------------------------------------------------
## likfit: likelihood maximisation using the function optim.
## likfit: Use control() to pass additional
## arguments for the maximisation function.
## For further details see documentation for optim.
## likfit: It is highly advisable to run this function several
## times with different initial values for the parameters.
## likfit: WARNING: This step can be time demanding!
## ---------------------------------------------------------------
## likfit: end of numerical maximisation.ml1## likfit: estimated model parameters:
## beta tausq sigmasq phi
## " 2.0374" " 0.1919" " 0.2245" "100.9552"
## Practical Range with cor=0.05 for asymptotic range: 302.4349
##
## likfit: maximised log-likelihood = -284.6summary(ml1)## Summary of the parameter estimation
## -----------------------------------
## Estimation method: maximum likelihood
##
## Parameters of the mean component (trend):
## beta
## 2.0374
##
## Parameters of the spatial component:
## correlation function: exponential
## (estimated) variance parameter sigmasq (partial sill) = 0.2245
## (estimated) cor. fct. parameter phi (range parameter) = 101
## anisotropy parameters:
## (fixed) anisotropy angle = 0 ( 0 degrees )
## (fixed) anisotropy ratio = 1
##
## Parameter of the error component:
## (estimated) nugget = 0.1919
##
## Transformation parameter:
## (fixed) Box-Cox parameter = 1 (no transformation)
##
## Practical Range with cor=0.05 for asymptotic range: 302.4349
##
## Maximised Likelihood:
## log.L n.params AIC BIC
## "-284.6" "4" "577.2" "592.7"
##
## non spatial model:
## log.L n.params AIC BIC
## "-318.7" "2" "641.4" "649.1"
##
## Call:
## likfit(geodata = geosuelo, ini.cov.pars = c(0.3, 80.13), nugget = 0.04,
## kappa = 0.5, cov.model = "matern")ml2 <- likfit(geosuelo, ini=c(0.2940,32.48), nug=0.0367,
cov.model= "matern",kappa=0.5)## ---------------------------------------------------------------
## likfit: likelihood maximisation using the function optim.
## likfit: Use control() to pass additional
## arguments for the maximisation function.
## For further details see documentation for optim.
## likfit: It is highly advisable to run this function several
## times with different initial values for the parameters.
## likfit: WARNING: This step can be time demanding!
## ---------------------------------------------------------------
## likfit: end of numerical maximisation.ml2## likfit: estimated model parameters:
## beta tausq sigmasq phi
## " 2.0373" " 0.1919" " 0.2244" "100.9269"
## Practical Range with cor=0.05 for asymptotic range: 302.35
##
## likfit: maximised log-likelihood = -284.6summary(ml1)## Summary of the parameter estimation
## -----------------------------------
## Estimation method: maximum likelihood
##
## Parameters of the mean component (trend):
## beta
## 2.0374
##
## Parameters of the spatial component:
## correlation function: exponential
## (estimated) variance parameter sigmasq (partial sill) = 0.2245
## (estimated) cor. fct. parameter phi (range parameter) = 101
## anisotropy parameters:
## (fixed) anisotropy angle = 0 ( 0 degrees )
## (fixed) anisotropy ratio = 1
##
## Parameter of the error component:
## (estimated) nugget = 0.1919
##
## Transformation parameter:
## (fixed) Box-Cox parameter = 1 (no transformation)
##
## Practical Range with cor=0.05 for asymptotic range: 302.4349
##
## Maximised Likelihood:
## log.L n.params AIC BIC
## "-284.6" "4" "577.2" "592.7"
##
## non spatial model:
## log.L n.params AIC BIC
## "-318.7" "2" "641.4" "649.1"
##
## Call:
## likfit(geodata = geosuelo, ini.cov.pars = c(0.3, 80.13), nugget = 0.04,
## kappa = 0.5, cov.model = "matern")names(ml1)## [1] "cov.model" "nugget"
## [3] "cov.pars" "sigmasq"
## [5] "phi" "kappa"
## [7] "beta" "beta.var"
## [9] "lambda" "aniso.pars"
## [11] "tausq" "practicalRange"
## [13] "method.lik" "trend"
## [15] "loglik" "npars"
## [17] "AIC" "BIC"
## [19] "parameters.summary" "info.minimisation.function"
## [21] "max.dist" "trend"
## [23] "trend.matrix" "transform.info"
## [25] "nospatial" "model.components"
## [27] "call"ml2 <- likfit(geosuelo, ini=c(0.2940,32.48), nug=0.0367,
cov.model= "matern",trend="1st", kappa=1)## ---------------------------------------------------------------
## likfit: likelihood maximisation using the function optim.
## likfit: Use control() to pass additional
## arguments for the maximisation function.
## For further details see documentation for optim.
## likfit: It is highly advisable to run this function several
## times with different initial values for the parameters.
## likfit: WARNING: This step can be time demanding!
## ---------------------------------------------------------------
## likfit: end of numerical maximisation.ml2## likfit: estimated model parameters:
## beta0 beta1 beta2 tausq sigmasq
## "-4192.6300" " -0.0004" " 0.0006" " 0.2169" " 0.1623"
## phi
## " 51.8679"
## Practical Range with cor=0.05 for asymptotic range: 207.3951
##
## likfit: maximised log-likelihood = -284.5summary(ml2)## Summary of the parameter estimation
## -----------------------------------
## Estimation method: maximum likelihood
##
## Parameters of the mean component (trend):
## beta0 beta1 beta2
## -4192.6300 -0.0004 0.0006
##
## Parameters of the spatial component:
## correlation function: matern
## (estimated) variance parameter sigmasq (partial sill) = 0.1623
## (estimated) cor. fct. parameter phi (range parameter) = 51.87
## (fixed) extra parameter kappa = 1
## anisotropy parameters:
## (fixed) anisotropy angle = 0 ( 0 degrees )
## (fixed) anisotropy ratio = 1
##
## Parameter of the error component:
## (estimated) nugget = 0.2169
##
## Transformation parameter:
## (fixed) Box-Cox parameter = 1 (no transformation)
##
## Practical Range with cor=0.05 for asymptotic range: 207.3951
##
## Maximised Likelihood:
## log.L n.params AIC BIC
## "-284.5" "6" "580.9" "604.2"
##
## non spatial model:
## log.L n.params AIC BIC
## "-313.7" "4" "635.3" "650.8"
##
## Call:
## likfit(geodata = geosuelo, trend = "1st", ini.cov.pars = c(0.294,
## 32.48), nugget = 0.0367, kappa = 1, cov.model = "matern")ml3 <- likfit(geosuelo, ini=c(0.2940,32.48), nug=0.0367,
cov.model="spherical",trend="1st")## kappa not used for the spherical correlation function
## ---------------------------------------------------------------
## likfit: likelihood maximisation using the function optim.
## likfit: Use control() to pass additional
## arguments for the maximisation function.
## For further details see documentation for optim.
## likfit: It is highly advisable to run this function several
## times with different initial values for the parameters.
## likfit: WARNING: This step can be time demanding!
## ---------------------------------------------------------------
## likfit: end of numerical maximisation.ml3## likfit: estimated model parameters:
## beta0 beta1 beta2 tausq sigmasq
## "-3762.8820" " -0.0003" " 0.0005" " 0.0717" " 0.2627"
## phi
## " 49.8717"
## Practical Range with cor=0.05 for asymptotic range: 49.8717
##
## likfit: maximised log-likelihood = -290.9summary(ml3)## Summary of the parameter estimation
## -----------------------------------
## Estimation method: maximum likelihood
##
## Parameters of the mean component (trend):
## beta0 beta1 beta2
## -3762.8820 -0.0003 0.0005
##
## Parameters of the spatial component:
## correlation function: spherical
## (estimated) variance parameter sigmasq (partial sill) = 0.2627
## (estimated) cor. fct. parameter phi (range parameter) = 49.87
## anisotropy parameters:
## (fixed) anisotropy angle = 0 ( 0 degrees )
## (fixed) anisotropy ratio = 1
##
## Parameter of the error component:
## (estimated) nugget = 0.0717
##
## Transformation parameter:
## (fixed) Box-Cox parameter = 1 (no transformation)
##
## Practical Range with cor=0.05 for asymptotic range: 49.8717
##
## Maximised Likelihood:
## log.L n.params AIC BIC
## "-290.9" "6" "593.8" "617.1"
##
## non spatial model:
## log.L n.params AIC BIC
## "-313.7" "4" "635.3" "650.8"
##
## Call:
## likfit(geodata = geosuelo, trend = "1st", ini.cov.pars = c(0.294,
## 32.48), nugget = 0.0367, cov.model = "spherical")ml4 <- likfit(geosuelo, ini=c(0.2940,32.48), nug=0.0367,
cov.model="powered.exponential",trend="1st")## ---------------------------------------------------------------
## likfit: likelihood maximisation using the function optim.
## likfit: Use control() to pass additional
## arguments for the maximisation function.
## For further details see documentation for optim.
## likfit: It is highly advisable to run this function several
## times with different initial values for the parameters.
## likfit: WARNING: This step can be time demanding!
## ---------------------------------------------------------------
## likfit: end of numerical maximisation.ml4## likfit: estimated model parameters:
## beta0 beta1 beta2 tausq sigmasq
## "-4177.7057" " -0.0004" " 0.0006" " 0.0000" " 0.3823"
## phi
## " 26.7915"
## Practical Range with cor=0.05 for asymptotic range: 240.4383
##
## likfit: maximised log-likelihood = -283.7summary(ml4)## Summary of the parameter estimation
## -----------------------------------
## Estimation method: maximum likelihood
##
## Parameters of the mean component (trend):
## beta0 beta1 beta2
## -4177.7057 -0.0004 0.0006
##
## Parameters of the spatial component:
## correlation function: powered.exponential
## (estimated) variance parameter sigmasq (partial sill) = 0.3823
## (estimated) cor. fct. parameter phi (range parameter) = 26.79
## (fixed) extra parameter kappa = 0.5
## anisotropy parameters:
## (fixed) anisotropy angle = 0 ( 0 degrees )
## (fixed) anisotropy ratio = 1
##
## Parameter of the error component:
## (estimated) nugget = 0
##
## Transformation parameter:
## (fixed) Box-Cox parameter = 1 (no transformation)
##
## Practical Range with cor=0.05 for asymptotic range: 240.4383
##
## Maximised Likelihood:
## log.L n.params AIC BIC
## "-283.7" "6" "579.4" "602.7"
##
## non spatial model:
## log.L n.params AIC BIC
## "-313.7" "4" "635.3" "650.8"
##
## Call:
## likfit(geodata = geosuelo, trend = "1st", ini.cov.pars = c(0.294,
## 32.48), nugget = 0.0367, cov.model = "powered.exponential")ml5 <- likfit(geosuelo, ini=c(0.2940,32.48), nug=0.0367,
cov.model="circular",trend="1st")## kappa not used for the circular correlation function
## ---------------------------------------------------------------
## likfit: likelihood maximisation using the function optim.
## likfit: Use control() to pass additional
## arguments for the maximisation function.
## For further details see documentation for optim.
## likfit: It is highly advisable to run this function several
## times with different initial values for the parameters.
## likfit: WARNING: This step can be time demanding!
## ---------------------------------------------------------------
## likfit: end of numerical maximisation.ml5## likfit: estimated model parameters:
## beta0 beta1 beta2 tausq sigmasq
## "-3762.8112" " -0.0003" " 0.0005" " 0.0946" " 0.2397"
## phi
## " 45.7679"
## Practical Range with cor=0.05 for asymptotic range: 45.76795
##
## likfit: maximised log-likelihood = -290.9summary(ml5)## Summary of the parameter estimation
## -----------------------------------
## Estimation method: maximum likelihood
##
## Parameters of the mean component (trend):
## beta0 beta1 beta2
## -3762.8112 -0.0003 0.0005
##
## Parameters of the spatial component:
## correlation function: circular
## (estimated) variance parameter sigmasq (partial sill) = 0.2397
## (estimated) cor. fct. parameter phi (range parameter) = 45.77
## anisotropy parameters:
## (fixed) anisotropy angle = 0 ( 0 degrees )
## (fixed) anisotropy ratio = 1
##
## Parameter of the error component:
## (estimated) nugget = 0.0946
##
## Transformation parameter:
## (fixed) Box-Cox parameter = 1 (no transformation)
##
## Practical Range with cor=0.05 for asymptotic range: 45.76795
##
## Maximised Likelihood:
## log.L n.params AIC BIC
## "-290.9" "6" "593.8" "617.1"
##
## non spatial model:
## log.L n.params AIC BIC
## "-313.7" "4" "635.3" "650.8"
##
## Call:
## likfit(geodata = geosuelo, trend = "1st", ini.cov.pars = c(0.294,
## 32.48), nugget = 0.0367, cov.model = "circular")plot(variog(geosuelo, max.dist=450, trend="1st",dir=pi/4,estimator="classical"))## variog: computing variogram for direction = 45 degrees (0.785 radians)
## tolerance angle = 22.5 degrees (0.393 radians)lines.variomodel(ml1,col="yellow")
lines.variomodel(ml2,col="red")
lines.variomodel(ml3,col="green")
lines.variomodel(ml4,col="black")
lines.variomodel(ml5,col="blue")
legend("bottomright", c("matérn k=0,5", "matérn k=1,0","esférico","gaussiano",
"circular"), lty=1, col=c("red","yellow","green","black","blue"))
#http://listas.inf.ufpr.br/pipermail/r-br/2014-October/014618.html
title("Ajuste de variograma con diferentes funciones de correlación")
Los modelos suponiendo una función de correlación espacial mattern tendrian un mejor ajuste comparado com los ajustes suponiendo una función de correlación esférica?
plot(variog(geosuelo, max.dist=450, trend="1st",dir=pi/4,estimator="classical"))## variog: computing variogram for direction = 45 degrees (0.785 radians)
## tolerance angle = 22.5 degrees (0.393 radians)lines.variomodel(ml1,col="black")
legend("bottomright", c("matérn k=0,5"), lty=1, col=c("black"))
title("Ajuste del variograma con función de correlación Matern")
ml1
Malla de puntos de predicción
gr <- pred_grid(geosuelo$borders, by=20)
points(geosuelo)
points(gr, col=2, pch=19, cex=0.3)
gr0 <- gr[.geoR_inout(gr, geosuelo$borders),]
points(gr0, col=3, pch=19, cex=0.3)
#title("Malla de puntos de predicción")
dim(gr)## [1] 1178 2dim(gr0)## [1] 543 2### Simulación!!
gr = pred_grid(geosuelo$borders, by=20)
s.out = output.control(n.predictive = 1000, n.post=1000, quant=0.95, thres=2.5)
geosuelo.kc = krige.conv(geosuelo, loc=gr, krige=krige.control(obj=ml1), output = s.out)## krige.conv: results will be returned only for prediction locations inside the borders
## krige.conv: model with constant mean
## krige.conv: sampling from the predictive distribution (conditional simulations)
## krige.conv: Kriging performed using global neighbourhoodimage(geosuelo.kc, col = terrain.colors(200), main="Mapa de P(Y > 2.5)",
val=(1-geosuelo.kc$probabilities.simulations),
x.leg=c(607500, 607890), y.leg=c(7250130,7250160))
Viendo el mapa del kriging
KC <- krige.control(type="OK", obj.model=ml1)
kc1 <- krige.conv(geosuelo, loc=gr, krige=KC)## krige.conv: results will be returned only for prediction locations inside the borders
## krige.conv: model with constant mean
## krige.conv: Kriging performed using global neighbourhoodimage(kc1,main="Mapa de interpolación")
image(kc1, col=terrain.colors(200),
x.leg=c(607500,607800), y.leg=c(7250150, 7250180), main="Mapa de interpolación")