Correlacion espacial variables cienaga

Se realiza la comparación entre las variables temperatura

Carga de datos

## GeoR - Clase 1
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)

datos <- read_excel("D:/Posgrado/1. Primer semestre/1. Clase tratamiento de datos espaciales/Clase 8/Material_Clase/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
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 los datos en geoR temperatura

geodatos <- as.geodata(datos,coords.col=1:2,data.col = 4)
class(geodatos) ##La clase es geodata
## [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 temperatura

summary(dist(geodatos$coords[1:5,]))
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    1624    1870    3282    3415    4368    6268
plot(variograma <- variog(geodatos,option = "bin"))
## variog: computing omnidirectional variogram

variograma <- variog(geodatos,option = "bin", uve=seq(0,20000,1500)) # El 1500 es el valor menor mas un poquito, el 15000 es el 3 cuartil
## variog: computing omnidirectional variogram
plot(variograma) #El variograma se debe estabilizar

geodatos2 <- as.geodata(datos,coords.col=1:2,data.col = 4)
geodatos2$data <- sample((geodatos2$data))
variograma2 <- variog(geodatos2,option = "bin", uve=seq(0,15000,1500)) # El 1500 es el valor menor mas un poquito, el 15000 es el 3 cuartil
## variog: computing omnidirectional variogram
lines(variograma2, ylim=c(0,8),type = "l", add=TRUE) 
## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a graphical
## parameter

Semivariograma con bandas de no correlación - Montecarlo - Temperatura

variograma <- variog(geodatos,option = "bin", uve=seq(0,15000,1500)) # El 1500 es el valor menor mas un poquito, el 15000 es el 3 cuartil
## variog: computing omnidirectional variogram
plot(variograma) #El variograma se debe estabilizar

variogram_mc <- variog.mc.env(geodatos ,obj=variograma)
## variog.env: generating 99 simulations by permutating data values
## variog.env: computing the empirical variogram for the 99 simulations
## variog.env: computing the envelops
variogram_mc
## $u
##  [1]  1500  3000  4500  6000  7500  9000 10500 12000 13500 15000
## 
## $v.lower
##  [1] 3.791020 3.676793 4.005769 4.019683 4.118401 4.008451 3.881208 3.976146
##  [9] 3.968136 3.113690
## 
## $v.upper
##  [1] 5.546393 5.370859 5.123543 5.076164 4.974208 5.231196 4.968643 5.173638
##  [9] 4.911755 5.573611
## 
## $call
## variog.mc.env(geodata = geodatos, obj.variog = variograma)
## 
## attr(,"class")
## [1] "variogram.envelope"

Semivariograma muestral profundidad

geodatos_prof <- as.geodata(datos,coords.col=1:2,data.col = 3)
plot(geodatos_prof)

variograma_prof <- variog(geodatos_prof,option = "bin", uve=seq(0,15000,1500)) # El 1500 es el valor menor mas un poquito, el 15000 es el 3 cuartil
## variog: computing omnidirectional variogram
plot(variograma_prof) #El variograma se debe estabilizar
variogram_mc_prof <- variog.mc.env(geodatos_prof ,obj=variograma_prof)
## variog.env: generating 99 simulations by permutating data values
## variog.env: computing the empirical variogram for the 99 simulations
## variog.env: computing the envelops
variogram_mc_prof
## $u
##  [1]  1500  3000  4500  6000  7500  9000 10500 12000 13500 15000
## 
## $v.lower
##  [1] 0.09944080 0.09986818 0.10305192 0.10779195 0.10711374 0.10872510
##  [7] 0.10545403 0.10837315 0.11060225 0.09637698
## 
## $v.upper
##  [1] 0.1520386 0.1601374 0.1464778 0.1530638 0.1467506 0.1513766 0.1447249
##  [8] 0.1545014 0.1469958 0.1671734
## 
## $call
## variog.mc.env(geodata = geodatos_prof, obj.variog = variograma_prof)
## 
## attr(,"class")
## [1] "variogram.envelope"
lines(variogram_mc_prof,col="red")

Semivariograma muestral salinidad

geodatos_sali <- as.geodata(datos,coords.col=1:2,data.col = 5)
plot(geodatos_sali) 

variograma_sali <- variog(geodatos_sali,option = "bin", uve=seq(0,15000,1500)) # El 1500 es el valor menor mas un poquito, el 15000 es el 3 cuartil
## variog: computing omnidirectional variogram
plot(variograma_sali,ylim=c(0,15)) #El variograma se debe estabilizar
variogram_mc_sali <- variog.mc.env(geodatos_sali ,obj=variograma_sali)
## variog.env: generating 99 simulations by permutating data values
## variog.env: computing the empirical variogram for the 99 simulations
## variog.env: computing the envelops
variogram_mc_sali
## $u
##  [1]  1500  3000  4500  6000  7500  9000 10500 12000 13500 15000
## 
## $v.lower
##  [1] 5.848804 5.549262 5.592008 5.878701 4.966265 5.765500 5.339856 6.026247
##  [9] 5.670300 3.110333
## 
## $v.upper
##  [1]  9.958134 11.728203 10.910269 10.990911 10.767095 11.977444 10.767600
##  [8] 11.907931 10.208169 12.523672
## 
## $call
## variog.mc.env(geodata = geodatos_sali, obj.variog = variograma_sali)
## 
## attr(,"class")
## [1] "variogram.envelope"
lines(variogram_mc_sali,col="red")

Semivariograma muestral oxigeno

geodatos_oxig <- as.geodata(datos,coords.col=1:2,data.col = 6)
plot(geodatos_oxig) 

variograma_oxig <- variog(geodatos_oxig,option = "bin", uve=seq(0,20000,1500)) # El 1500 es el valor menor mas un poquito, el 15000 es el 3 cuartil
## variog: computing omnidirectional variogram
plot(variograma_oxig, ylim=c(0,17)) #El variograma se debe estabilizar
variogram_mc_oxig <- variog.mc.env(geodatos_oxig ,obj=variograma_oxig)
## variog.env: generating 99 simulations by permutating data values
## variog.env: computing the empirical variogram for the 99 simulations
## variog.env: computing the envelops
variogram_mc_oxig
## $u
##  [1]  1500  3000  4500  6000  7500  9000 10500 12000 13500 15000 16500 18000
## [13] 19500
## 
## $v.lower
##  [1] 8.832281 8.178658 9.163914 9.059050 9.562933 9.233501 9.643862 9.445938
##  [9] 9.428919 9.224932 8.745809 8.943371 7.763669
## 
## $v.upper
##  [1] 12.26846 12.48240 11.50115 12.26018 11.83460 12.27318 11.84315 12.03275
##  [9] 12.13953 11.79883 12.07975 12.44993 14.42844
## 
## $call
## variog.mc.env(geodata = geodatos_oxig, obj.variog = variograma_oxig)
## 
## attr(,"class")
## [1] "variogram.envelope"
lines(variogram_mc_oxig,col="red")

Ajuste del Modelo de Semivarianza

plot(variograma) #El variograma se debe estabilizar

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: 2739.11345103064
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: 1218.30517924673
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: 1287.33140895057
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)

require(rasterVis)
## Loading required package: rasterVis
## Loading required package: raster
## Loading required package: sp
## Loading required package: lattice
## Loading required package: latticeExtra
require(raster)
pred=cbind(geodatos_grid,geodatos_ko$predict)
temp_predict=rasterFromXYZ(cbind(geodatos_grid,geodatos_ko$predict))
plot(temp_predict)

temp_error=rasterFromXYZ(cbind(geodatos_grid,geodatos_ko$krige.var))
plot(temp_error)