Geoestadística con sgeostat

Data Load

rm(list=ls())
aquifer=read.table("data/aquifer.txt",head=T,dec=",")
head(aquifer)

       Este     Norte Profundidad
1  42.78275 127.62282        1464
2 -27.39691  90.78732        2553
3  -1.16289  84.89600        2158
4 -18.61823  76.45199        2455
5  96.46549  64.58058        1756
6 108.56243  82.92325        1702

Libraries

library(scatterplot3d)
library(ggplot2)
library(cowplot)
library(sgeostat)

Including Plots

g1=ggplot(aquifer, aes(Profundidad, Este)) + 
  geom_point() + 
  geom_line() +
  xlab("Este") + 
  ylab("Profundidad")

g2=ggplot(aquifer, aes(Profundidad, Norte)) + 
  geom_point() + 
  geom_line() +
  xlab("Norte") + 
  ylab("Profundidad")

g3=ggplot(aquifer, aes(Profundidad, Este*Norte)) + 
  geom_point() + 
  geom_line() +
  xlab("Interacción este,norte") + 
  ylab("Profundidad")

plot_grid(g1,g2,g3)

cor(aquifer)

                  Este      Norte Profundidad
Este         1.0000000  0.1147565  -0.7788885
Norte        0.1147565  1.0000000  -0.6200923
Profundidad -0.7788885 -0.6200923   1.0000000

scatterplot3d(aquifer, highlight.3d=TRUE, col.axis="blue",
col.grid="lightblue", main="Tendencia de Profundidad", pch=20)

reg1 <- lm(Profundidad ~ Este + Norte, data = aquifer)
residuales1  <-  residuals(reg1)
summary(reg1)


Call:
lm(formula = Profundidad ~ Este + Norte, data = aquifer)

Residuals:
    Min      1Q  Median      3Q     Max 
-366.96 -161.53  -30.71  148.15  651.20 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 2591.4302    38.9599   66.52   <2e-16 ***
Este          -6.7514     0.3438  -19.64   <2e-16 ***
Norte         -5.9872     0.4066  -14.73   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 203.3 on 82 degrees of freedom
Multiple R-squared:  0.8921,    Adjusted R-squared:  0.8894 
F-statistic: 338.9 on 2 and 82 DF,  p-value: < 2.2e-16

anova(reg1)

Analysis of Variance Table

Response: Profundidad
          Df   Sum Sq  Mean Sq F value    Pr(>F)    
Este       1 19045642 19045642  460.95 < 2.2e-16 ***
Norte      1  8960172  8960172  216.86 < 2.2e-16 ***
Residuals 82  3388069    41318                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

reg2 <- lm(Profundidad ~ Este*Norte, data = aquifer)
residuales2  <-  residuals(reg2)
summary(reg2)


Call:
lm(formula = Profundidad ~ Este * Norte, data = aquifer)

Residuals:
    Min      1Q  Median      3Q     Max 
-406.30 -138.88  -13.04  129.36  722.48 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  2.627e+03  3.833e+01  68.546  < 2e-16 ***
Este        -8.287e+00  5.658e-01 -14.646  < 2e-16 ***
Norte       -6.649e+00  4.327e-01 -15.366  < 2e-16 ***
Este:Norte   2.452e-02  7.401e-03   3.314  0.00138 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 191.9 on 81 degrees of freedom
Multiple R-squared:  0.905, Adjusted R-squared:  0.9014 
F-statistic: 257.1 on 3 and 81 DF,  p-value: < 2.2e-16

anova(reg2)

Analysis of Variance Table

Response: Profundidad
           Df   Sum Sq  Mean Sq F value    Pr(>F)    
Este        1 19045642 19045642  517.06 < 2.2e-16 ***
Norte       1  8960172  8960172  243.25 < 2.2e-16 ***
Este:Norte  1   404448   404448   10.98  0.001379 ** 
Residuals  81  2983621    36835                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

reg3 <- lm(Profundidad ~ Este*Norte+I(Este^2)*I(Norte^2), data = aquifer)
residuales3  <-  residuals(reg3)
summary(reg3)


Call:
lm(formula = Profundidad ~ Este * Norte + I(Este^2) * I(Norte^2), 
    data = aquifer)

Residuals:
   Min     1Q Median     3Q    Max 
-372.7 -133.6  -20.3  129.9  505.1 

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)    
(Intercept)           2.538e+03  7.038e+01  36.055   <2e-16 ***
Este                 -7.728e+00  6.028e-01 -12.822   <2e-16 ***
Norte                -3.075e+00  1.770e+00  -1.737   0.0863 .  
I(Este^2)            -6.792e-03  5.967e-03  -1.138   0.2585    
I(Norte^2)           -2.372e-02  9.049e-03  -2.622   0.0105 *  
Este:Norte            1.155e-02  9.680e-03   1.193   0.2365    
I(Este^2):I(Norte^2)  2.251e-06  9.541e-07   2.360   0.0208 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 180.7 on 78 degrees of freedom
Multiple R-squared:  0.9189,    Adjusted R-squared:  0.9126 
F-statistic: 147.2 on 6 and 78 DF,  p-value: < 2.2e-16

anova(reg3)

Analysis of Variance Table

Response: Profundidad
                     Df   Sum Sq  Mean Sq  F value    Pr(>F)    
Este                  1 19045642 19045642 583.2335 < 2.2e-16 ***
Norte                 1  8960172  8960172 274.3868 < 2.2e-16 ***
I(Este^2)             1    55368    55368   1.6955 0.1967061    
I(Norte^2)            1   152170   152170   4.6599 0.0339500 *  
Este:Norte            1   451567   451567  13.8283 0.0003755 ***
I(Este^2):I(Norte^2)  1   181854   181854   5.5689 0.0207829 *  
Residuals            78  2547110    32655                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

aquifer=data.frame(aquifer,resi=residuales2)
aquifer_points=point(aquifer, x="Este", y="Norte")
aquifer_pair=pair(aquifer_points,num.lags=10)

....................................................................................

str(aquifer_pair)

List of 5
 $ from: num [1:3570] 1 1 1 1 1 1 1 1 1 1 ...
 $ to  : num [1:3570] 2 3 4 5 6 7 8 9 10 11 ...
 $ lags: Factor w/ 10 levels "1","2","3","4",..: 3 3 3 4 3 4 4 4 4 4 ...
 $ dist: num [1:3570] 79.3 61.3 79.9 82.8 79.5 ...
 $ bins: num [1:10] 13.6 40.7 67.8 94.9 122 ...
 - attr(*, "type")= chr "isotropic"
 - attr(*, "class")= chr "pair"

aquifer.v<-est.variogram(aquifer_points,aquifer_pair,'resi')

g4=ggplot(aquifer, aes(resi, Este)) + 
  geom_point() + 
  geom_line() +
  xlab("Este") + 
  ylab("residuales2")

g5=ggplot(aquifer, aes(resi, Norte)) + 
  geom_point() + 
  geom_line() +
  xlab("Norte") + 
  ylab("residuales2")

plot_grid(g4,g5)

aquifer_points=point(aquifer, x="Este", y="Norte")
fit.trend(aquifer_points,at="Profundidad", np=2, plot.it=TRUE)

$beta
      x^0 y^0       x^1 y^0       x^2 y^0       x^0 y^1       x^1 y^1 
 2.481430e+03 -8.373708e+00  1.416675e-03 -2.043419e+00  2.680056e-02 
      x^0 y^2 
-2.464371e-02 

$R
       x^0 y^0   x^1 y^0    x^2 y^0    x^0 y^1    x^1 y^1   x^0 y^2
[1,] -9.219544 -155.6739 -41051.636 -731.67314 -16082.944 -85540.31
[2,]  0.000000  595.1832   3500.219   57.75539  38829.771  12491.66
[3,]  0.000000    0.0000  39397.313 -117.36878   1909.315 -23722.80
[4,]  0.000000    0.0000      0.000  485.98967  14332.040  91118.22
[5,]  0.000000    0.0000      0.000    0.00000  25401.055   3240.90
[6,]  0.000000    0.0000      0.000    0.00000      0.000  19989.20

$np
[1] 2

$x
 [1]   42.78275  -27.39691   -1.16289  -18.61823   96.46549  108.56243
 [7]   88.36356   90.04213   93.17269   97.61099   90.62946   92.55262
[13]   99.48996  -24.06744  -26.06285   56.27842   73.03881   80.26679
[19]   80.23009   68.83845   76.39921   64.46148   43.39657   39.07769
[25]  112.80450   54.25899    6.13202   -3.80469   -2.23054   -2.36177
[31]   -2.18890   63.22428  -10.77860  -18.98889  -38.57884   83.14496
[37]  -21.80248  -23.56457  -20.11299  -16.62654   29.90748  100.91568
[43]  101.29544  103.26625  -14.31073  -18.13447  -18.12151   -9.88796
[49]  -12.16336   11.65754   61.69122   69.57896   66.72205  -36.65446
[55]  -19.55102  -21.29791  -22.36166   21.14719    7.68461   -8.33227
[61]   56.70724   59.00052   68.96893   70.90225   73.00243   59.66237
[67]   61.87249   63.70810    5.62706   18.24739   85.68824  105.07646
[73] -101.64278 -145.23654  -73.99313  -94.48182  -88.84983 -120.25898
[79]  -86.02454  -72.79097 -100.17372  -78.83539  -83.69063  -95.61661
[85]  -87.55480

$y
 [1] 127.62282  90.78732  84.89600  76.45199  64.58058  82.92325  56.45348
 [8]  39.25820  33.05852  56.27887  35.08169  41.75238  59.15785 184.76636
[15] 114.07479  26.84826  18.88140  12.61593  14.61795 107.77423  95.99380
[22] 110.39641  53.61499  61.99805  45.54766 147.81987  48.32772  40.40450
[29]  29.91113  33.82002  33.68207  79.49924 175.11346 171.91695 158.52742
[36] 159.11559  15.02551   9.41441  22.09269  17.25621 175.12875  22.97808
[43]  22.96385  20.34239  31.26545  30.18118  29.53241  38.14483  39.11081
[50]  18.73347  32.94906  33.80841  33.93264 150.91457 137.78404 131.82542
[57] 137.13680 139.26199 126.83751 107.77691 171.26443 164.54863 177.24820
[64] 161.38136 162.98959 170.10544 174.30177 173.91454  79.08730  77.39191
[71] 139.81702 132.03181  10.65106  28.02333  87.97270  86.62606  76.70991
[78]  80.76485  54.36334  43.09215  42.89881  40.82141  46.50482  35.82183
[85]  29.39267

$z
 [1] 1464 2553 2158 2455 1756 1702 1805 1797 1714 1466 1729 1638 1736 1476 2200
[16] 1999 1680 1806 1682 1306 1722 1437 1828 2118 1725 1606 2648 2560 2544 2386
[31] 2400 1757 1402 1364 1735 1376 2729 2766 2736 2432 1024 1611 1548 1591 2540
[46] 2352 2528 2575 2468 2646 1739 1674 1868 1865 1777 1579 1771 1408 1527 2003
[61] 1386 1089 1384 1030 1092 1161 1415 1231 2300 2238 1038 1332 3510 3490 2594
[76] 2650 2533 3571 2811 2728 3136 2553 2798 2691 2946

$residuals
 [1] -145.932017  296.391955   20.569629  155.586776  136.944207  210.578982
 [7]  112.643763   81.535500   12.407325 -165.733666   11.643984  -55.843867
[13]  123.038140  130.250727  132.838620   16.473072 -186.973641   -9.864104
[19] -133.020821 -298.072286   98.737035 -175.328351 -174.667016  118.113364
[25]  176.632628  200.333264  366.232978  173.604750  128.842139  -15.778284
[31]   -1.005758  -17.176812   -5.743382 -109.803640   35.578021  175.509274
[37]  109.375693  113.827801  154.658230 -138.758151 -234.947039  -41.999962
[43] -102.169175  -45.349545   38.415648 -182.959426   -9.456222  134.544149
[49]   14.873572  303.070200 -191.631118 -197.446346  -23.989926   92.632496
[55]  -47.092725 -308.538280  -72.511843 -213.402614 -260.643390  -17.741523
[61]  187.380986 -159.999448  282.152142 -199.908135 -116.838018  -37.190026
[67]  262.093246   81.109636  169.467368  176.796541 -289.932780   42.387375
[73]  216.381585  -51.786437   30.159248  -53.946573 -219.188525  648.160187
[79]  -92.004756 -152.583829   49.711612 -386.649271 -141.519561 -407.429504
[85] -129.126052

attr(,"class")
[1] "trend.surface"

g6=ggplot(aquifer.v, aes(resi, Norte)) + 
  geom_point() + 
  geom_line() +
  xlab("Norte") + 
  ylab("residuales2")

g6=ggplot(aquifer.v, aes(bins, classic)) + 
  geom_point() + 
  geom_line() +
  xlab("Rezago espacial, h") + 
  ylab("Estimador clásico del variograma")

g7=ggplot(aquifer.v, aes(bins, robust)) + 
  geom_point() + 
  geom_line() +
  xlab("Rezago espacial, h") + 
  ylab("Estimador robusto 1 del variograma")

g8=ggplot(aquifer.v, aes(bins, med)) + 
  geom_point() + 
  geom_line() +
  xlab("Rezago espacial, h") + 
  ylab("Estimador robusto 2 del variograma")

plot_grid(g6,g7,g8,nrow=1,ncol=3)

#par(mfrow=c(1,3))
print(aquifer.v)

   lags      bins   classic    robust       med   n
1     1  13.55308  43779.20  44355.34  47948.45 285
2     2  40.65923  71039.50  71176.29  73188.30 350
3     3  67.76539  80041.91  85367.59  93223.52 492
4     4  94.87154  67197.27  68067.40  73056.46 719
5     5 121.97770  73572.25  68052.99  66133.91 612
6     6 149.08385  57650.90  58608.95  58819.91 521
7     7 176.19001  65498.82  62167.57  68112.31 356
8     8 203.29616 130414.72 107613.55  77805.71 173
9     9 230.40231 161738.13 134102.60 123952.77  43
10   10 257.50847  35525.99  45217.14  58333.98  19

plot(aquifer.v$robust)

plot(aquifer.v$med,col='blue')
points(aquifer.v$robust,col="red")

aquifer.vmodExp<-fit.exponential(aquifer.v,c0=0,ce=40000,ae=20,plot.it=TRUE,iterations=30)

Initial parameter estimates:  0 40000 20

Iteration: 1 
Gradient vector:  -4432.441 977.0988 -8.943538 
New parameter estimates:  1e-06 40977.1 11.05646 

rse.dif =  3232643827 (rse = 3232643827 )  ;  parm.dist =  977.1397

Iteration: 2 
Gradient vector:  -26700.7 22493.46 -2.800242 
New parameter estimates:  1e-06 63470.56 8.256219 

rse.dif =  -17644208 (rse = 3.215e+09 )  ;  parm.dist =  22493.46

Iteration: 3 
Gradient vector:  -11057.27 -15597.73 2.315183 
New parameter estimates:  1e-06 47872.83 10.5714 

rse.dif =  -3772568 (rse = 3211227051 )  ;  parm.dist =  15597.73

Iteration: 4 
Gradient vector:  -27525.12 16431.58 -1.824505 
New parameter estimates:  1e-06 64304.41 8.746897 

rse.dif =  3032851 (rse = 3214259902 )  ;  parm.dist =  16431.58

Iteration: 5 
Gradient vector:  -20442.22 -7053.019 1.144197 
New parameter estimates:  1e-06 57251.39 9.891094 

rse.dif =  -2468665 (rse = 3211791237 )  ;  parm.dist =  7053.019

Iteration: 6 
Gradient vector:  -27557.41 7097.539 -0.7122805 
New parameter estimates:  1e-06 64348.93 9.178813 

rse.dif =  1486180 (rse = 3213277417 )  ;  parm.dist =  7097.539

Iteration: 7 
Gradient vector:  -24787.06 -2758.919 0.3605893 
New parameter estimates:  1e-06 61590.01 9.539403 

rse.dif =  -951749.7 (rse = 3212325667 )  ;  parm.dist =  2758.919

Iteration: 8 
Gradient vector:  -26691.4 1898.737 -0.1885371 
New parameter estimates:  1e-06 63488.75 9.350866 

rse.dif =  471370.4 (rse = 3212797038 )  ;  parm.dist =  1898.737

Iteration: 9 
Gradient vector:  -25850.35 -838.0686 0.09276125 
New parameter estimates:  1e-06 62650.68 9.443627 

rse.dif =  -249219.6 (rse = 3212547818 )  ;  parm.dist =  838.0686

Iteration: 10 
Gradient vector:  -26302.53 450.7265 -0.04631475 
New parameter estimates:  1e-06 63101.41 9.397312 

rse.dif =  121873.4 (rse = 3212669692 )  ;  parm.dist =  450.7265

Iteration: 11 
Gradient vector:  -26086.54 -215.2624 0.02285916 
New parameter estimates:  1e-06 62886.14 9.420171 

rse.dif =  -61031.79 (rse = 3212608660 )  ;  parm.dist =  215.2624

Iteration: 12 
Gradient vector:  -26195.52 108.6221 -0.01133309 
New parameter estimates:  1e-06 62994.77 9.408838 

rse.dif =  30077.83 (rse = 3212638738 )  ;  parm.dist =  108.6221

Iteration: 13 
Gradient vector:  -26142.08 -53.26613 0.005604603 
New parameter estimates:  1e-06 62941.5 9.414443 

rse.dif =  -14922.96 (rse = 3212623815 )  ;  parm.dist =  53.26613

Iteration: 14 
Gradient vector:  -26168.65 26.48517 -0.002774911 
New parameter estimates:  1e-06 62967.99 9.411668 

rse.dif =  7377.216 (rse = 3212631192 )  ;  parm.dist =  26.48517

Iteration: 15 
Gradient vector:  -26155.53 -13.07801 0.001373075 
New parameter estimates:  1e-06 62954.91 9.413041 

rse.dif =  -3653.216 (rse = 3212627539 )  ;  parm.dist =  13.07801

Iteration: 16 
Gradient vector:  -26162.03 6.479831 -0.0006796194 
New parameter estimates:  1e-06 62961.39 9.412361 

rse.dif =  1807.514 (rse = 3212629346 )  ;  parm.dist =  6.479831

Iteration: 17 
Gradient vector:  -26158.82 -3.20516 0.0003363367 
New parameter estimates:  1e-06 62958.18 9.412698 

rse.dif =  -894.6895 (rse = 3212628451 )  ;  parm.dist =  3.20516

Iteration: 18 
Gradient vector:  -26160.41 1.586717 -0.0001664615 
New parameter estimates:  1e-06 62959.77 9.412531 

rse.dif =  442.763 (rse = 3212628894 )  ;  parm.dist =  1.586717

Iteration: 19 
Gradient vector:  -26159.62 -0.7851797 8.238305e-05 
New parameter estimates:  1e-06 62958.98 9.412613 

rse.dif =  -219.1369 (rse = 3212628675 )  ;  parm.dist =  0.7851797

Iteration: 20 
Gradient vector:  -26160.01 0.3886224 -4.077271e-05 
New parameter estimates:  1e-06 62959.37 9.412573 

rse.dif =  108.4519 (rse = 3212628784 )  ;  parm.dist =  0.3886224

Iteration: 21 
Gradient vector:  -26159.82 -0.192328 2.01789e-05 
New parameter estimates:  1e-06 62959.18 9.412593 

rse.dif =  -53.67477 (rse = 3212628730 )  ;  parm.dist =  0.192328

Iteration: 22 
Gradient vector:  -26159.91 0.09518727 -9.986825e-06 
New parameter estimates:  1e-06 62959.28 9.412583 

rse.dif =  26.56425 (rse = 3212628756 )  ;  parm.dist =  0.09518727

Iteration: 23 
Gradient vector:  -26159.86 -0.04710907 4.94261e-06 
New parameter estimates:  1e-06 62959.23 9.412588 

rse.dif =  -13.14703 (rse = 3212628743 )  ;  parm.dist =  0.04710907

Iteration: 24 
Gradient vector:  -26159.89 0.023315 -2.446165e-06 
New parameter estimates:  1e-06 62959.25 9.412585 

rse.dif =  6.506634 (rse = 3212628750 )  ;  parm.dist =  0.023315

Iteration: 25 
Gradient vector:  -26159.88 -0.01153889 1.21064e-06 
New parameter estimates:  1e-06 62959.24 9.412587 

rse.dif =  -3.220222 (rse = 3212628747 )  ;  parm.dist =  0.01153889

Iteration: 26 
Gradient vector:  -26159.88 0.005710764 -5.991627e-07 
New parameter estimates:  1e-06 62959.25 9.412586 

rse.dif =  1.593734 (rse = 3212628748 )  ;  parm.dist =  0.005710764

Iteration: 27 
Gradient vector:  -26159.88 -0.002826342 2.965346e-07 
New parameter estimates:  1e-06 62959.24 9.412586 

rse.dif =  -0.7887611 (rse = 3212628747 )  ;  parm.dist =  0.002826342

Iteration: 28 
Gradient vector:  -26159.88 0.001398795 -1.467591e-07 
New parameter estimates:  1e-06 62959.24 9.412586 

rse.dif =  0.3903689 (rse = 3212628748 )  ;  parm.dist =  0.001398795

Iteration: 29 
Gradient vector:  -26159.88 -0.0006922812 7.263288e-08 
New parameter estimates:  1e-06 62959.24 9.412586 

rse.dif =  -0.1932006 (rse = 3212628748 )  ;  parm.dist =  0.0006922812

Iteration: 30 
Gradient vector:  -26159.88 0.000342624 -3.594748e-08 
New parameter estimates:  1e-06 62959.24 9.412586 

rse.dif =  0.09561825 (rse = 3212628748 )  ;  parm.dist =  0.000342624

Convergence not achieved!

aquifer.vmodGau<-fit.gaussian(aquifer.v,c0=0,cg=50000,ag=50,plot.it=TRUE,iterations=30)

Initial parameter estimates:  0 50000 50

Iteration: 1 
Gradient vector:  19162.34 -33401.14 -11.41191 
New parameter estimates:  19162.34 16598.86 38.58809 

rse.dif =  3299750048 (rse = 3299750048 )  ;  parm.dist =  38507.55

Iteration: 2 
Gradient vector:  -1294.927 2010.017 -18.77473 
New parameter estimates:  17867.41 18608.87 19.81336 

rse.dif =  -66430135 (rse = 3233319913 )  ;  parm.dist =  2391.1

Iteration: 3 
Gradient vector:  3201.043 -2835.169 9.216254 
New parameter estimates:  21068.46 15773.71 29.02961 

rse.dif =  -24694350 (rse = 3208625564 )  ;  parm.dist =  4276.09

Iteration: 4 
Gradient vector:  -4345.272 4292.413 -6.361973 
New parameter estimates:  16723.18 20066.12 22.66764 

rse.dif =  4004881 (rse = 3212630445 )  ;  parm.dist =  6107.884

Iteration: 5 
Gradient vector:  53.88685 -4.270081 2.074271 
New parameter estimates:  16777.07 20061.85 24.74191 

rse.dif =  -3703977 (rse = 3208926468 )  ;  parm.dist =  54.09555

Iteration: 6 
Gradient vector:  -391.4471 384.4526 -0.5571294 
New parameter estimates:  16385.62 20446.3 24.18478 

rse.dif =  588163 (rse = 3209514631 )  ;  parm.dist =  548.6666

Iteration: 7 
Gradient vector:  29.55911 -27.0943 0.07968918 
New parameter estimates:  16415.18 20419.21 24.26447 

rse.dif =  -201438.9 (rse = 3209313192 )  ;  parm.dist =  40.09799

Iteration: 8 
Gradient vector:  -6.581211 6.259206 -0.01207028 
New parameter estimates:  16408.6 20425.47 24.2524 

rse.dif =  26607.8 (rse = 3209339800 )  ;  parm.dist =  9.082408

Iteration: 9 
Gradient vector:  0.9423146 -0.8928955 0.001794561 
New parameter estimates:  16409.54 20424.57 24.25419 

rse.dif =  -4077.43 (rse = 3209335722 )  ;  parm.dist =  1.298161

Iteration: 10 
Gradient vector:  -0.1413215 0.1339887 -0.0002673761 
New parameter estimates:  16409.4 20424.71 24.25393 

rse.dif =  605.1536 (rse = 3209336327 )  ;  parm.dist =  0.194743

Iteration: 11 
Gradient vector:  0.02102884 -0.01993597 3.982407e-05 
New parameter estimates:  16409.42 20424.69 24.25397 

rse.dif =  -90.18701 (rse = 3209336237 )  ;  parm.dist =  0.02897682

Iteration: 12 
Gradient vector:  -0.003132718 0.00296995 -5.931842e-06 
New parameter estimates:  16409.42 20424.69 24.25396 

rse.dif =  13.43229 (rse = 3209336251 )  ;  parm.dist =  0.004316777

Iteration: 13 
Gradient vector:  0.0004666088 -0.0004423641 8.835486e-07 
New parameter estimates:  16409.42 20424.69 24.25396 

rse.dif =  -2.000767 (rse = 3209336249 )  ;  parm.dist =  0.0006429701

Iteration: 14 
Gradient vector:  -6.950171e-05 6.589045e-05 -1.316048e-07 
New parameter estimates:  16409.42 20424.69 24.25396 

rse.dif =  0.2980142 (rse = 3209336249 )  ;  parm.dist =  9.577086e-05

Iteration: 15 
Gradient vector:  1.035229e-05 -9.814388e-06 1.960254e-08 
New parameter estimates:  16409.42 20424.69 24.25396 

rse.dif =  -0.04438925 (rse = 3209336249 )  ;  parm.dist =  1.426508e-05

Iteration: 16 
Gradient vector:  -1.542e-06 1.46188e-06 -2.919839e-09 
New parameter estimates:  16409.42 20424.69 24.25396 

rse.dif =  0.006612301 (rse = 3209336249 )  ;  parm.dist =  2.124821e-06

Iteration: 17 
Gradient vector:  2.296994e-07 -2.177628e-07 4.349363e-10 
New parameter estimates:  16409.42 20424.69 24.25396 

rse.dif =  -0.0009851456 (rse = 3209336249 )  ;  parm.dist =  3.165152e-07

Iteration: 18 
Gradient vector:  -3.420782e-08 3.242781e-08 -6.477718e-11 
New parameter estimates:  16409.42 20424.69 24.25396 

rse.dif =  0.0001473427 (rse = 3209336249 )  ;  parm.dist =  4.713621e-08

Iteration: 19 
Gradient vector:  5.086605e-09 -4.821087e-09 9.637383e-12 
New parameter estimates:  16409.42 20424.69 24.25396 

rse.dif =  -2.193451e-05 (rse = 3209336249 )  ;  parm.dist =  7.007276e-09

Iteration: 20 
Gradient vector:  -7.583935e-10 7.161523e-10 -1.439258e-12 
New parameter estimates:  16409.42 20424.69 24.25396 

rse.dif =  2.861023e-06 (rse = 3209336249 )  ;  parm.dist =  1.042223e-09

Iteration: 21 
Gradient vector:  1.019258e-10 -9.775917e-11 1.99309e-13 
New parameter estimates:  16409.42 20424.69 24.25396 

rse.dif =  -4.768372e-07 (rse = 3209336249 )  ;  parm.dist =  1.415077e-10 

Convergence achieved by sums of squares.

Final parameter estimates:  16409.42 20424.69 24.25396

aquifer.vmodWave<-fit.wave(aquifer.v,c0=0,cw=40000,aw=10,plot.it=TRUE,iterations=30,weighted=T)

Initial parameter estimates:  0 40000 10

Iteration: 1 
Gradient vector:  18650.32 -21981.27 -0.7942028 
New parameter estimates:  18650.32 18018.73 9.205797 

rse.dif =  3409704989 (rse = 3409704989 )  ;  parm.dist =  28827.26

Iteration: 2 
Gradient vector:  812.9227 -1109.399 -1.187299 
New parameter estimates:  19463.25 16909.33 8.018498 

rse.dif =  -289093760 (rse = 3120611230 )  ;  parm.dist =  1375.358

Iteration: 3 
Gradient vector:  -6990.158 6973.566 0.9858099 
New parameter estimates:  12473.09 23882.9 9.004308 

rse.dif =  24044562 (rse = 3144655792 )  ;  parm.dist =  9873.851

Iteration: 4 
Gradient vector:  7025.438 -6960.473 -1.260353 
New parameter estimates:  19498.53 16922.43 7.743955 

rse.dif =  -56767551 (rse = 3087888241 )  ;  parm.dist =  9889.639

Iteration: 5 
Gradient vector:  -9210.154 9213.61 1.066674 
New parameter estimates:  10288.37 26136.04 8.810629 

rse.dif =  175986924 (rse = 3263875165 )  ;  parm.dist =  13027.57

Iteration: 6 
Gradient vector:  11994.7 -11983.26 -2.255679 
New parameter estimates:  22283.07 14152.77 6.55495 

rse.dif =  -196728543 (rse = 3067146622 )  ;  parm.dist =  16954.98

Iteration: 7 
Gradient vector:  -14060.45 14195.04 -1.578095 
New parameter estimates:  8222.625 28347.81 4.976855 

rse.dif =  147278852 (rse = 3214425474 )  ;  parm.dist =  19979.87

Iteration: 8 
Gradient vector:  -15826.64 16212.91 0.3854677 
New parameter estimates:  1e-06 44560.72 5.362323 

rse.dif =  -46983778 (rse = 3167441696 )  ;  parm.dist =  18178.84

Iteration: 9 
Gradient vector:  13145.08 -21444.98 -0.8756698 
New parameter estimates:  13145.08 23115.75 4.486653 

rse.dif =  -757940879 (rse = 2409500817 )  ;  parm.dist =  25153.13

Iteration: 10 
Gradient vector:  -9434763 9682459 25.73116 
New parameter estimates:  1e-06 9705575 30.21781 

rse.dif =  1636307005 (rse = 4045807822 )  ;  parm.dist =  9682468

Iteration: 11 
Gradient vector:  20962.2 -9688482 0.02156687 
New parameter estimates:  20962.2 17093.21 30.23938 

rse.dif =  83628062 (rse = 4129435883 )  ;  parm.dist =  9688504

Iteration: 12 
Gradient vector:  7173.136 -8587.116 1.22582 
New parameter estimates:  28135.34 8506.099 31.4652 

rse.dif =  -628497356 (rse = 3500938527 )  ;  parm.dist =  11188.94

Iteration: 13 
Gradient vector:  2974.651 -2890.861 -4.19572 
New parameter estimates:  31109.99 5615.237 27.26947 

rse.dif =  -192443200 (rse = 3308495327 )  ;  parm.dist =  4147.969

Iteration: 14 
Gradient vector:  -2399.351 1443.698 15.69929 
New parameter estimates:  28710.64 7058.936 42.96876 

rse.dif =  147479203 (rse = 3455974530 )  ;  parm.dist =  2800.25

Iteration: 15 
Gradient vector:  4786.661 2165.107 -43.14322 
New parameter estimates:  33497.3 9224.042 1e-06 

rse.dif =  -686128323 (rse = 2769846206 )  ;  parm.dist =  5253.728

Iteration: 16

Warning in lsfit(xmat, y, wt = w, intercept = FALSE): 'X' matrix was collinear

Gradient vector:  -7188.309 -5.926894e-07 0 
New parameter estimates:  26308.99 9224.042 1e-06 

rse.dif =  686457465 (rse = 3456303671 )  ;  parm.dist =  7188.309

Iteration: 17

Warning in lsfit(xmat, y, wt = w, intercept = FALSE): 'X' matrix was collinear

Gradient vector:  -5.339328e-06 -5.926894e-07 0 
New parameter estimates:  26308.99 9224.042 1e-06 

rse.dif =  0.4889331 (rse = 3456303672 )  ;  parm.dist =  5.372122e-06

Iteration: 18

Warning in lsfit(xmat, y, wt = w, intercept = FALSE): 'X' matrix was collinear

Gradient vector:  5.926852e-07 -5.926894e-07 0 
New parameter estimates:  26308.99 9224.042 1e-06 

rse.dif =  1.907349e-06 (rse = 3456303672 )  ;  parm.dist =  8.381857e-07

Iteration: 19

Warning in lsfit(xmat, y, wt = w, intercept = FALSE): 'X' matrix was collinear

Gradient vector:  5.926931e-07 -5.926894e-07 0 
New parameter estimates:  26308.99 9224.042 1e-06 

rse.dif =  -9.536743e-07 (rse = 3456303672 )  ;  parm.dist =  8.381908e-07 

Convergence achieved by sums of squares.

Final parameter estimates:  26308.99 9224.042 1e-06

curve(65000*(1-(14/x)*sin(x/14)),0,300,ylim=c(0,200000))
points(aquifer.v$bins,aquifer.v$classic,col=3)
text(aquifer.v$bins,aquifer.v$classic,aquifer.v$n,col=2)

curve(200000*(1-exp(-x/170)),0,300)
points(aquifer.v$bins,aquifer.v$classic,col=2)

curve(65000*(1-(14/x)*sin(x/14)),0,300,ylim=c(0,200000))
points(aquifer.v$bins,aquifer.v$classic,col=3)
text(aquifer.v$bins,aquifer.v$classic,aquifer.v$n,col=2)

aquifer.vmodExp<-fit.exponential(aquifer.v,c0=0,ce=200000,ae=170,plot.it=TRUE,iterations=30,weighted=T)

Initial parameter estimates:  0 2e+05 170

Iteration: 1 
Gradient vector:  16365.66 -238859.4 -103.7436 
New parameter estimates:  16365.66 1e-06 66.25643 

rse.dif =  3826411368 (rse = 3826411368 )  ;  parm.dist =  200668.5

Iteration: 2 
Gradient vector:  7737.246 16547.95 166070861252 
New parameter estimates:  24102.91 16547.95 166070861318 

rse.dif =  -767474321 (rse = 3058937047 )  ;  parm.dist =  166070861252

Iteration: 3

Warning in lsfit(xmat, y, wt = w, intercept = FALSE): 'X' matrix was collinear

Gradient vector:  3355.03 1.242479e+13 0 
New parameter estimates:  27457.94 1.242479e+13 166070861318 

rse.dif =  -120011141 (rse = 2938925906 )  ;  parm.dist =  1.242479e+13

Iteration: 4

Warning in lsfit(xmat, y, wt = w, intercept = FALSE): 'X' matrix was collinear

Gradient vector:  423.3165 -663474885968 0 
New parameter estimates:  27881.25 1.176131e+13 166070861318 

rse.dif =  11801483 (rse = 2950727388 )  ;  parm.dist =  663474885968

Iteration: 5

Warning in lsfit(xmat, y, wt = w, intercept = FALSE): 'X' matrix was collinear

Gradient vector:  3.873181 -6320523090 0 
New parameter estimates:  27885.12 1.175499e+13 166070861318 

rse.dif =  128956.4 (rse = 2950856345 )  ;  parm.dist =  6320523090

Iteration: 6

Warning in lsfit(xmat, y, wt = w, intercept = FALSE): 'X' matrix was collinear

Gradient vector:  0.02266712 -36921321 0 
New parameter estimates:  27885.15 1.175495e+13 166070861318 

rse.dif =  752.3639 (rse = 2950857097 )  ;  parm.dist =  36921321

Iteration: 7

Warning in lsfit(xmat, y, wt = w, intercept = FALSE): 'X' matrix was collinear

Gradient vector:  0.0001316946 -214507.3 0 
New parameter estimates:  27885.15 1.175495e+13 166070861318 

rse.dif =  4.371067 (rse = 2950857102 )  ;  parm.dist =  214507.3

Iteration: 8

Warning in lsfit(xmat, y, wt = w, intercept = FALSE): 'X' matrix was collinear

Gradient vector:  7.650992e-07 -1246.212 0 
New parameter estimates:  27885.15 1.175495e+13 166070861318 

rse.dif =  0.02539396 (rse = 2950857102 )  ;  parm.dist =  1246.213

Iteration: 9

Warning in lsfit(xmat, y, wt = w, intercept = FALSE): 'X' matrix was collinear

Gradient vector:  4.447233e-09 -7.24441 0 
New parameter estimates:  27885.15 1.175495e+13 166070861318 

rse.dif =  0.0001473427 (rse = 2950857102 )  ;  parm.dist =  7.244141

Iteration: 10

Warning in lsfit(xmat, y, wt = w, intercept = FALSE): 'X' matrix was collinear

Gradient vector:  2.31966e-11 -0.03646515 0 
New parameter estimates:  27885.15 1.175495e+13 166070861318 

rse.dif =  9.536743e-07 (rse = 2950857102 )  ;  parm.dist =  0.03710938 

Convergence achieved by sums of squares.

Final parameter estimates:  27885.15 1.175495e+13 166070861318

aquifer.vmodwave<-fit.wave(aquifer.v,c0=4000,cw=30000,aw=15,plot.it=TRUE,iterations=0,weighted=T)

Convergence not achieved!

aquifer.vmodExp_0<-fit.exponential(aquifer.v,c0=0,ce=200000,ae=170,plot.it=TRUE,iterations=0,weighted=T)

Convergence not achieved!

aquifer.vmodwave_0<-fit.wave(aquifer.v,c0=4000,cw=30000,aw=15,plot.it=TRUE,iterations=0,weighted=T)

Convergence not achieved!

aquifer.spherical<-fit.spherical(aquifer.v,c0=0,cs=35000,as=70,plot.it=TRUE,iterations=0,weighted=T)

Convergence not achieved!

ggplot(aquifer.v, aes(bins, classic)) + 
  geom_point() + 
  geom_line() +
  xlab("Rezago espacial, h") + 
  ylab("Estimador clásico del variograma")+
  xlim(0, 300) +
  geom_function(aes(color = "Exponencial"),
    fun =~4000+150000*(1-exp(-.x/100)) 
    ) +
  geom_function(aes(color = "Seno cardinal"),
    fun =~4000+30000*(1-((15/.x)*sin(.x/15)))             
    ) + xlab("Rezago espacial") + ylab("Modelos teóricos de  semivariogramas")

Warning: Removed 1 row containing missing values or values outside the scale range
(`geom_function()`).

Kriging_aquifer <- point(data.frame(list(x=10,y=80)))
Kriging_aquifer <- krige(Kriging_aquifer, aquifer_points, 'resi', aquifer.vmodExp_0)


Using all points.
  Preparing the kriging system matrix...
  Inverting the matrix...
  Predicting.

Kriging_aquifer


Point object: x 

   Locations: 1

   Attributes:
      x
      y
      do
      zhat
      sigma2hat

Kriging_aquifer$sigma2hat

[1] 7010.452

Kriging_aquifer <- point(data.frame(list(x=10,y=80)))
Kriging_aquifer <- krige(Kriging_aquifer, aquifer_points, 'resi', aquifer.vmodwave_0)


Using all points.
  Preparing the kriging system matrix...
  Inverting the matrix...
  Predicting.

Kriging_aquifer


Point object: x 

   Locations: 1

   Attributes:
      x
      y
      do
      zhat
      sigma2hat

Kriging_aquifer$zhat

[1] 196.2781

Kriging_aquifer$sigma2hat

[1] 5169.927

grid <- list(x=seq(min(aquifer$Este),max(aquifer$Este),by=20),y=seq(min(aquifer$Norte),max(aquifer$Norte),by=10))
grid$xr <- range(grid$x)
grid$xs <- grid$xr[2] - grid$xr[1]
grid$yr <- range(grid$y)
grid$ys <- grid$yr[2] - grid$yr[1]
grid$max <- max(grid$xs, grid$ys)
grid$xy <- data.frame(cbind(c(matrix(grid$x, length(grid$x), length(grid$y))),
c(matrix(grid$y, length(grid$x), length(grid$y), byrow=TRUE))))
colnames(grid$xy) <- c("x", "y")
grid$point <- point(grid$xy)
grid$krige <- krige(grid$point,aquifer_points,'resi',aquifer.vmodwave_0,maxdist=180,extrap=FALSE)


Using points within 180 units of prediction points.
  Predicting..........................................................................................................................................................................................................................................

op <- par(no.readonly = TRUE)
par(pty="s")
plot(grid$xy, type="n", xlim=c(grid$xr[1], grid$xr[1]+grid$max),ylim=c(grid$yr[1], grid$yr[1]+grid$max))
image(grid$x,grid$y,matrix(grid$krige$zhat,length(grid$x),length(grid$y)),add=TRUE)
contour(grid$x,grid$y,matrix(grid$krige$zhat,length(grid$x),length(grid$y)),add=TRUE)

x11()
op <- par(no.readonly = TRUE)
par(pty="s")
plot(grid$xy, type="n", xlim=c(grid$xr[1], grid$xr[1]+grid$max),ylim=c(grid$yr[1], grid$yr[1]+grid$max))
image(grid$x,grid$y,matrix(grid$krige$sigma2hat,length(grid$x),length(grid$y)), add=TRUE)
contour(grid$x,grid$y,matrix(grid$krige$sigma2hat,length(grid$x),length(grid$y)),add=TRUE)