AUTOCORRELACION ESPACIAL VARIABLE DE SALINIDAD
library(geoR)
## --------------------------------------------------------------
## Analysis of Geostatistical Data
## For an Introduction to geoR go to http://www.leg.ufpr.br/geoR
## geoR version 1.8-1 (built on 2020-02-08) is now loaded
## --------------------------------------------------------------
library(readxl)
#----- 1. Cargar los datos desde tabla de excel
datos=read_excel("D:/ESPECIALIZACION/SEMESTRE_1/1. Tratamiento de datos/Clase_9_Geoestadistica3/base_cienaga.xls")
datos
## # A tibble: 114 x 6
## Este Norte prof temp sali oxig
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 976952 1706444 1.75 26.4 29.0 6.34
## 2 970883. 1704880. 1.4 30.3 13.0 9.42
## 3 972704. 1704827. 1.1 30.4 19.2 7.88
## 4 974718 1704874 1.5 28.3 35.0 5.67
## 5 976538 1704874 0.5 27 34.0 5.49
## 6 955344. 1703160. 1.62 26 18.0 5.07
## 7 957298. 1703153. 1.55 26.5 18.8 4.53
## 8 959187. 1703146. 1.41 26.9 15.5 5.51
## 9 961189. 1703130. 0.25 27.5 18.9 6.69
## 10 963090. 1703199. 1.54 29 17.2 8.32
## # ... with 104 more rows
plot(datos[,1:2]) # visualizar coordenadas col 1 y 2 (distribucion espacial)
#----- 2. Convertir datos a geodatos,
geodatos5=as.geodata(datos,coords.col=1:2,data.col=5)
class(geodatos5)
## [1] "geodata"
geodatos5
## $coords
## Este Norte
## [1,] 976952.0 1706444
## [2,] 970882.6 1704880
## [3,] 972704.2 1704827
## [4,] 974718.0 1704874
## [5,] 976538.0 1704874
## [6,] 955343.6 1703160
## [7,] 957297.9 1703153
## [8,] 959186.6 1703146
## [9,] 961189.0 1703130
## [10,] 963090.1 1703199
## [11,] 964999.8 1703147
## [12,] 966967.1 1703160
## [13,] 968960.1 1703126
## [14,] 970905.6 1703116
## [15,] 972615.4 1703029
## [16,] 974694.0 1703070
## [17,] 955374.4 1701239
## [18,] 957299.0 1701231
## [19,] 959186.0 1701241
## [20,] 961173.9 1701221
## [21,] 963061.0 1701237
## [22,] 964966.7 1701202
## [23,] 966985.1 1701245
## [24,] 968980.6 1701206
## [25,] 970818.9 1701205
## [26,] 972695.9 1701130
## [27,] 974699.0 1701133
## [28,] 953180.2 1699228
## [29,] 955353.5 1699227
## [30,] 957270.0 1699205
## [31,] 959183.6 1699265
## [32,] 961148.7 1699222
## [33,] 963078.2 1699297
## [34,] 964978.1 1699247
## [35,] 966991.4 1699272
## [36,] 968883.2 1699239
## [37,] 970785.0 1699269
## [38,] 972755.3 1699166
## [39,] 974763.8 1699209
## [40,] 953280.1 1697307
## [41,] 955393.9 1697383
## [42,] 957294.0 1697305
## [43,] 959144.4 1697424
## [44,] 961172.9 1697311
## [45,] 963052.6 1697384
## [46,] 964967.9 1697350
## [47,] 966960.2 1697405
## [48,] 969019.2 1697291
## [49,] 970820.8 1697331
## [50,] 972692.2 1697312
## [51,] 955389.4 1695340
## [52,] 957272.0 1695313
## [53,] 959190.6 1695330
## [54,] 961161.8 1695337
## [55,] 963083.2 1695344
## [56,] 964929.4 1695307
## [57,] 966984.6 1695356
## [58,] 968947.0 1695369
## [59,] 970923.0 1695411
## [60,] 972704.7 1695314
## [61,] 957305.4 1693421
## [62,] 959249.3 1693418
## [63,] 961193.3 1693539
## [64,] 963076.3 1693475
## [65,] 964959.5 1693473
## [66,] 966964.0 1693471
## [67,] 968938.3 1693469
## [68,] 970851.7 1693468
## [69,] 972492.2 1693958
## [70,] 957212.0 1691577
## [71,] 959155.9 1691575
## [72,] 961160.5 1691450
## [73,] 963043.7 1691447
## [74,] 964927.0 1691476
## [75,] 966840.6 1691505
## [76,] 968906.0 1691503
## [77,] 970819.6 1691501
## [78,] 957179.0 1689549
## [79,] 959183.8 1689547
## [80,] 961158.2 1689514
## [81,] 963041.6 1689573
## [82,] 964955.4 1689663
## [83,] 967020.9 1689538
## [84,] 968934.6 1689567
## [85,] 970939.5 1689627
## [86,] 959181.4 1687550
## [87,] 961186.4 1687609
## [88,] 963009.1 1687576
## [89,] 964953.3 1687604
## [90,] 966988.6 1687602
## [91,] 968902.4 1687601
## [92,] 970907.7 1687998
## [93,] 959209.5 1685706
## [94,] 961093.1 1685704
## [95,] 963067.7 1685701
## [96,] 965012.1 1685730
## [97,] 966956.3 1685697
## [98,] 969022.2 1685695
## [99,] 959207.2 1683801
## [100,] 961121.1 1683706
## [101,] 963035.2 1683704
## [102,] 964858.1 1683764
## [103,] 966954.3 1683639
## [104,] 968443.1 1683699
## [105,] 959174.4 1681804
## [106,] 961149.3 1681801
## [107,] 963063.5 1681861
## [108,] 964916.8 1681797
## [109,] 966922.1 1681826
## [110,] 959172.1 1679929
## [111,] 961147.1 1679927
## [112,] 963061.4 1679925
## [113,] 964945.2 1679923
## [114,] 966950.6 1679921
##
## $data
## [1] 28.95 13.03 19.19 34.95 33.99 18.03 18.85 15.47 18.91 17.15 16.41 16.27
## [13] 14.47 13.88 17.89 16.74 18.78 17.28 18.64 15.20 16.27 16.61 15.94 17.55
## [25] 16.95 18.50 15.34 18.44 18.85 17.69 16.21 15.54 15.34 15.60 16.01 16.68
## [37] 16.48 16.21 15.07 18.10 18.16 17.69 16.74 15.40 14.74 14.94 15.07 16.48
## [49] 16.34 16.01 18.71 17.49 16.68 16.34 15.01 15.81 15.47 17.08 16.74 15.20
## [61] 19.25 17.96 17.42 16.21 16.74 16.34 16.54 17.28 16.21 19.32 18.23 17.89
## [73] 17.69 17.62 16.88 17.01 17.08 18.64 18.37 18.50 18.44 18.23 17.15 17.35
## [85] 16.74 19.19 18.98 18.71 18.16 17.55 16.74 16.95 19.46 19.12 19.05 17.96
## [97] 17.76 16.88 19.53 19.25 18.91 18.23 16.88 16.74 18.98 18.91 18.71 18.23
## [109] 17.15 18.16 18.50 17.89 17.55 16.61
##
## attr(,"class")
## [1] "geodata"
plot(geodatos5)
#----- 3. Cálculo del semivariagram muestral
summary(dist(geodatos5$coords))
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1373 6994 11204 11531 15606 31924
variograma5=variog(geodatos5,option= "bin",uvec=seq(0,20000,1500))
## variog: computing omnidirectional variogram
variograma5
## $u
## [1] 1500 3000 4500 6000 7500 9000 10500 12000 13500 15000 16500 18000
## [13] 19500
##
## $v
## [1] 2.139230 3.323165 3.205625 4.218861 4.155430 5.149708 6.334606
## [8] 5.290376 6.869555 9.191472 11.714082 11.160215 14.239824
##
## $n
## [1] 201 198 494 584 688 510 571 650 601 438 316 396 298
##
## $sd
## [1] 14.64185 20.04953 19.55825 22.52982 19.85200 22.48796 26.39760 22.22508
## [9] 26.67721 31.63044 35.42439 32.96839 36.85020
##
## $bins.lim
## [1] 1.000e-12 7.500e+02 2.250e+03 3.750e+03 5.250e+03 6.750e+03 8.250e+03
## [8] 9.750e+03 1.125e+04 1.275e+04 1.425e+04 1.575e+04 1.725e+04 1.875e+04
## [15] 2.025e+04
##
## $ind.bin
## [1] FALSE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [13] TRUE TRUE
##
## $var.mark
## [1] 8.13333
##
## $beta.ols
## [1] 17.62623
##
## $output.type
## [1] "bin"
##
## $max.dist
## [1] 20250
##
## $estimator.type
## [1] "classical"
##
## $n.data
## [1] 114
##
## $lambda
## [1] 1
##
## $trend
## [1] "cte"
##
## $pairs.min
## [1] 2
##
## $nugget.tolerance
## [1] 1e-12
##
## $direction
## [1] "omnidirectional"
##
## $tolerance
## [1] "none"
##
## $uvec
## [1] 0 1500 3000 4500 6000 7500 9000 10500 12000 13500 15000 16500
## [13] 18000 19500
##
## $call
## variog(geodata = geodatos5, uvec = seq(0, 20000, 1500), option = "bin")
##
## attr(,"class")
## [1] "variogram"
plot(variograma5, pch=16)
#----- 4. Permuto datos originales para hacer el calculo del semivariograma
geodatos5.env=variog.mc.env(geodatos5,obj=variograma5) ## metodo montecarlo (variog.mc) calcula el semivariagrama sin autocorrelacion
## variog.env: generating 99 simulations by permutating data values
## variog.env: computing the empirical variogram for the 99 simulations
## variog.env: computing the envelops
plot(variograma5,pch=16,main=names(datos)[5],envelope =geodatos5.env)
#----- 5. Encuentro la curva que mas se ajuste al semivariograma usuando los modelo gaus, exp o spe (meseta y rango)
ini.vals = expand.grid(seq(2,8,l=10), seq(10000,15000,l=10)) ## rango para encontrar la meseta (valores maximo donde esta la meseta (entre 5,7)
model_mco_exp=variofit(variograma5, ini=ini.vals, cov.model="exponential",wei="npair", min="optim")
## variofit: covariance model used is exponential
## variofit: weights used: npairs
## variofit: minimisation function used: optim
## variofit: searching for best initial value ... selected values:
## sigmasq phi tausq kappa
## initial.value "8" "10000" "0" "0.5"
## status "est" "est" "est" "fix"
## loss value: 38967.7442784068
model_mco_gaus=variofit(variograma5, ini=ini.vals, cov.model="gaussian", wei="npair", min="optim",nugget = 0)
## variofit: covariance model used is gaussian
## variofit: weights used: npairs
## variofit: minimisation function used: optim
## variofit: searching for best initial value ... selected values:
## sigmasq phi tausq kappa
## initial.value "8" "10000" "0" "0.5"
## status "est" "est" "est" "fix"
## loss value: 31695.8923874857
model_mco_spe=variofit(variograma5, ini=ini.vals, cov.model="spheric",fix.nug=TRUE, wei="npair", min="optim")
## variofit: covariance model used is spherical
## variofit: weights used: npairs
## variofit: minimisation function used: optim
## variofit: searching for best initial value ... selected values:
## sigmasq phi tausq kappa
## initial.value "8" "15000" "0" "0.5"
## status "est" "est" "fix" "fix"
## loss value: 27179.3901649479
lines(model_mco_exp,col="blue")
lines(model_mco_gaus,col="red")
lines(model_mco_spe,col="purple")
## Suma mÃnima ponderada de cuadrados:
#--gaussino: 427331.7
#--esferico:7511.132
#--exponential:1945.665 *
#----- 6. Prediccion espacial con el metodo de Kriging
plot(geodatos5$coords)
loc0=cbind( 966988.6, 1687602)
geodatos_ko=krige.conv(geodatos5, loc=loc0,
krige= krige.control(nugget=0,trend.d="cte",
trend.l="cte",cov.pars=c(sigmasq=9.04, phi=13858))) # sigmasq(9.04), phi(13858) del exponencial.
## krige.conv: model with constant mean
## krige.conv: Kriging performed using global neighbourhood
#----- 7. Generacion de la grilla
geodatos_grid=expand.grid(Este=seq(952023,979916,l=100), Norte=seq(1677586,1708649,l=100))
plot(geodatos_grid)
points(geodatos5$coord,col="red",pch=16)
geodatos_ko=krige.conv(geodatos5, loc=geodatos_grid,
krige= krige.control(nugget=0,trend.d="cte",
trend.l="cte",cov.pars=c(sigmasq=10.0246, phi=8804.8018)))
## krige.conv: model with constant mean
## krige.conv: Kriging performed using global neighbourhood
image(geodatos_ko, main="kriging Predict", xlab="East", ylab="North")
image(geodatos_ko, main="kriging StDv Predicted",val=sqrt(geodatos_ko$krige.var), xlab="East", ylab="North")
