Trabajo Final Estadística espacial: análisis de camg - mg020
Rattia, D., Rubio, I., Ruffa, F.
2023-04-13
Análisis geoestadístico
La presente publicación corresponde al trabajo final de la matería de geoestadística para una universidad argentina. Se desarrolla según el siguiente planteamiento:
En este informe se analizan los datos de “camg” de la Librería de geoR, de acuerdo a la evolución de la variable “contenido de magnesio en la capa de suelo de 0-20 cm (mg020)”. El dataset contiene 178 observaciones y, a los fines del presente estudio, se consideran 3 de las 10 variables que lo describen (“east”: coordenadas este, en metros, “north”: coordenadas norte, en metros y “mg020”: concentración de magnesio en ese rango de profundidad del suelo).
Cabe destacar que el magnesio es un nutriente esencial para el crecimiento de las plantas, por lo que la medición de la cantidad disponible (y su interacción con otros minerales, como el calcio) en el suelo sirve para detectar carencias del mismo y, en consecuencia, es fundamental a la hora de evaluar la necesidad de incorporar fertilizantes u otras prácticas agrícolas que contribuyan a optimizar crecimiento y rendimiento.
Se realiza un análisis exploratorio y descriptivo del conjunto de datos y se implementan distintas técnicas de análisis geoestadístico. Para identificar la estructura espacial de los datos se estudia la independencia o dependencia de esos datos, la estacionariedad, la isotropía y la existencia de posibles outliers. Se finaliza con la predicción de valores no observables mediante kriging.
Cabe mencionar que los datos que se analizarán a continuación no están proyectados, por lo que se trabajará con el plano sin proyección geográfica.
Librerías necesarias
library(gstat)
library(geoR)
library(spdep)
library(sp) # para trabajar con data frame con el tipo dato "sp". Alternativa a usar GeoR.
library(dplyr) #para seleccionar variables
library(skimr)
library(sm)
library(tidyverse)
options(scipen=999)Estadísticas descriptivas
Data: El set de datos es un subset de “camg” disponible en R se analizan los datos y se observan las posibles tendencias en base al grafico sobre la variable “datosg” #Variable de estudio: Contenido de magnesio en la capa de suelo de 0-20 cm (camg de la librería geoR)
Los datos muestran 178 observaciones y 10 variables.
Descripción del dataset
Los datos presentes en camg son datos con distancias planas (en metros) y también de tipo discreto (puntuales). Cada punto significa la toma de una muestra en terreno de la concentración de Calcio o Magnesio en distintas facies de suelo. En este informe, se analizará la distribución de la variable Magnesio en la facie de 0 a 20 centímetros de profundidad. Sin embargo, antes se realiza una exploración sencilla sobre el dataset completo
skim(camg)| Name | camg |
| Number of rows | 178 |
| Number of columns | 10 |
| _______________________ | |
| Column type frequency: | |
| factor | 1 |
| numeric | 9 |
| ________________________ | |
| Group variables | None |
Variable type: factor
| skim_variable | n_missing | complete_rate | ordered | n_unique | top_counts |
|---|---|---|---|---|---|
| region | 0 | 1 | FALSE | 3 | 3: 116, 2: 48, 1: 14 |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|
| east | 0 | 1 | 5489.90 | 248.81 | 4957.0 | 5291.00 | 5501.00 | 5695.50 | 5961.0 | ▃▆▇▇▅ |
| north | 0 | 1 | 5225.74 | 196.76 | 4829.0 | 5068.75 | 5198.50 | 5363.75 | 5720.0 | ▃▇▆▅▂ |
| elevation | 0 | 1 | 5.52 | 0.67 | 3.3 | 5.20 | 5.65 | 6.00 | 6.6 | ▁▁▃▇▆ |
| ca020 | 0 | 1 | 50.68 | 11.08 | 21.0 | 43.00 | 50.50 | 58.00 | 78.0 | ▁▅▇▆▂ |
| mg020 | 0 | 1 | 27.34 | 6.28 | 11.0 | 23.00 | 27.00 | 32.00 | 46.0 | ▂▇▇▅▁ |
| ctc020 | 0 | 1 | 132.19 | 18.36 | 72.0 | 122.20 | 132.25 | 141.08 | 186.2 | ▁▂▇▃▁ |
| ca2040 | 0 | 1 | 45.01 | 13.54 | 21.0 | 34.00 | 44.00 | 56.00 | 79.0 | ▆▇▆▆▂ |
| mg2040 | 0 | 1 | 25.58 | 6.60 | 10.0 | 21.00 | 25.50 | 30.00 | 43.0 | ▂▇▇▅▂ |
| ctc2040 | 0 | 1 | 126.98 | 22.05 | 73.1 | 114.40 | 128.00 | 138.67 | 212.4 | ▂▇▇▂▁ |
El dataset consta de 178 observaciones (puntos de toma de muestras), una variable de altimetría (elevation), un factor (categoría) de regionalización y distintas mediciones de las concentraciones de los minerales calcio y magnesio a distintas profundidades. Las distrubuciones se asemejan a la normal, salvo por la medida de calcio en profundidades de 20 a 40 cm que se asemejan una distribución uniforme.
# Observación de la Región mediante la variable de estudio mg020
ggplot(camg) +
aes(x = east, y = north, colour = region, size = mg020) +
geom_point(shape = "circle") +
scale_color_hue(direction = 1) +
labs(
x = "Long (este)",
y = "Lat (norte)",
title = "Scatter plot: concentracion de mg 0 a 20 cm segun Region",
subtitle = "dataset: camg (geoR)"
) +
theme_minimal() +
theme(
plot.title = element_text(face = "bold"),
plot.subtitle = element_text(face = "italic")
)
with(camg, plot(mg020 ~ region))
Analizamos normalidad en los datos
qqnorm(camg$mg20)
qqline(camg$mg20, col=2)
#test de normalidad
shapiro.test(camg$mg20)##
## Shapiro-Wilk normality test
##
## data: camg$mg20
## W = 0.99322, p-value = 0.5817mg20 <- as.geodata(camg, data.col=6)
plot(mg20)
points(mg20)
Generar objeto espacial (Georeferencia de los datos)
coordinates(camg)=~east + north# En grilla se define un vecindario; los vecinos de cada punto,
# son todas las posiciones que estan a una distancia del punto
# mayor a 0 y menor que 100,
grilla <- dnearneigh(camg,0,100)
pixel <-coordinates(camg)
plot(grilla ,pixel, pch=15, cex=0.5, col="Blue")
Indices de Moran Global y Local, Indice de Geary
moran.test(camg$mg020, nb2listw(grilla, style = "W"),randomisation=FALSE)##
## Moran I test under normality
##
## data: camg$mg020
## weights: nb2listw(grilla, style = "W")
##
## Moran I statistic standard deviate = 10.207, p-value <
## 0.00000000000000022
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.399796322 -0.005649718 0.001577812geary.test(camg$mg020, nb2listw(grilla, style="W"),randomisation=FALSE)##
## Geary C test under normality
##
## data: camg$mg020
## weights: nb2listw(grilla, style = "W")
##
## Geary C statistic standard deviate = 9.934, p-value <
## 0.00000000000000022
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic Expectation Variance
## 0.593715071 1.000000000 0.001672695pesos <- nb2listw(grilla, style = "W")
ML <- localmoran(camg$mg020,pesos,alternative ="less")IML<-printCoefmat(data.frame(ML,row.names=mg20$.),check.names=FALSE)## Ii E.Ii Var.Ii Z.Ii Pr.z...E.Ii..
## [1,] 1.430821086 -0.012571897 0.728184839 1.691467950 0.9546
## [2,] 1.111854685 -0.007762920 0.267982458 2.162800872 0.9847
## [3,] 0.882726896 -0.004107580 0.100474681 2.797784435 0.9974
## [4,] -0.276432091 -0.038488199 1.066680609 -0.230386677 0.4089
## [5,] -0.794276900 -0.006401889 0.183345862 -1.840016723 0.0329
## [6,] 0.023794716 -0.000016385 0.000402432 1.186952361 0.8824
## [7,] 0.018673489 -0.000016385 0.000472266 0.860029203 0.8051
## [8,] -0.034239027 -0.001022576 0.029444442 -0.193576025 0.4233
## [9,] -0.111127080 -0.000787636 0.027380982 -0.666816412 0.2524
## [10,] -0.045375989 -0.000257806 0.007429044 -0.523461979 0.3003
## [11,] 0.297309204 -0.008467740 0.206220857 0.673346244 0.7496
## [12,] -0.056751941 -0.000257806 0.008966986 -0.596595694 0.2754
## [13,] -0.023312377 -0.006401889 0.221301672 -0.035947095 0.4857
## [14,] 0.885970339 -0.004107580 0.117909993 2.592105777 0.9952
## [15,] 0.721736006 -0.005791046 0.123009514 2.074339294 0.9810
## [16,] -0.150352605 -0.001022576 0.021825047 -1.010810058 0.1561
## [17,] 0.151743450 -0.001605875 0.030268284 0.881430287 0.8110
## [18,] 0.180326724 -0.005791046 0.108694798 0.564525365 0.7138
## [19,] -0.056751941 -0.000257806 0.004865806 -0.809889838 0.2090
## [20,] 0.190705744 -0.001022576 0.025090502 1.210407452 0.8869
## [21,] 0.178168981 -0.010822000 0.262930032 0.368570848 0.6438
## [22,] -0.141641398 -0.001605875 0.039379655 -0.705670994 0.2402
## [23,] -0.016458127 -0.000016385 0.000350057 -0.878777099 0.1898
## [24,] -0.056751941 -0.000257806 0.005506616 -0.761309023 0.2232
## [25,] 0.075266461 -0.001022576 0.019285249 0.549350480 0.7086
## [26,] 0.094337910 -0.007762920 0.222020350 0.216687013 0.5858
## [27,] 0.247735727 -0.001605875 0.039379655 1.256489295 0.8955
## [28,] 0.250398765 -0.007762920 0.189190273 0.593529604 0.7236
## [29,] -0.662003575 -0.006401889 0.120086178 -1.891879841 0.0293
## [30,] -0.008571438 -0.004107580 0.069091119 -0.016982408 0.4932
## [31,] 0.240820022 -0.000787636 0.013292512 2.095594121 0.9819
## [32,] 0.228524608 -0.002712523 0.041287300 1.138018519 0.8724
## [33,] 0.043193924 -0.003135415 0.052790399 0.201640901 0.5799
## [34,] 0.376406554 -0.001022576 0.021825047 2.554805364 0.9947
## [35,] -0.211044265 -0.000787636 0.016814629 -1.621459295 0.0525
## [36,] 0.173125141 -0.008467740 0.158507012 0.456115237 0.6758
## [37,] 0.348767292 -0.004624447 0.098344679 1.126888532 0.8701
## [38,] 0.020107433 -0.000787636 0.019330428 0.150287575 0.5597
## [39,] -0.113536105 -0.001022576 0.019285249 -0.810199785 0.2089
## [40,] -0.044018522 -0.000398770 0.006732445 -0.531614477 0.2975
## [41,] 0.059684275 -0.001605875 0.024470138 0.391807249 0.6524
## [42,] 0.475069039 -0.005791046 0.097243024 1.542018199 0.9385
## [43,] 0.263877834 -0.004107580 0.069091119 1.019530081 0.8460
## [44,] 0.316933745 -0.003135415 0.047703898 1.465436127 0.9286
## [45,] 0.753330562 -0.003135415 0.059007232 3.114130291 0.9991
## [46,] 1.067842178 -0.013464669 0.250773953 2.159273922 0.9846
## [47,] 0.164903865 -0.000398770 0.007525287 1.905540699 0.9716
## [48,] 0.180995898 -0.001022576 0.021825047 1.232077062 0.8910
## [49,] -0.073887566 -0.001022576 0.021825047 -0.493220723 0.3109
## [50,] 0.475126396 -0.007762920 0.267982458 0.932812623 0.8245
## [51,] 0.142279422 -0.001605875 0.046213183 0.669319497 0.7484
## [52,] 0.100572150 -0.002712523 0.066443418 0.400691285 0.6557
## [53,] -0.210384651 -0.006401889 0.107434242 -0.622331962 0.2669
## [54,] -0.171206707 -0.001022576 0.015590997 -1.362956901 0.0864
## [55,] 0.450035571 -0.010023204 0.187329298 1.062944199 0.8561
## [56,] 0.035393271 -0.000257806 0.004353159 0.540344072 0.7055
## [57,] 0.384010438 -0.001934791 0.036455831 2.021353062 0.9784
## [58,] 1.677315980 -0.016395747 0.344552693 2.885437488 0.9980
## [59,] 1.005585130 -0.006401889 0.135901099 2.745133418 0.9970
## [60,] 0.150166112 -0.000398770 0.008516341 1.631537904 0.9486
## [61,] 0.176514824 -0.006401889 0.135901099 0.496183028 0.6901
## [62,] 0.013224504 -0.000016385 0.000570033 0.554583783 0.7104
## [63,] 0.318367689 -0.004107580 0.142319430 0.854799608 0.8037
## [64,] 0.655303451 -0.005791046 0.141414150 1.757992862 0.9606
## [65,] 0.277500298 -0.004107580 0.069091119 1.071355705 0.8580
## [66,] 0.089701493 -0.000787636 0.010944433 0.864967456 0.8065
## [67,] 0.005246564 -0.000016385 0.000250067 0.332813289 0.6304
## [68,] -0.006965421 -0.002712523 0.045689621 -0.019896479 0.4921
## [69,] 0.070582246 -0.000398770 0.006732445 0.865079082 0.8065
## [70,] 0.804073044 -0.004624447 0.077744670 2.900352506 0.9981
## [71,] 1.553471056 -0.019615233 0.363048078 2.610781189 0.9955
## [72,] 0.157766013 -0.004624447 0.113058972 0.482956478 0.6854
## [73,] 0.003063990 -0.000787636 0.019330428 0.027702779 0.5111
## [74,] 0.314065858 -0.002712523 0.077973305 1.134442696 0.8717
## [75,] 0.585982523 -0.005791046 0.141414150 1.573653563 0.9422
## [76,] 0.205028792 -0.000787636 0.013292512 1.785157386 0.9629
## [77,] 0.004334054 -0.000016385 0.000227847 0.288211278 0.6134
## [78,] -0.174022450 -0.002712523 0.041287300 -0.843090680 0.1996
## [79,] -0.027212704 -0.000016385 0.000276731 -1.634862151 0.0510
## [80,] 0.868555892 -0.008467740 0.158507012 2.202860809 0.9862
## [81,] 1.259735689 -0.008467740 0.158507012 3.185405194 0.9993
## [82,] -0.054099145 -0.000016385 0.000350057 -2.890611589 0.0019
## [83,] 0.228234095 -0.000398770 0.009790553 2.310654736 0.9896
## [84,] 0.071729400 -0.010023204 0.345221135 0.139140292 0.5553
## [85,] 0.142661807 -0.010023204 0.582058891 0.200130344 0.5793
## [86,] 0.199063587 -0.001605875 0.094047882 0.654344982 0.7436
## [87,] 0.341740969 -0.002712523 0.057796002 1.432787192 0.9240
## [88,] 0.081594932 -0.002712523 0.045689621 0.394418430 0.6534
## [89,] -0.076659857 -0.000257806 0.003584188 -1.276173263 0.1009
## [90,] 0.529444178 -0.004624447 0.070253757 2.014940788 0.9780
## [91,] 0.659531535 -0.004624447 0.077744670 2.381961722 0.9914
## [92,] 1.116127837 -0.006401889 0.107434242 3.424731189 0.9997
## [93,] 1.421357058 -0.013464669 0.224353126 3.029227437 0.9988
## [94,] 1.065108126 -0.006401889 0.107434242 3.269074914 0.9995
## [95,] -0.040118196 -0.000016385 0.000350057 -2.143358778 0.0160
## [96,] 0.424766299 -0.002712523 0.118327910 1.242713472 0.8930
## [97,] -0.843700154 -0.004624447 0.407344826 -1.314679968 0.0943
## [98,] 0.093015496 -0.000787636 0.014857897 0.769553223 0.7792
## [99,] 0.009524930 -0.000063373 0.001070284 0.293084051 0.6153
## [100,] 0.883628232 -0.010822000 0.148864798 2.318249502 0.9898
## [101,] 2.017330065 -0.031004148 0.458526513 3.024952730 0.9988
## [102,] 1.943156076 -0.026919434 0.399795094 3.115761038 0.9991
## [103,] 1.322759101 -0.008467740 0.141807166 3.535110364 0.9998
## [104,] 1.495377224 -0.010822000 0.180803063 3.542252596 0.9998
## [105,] 0.977216947 -0.008467740 0.158507012 2.475789805 0.9934
## [106,] 0.738162626 -0.008467740 0.158507012 1.875346014 0.9696
## [107,] 0.015928512 -0.001022576 0.025090502 0.107014564 0.5426
## [108,] 1.036945475 -0.007762920 0.267982458 2.018096372 0.9782
## [109,] 0.159438947 -0.004107580 0.117909993 0.476284027 0.6831
## [110,] 0.087215991 -0.000787636 0.016814629 0.678667295 0.7513
## [111,] 0.113026974 -0.001022576 0.017253410 0.868272614 0.8074
## [112,] 0.952744310 -0.013464669 0.224353126 2.039881815 0.9793
## [113,] 1.230040312 -0.010822000 0.180803063 2.918237956 0.9982
## [114,] 0.561564513 -0.004624447 0.086900229 1.920661962 0.9726
## [115,] 0.311293567 -0.001934791 0.036455831 1.640505063 0.9495
## [116,] 1.252151721 -0.010822000 0.308556126 2.273668932 0.9885
## [117,] 0.070582246 -0.000398770 0.009790553 0.717362393 0.7634
## [118,] 0.164993088 -0.001605875 0.055780123 0.705395298 0.7597
## [119,] 0.079759484 -0.000257806 0.005506616 1.078304589 0.8596
## [120,] -0.056465152 -0.001605875 0.024470138 -0.350696852 0.3629
## [121,] 0.055147799 -0.000398770 0.006083755 0.712149826 0.7618
## [122,] 0.072746197 -0.000063373 0.000967159 2.341205078 0.9904
## [123,] -0.027212704 -0.000016385 0.000276731 -1.634862151 0.0510
## [124,] -0.393027633 -0.004107580 0.077227598 -1.399503772 0.0808
## [125,] -0.119788099 -0.010023204 0.167592819 -0.268124044 0.3943
## [126,] 0.063508124 -0.000398770 0.011489501 0.596207066 0.7245
## [127,] 1.770092127 -0.029641496 1.258130773 1.604520750 0.9457
## [128,] 0.335718406 -0.000257806 0.022808468 2.224643759 0.9869
## [129,] 0.072245620 -0.001022576 0.035539959 0.388648538 0.6512
## [130,] -0.011726114 -0.000016385 0.000276731 -0.703911157 0.2407
## [131,] -1.723510552 -0.025650758 0.422122663 -2.613259209 0.0045
## [132,] 0.689315846 -0.010822000 0.163382191 1.732132842 0.9584
## [133,] 0.133417651 -0.001022576 0.017253410 1.023509233 0.8470
## [134,] 0.058651835 -0.001022576 0.017253410 0.454308297 0.6752
## [135,] -0.036160511 -0.000063373 0.001070284 -1.103375250 0.1349
## [136,] 0.032152559 -0.000016385 0.000276731 1.933783327 0.9734
## [137,] 0.511376489 -0.001605875 0.046213183 2.386269512 0.9915
## [138,] 1.062106404 -0.018534510 1.067066457 1.046130038 0.8522
## [139,] 0.692985438 -0.004107580 0.178933586 1.647952918 0.9503
## [140,] 0.040264582 -0.000398770 0.009790553 0.410960020 0.6594
## [141,] 0.146982757 -0.000398770 0.006732445 1.796208121 0.9638
## [142,] 1.658244531 -0.050227092 0.728082434 2.002246270 0.9774
## [143,] 1.076063454 -0.013464669 0.224353126 2.300235926 0.9893
## [144,] 0.493308801 -0.004624447 0.064011330 1.968079789 0.9755
## [145,] 0.127653198 -0.004107580 0.069091119 0.501273836 0.6919
## [146,] -0.166400823 -0.001022576 0.019285249 -1.190873859 0.1169
## [147,] 0.578516912 -0.002712523 0.051070235 2.571956814 0.9949
## [148,] 1.646629588 -0.025650758 1.093221866 1.599392043 0.9451
## [149,] 1.436461264 -0.007762920 0.451830888 2.148555887 0.9842
## [150,] 1.891633181 -0.018534510 0.632880795 2.401101596 0.9918
## [151,] -0.566662263 -0.001022576 0.035539959 -3.000415601 0.0013
## [152,] 0.022401743 -0.001022576 0.029444442 0.136510260 0.5543
## [153,] 0.277012757 -0.001934791 0.032614948 1.544593488 0.9388
## [154,] 0.401953561 -0.001934791 0.029472409 2.352629368 0.9907
## [155,] 1.448341269 -0.023123129 0.344753886 2.506081402 0.9939
## [156,] 1.322701743 -0.013464669 0.224353126 2.820944148 0.9976
## [157,] 0.020107433 -0.000787636 0.014857897 0.171421436 0.5681
## [158,] 0.110095356 -0.000257806 0.004865806 1.582003244 0.9432
## [159,] 0.763320369 -0.010023204 0.211999898 1.679595316 0.9535
## [160,] 1.206791307 -0.018534510 0.524334380 1.692183346 0.9547
## [161,] 2.701868638 -0.012571897 2.209664340 1.826068978 0.9661
## [162,] 1.777548633 -0.015408999 0.527832427 2.467868379 0.9932
## [163,] 0.718581330 -0.007762920 0.681641426 0.879760614 0.8105
## [164,] -0.006932645 -0.000016385 0.000402432 -0.344766565 0.3651
## [165,] -0.592329852 -0.007762920 0.130096135 -1.620697816 0.0525
## [166,] 0.554911015 -0.006401889 0.107434242 1.712512164 0.9566
## [167,] 0.520094866 -0.001934791 0.032614948 2.890592207 0.9981
## [168,] -0.039042738 -0.000016385 0.000350057 -2.085877793 0.0185
## [169,] -0.129227544 -0.004107580 0.100474681 -0.394728327 0.3465
## [170,] 0.849672449 -0.007762920 0.267982458 1.656335123 0.9512
## [171,] 2.940190054 -0.012571897 2.209664340 1.986393486 0.9765
## [172,] 3.255084055 -0.021948430 1.259216880 2.920320237 0.9983
## [173,] -0.322891861 -0.005791046 0.337730265 -0.545647751 0.2927
## [174,] -0.005703551 -0.000063373 0.001826533 -0.131971171 0.4475
## [175,] 0.113887340 -0.004624447 0.098344679 0.377908026 0.6473
## [176,] 0.248964822 -0.001934791 0.041256933 1.235240489 0.8916
## [177,] 0.715283260 -0.004624447 0.113058972 2.141037663 0.9839
## [178,] 0.192324052 -0.001022576 0.029444442 1.126769115 0.8701Veamos los 10 primeros
head(IML,10)Plot de Moran Local
M<-moran.plot(camg$mg020,pesos,zero.policy=F,col=3, quiet=T,labels=T,xlab="mg20", ylab="lag(mg20)", main= "Indice de Moran Local")
Muestra los que cumplen la condición para rechazar la hipótesis nula con IML negativo
local_tabla_resultado <- IML %>%
filter(Pr.z...E.Ii..<= 0.05)
local_tabla_resultadoEn base al IML.”Les” H0=No hay correlacion negativa, si rechaza H0 hay correlacion negativa. Si los IML son negativos y el p-value es < a 0.05 (correlacion negativa), esos son outliers.
Discusión sobre el supuesto de estacionalidad
Análisis de isotropia
v_1 <- variogram(camg$mg020~1, camg, cutoff=800, width=100, map=T)
#cutoff: maxima distancia para la cual tiene sentido considerar la maxima dist entre punto. #Width:Ancho de los intervalos.
plot(v_1) # me dice que es un proc.anisotropico.
¿El proceso es isotropico? NO, es anisotropico. Se requiere definir una distancia angular.
El grafico muestra la concentracion de magnesio en zonas, en particular al norte/noroeste/noreste, sur/suroeste/sureste y, en la diagonal. Esta dispersion muestra que al variar la direccion la concentracion de magnesio varia. Por ello, a los fines del analisis, el modelo se aplica en relacion a una distancia fija y a un angulo determinado.
TEST DE INDEPENDENCIA Y DE ISOTROPIA
# Vemos la dependencia o independencia
sm.variogram(mg20, mg20$data, model = "independent")## Test of spatial independence: p = 0.019
# Vemos con parámetro isotropia
#sm.variogram(mg020, mg020$data, model = "isotropic") #OJO ... CORRER ESTO TARDA MUCHO Y CONSUME RECURSOS! el resultado es 0.023Eliminación de Outliers
data(camg)
data_camg <- camg
data_camg <- slice(data_camg, c(-5, -29, -131, -151))Datos sin outliers (variable de estudio)
mg20 <- as.geodata(data_camg, data.col=6)Análisis estructural
variograma - Variograma nube
nube_clasica <- variog(mg20, option = "cloud")## variog: computing omnidirectional variogram#tengo que tranformar los datos a geodata y luego hago el varog.nube
plot(nube_clasica, main = "classical estimator")
Paso data_camg a un objeto espacial. Lo transformo a objeto espacial para construir el variogram empirico (funcion variogram).
coordinates(data_camg)=~east + northDatos direccionales (Anisotropia)
plot(variogram(data_camg$mg020 ~ 1, data_camg, cutoff = 800, alpha = c(0, 45, 90, 135, 180, 225)))
Alpha es el angulo: se eleigio el angulo perpendicular en 135. El 0 es el norte, en el sentido de las agujas del reloj y luego hacia 90 es hacia el este, 180 es el sur. El angulo de 135 es una orientacion sureste.
Generamos un variograma empirico sin tendencia
v <- variogram(data_camg$mg020 ~ 1, data_camg, cutoff = 800, alpha = 135)plot(v, main="variograma empirico sin tendencia")
Ahora ajusto al variograma empirico, un variograma teorico
un modelo exponencial, sin tendencia
vt_exp = fit.variogram(v, vgm(50, "Exp", 500, 15))#valores inciales de sill (meseta), rango y nuggetplot(v , vt_exp, main="Varigrama teorico exponencial")
Ahora con el modelo esferico, sin tendencia
vt_sph = fit.variogram(v, vgm(50, "Sph", 500, 15)) #cambie a un modelo esferico (antes era exponencial)plot(v , vt_sph,main="Varigrama teorico esferico")
Comparo los modelos. Calculo la Suma de cuadrado del error para cada uno de los modelos ajustados
attr(vt_exp, 'SSErr')#es mejor (exp) el que me de menor ## [1] 0.5814336attr(vt_sph, 'SSErr')# el esferico me da mayor error ## [1] 0.8332345Comparo distintivos tipos de modelos teoricos
vt_sel = fit.variogram(v, vgm(50, c("Exp", "Sph" ,"Gau" ,"Mat"), 500, 15))
vt_sel # se quedo con un modelo exponencial con nugget = 0, rango=108,1759*3, sill=43.7959plot(v , vt_sel, main="Varigrama teorico segun R")
attr(vt_sel, 'SSErr')## [1] 0.5814335Prediccion (Kriging)
Vamos a quedarnos con el modelo exponencial
los usamos para predecir. #Generamos una grilla de prediccion que llamamos g2. Elegimos posiciones en el espacio donde no tenemos datos
points(mg20)#grafico con el dato de tipo geodata
g2 <- expand.grid(x=seq(5600, 5800, by=15), y=seq(4900,5100, by=15))
#View(g2)
class(g2)## [1] "data.frame"gridded(g2) = ~x+y #el df lo trabajo como objeto espacialplot(g2)#grilla de prediccion
Realizamos la prediccion sobre la grilla g2 con
modelos Exp
k1<- krige(data_camg$mg020~1, data_camg, g2, model = vt_exp, nmax=174)#nmax=maximo de obervaciones mas cercadas, son todas en este caso## [using ordinary kriging]#el nmax toma la distancia euclidia. Habitualmente se usan todas las observaciones.
#sin tendencia ~1, datos g2: la grilla que yo cree arriba.
#con el modelo exponencial vtMapa de predicciones y mapa de variaciones
spplot(k1["var1.pred"], main = "Kriging ordinario: Valores Predichos (Exp)", col.regions=terrain.colors)
spplot(k1["var1.var"], main = "Kriging ordinario: Varianza de las Predicciones (Exp)", col.regions=terrain.colors)
Mapa con superposición de grilla y curvas de interpolación
image(k1, col=terrain.colors(21), x.leg=c(5000, 5600), y.leg=c(4700, 4800))
points(mg20, add=T, col=2)
contour(k1, add=T, nlev=20,col="blue")
plot(g2, add=T, col="white")