require(geoR)
## Loading required package: 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)
datos= read_excel("~/Especializacion_geomatica/Tratamiento_de_datos_espaciales/semivariograma_clas8/base_cienaga.xls")
head(datos)
## # A tibble: 6 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
plot(datos[,1:2])
#convertir en dato GeoR
geodatos=as.geodata(datos,coords.col=1:2,data.col=4)
class(geodatos)
## [1] "geodata"
geodatos
## $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] 26.4 30.3 30.4 28.3 27.0 26.0 26.5 26.9 27.5 29.0 29.2 31.3 31.3 29.4 30.5
## [16] 31.7 26.2 26.8 27.2 27.8 29.0 29.5 32.5 31.0 31.3 31.0 31.7 26.1 26.2 26.7
## [31] 27.0 27.5 28.6 29.5 30.4 33.1 32.1 32.3 32.4 26.0 26.0 26.6 26.5 27.3 28.1
## [46] 29.2 30.9 33.0 32.3 31.5 26.4 26.3 26.7 27.8 28.1 30.4 30.3 32.4 32.3 33.0
## [61] 29.6 32.1 31.5 30.6 31.3 32.3 31.7 32.2 33.2 31.4 30.8 29.9 31.6 31.7 31.8
## [76] 31.4 31.2 30.2 30.1 31.2 30.8 31.8 31.7 32.0 31.2 30.2 28.2 30.1 30.6 30.4
## [91] 30.8 30.5 28.0 27.3 27.5 28.4 29.0 30.3 27.4 27.5 27.6 27.4 28.4 28.8 27.3
## [106] 27.2 27.0 27.4 28.3 27.1 26.8 27.1 27.9 29.3
##
## attr(,"class")
## [1] "geodata"
plot(geodatos)
#semivariograma muestral
summary(dist(geodatos$coords))
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1373 6994 11204 11531 15606 31924
plot(variog(geodatos,option = "cloud"))
## variog: computing omnidirectional variogram
variograma=variog(geodatos,option= "bin",uvec=seq(0,20000,1500))
## variog: computing omnidirectional variogram
variograma
## $u
## [1] 1500 3000 4500 6000 7500 9000 10500 12000 13500 15000 16500 18000
## [13] 19500
##
## $v
## [1] 0.6501244 1.0682828 1.6100709 2.5685702 3.6511337 4.4781569 5.3035552
## [8] 5.8860615 6.2235025 6.2811187 6.2684335 6.1047854 5.6110906
##
## $n
## [1] 201 198 494 584 688 510 571 650 601 438 316 396 298
##
## $sd
## [1] 1.447806 2.224779 2.572577 3.141565 3.829917 4.250190 4.783688 5.363915
## [9] 5.957448 6.290097 6.294742 6.538931 6.275349
##
## $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] 4.503763
##
## $beta.ols
## [1] 29.43684
##
## $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 = geodatos, uvec = seq(0, 20000, 1500), option = "bin")
##
## attr(,"class")
## [1] "variogram"
plot(variograma, pch=16)
# no correlacion espacial con monte carlo
geodatos2=as.geodata(datos,coords.col=1:2,data.col=4)
geodatos2$data=sample(geodatos2$data)
plot(geodatos2)
#semivariograma muestral
variograma2=variog(geodatos2,option= "bin",uvec=seq(0,20000,1500))
## variog: computing omnidirectional variogram
plot(variograma2,ylim=c(0,8),type="l")
geodatos2=as.geodata(datos,coords.col=1:2,data.col=4)
geodatos2$data=sample(geodatos2$data)
variograma2=variog(geodatos2,option = "bin",uvec=seq(1000,1500))
## variog: computing omnidirectional variogram
lines(variograma2,ylim=c(0,8),type="l",add=TRUE)
## Warning in max(x$u): ningun argumento finito para max; retornando -Inf
## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a graphical
## parameter
##semivariograma con bandas de no correlacion
variograma_mc=variog.mc.env(geodatos,obj=variograma,nsim = 99)
## variog.env: generating 99 simulations by permutating data values
## variog.env: computing the empirical variogram for the 99 simulations
## variog.env: computing the envelops
##semivariograma para profundidad
geodatos_pro=as.geodata(datos,coords.col=1:2,data.col=3)
class(geodatos_pro)
## [1] "geodata"
geodatos_pro
## $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] 1.75 1.40 1.10 1.50 0.50 1.62 1.55 1.41 0.25 1.54 1.58 2.50 1.50 1.60 1.05
## [16] 1.30 1.65 1.80 1.41 1.85 1.73 1.65 1.21 1.28 1.35 1.15 0.95 1.80 1.80 1.92
## [31] 1.90 1.90 2.00 1.78 1.76 0.90 1.35 0.85 0.95 1.90 1.50 1.70 1.60 1.71 1.80
## [46] 1.63 1.40 1.80 1.25 1.10 2.00 1.90 1.90 1.60 1.75 1.67 1.54 1.87 1.45 1.40
## [61] 1.90 1.80 1.60 1.80 1.30 1.70 1.60 1.40 0.90 1.80 1.70 1.40 1.40 1.30 1.60
## [76] 1.70 1.90 1.70 1.72 1.51 1.50 1.40 1.42 1.60 1.70 1.70 1.35 1.30 1.40 1.60
## [91] 1.78 0.40 1.50 1.40 0.99 1.40 1.47 1.80 1.38 1.50 1.20 1.40 1.55 1.30 1.20
## [106] 1.15 0.60 1.10 1.50 0.90 0.90 1.70 1.30 0.80
##
## attr(,"class")
## [1] "geodata"
plot(geodatos_pro)
## Clase 19/10/2020 continuacion
plot(datos[,1:2])
geodatos=as.geodata(datos,coords.col=1:2,data.col=4)
plot(geodatos)
summary(dist(geodatos$coords))
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1373 6994 11204 11531 15606 31924
variograma=variog(geodatos,option= "bin",uvec=seq(0,20000,1500))
## variog: computing omnidirectional variogram
variograma
## $u
## [1] 1500 3000 4500 6000 7500 9000 10500 12000 13500 15000 16500 18000
## [13] 19500
##
## $v
## [1] 0.6501244 1.0682828 1.6100709 2.5685702 3.6511337 4.4781569 5.3035552
## [8] 5.8860615 6.2235025 6.2811187 6.2684335 6.1047854 5.6110906
##
## $n
## [1] 201 198 494 584 688 510 571 650 601 438 316 396 298
##
## $sd
## [1] 1.447806 2.224779 2.572577 3.141565 3.829917 4.250190 4.783688 5.363915
## [9] 5.957448 6.290097 6.294742 6.538931 6.275349
##
## $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] 4.503763
##
## $beta.ols
## [1] 29.43684
##
## $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 = geodatos, uvec = seq(0, 20000, 1500), option = "bin")
##
## attr(,"class")
## [1] "variogram"
plot(variograma, pch=16)
C1=5
a=15000
h=seq(0,20000,100)
semi=C1*(1-exp(-h^2/a^2))
ini.vals = expand.grid(seq(5,7,l=10), seq(10000,15000,l=10))
model_mco_exp=variofit(variograma, 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 "7" "10000" "0" "0.5"
## status "est" "est" "est" "fix"
## loss value: 2930.95365850542
model_mco_gaus=variofit(variograma, 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 "7" "10000" "0" "0.5"
## status "est" "est" "est" "fix"
## loss value: 1847.1908097395
model_mco_spe=variofit(variograma, 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 "5.89" "15000" "0" "0.5"
## status "est" "est" "fix" "fix"
## loss value: 1374.30780950654
lines(model_mco_exp,col="blue")
lines(model_mco_gaus,col="red")
lines(model_mco_spe,col="purple")
##Predicción Espacial Kriging
geodatos_grid=expand.grid(Este=seq(952023,979916,l=100), Norte=seq(1677586,1708649,l=100))
plot(geodatos_grid)
geodatos_ko=krige.conv(geodatos, loc=geodatos_grid,
krige= krige.control(nugget=0,trend.d="cte",
trend.l="cte",cov.pars=c(sigmasq=6.5907, phi=8804.8018)))
## krige.conv: model with constant mean
## krige.conv: Kriging performed using global neighbourhood
par(mfrow=c(1,1))
image(geodatos_ko, main="kriging Predict", xlab="East", ylab="North")
contour(geodatos_ko,main="kriging Predict", add=TRUE, drawlabels=TRUE)
image(geodatos_ko, main="kriging StDv Predicted",val=sqrt(geodatos_ko$krige.var), xlab="East", ylab="North")
contour(geodatos_ko,main="kriging StDv Predict",val=sqrt(geodatos_ko$krige.var), add=TRUE, drawlabels=TRUE)
