Geoestadistica Con R

Geoestad stica con R
Jorge Gaspar Sanz Salinas Septiembre de 2005

Resumen: A lo largo de la asignatura de doctorado Prediccion y analisis de modelos superciales mediante sistemas de informacion geograca se ha cubierto el desarrollo del estudio de la distribucion espacial de una o varias variables, as como su modelizacion mediante m etodos geoestad siticos (krigeado ). En este trabajo se presenta un resumen de dicho desarrollo utilizando los mismos datos de partida pero empleando para el mismo herramientas de Software Libre, principalmente una herramienta estad stica R y un Sistema de Informacion Geograca, GRASS, ambos funcionando bajo el Sistema Operativo Linux.
Indice
0. Introduccion 0.1. R . . . . . . . . . 0.2. gstat . . . . . . . 0.3. GRASS . . . . . . 0.4. Datos de trabajo 4 4 4 5 5
3.5. Diagramas de dispersion cruzados . . . . . . . . . . . 20 3.6. Scripts . . . . . . . . . . . . . 20 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4. Estimacion. M etodos deterministas 24 4.1. Scripts . . . . . . . . . . . . . 27 5. Continuidad espacial de V 5.1. Variograma omnidireccional 5.2. Variograma supercial . . . 5.3. Variogramas direccionales . 5.4. Variogramas cruzados . . . . 5.5. Scripts . . . . . . . . . . . . . 30 30 30 30 33 36
univariada 1. Descripcion 6 1.1. Carga y visualizacion de los datos . . . . . . . . . . . . . . 6 1.2. M etodos gracos para la descripcion univariada . . . 6 1.3. M etodos num ericos . . . . . 9 1.4. Scripts . . . . . . . . . . . . . 10 bivariada 2. Descripcion 12 2.1. M etodos gracos . . . . . . . 12 2.2. M etodos num ericos . . . . . 13 2.3. Scripts . . . . . . . . . . . . . 14 espacial 3. Descripcion 3.1. Visualizacion espacial de datos . . . . . . . . . . . . . . 3.2. Ventanas moviles y el efecto proporcional . . . . . . . . . 3.3. Continuidad espacial . . . . 3.4. Variograma . . . . . . . . . . 15 15 18 18 20
del variograma ex6. Modelizacion perimental 39 6.1. Estimacion automatizada del modelo . . . . . . . . . . 41 6.2. Scripts . . . . . . . . . . . . . 44 8. Kriging 8.1. wlc . . . . . . . . . . . . . . . 8.2. Krigeado Ordinario (KO) . . . 8.3. Krigeado Universal (KU) . . . 8.4. Krigeado por bloques (KUB) 8.5. Krigeado Local (KUL) . . . . 8.6. Cokrigeado (CKO) . . . . . . 8.7. Resultados . . . . . . . . . . 8.8. Scripts . . . . . . . . . . . . . 45 45 45 46 46 46 48 51 56
Modelizacion geoestad stica con R
Indice de guras
1. 2. 3. 4. 5. 6. 7. 8. 9. 10. 11. 12. 13. 14. 15. 16. 17. 18. 19. 20. 21. 22. 23. 24. 25. 26. Distribucion de U y V . . . . Histograma de V . . . . . . . Histograma acumulado de V Graco de probabilidad uniforme de V . . . . . . . . Graco de probabilidad normal . . . . . . . . . . . . . Graco de probabilidad lognormal . . . . . . . . . . . . . Graco de caja y bigotes de V Graco de caja y bigotes de VyU. . . . . . . . . . . . . . Graco de cuantiles de V y U Graco de dispersion . . . . Distribucion de V . . . . . . Mapa graduado de color de V Mapa graduado de tamano de V . . . . . . . . . . . . . . Mapas de indicadores . . . . Mapa de supercie interpolada . . . . . . . . . . . . . . Media y varianza en ventana de 3x3 . . . . . . . . . . . Graco de dispersion de media y varianza . . . . . . . h-Scatterplots de direccion N-S . . . . . . . . . . . . . . . h-Scatterplots de direccion E-W . . . . . . . . . . . . . . h-Scatterplots cruzado de U y V en direccion N-S . . . . . Distribucion de wlm . . . . . Mapas generados por GRASS Histogramas de los conjuntos de datos . . . . . . . . . . Variogramas omnidireccionales (i)) . . . . . . . . . . . . Mapa del variograma supercial . . . . . . . . . . . . Isol neas del variograma supercial . . . . . . . . . . . 6 7 7 7 8 8 9 12 12 13 15 16 16 17 18 19 19 20 21 21 24 25 26 31 31 32
27. Variogramas direccionales . 28. Deteccion de ejes de anisotrop a . . . . . . . . . . . . . 29. Variogramas por tolerancias (i) . . . . . . . . . . . . . 30. Variogramas por tolerancias (ii) . . . . . . . . . . . . . 31. Variogramas cruzados . . . . 32. Modelos de variogramas disponibles . . . . . . . . . . 33. Modelos esf erico de rango 30 y meseta parcial de 92000ppm . . . . . . . . . . 34. Modelos de variograma combinados . . . . . . . . . . 35. Modelo de variograma de V ajustado . . . . . . . . . . . . 36. Modelado interactivo del variograma . . . . . . . . . . . . 37. Descripcion de geoR del conjunto de datos . . . . . . 38. Modelo ajustado por geoR . 39. Conjunto de datos wlc . . . . 40. Error del Krigeado ordinario 41. Error del Krigeado universal 42. Error del Krigeado universal por bloques . . . . . . . . 43. Error del Krigeado local universal . . . . . . . . . . . 44. Diferencias con wlc de la modelizacion de U . . . . . . 45. Prediccion en la modelizacion de V (wlm) . . . . . . . . 46. Desviacion t pica en la modelizacion de V (wlm) . . . . . 47. Diagramas de caja y bigote de las diferencias . . . . . . . 48. Prediccion en la modelizacion de U (wlm) . . . . . . . . 49. Desviacion t pica en la modelizacion de U (wlm) . . . . . 2.
32 33 34 35 35 39
40 40 41 42 43 43 45 46 47 47 48 50 51 52 53 54 54
Indice de cuadros
1. Resumen de estad sticos de datos de validacion y estimados . . . . . . . . . . . . . 26
3.
Estad sticos de los errores en los m etodos de krigeado de V . . . . . . . . . . . . . . 52 Estad sticos de los errores en los m etodos de krigeado de U . . . . . . . . . . . . . . 55
Indice de listados
1. 2. 3. 4. 5. 6. 7. R-Script del tema 1 . . . . . . R-Script del tema 2 . . . . . . R-Script en Linux del tema 3 . R-Script en Windows del tema 3 . . . . . . . . . . . . . . R-Script del tema 4 . . . . . . Script para GRASS . . . . . . Script para ps.map del m etodo IDW . . . . . . . . . 10 14 22 23 27 27 28
8. 9.
10. 11. 12. 13.
Script para ps.map del m etodo RST . . . . . . . . . . Script para ps.map del m etodo Pol gonos de inuencia . . . . . . . . . . . . Funciones para imprimir variogramas . . . . . . . . . . R-Script del tema 5 . . . . . . R-Script del tema 6 . . . . . . R-Script del tema 8 . . . . . .
28
29 36 36 44 56
Tema 0 Introduccion
0.1. R
R [5] es un conjunto integrado de herramientas para manipular datos, realizar todo
tipo de calculos con los mismos y tambi en es capaz de realizar toda clase de gracos estad sticos. En [4] se citan las siguientes caracter sticas: es multiplataforma, almacenamiento y manipulacion efectiva de datos, operadores para calculo sobre variables indexadas (Arrays), en particular matrices, una amplia, coherente e integrada coleccion de herramientas para analisis de datos, posibilidades gracas para analisis de datos, que funcionan directamente sobre pantalla o impresora, y un lenguaje de programacion bien desarrollado, simple y efectivo, que incluye condicionales, ciclos, funciones recursivas y posibilidad de entradas y salidas. (Debe destacarse que muchas de las funciones suministradas con el sistema estan escritas en el lenguaje R).
R puede extenderse mediante paquetes. En Linux , basta con ejecutar el comando install.packages(paquete) para conectar a la red de servidores CRAN (Comprehensive R Archive Network ) descarga el codigo fuente y si se dispone de los compiladores pertinentes (C++, Fortran, ...) genera los binarios adaptados perfectamente a la maqui na. En Windows, al ejecutar dicho comando se descargan directamente los binarios. En esta ultima plataforma se dispone de una interfaz graca un poco mas elaborada y permite ademas exportar al formato Windows MetaFile. Los paquetes empleados en el trabajo, ademas de los que se incluyen por defecto en R son los paquetes de geoestad stica gstat [3], y geoR [6], y el paquete para presentacion de gracos lattice [7]. En denitiva, se dispone de un sistema ampliable que se maneja como una consola de entrada de comandos que permite adquirir datos desde cheros, manipularlos, crear nuevos datos y por ultimo o bien ver las gracas por pantalla o mandarlas a cheros PostScript o raster . Otra caracter stica importante es la posibilidad de ejecutar secuencias de comandos en forma de scripts.
0.2. gstat
gstat [2] es un software para llevar a cabo modelizacion, prediccion y simulacion de
datos geoestad sticos. Al igual que el anterior, es Software Libre bajo licencia (GNU1 ). Puede usarse de muy diversas formas, directamente tanto de forma no interactiva (mediante cheros de parametros) como interactiva mostrando los resultados utilizando el programa para presentacion de gracos gnuplot. Pero su uso mas interesante es integrado con otras herramientas. En este sentido se ha conseguido que gstat funcione con GRASS, Idrisi, PCRaster y con R.
1
http://www.gnu.org
En este trabajo se ha usado con R porque este ultimo ofrece caracter sticas muy interesantes para la manipulacion de datos, presentacion de todo tipo de gracas y repeticion de tareas mediante sentencias de control (bucles, condicionales, etc). Por otro lado, no ha sido posible compilar gstat para que trabaje conjuntamente con GRASS en su version 6.
0.3. GRASS
Este ya veterano software para la gestion de informacion geograca dispone de herramientas para la modelizacion de variables espaciales mediante m etodos determin sticos. Se ha usado en este trabajo para la obtencion de la modelizacion por pol gonos de inuencia (Voronoi), Splines de tension y por el m etodo de pesos inversos a la distancia. Ademas se ha utilizado para la presentacion de la cartograf a, maquetando sencillos mapas con salida PostScript.
0.4. Datos de trabajo

Los datos con los que se va a trabajar durante todo el proyecto son los utilizados en el libro Applied Geostatistics de Issaks y Srivastava. El conjunto de datos walker esta a su vez dividido en tres grupos, wlc que es una malla de 78000 puntos que sirven para validacion, wlm es la malla irregular de 470 puntos y wle una malla de 100 puntos para algunos calculos estad sticos. gstat dispone del conjunto wlm, los otros dos seran cargados desde cheros de texto separados por comas (CSV) para poder operar con ellos.
univariada Tema 1 Descripcion

de los datos 1.1. Carga y visualizacion
En este cap tulo se van a usar los datos wle redondeados a valores enteros. El primer paso sera por tanto cargar el chero walker10.asc, que es un chero de texto separado por tabuladores que se puede importar directamente con la orden read.delim2 para a continuacion redondearlo. En la gura 1 se muestra la distribucion de los datos, as como los valores que toman las variables U y V.
250
15 81 16 82
12 77 7 61 9 74 8 70 18 88 16 82 15 80 15 80 17 84 28 100 12
24 103 34 110 22 97 27 103 20 94 16 86 15 85 15 83 11 74 4 47
27 112 36 121 24 105 27 111 27 110 23 101 16 90 15 87 29 108 32 111 14
30 123 29 119 25 112 32 122 29 116 24 109 17 97 16 94 37 121 38 124
0 19 7 77 10 91 4 64 19 108 25 113 18 101 17 99 55 143 20 109 16 X
2 40 4 52 7 73 10 84 7 73 7 79 14 96 13 95 11 91 0 0
18 111 18 111 19 115 15 105 16 107 15 102 6 72 2 48 3 52 14 98 18
18 114 18 117 19 118 17 113 19 118 21 120 28 128 40 139 34 136 31 134
18 120 20 124 22 129 19 123 22 127 20 121 25 130 38 145 35 144 34 144 20
248
16 82 21 88
246
21 89 15 77
244
14 74 14 75
242
16 77 22 87
Figura 1: Distribucion de U y V
univariada 1.2. M etodos gr acos para la descripcion

El m etodo graco mas utilizado es el histograma, en el que debemos integrar la variable en clases. La variable V se var a entre 0ppm y 145ppm por lo que dividirla en clases de 10 unidades es conveniente (g. 2 en la pagina siguiente). Otro graco interesante es el histograma acumulado en el que a partir de las mismas clases del histograma anterior se muestra la suma acumulada (g 3 en la pagina siguiente). El graco de probabilidad acumulada muestra la proporcion de datos para cada punto que son menores que e l (g 4 en la pagina siguiente) Las guras 5 en la pagina 8 y 6 en la pagina 8 muestran la similitud de nuestra muestra con la distribucion normal y lognormal. Las l neas trazadas pasan por el primer y tercer cuartil.
17 16 15 14 14 12 11
Frequency
10
12
14
4 3 3 2 1 1 0 1 2
0 0
50 V
100
150
Figura 2: Histograma de V
Frecuencia acumulada
20
40
60
80
100
10
20
30
40
50
60
70
80
90 100
120
140
Variable
Figura 3: Histograma acumulado de V
Probabilidad acumulada
0.0
0.2
0.4
0.6
0.8
1.0
50 V
100
150
Figura 4: Graco de probabilidad uniforme de V
Theoretical Quantiles
2 0
50 Sample Quantiles
100
150
Figura 5: Graco de probabilidad normal
Theoretical Quantiles
3.0
3.5
4.0 Sample Quantiles
4.5
5.0
Figura 6: Graco de probabilidad lognormal
1.3. M etodos num ericos

1.3.1. Medidas de localizacion Se puede solicitar una descripcion sencilla de nuestros datos con el comando summary(V) que devuelve tanto los valores maximos y m nimos, la media, la mediana y el segundo y tercer cuartil. En cualquier caso estan disponibles comandos como min, max, mean y median. Para calcular la moda no hay un comando denido, pero a partir de la tabla denida del corte de V (tcutV) donde se almacenan las frecuencias relativas en las clases denidas previamente (secV), podemos solicitar aquella clase que almacene el valor maximo con el comando tcutV[tcutV==max(tcutV)]. En resumen: M nimo Maximo Media Mediana Moda Rango 0 145 100.5 97.55 110-120 145
R puede calcular cualquier cuantil de una muestra, por ejemplo los cuartiles con
el comando quantile y pasando un vector con los valores de los cuantiles a obtener, en este caso una secuencia de 0 a 1 cada 0.25 unidades:
> print(cuantiles<-quantile(V,seq(0,1,.25))) 0% 25 % 50 % 75 % 100 % 0.00 81.75 100.50 116.25 145.00
Una forma graca de ver tanto los cuantiles como la distribucion de la muestra y si existen valores alejados de la media (outliers) es el diagrama de caja y bigotes (box and whisker ). La gura 7 muestra el de la variable V.
50
100
150
Figura 7: Graco de caja y bigotes de V
1.3.2. Medidas de dispersion Las medidas de dispersion como la varianza, la desviacion t pica y el rango intercuantil son sencillos de calcular:
> dt<-sqrt(var(V)) > dt2;dt;as.numeric(cuantiles["75 %"]-cuantiles["25 %"])
10
[1] 695.3409 [1] 26.36932 [1] 34.5
1.3.3. Medidas de forma El coeciente de sesgo o asimetr a (skewness) se calcula a partir de la formula CS =
1 n n i=1 (xi 3
m)3
(1)
El coeciente de curtosis o apuntalamiento se calcula como K=

(xi m)4 n i=1 n 4
(2)
El coeciente de variacion no es mas que el cociente entre la desviacion t pica y la media, siendo trivial su calculo. En R estos tres coecientes se calculan como:
> media<-mean(V); > print(CS<-sum((V-media)3)/(length(V)*dt3)) [1] -0.7665234 > print(K<-sum((V-media)4/length(V))/dt4-3) [1] 1.187891 > print(CV<-dt/media) [1] 0.2703159
1.4. Scripts
Listado 1: R-Script del tema 1
N UNIVARIADA 1 #TEMA 1 - DESCRIPCIO 2 ps.options(family="Bookman",pointsize
3 4 5 6 7 8 9 25 26 27 28 29 30 31 32 33 34 35 36
dev.off() #Discretizar V para hallar frecuencias lsecV<-length(secV) cutV<-cut(V,secV) tcutV<-table(cutV) #Obtener las frecuencias acumuladas sumfreq<-rep(0,length(secV)) for (i in 1:lsecV){ sumfreq[i]<-length(V [V<secV[i]])} #Imprimir el histograma postscript("imgs/01/histAcumV.ps") barplot(sumfreq,xlab="Variable",ylab=" Frecuencia acumulada",names.arg=secV ) dev.off() #Gr afico de probabilidad acumulada postscript("imgs/01/probUnif.ps") plot.ecdf(V,pch=1,xlab="V",ylab=" Probabilidad acumulada",main="") dev.off()
=15) rm(list=ls()) #Cargar wle y redondearlo wle<-read.delim2("walker10.asc") wle<-round(wle,0)
#Adjuntar los datos de wle para acceder directamente 10 attach(wle)
11 12 #Ver los datos U y V 37 13 postscript("imgs/01/UyV.ps") 38 14 plot(X,Y,xlim=c(min(X),max(X)*1.01),ylim 39 =c(min(Y)*.999,max(Y)*1.001),pch=3) 40 15 text(X+.3,Y+.2,U) 41 16 text(X+.3,Y-.2,V) 17 dev.off() 42 18 43 19 #Mostrar el histograma de V 44 20 secV<-seq(0,150,10) 21 postscript("imgs/01/histV.ps") 45 22 hist(V,breaks=secV, labels=TRUE, col=" 46
lightgray", axes=FALSE,main="")
23 axis(2,at=seq(0,18,2)) 24 axis(1,at=seq(0,150,50))
#Gr afico de probabilidad normal y lognormal postscript("imgs/01/probNormal.ps") qqnorm(V,datax=TRUE,pch=3,main=""); qqline(V,datax=TRUE) 47 dev.off()
11
48 49 logV<-log(V)[is.finite(log(V))] 50 postscript("imgs/01/probLogNormal.ps") 51 qqnorm(logV,datax=TRUE,pch=3,main="");
62 print(as.numeric(cuantiles["75 %"]-
cuantiles["25 %"]))
63 64 #Diagrama de caja y bigotes de V 65 postscript("imgs/01/cajaybig.ps") 66 boxplot(V,horizontal=TRUE,col="lightgray
qqline(logV,datax=TRUE)
52 dev.off() 53 54 #Medidas de localizaci on, dispersi on,
",ylab="V")
67 dev.off() 68 69 #Medidas de forma: sesgo, apuntalamiento
etc media<-mean(V) dt<-var(V) print("Cuartiles") print(cuantiles<-quantile(V,seq(0,1,.25) )) on t pica y 59 print("Varianza, desviaci rango intercuart lico") 60 print(dt2) 61 print(dt)
55 56 57 58
y variaci on
70 print("Sesgo, apuntalamiento y variaci on
")
71 print(CS<-sum((V-media)3)/(length(V)*dt
3))
72 print(K<-sum((V-media)4/length(V))/dt
4-3)
73 print(CV<-dt/media)
12
bivariada Tema 2 Descripcion

2.1. M etodos gr acos
La visualizacion de pares de histogramas y sobre todo de gracos de caja y bigotes ( 8)pueden aportar informacion de como son dos variables.
Variables V 0
50
100
150
Figura 8: Graco de caja y bigotes de V y U El diagrama de cuantiles muestra cada cuantil de una variable contra el mismo cuantil de la otra formando un graco de puntos. Si estos puntos adoptan la forma de una l nea signica que ambas variables son similares pero su localizacion y dispersion son diferentes.
100% 50 U=V
40
30 75% 20 25% 10 Cuartiles 0 0% 50%
50 V
100
150
Figura 9: Graco de cuantiles de V y U Por ultimo, el diagrama de dispersion ( 10 en la pagina siguiente)puede mostrar informacion sobre las tendencias de ambas variables y la existencia de outliers que
13
pueden ser susceptibles de ser eliminados.
U 0 0 10 20
30
40
50
50 V
100
150
Figura 10: Graco de dispersion
2.2. M etodos num ericos

El coeciente de correlacion o de Pearson es el indicador mas utilizado para comprobar la relacion entre variables, en R se corresponde con el comando cor(U,V). =
1 n n i=1 (xi
mx ) (yi my ) x y
(3)
El numerador de esta ultima ecuacion es otro indicador, la covarianza, que se obtiene con el comando cov(U,V). El coeciente de correlacion de orden utiliza el orden en que aparece un valor determinado en lugar de su propio valor. Por esta razon este coeciente es menos sensible a valores extremos y por tanto, si diere mucho del coeciente de Pearson indica la existencia de estos valores extremos. En R se obtiene con el mismo comando que el de correlacion, pero cambiando el m etodo. Orden =
1 n n i=1 (Rxi
mRx ) (Ryi mRy ) Rx Ry
(4)
> cov(U,V);cor(U,V,method="pearson");cor(U,V,method="spearman") [1] 218.2778 [1] 0.8393395 [1] 0.8576604
14
2.3. Scripts
N BIVARIADA 1 #TEMA 2 - DESCRIPCIO 2 ps.options(family="Bookman",pointsize
3 4 5 6 7 8 9 23 24 25
=15) rm(list=ls()) #Cargar wle y redondearlo wle<-read.delim2("walker10.asc") wle<-round(wle,0)

11 12 #Diagrama de caja y bigotes de V y U 13 postscript("imgs/02/cajaybig.ps") 14 boxplot(data.frame(V,U),horizontal=TRUE,
26 27 28 29 30 31 32 33 34 35 36 37 38 39
(U,seq(0,1,.05)),pch=21,bg=" lightgray") points(quantile(V,seq(0,1,.25)),quantile (U,seq(0,1,.25)),pch=21,bg="white") # Texto en los cuartiles text(quantile(V,seq(0,1,.25))-2,quantile (U,seq(0,1,.25))+2,c("0 %","25 %","50 % ","75 %","100 %")) # Linea con U=V y su texto lines(c(0,150),c(0,150),type="l") text(60,50,"U=V") # Leyenda legend(100,8,c("Cuartiles"),bg="white", pch=21) dev.off() #Diagramas de dispersi on postscript("imgs/02/dispers.ps") plot(V,U,pch=21,bg="lightgray") dev.off() #Covarianza y coeficientes de correlaci on print("Covarianza y coeficiente de correlaci on de Pearson y de Spearman ") print(covar<-cov(U,V)) print(rho<-cor(U,V,method="pearson")) print(spear<-cor(U,V,method="spearman"))
col="lightgray",ylab="Variables", boxwex=0.5) 15 dev.off() # Diagrama qqplot postscript("imgs/02/qqnorm.ps") # Dibuja s olo los ejes qqplot(V,U,xlab="V",ylab="U",las=1,xlim= c(-1,148),ylim=c(0,57),type="n") 21 # Dibuja los puntos por cuantiles del 5 % y los cuartiles 22 points(quantile(V,seq(0,1,.05)),quantile
16 17 18 19 20
40 41 42
15
espacial Tema 3 Descripcion

espacial de datos 3.1. Visualizacion
3.1.1. Mapas de localizacion Se trata de trazar un mapa de puntos e indicar por ejemplo los 10 valores maximos y los 10 valores m nimos. En R se trata de obtener los conjuntos de datos de wle con los valores maximos y m nimos y pintarlos sobre un mapa de localizacion de V.
250
81 82
77 61 74 70 88 82 80 80 84 100 Mximos Mnimos 12
103 110 97 103 94 86 85 83 74 47
112 121 105 111 110 101 90 87 108 111
123 119 112 122 116 109 97 94 121 124
19 77 91 64 108 113 101 99 143 109
40 52 73 84 73 79 96 95 91 0
111 111 115 105 107 102 72 48 52 98
114 117 118 113 118 120 128 139 136 134
120 124 129 123 127 121 130 145 144 144
248
82 88
246
89 77
Y 244
74 75
242
77 87
240
14 X
16
18
20
Figura 11: Distribucion de V
3.1.2. Mapas de s mbolos graduados Utilizando la biblioteca lattice, se pueden generar mapas de s mbolos graduados por color con el comando levelplot. R dispone de diferentes paletas de color, pero en este caso se usara una escala de grises. Otro tipo de mapa es el de gradacion de s mbolos por tamano. En este caso el paquete gstat proporciona el comando bubble. El resultado es mejorable, pero sirve a modo de ejemplo. 3.1.3. Mapas de indicadores Se trata de mapas graduados de color con solo dos niveles. El umbral entre ambos niveles se va variando y se observa el conjunto de mapas generados. Estos mapas pueden mostrar alineaciones en la distribucion espacial de la variable y la ubicacion de maximos y m nimos. Las guras 14(c) y 14(d) muestran la alineacion norte-sur de los datos.
16
140
120 248
100
246
80
60
244 40
20 242
12
14
16
18
Figura 12: Mapa graduado de color de V
250
248
246
0 81.75 100.5 116.25 145
y
244 242
12
14
16
18
20
Figura 13: Mapa graduado de tamano de V
17
Umbral = 73.375
Umbral = 81.75
140
140
120
120
248 100
248 100
246
80
246
80
60 244 244
60
40
40
242
20
242
20
0 12 14 16 18 12 14 16 18
(a) Cuantil 12.5 %

Umbral = 89.125
(b) Cuantil 25 %
Umbral = 100.5
140
140
120
120
248 100
248 100
246
80
246
80
60 244 244
60
40
40
242
20
242
20
0 12 14 16 18 12 14 16 18
(c) Cuantil 37.5 %

Umbral = 109.875
(d) Cuantil 50 %
Umbral = 116.25
140
140
120
120
248 100
248 100
246
80
246
80
60 244 244
60
40
40
242
20
242
20
0 12 14 16 18 12 14 16 18
(e) Cuantil 62.5 %
(f) Cuantil 75 %
Figura 14: Mapas de indicadores
18
3.1.4. Mapas de supercies interpoladas Es posible interpolar polinomios de un orden determinado, que pasen por todos los puntos de la muestra. Estos polinomios suelen ofrecer supercies suavizadas pero que muestran las tendencias y la distribucion de la variable. Para poder obtener estos polinomios en R, se ha de cargar el paquete spatial que dispone de la funcion surf.ls que obtiene el mejor polinomio de hasta grado 6 que se ajusta a los datos mediante m nimos cuadrados (g. 15). Este mapa ha sido generado en Windows ya que el paquete spatial no ha sido posible instalarlo en Linux, por lo que el script que genera este graco esta separado del resto del tema.
242
244
246
248
250
12
14
16
18
20
Figura 15: Mapa de supercie interpolada
3.2. Ventanas moviles y el efecto proporcional

Se pueden trazar mapas en las que el valor de cada punto es la media o la varianza de sus vecinos. Por lo tanto se crea una ventana movil de tamano impar que va recorriendo los valores. En R esta operacion requiere algo de programacion que se incluye en el script del tema y que se presenta en la gura 16 en la pagina siguiente. Un graco de dispersion entre las medias y las varianzas demuestra la existencia de una relacion entre ambas. En este caso, en la gura 17 en la pagina siguiente, se puede ver que ambas medidas no estan correladas, siendo el coeciente de correlacion () bastante bajo.
3.3. Continuidad espacial

de tipo h 3.3.1. Diagramas de dispersion Este tipo de diagramas muestran contrapuestas la variable contra e sta a una distancia determinada. Cuando los datos estan en forma de malla, se pueden obtener de forma sencilla los mapas de dispersion en las direcciones norte-sur y este-oeste a diferentes distancias. Obtener estos diagramas en R es relativamente sencillo debido a la exibilidad en la manipulacion de variables indexadas (arrays). Se han obtenido los diagramas con
19
250
85.22 235.94 85.22 255.19 87.22 111.69 86.33 90.25 83.89 39.36 80.22 17.94 79.11 18.61 78.56 201.78
95.56 409.53 94.67 443.25 94.67 220.5 93.89 186.86 90.67 92.75 86 42.5 85.67 91.25 86 377.5
111.33 75.75 111.11 72.86 107.78 79.94 105.78 122.94 98.67 123.5 92.44 74.03 93.22 199.44 94.33 611.5
97.67 1095.75 102.44 428.53 104.33 300.25 106 281 105 69.5 99 69.75 104.44 311.53 110.67 289.25
78.44 1332.03 88.22 616.94 93.67 459.75 96.44 457.78 99.11 219.86 98.11 93.36 104.11 288.36 97.33 1628.25
76.56 1165.53 85.78 468.19 91.11 342.86 92.78 320.44 94.56 248.78 89.44 392.78 88.56 823.28 81.67 1753
94.56 955.28 98.67 553.25 100.67 353.25 100.11 299.61 99.44 441.03 97.67 828.75 95.22 1163.69 88.11 2227.36
117.67 35.5 117.22 53.69 117.22 68.19 115.11 76.86 113.89 334.86 111.67 1042.25 110.44 1659.03 115.56 1588.03
Y 242 244
246
248
12
14 X
16
18
20
Figura 16: Media y varianza en ventana de 3x3
2000
= 0.11
varianza
500
1000
1500
80
90 media
100
110
Figura 17: Graco de dispersion de media y varianza
20
distancia (h) variando de 1 a 4 puntos y en las direcciones norte-sur (g. 18) y esteoeste (g:03:vew). Se aprecia la alta correlacion en las direcciones N-S (0,5 < < 0,7) y la baja correlacion en direccion E-W (0,4 < < 0,3).
h=1
140 140
h=2
120
= 0.74
120
= 0.59
100
V(t)
80
V(t) x=y
60
40
40
60
80
100
x=y 20
20 0
50 V(t+h)
100
150
50 V(t+h)
100
150
h=3
h=4
120
120
= 0.56
= 0.48
100
V(t)
80
V(t)
60
40
40
60
80
100
x=y 20 20
x=y
50 V(t+h)
100
150
50 V(t+h)
100
150
Figura 18: h-Scatterplots de direccion N-S
3.4. Variograma
La funcion del variograma muestra la variacion de la variable agrupando los datos segun sus distancias relativas. Ha de establecerse por tanto el numero de ((cajas)) o lags en los que queremos dividir los datos. Por otro lado, como en el apartado anterior, podemos obtener el variograma en una direccion determinada o en todas direcciones (variograma omnidireccional). Formalmente la funcion del variograma se expresa como: 1 (h) = 2N (h)
N (h)
(vi vi+h )2
i=1
(5)
En R existen diversos paquetes que calculan el variograma pero probablemente el mas completo es el paquete gstat. El uso de esta funcion se hara en el tema 5.
cruzados 3.5. Diagramas de dispersion

El ultimo graco que se va a mostrar en este tema es el de dispersion cruzada entre dos variables. Es un h-Scatterplot en el que en lugar de ccontrastarla misma variable se utiliza el valor de otra. Se ha calculado solo la variacion en la direccion N-S de las variables U y V.
3.6. Scripts
21
h=1
140 140
h=2
120
100
= 0.21 80 V(t) V(t) 80
100
120
= 0.41
60
40
x=y 20 20
40
60
x=y
50 V(t+h)
100
150
0 0
50 V(t+h)
100
150
h=3
140 140
h=4
120
100
= 0.35 V(t) 80 V(t)
100
120
= 0.3 80 40 60
40
60
x=y 20 x=y 20 0 50 V(t+h) 100 150 0 0
50 V(t+h)
100
150
Figura 19: h-Scatterplots de direccion E-W
h=1
40
h=2
50
30
30
40
x=y
U(t)
U(t)
x=y
20
= 0.6 10
10
50 V(t+h)
100
150
20
= 0.45
50 V(t+h)
100
150
h=3
35 35
h=4
30
x=y
30
x=y
25
20
U(t)
15
U(t)
10
50 V(t+h)
100
150
10
15
20
= 0.36
25
= 0.28
50 V(t+h)
100
150
Figura 20: h-Scatterplots cruzado de U y V en direccion N-S
22
Listado 3: R-Script en Linux del tema 3

N ESPACIAL 1 #TEMA 3 : DESCRIPCIO 2 ps.options(family="Bookman",pointsize =15) 3 rm(list=ls())
4 5 6 7 8 9 10 11 12 13
48 postscript("imgs/03/indics %02d.ps",
onefile=FALSE)
49 for (i in 1:length(sec2)){ 50 tit<-as.expression(substitute(
#Cargar paquetes library(lattice) library(gstat) #Cargar wle y redondearlo wle<-read.delim2("walker10.asc") wle<-round(wle,0)
Umbral==p,list(p=sec2[i]))) lv<-levelplot(VX*Y,wle,col. regions=c("lightgray","white "),at=c(min(V),sec2[i],max(V )),main=tit); 52 print(lv) 53 } 54 dev.off()

51 55 56 57 # Mapas de isol neas y de superficies
interpoladas en Windows
58 59 #Conseguir una matriz con los valores de 60 61 62 63 64 65

15 16 #Gr afico de V con m aximos y m nimos 17 #Obtener los 10 valores m aximos y 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33
68 69 70 71 72 73 74 75 76 77 78 #Gr afico de s mbolos graduados 79 postscript("imgs/03/simbs.ps") 80 lv<-levelplot(VX*Y,wle,col.regions=gray 81 (seq(1,0.3,len=18)),aspect=mapasp( 82
m nimos Vmax=sort(V)[90:100][1] Vmin=sort(V)[1:10][10] dfVmax= subset(wle,V>=Vmax) dfVmin= subset(wle,V<=Vmin) #Ya se puede hacer el dibujo postscript("imgs/03/maxmin.ps") plot(X,Y,pch=3,xlim=c(11,20.2),ylim=c (239.5,250.5)) text(X+0.2,Y+0.2,V) points(dfVmax$X,dfVmax$Y,bg="black",pch =21) points(dfVmin$X,dfVmin$Y,bg="lightgray", pch=21) legend(11,240.8,c("M aximos","M nimos"), pch=21,pt.bg=c("black","lightgray")) dev.off()
V imgV<-xyz2img(wle,"V",2,1) matVtemp<-imgV$z matV<-matVtemp*0 m<-dim(matVtemp)[1] n<-dim(matVtemp)[2] for(i in 1:m){for(j in 1:n){matV[i,j]<matVtemp[m+1-i,j] }}
66 67 imgV2<-list(sort(imgV$y),sort(imgV$x,
decreasing=TRUE),matV) names(imgV2)<-c("x","y","z") #Calcular una media y varianza 3x3 media3<-matV*0 var3<-media3 seqm<-seq(2,m-1) seqn<-seq(2,n-1) for(i in seqm){ for(j in seqn){ i1<-i-1; i2<-i+1; j1<-j-1; j2<-j+1; matVtemp<-array(matV[i1: i2,j1:j2]); media3[i,j]<-mean( matVtemp); var3[i,j]<-var(matVtemp) } } media3<-round(media3[seqm,seqn],2) var3<-round(var3[seqm,seqn],2)
#Gr afico de burbujas postscript("imgs/03/bubbles.ps") # bub<-bubble(wle,1,2,"V",key.entries= quantile(V,seq(0,1,.25))) 40 bub<-bubble(wle,1,2,"V",fill=TRUE,pch =21,key.entries=quantile(V,seq (0,1,.25))) 41 print(bub) 42 dev.off()
43 44 #Mapas de indicadores (por 12.5 %, 25 %,
34 35 36 37 38 39
wle)) print(lv) dev.off()
83 84 85 86 87 88 89 90 91 92 93 94 95
37.5 %, 50 %, 62.5 %, 75 %)
45 # sec2<-quantile(V,c
(12.5,25,37.5,50,62.5,75)/100)
46 sec1<-seq(0.125,.75,.125) 47 sec2<-quantile(V,sec1)
#Mostrar resultados postscript("imgs/03/mediavar.ps") plot(X,Y,pch=3) for(i in seqm){ for(j in seqn){ text(imgV2$x[i],imgV2$y[ j]+0.25,media3[j-1,i -1]); 96 text(imgV2$x[i],imgV2$y[ j]-0.25,var3[j-1,i
23
-1])
97 98 99 100 101 102 103 104 105
139 140
} } dev.off()
#Gr afico de media contra varianza 3x3 141 postscript("imgs/03/mediavar2.ps") media<-array(media3,dim=64) 142 varianza<-array(var3,dim=64) plot(media,varianza,pch=21,bg="lightgray 143 144 ",xlab="media",ylab="varianza") 106 text(105,2000,as.expression(substitute( 145 146 rho==ro,list(ro=round(cor(media, 147 varianza),2))))) 148 107 dev.off()
lines(c(0,150),c(0,150));text (40,30,"x=y"); text(10,90,as.expression( substitute(rho==ro,list(ro= round(cor(matV1,matV2),2)))) ) title(as.expression(substitute(h ==hh,list(hh=h)))) } dev.off()
108 109 110 #Continuidad espacial en V 111 #NS 112 postscript("imgs/03/V-NS %02d.ps",onefile 152 153 #NS =FALSE) 113 for(i in 1:4){ 154 postscript("imgs/03/VU-NS %02d.ps", onefile=FALSE) 114 h<-i; 155 for(i in 1:4){ 115 fil1<-seq(1+h,m); 156 h<-i; 116 col1<-seq(1,n); fil1<-seq(1+h,m); 117 matV1<-array(matV[fil1,col1],dim 157 col1<-seq(1,n); =length(fil1)*length(col1)); 158 159 matV1<-array(matV[fil1,col1],dim 118 fil2<-seq(1,m-h); =length(fil1)*length(col1)); 119 col2<-seq(1,n); fil2<-seq(1,m-h); 120 matV2<-array(matV[fil2,col2],dim 160 col2<-seq(1,n); =length(fil2)*length(col2)); 161 matU2<-array(matU[fil2,col2],dim 121 plot(matV1,matV2,xlab="V(t+h)", 162 122 123
#Continuidad espacial en V cruzado con U imgU<-xyz2img(wle,"U",2,1) matUtemp<-imgU$z matU<-matUtemp*0 149 m<-dim(matUtemp)[1] 150 n<-dim(matUtemp)[2] 151 for(i in 1:m){for(j in 1:n){matU[i,j]<matUtemp[m+1-i,j] }}
124 125 126 127 128 129 130 131 132 133 134 135 136 137 138
ylab="V(t)"); 163 lines(c(0,150),c(0,150));text (40,30,"x=y"); 164 text(40,120,as.expression( substitute(rho==ro,list(ro= round(cor(matV1,matV2),2)))) 165 ) title(as.expression(substitute(h ==hh,list(hh=h)))) } dev.off()
166 167 } 168 #EW postscript("imgs/03/V-EW %02d.ps",onefile =FALSE) for(i in 1:4){ h<-i; 1 fil1<-seq(1,m); 2 col1<-seq(1+h,n); 3 matV1<-array(matV[fil1,col1],dim 4 =length(fil1)*length(col1)); 5 fil2<-seq(1,m); col2<-seq(1,n-h); 6 matV2<-array(matV[fil2,col2],dim 7 =length(fil2)*length(col2)); 8 plot(matV1,matV2,xlab="V(t+h)", 9 ylab="V(t)");
=length(fil2)*length(col2)); plot(matV1,matU2,xlab="V(t+h)", ylab="U(t)"); lines(c(0,150),c(0,150));text (40,30,"x=y"); text(10,20,as.expression( substitute(rho==ro,list(ro= round(cor(matV1,matU2),2)))) ) title(as.expression(substitute(h ==hh,list(hh=h)))) dev.off()
Listado 4: R-Script en Windows del tema 3

orden<-6 V.ls<-surf.ls(orden,wle$X,wle$Y,wle$V) V.tr<-trmat(V.ls,11,20,241,250,150) ps.options(family="Bookman",pointsize =15) postscript("imgs/03/interp.ps") image(V.tr,col=gray(seq(0.93,0.3,l=50))) contour(V.tr,add=TRUE,labcex=.8) dev.off()
24
M Tema 4 Estimacion. etodos deterministas

En esta seccion se va a estimar el valor de la variable V en toda la extension de trabajo a partir de los valores en el conjunto de datos wlm (g. 21). Este conjunto de datos se obtiene facilmente al estar presente en el paquete gstat y cargandose con ejecutar el comando data(walker).
250
200
150
100
50
50
100
150
200
250
Figura 21: Distribucion de wlm La estimacion de estos datos se ha realizado desde GRASS, para ello primero se ha exportado este conjunto de datos a un chero csv para a continuacion importarlo en GRASS como una cobertura vectorial (GRASS ya no utiliza sites). En GRASS se han utilizado las funciones v.voronoi, v.surf.idw y v.surf.rst que implementan el m etodo de pol gonos de inuencia, el de pesos inversos a la distancia y el de splines de tension respectivamente. Ademas de las coberturas raster con los valores estimados se han obtenido las curvas de nivel cada 250ppm. Finalmente se han maquetado tres sencillos mapas con los resultados de estos tres m etodos. Al nal del tema se presenta el script en GRASS que genera los mapas y exporta a cheros csv los valores de los tres mapas para poder estudiarlos en R. Igualmente se presentan los tres cheros que indican la maquetacion de los mapas. A continuacion se importa el conjunto de datos wlc que consiste en 78000 puntos de validacion y que se entiende como valores correctos.. En el cuadro 1 se muestran los estad sticos de estos cuatro conjunto de datos y en la gura 23 en la pagina 26 se pueden ver los histogramas. Se pueden observar las siguientes caracter sticas: Se observa como el m etodo de Pol gonos de inuencia es el que mejor mantiene las caracter sticas estad sticas.
Leyenda Modelizacion geoestad stica con R 1400

1400 1200 1000 800 800 600 600 400 400 200 200 0 0 100 200
25
Inverso la distancia de Splines de a Tensin de V (T=100) Polgonos de influencia deV V

Escala
0 50
Equidistancia de curvas = 250ppm
200
100
Figura 22: Mapas generados por GRASS
26
Cuadro 1: Resumen de estad sticos de datos de validacion y estimados

wlc
M n. Q1 Median. Media Q3 Max. 2 Simetr. Apunt. CV
0.00 67.79 221.30 278.00 428.30 1631.00 62423.16 249.84 0.90 0.77 0.90
Voronoi 0.00 70.70 224.40 275.10 425.90 1528.00 60039.41 245.03 0.89 1.25 0.89
IDW 0.46 163.00 297.20 316.80 428.50 1498.00 39458.06 198.64 0.63 1.39 0.63
RST -47.51 121.20 277.90 277.90 379.50 1530.00 44820.41 211.71 0.76 2.44 0.76
El m etodo de Splines de tension ofrece resultados negativos pese a calcularse con un valor de tension alto (100). Cabe destacar la diferencia entre los m etodos estad sticos y wlc en el apuntalamiento, claramente observable en los histogramas.
WLC
25000
Vor
20000
Frecuencia
15000
Frecuencia 0 500 V 1000 1500
10000
5000
0 0
5000
10000
15000
20000
500 V
1000
1500
IDW
15000
RST
Frecuencia
10000
Frecuencia 0 500 V 1000 1500
5000
5000
10000
15000
500 V
1000
1500
Figura 23: Histogramas de los conjuntos de datos
27
4.1. Scripts
N 1 #TEMA 4 - ESTIMACIO 2 ps.options(family="Bookman",pointsize
3 4 5 6 7 8 9 10 11 12 13 14
="Frecuencia")
49 hist(datos[[2]],main="Vor",xlab="V",ylab
="Frecuencia")
50 hist(datos[[3]],main="IDW",xlab="V",ylab
=15) rm(list=ls()) #Cargar los datos wlm library(lattice) library(gstat) data(walker) wlm<-walker attach(wlm)
="Frecuencia")
51 hist(datos[[4]],main="RST",xlab="V",ylab
="Frecuencia")
52 dev.off()
Listado 6: Script para GRASS

1 # CREAR MAPA DE SPLINES 2 g.remove rast=splines vect=spcontour 3 v.surf.rst input=wle elev=splines
#Ver wlm postscript("imgs/04/wlm.ps") xy<-xyplot(YX,wlm,pch=3,aspect=mapasp( wlm)) 15 print(xy) 16 dev.off()

17 18 #Importar valores reales wlc 19 wlc<-read.csv("wlc.csv",sep=";",dec=".",
zcolumn=u tension=100
4 r.contour input=splines output=spcontour
step=250
5 r.colors map=splines rules=grey 6 7 8 # CREAR MAPA DE VORONOI 9 g.remove rast=rvorwle vect=clvorwle,
vorwle
10 v.voronoi input=wle ouput=vorwle 11 v.clean input=vorwle output=clvorwle
header=TRUE)
20 21 #Importar los valores estimados en GRASS 22 idw<-read.csv(file="/geo/r/mapasgrass/
idw.values.asc",sep="|")
23 vor<-read.csv(file="/geo/r/mapasgrass/
type=boundary,line,centroid,area tool=break,rmdupl,snap 12 g.region vect=clvorwle 13 r.colors map=rvorwle rast=splines # CREAR MAPA IDW g.remove rast=idw vect=idwcontour v.surf.idw input=wle output=idw col=u r.colors map=idw rast=splines r.contour input=idw output=idwcontour step=250
14 15 mapasgrass/splines.values.asc",sep=" 16 17 |") 18 25 26 # #Mostrar medias, varianzas y demas.... 19 27 Desv<-function(datos){sqrt(var(datos))} 20 28 CS<-function(V){sum((V-mean(V))3)/( 21 length(V)*Desv(V)3)} 29 K<-function(V){sum((V-mean(V))4/length( 22 24 splines<-read.csv(file="/geo/r/
voronoi.values.asc",sep="|")
V))/Desv(V)4-3}
30 CV<-function(V){Desv(V)/mean(V)} 31 32 datos<-list(wlc$V,vor$value,idw$value, 33 34 35 36 37 38 39 40 41 42 43 44 45 46
# GENERAR MAPAS PS ps.map input=splines.psmap output= splines.ps 23 ps.map input=voronoi.psmap output= voronoi.ps 24 ps.map input=idw.psmap output=idw.ps
25 26 # EXPORTAR A CSV 27 r.to.vect input=idw output=idw feature=
splines$value) print("Resumen") print(lapply(datos,summary)) print("Varianza") print(lapply(datos,var)) pica") print("Desviaci on T print(lapply(datos,Desv)) print("Simetr a") print(lapply(datos,CS)) print("Apuntalamiento") print(lapply(datos,K)) print("Coeficiente de Variaci on") print(lapply(datos,CV))
point #convertir a puntos

28 r.to.vect input=splines output=splines
feature=point #convertir a puntos

29 r.to.vect input=rvorwle output=voronoi
feature=point #convertir a puntos

30 31 echo "SELECT * FROM idw" | db.select >
idw.values.asc #valores
32 echo "SELECT * FROM splines" | db.select
> splines.values.asc #valores

33 echo "SELECT * FROM voronoi" | db.select
postscript("imgs/04/hist %01d.ps",onefile > voronoi.values.asc #valores =FALSE) 34 47 # layout(matrix(1:4,2,2,byrow=TRUE)) 35 #v.out.ascii input=idw output=idw.coor. 48 hist(datos[[1]],main="WLC",xlab="V",ylab asc format=point # si quisieramos
28
coordenadas
Listado 7: Script para ps.map del m etodo IDW

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57
grid 50 color black numbers 2 black end vpoints wle color black fcolor black size 1 label u end vlines idwcontour color white width 0.5 end raster idw maploc 0.7 3 7 10
58 59 60 61 62 63 end
left 0.5 right 0.5 bottom 0.5 top 0.5 end
Listado 8: Script para ps.map del m etodo RST

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52
grid 50 color black numbers 2 black end vpoints wle color black fcolor black size 1 label u end vlines spcontour color white width 0.5 end raster splines maploc 0.7 3 7 10
rectangle 0 390 260 305 color black fcolor white end text 160 365 Inverso a la distancia de V font Bookman size 10 end text 20 385 Leyenda font Bookman size 6 end colortable y where 2 0.7 raster idw width 0.3 height 2 cols 6 font Bookman fontsize 8 color black end text 160 320 Equidistancia de curvas = 250ppm font Bookman size 6 end text 160 345 Escala font Bookman size 6 end scalebar f where 5 2.1 length 50 height 0.05 segment 5 numbers 5 fontsize 8 end paper a4
rectangle 0 390 260 305 color black fcolor white end text 160 365 Splines de Tensi on de V (T=100) font Bookman size 10 end text 20 385 Leyenda font Bookman size 6 end colortable y where 2 0.7 raster splines width 0.3 height 2 cols 6 font Bookman fontsize 8 color black end text 160 320 Equidistancia de curvas = 250ppm font Bookman size 6 end text 160 345 Escala font Bookman size 6 end scalebar f where 5 2.1 length 50 height 0.05
29
53 segment 5 54 numbers 5 55 fontsize 8 56 end 57 paper a4 58 left 0.5 59 right 0.5 60 bottom 0.5 61 top 0.5 62 end 63 end
Listado 9: Script para ps.map del m etodo Pol gonos de inuencia

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22
grid 50 color black numbers 2 black end vpoints wle color black fcolor black size 1 end vareas clvorwle color black fcolor none end raster rvorwle maploc 0.7 3 7 10
rectangle 0 390 260 305 color black fcolor white end text 160 365 Pol gonos de influencia de V
23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58
font Bookman size 10 end text 20 385 Leyenda font Bookman size 6 end colortable y where 2 0.7 raster rvorwle width 0.3 height 2 cols 6 font Bookman fontsize 8 color black end text 160 345 Escala font Bookman size 6 end scalebar f where 5 2.1 length 50 height 0.05 segment 5 numbers 5 fontsize 8 end paper a4 left 0.5 right 0.5 bottom 0.5 top 0.5 end end
30
Tema 5 Continuidad espacial de V

5.1. Variograma omnidireccional
El calculo del variograma omnidireccional con gstat es sencillo. En primer lugar se ha de crear un objeto de tipo gstat anadi endole los datos de trabajo. A continuacion se ejecuta el comando variogram con las opciones pertinentes. La opcion principal es el tamano de los lags. Ningun paquete de los consultados permite elegir una tolerancia para incluir pares en cada lag, siendo unicamente GSLIB el que acepta este parametro, pero no en R. Otros parametros permiten obtener el mapa supercial o la nube de puntos o establecer restricciones de direccion como se vera mas adelante. En el script del tema se han creado diversas funciones para realizar gracas de uno o varios variogramas presentando el numero de pares y una leyenda. Si se muestra el variograma en R, aparece el listado con los lags, el numero de pares, el valor del variograma, las direcciones horizontales y verticales (si las hubiera) y la variable utilizada.
> gstatV<-gstat(id="V",formula=V1,locations=X+Y,data=wlm) > varOmni<-variogram(gstatV,width=10,cutoff=101) > varOmni np dist gamma dir.hor dir.ver id 1 565 7.291342 42743.67 0 0 V 2 2072 15.022197 67877.29 0 0 V 3 2948 24.783924 79062.05 0 0 V 4 3210 34.757173 94338.18 0 0 V 5 4044 44.673417 88377.42 0 0 V 6 4265 54.887742 94888.71 0 0 V 7 4926 64.548384 92944.57 0 0 V 8 5196 74.614543 94322.57 0 0 V 9 5533 84.724877 89014.25 0 0 V 10 5167 94.880575 98948.24 0 0 V 11 707 100.470981 86139.86 0 0 V
La gura 24 en la pagina siguiente muestra los variogramas omnidireccionales con lags de tamano 5, 10, 15 y 20. Se descartan los extremos y entre el de 10 y el de 15 se opta por el de 10 por presentar equilibrio entre el numero de lags y la continuidad buscada. El parametro que habra que modicar sera tal vez la distancia maxima computada, para evitar el salto que aparece, dejandola en 80 metros en sucesivos calculos.
5.2. Variograma supercial

Estableciendo la opcion map=TRUE se obtiene el variograma supercial, que puede ser graado con los comandos levelplot y contourplot (gs. 25 y 26). Estas image nes pueden servir para conocer la existencia de anisotrop a geom etrica, mostrando los ejes de maxima y m nima continuidad. Este m etodo puede ser util para una primera aproximacion, pero es dif cil de cuanticar. Por esta razon se usaran los variogramas direccionales.
5.3. Variogramas direccionales

5.3.1. Busqueda de los ejes de anisotrop a En lugar de tomar todos los pares para calcular el variograma, se puede establecer una direccion determinada y una tolerancia angular para que solo se computen los
31
1e+05
Semivariograma
8e+04
Semivariograma
6e+04
4e+04
2e+04
2e+04
4e+04
6e+04
8e+04
1e+05
t=5
t = 10
0e+00
20
40
60 Distancia
80
100
0e+00
20
40
60 Distancia
80
100
1e+05
Semivariograma
8e+04
Semivariograma
6e+04
4e+04
2e+04
2e+04
4e+04
6e+04
8e+04
1e+05
t = 15
t = 20
0e+00
20
40
60 Distancia
80
100
0e+00
20
40
60 Distancia
80
100
Figura 24: Variogramas omnidireccionales (i))
200000
180000 100 160000
50
140000
map.dy
120000 0
100000
50 80000
60000 100
40000
100
50
50
100
map.dx
Figura 25: Mapa del variograma supercial
32
100
50
map.dy
50
100
100
50
50
100
map.dx
Figura 26: Isol neas del variograma supercial pares de vectores en ese rango de direcciones. Comparando diferentes variogramas direccionales se pueden obtener las direcciones de maxima y m nima continuidad . En este trabajo se van a graar los variogramas cada 15 grados sexagesimales, agrupados en dos gracos para que facilitar la lectura (g 27). Para obtener num ericamente las direcciones de los ejes se ha procedido del siguiente modo: se establece como valor de meseta 90.000ppm, a continuacion se halla la distancia a la que se corresponde dicho valor de semivarianza interpolando linealmente los valores alrededor de dicho valor. Por ultimo se ha creado un graco donde en el eje de las ordenadas se presentan las direcciones y en las abscisas la distancia obtenida. Se comprueba el gran salto entre los grados 135 y 140 debido a que el resalte que se observa en los variogramas desciende y por tanto la distancia avanza de forma abrupta. En cualquier caso, los valores obtenidos, unos 160 grados de maximo y en torno a los 90 grados de m nimo, se asemejan a los obtenidos en los apuntes de la asignatura.
80000 100000
Semivariograma
Semivariograma
60000
6e+04
8e+04
1e+05
20000
20
40
60 Distancia
80
100
0e+00
2e+04
0 30 60 90 120 150 180
15 45 75 105 135 165
40000
4e+04
20
40
60 Distancia
80
100
Figura 27: Variogramas direccionales
33
20
40
60
80
100
120
140
160
180
Distancia en lag = 4
30 0
40
50
60
20
40
60
80
100
120
140
160
180
Direccin
Figura 28: Deteccion de ejes de anisotrop a de la tolerancia angular 5.3.2. Obtencion Se pretende en este caso encontrar la tolerancia angular mas pequena posible que haga el variograma representativo y mas adaptado por tanto a la direccion elegida. Se trata en denitiva de un proceso repetitivo como el anterior, pero esta vez variaran las tolerancias y se buscaran los numeros de pares en el primer lag (el que menos pares suele presentar) mayor a 30 que a su vez generen un variograma que sea continuo (gs 29 y 30). Se observa que en los primeros variogramas la direccion de m nima continuidad presenta un ((efecto hueco)) muy importante en la distancia 60, en cualquier caso parece que 40 es una tolerancia suciente.
5.4. Variogramas cruzados

Estos variogramas presentan la continuidad espacial entre variables. Se pueden calcular por tanto los variogramas omnidireccionales y en las direcciones obtenidas anteriormente. La gura 31 presenta estos tres variogramas y se puede apreciar como tanto el omnidireccional como el de maxima continuidad no presentan grandes discontinuidades, el de m nima continuidad ofrece un aspecto poco claricador debido seguramente a la existencia de datos anomalos en alguna o ambas variables.
34
tol = 15
tol = 20
1e+05
80000 100000
424 396
613
696 729 537 1174 483 704 280 115 458 90 160 828 1065 694 1473 982 799 1598 1637 903 1481 902
608 876 1092 530 588 1165 1254 754
1121 640
8e+04
587 464 487 303 234
829
Semivariograma
Semivariograma 90 160 40 60 Distancia 80 100
6e+04
4e+04
2e+04
73
0e+00
20
20000
40000
60000
20
40
60 Distancia
80
100
tol = 25
120000
tol = 30
1e+05
968 783 822 891 1771 1473 503 883 546 305 118 90 160 985 1284 1089 1071 1889 1115 1940 Semivariograma 1090 1806 80000 1006 911 1296 1151
1062 2258 1193 2100 1295 2302 1464
8e+04
1773 1562
2052 1392
Semivariograma
585 1094 723 320 151
6e+04
4e+04
40000
60000
2e+04
0e+00
20
40
60 Distancia
80
100
20000
90 160
20
40
60 Distancia
80
100
Figura 29: Variogramas por tolerancias (i)
35
tol = 35
tol = 40
1e+05
1e+05
1171 1076 1361 744 1205 856 329 161 1752 1383
1368 2365 1993 1513 2485 1624
2415 1539 2677 1624 Semivariograma
1368 1146 1582 895 1309 980 334 179 1904 1563
1520 2656 1822 2855 1769 2941 1884
2728 1816
8e+04
8e+04
2271
Semivariograma
6e+04
4e+04
2e+04
90 160
4e+04
6e+04
0e+00
20
40
60 Distancia
80
100
0e+00
2e+04
90 160
20
40
60 Distancia
80
100
tol = 45
tol = 50
80000 100000
80000 100000
1502 1812 1362 1708 2255 1072 1495 1058 343 203 90 160 1847 2439 2949 1998 3103 2072 3284 2181
3002 2011
1626 2146 1627 1209 1114 359 230 90 160 2050 1856 2528 1716 2710 3143 2228 3354 2284 3564 2476
3230 2293
Semivariograma
Semivariograma 60 80 100
60000
40000
20000
20
40
20000
40000
60000
20
40
60 Distancia
80
100
Distancia
Figura 30: Variogramas por tolerancias (ii)
Semivariograma
2e+05
3e+05
4e+05
5e+05
Cruzado 85 160
0e+00
1e+05
20
40
60 Distancia
80
100
120
Figura 31: Variogramas cruzados
36
5.5. Scripts
Listado 10: Funciones para imprimir variogramas
1 #FUNCIONES PARA MOSTRAR VARIOGRAMAS 2 #-------------Plotear el primer 3 4 41 42
6 7 8 9 #-------------A nadir nuevos variogramas 10 addvar<-function(var,color,sepX,sepY,nps 11 12
variograma plotvar<-function(var,titulo,color,sepX, sepY,limY,limX,nps){ plot(var$dist,var$gamma,type="b", ylim=limY,xlim=limX,xlab=" Distancia",ylab="Semivariograma" ,col=color,main=titulo) if (nps) text(var$dist+sepX,var$ gamma-sepY,var$np,col=color) title(titulo) }
43 44 45 46 47 }
for (i in 2:nvars){ addvar(listvars[[i]],colores [i],seps2[i,1],seps2[i ,2],nps) } } #A nadir la leyenda legend(posley[1],posley[2],etiqs, colores)

1 #TEMA 5. CONTINUIDAD ESPACIAL 2 ps.options(family="Bookman",pointsize 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
=15) rm(list=ls()) library(lattice) library(gstat) #A nadir las funciones programadas source("Variogramas.R") #------------#Cargar los datos wlm data(walker) wlm<-walker attach(wlm)
){ points(var$dist,var$gamma,type="b", col=color) if (nps) text(var$dist+sepX,var$ gamma-sepY,var$np,col=color)
13 } 14 15 #-------------Funci on principal de
grafiado
16 pintavars<-function(listvars,colores, 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37
etiqs,seps,titulo,posley,nps=FALSE){ #Obtener los l mites nvars<-length(listvars) # Crear el objeto gstat maxtemp<-array(0,nvars) gstatV<-gstat(id="V",formula=V1, for (i in 1:nvars){ locations=X+Y,data=wlm) maxtemp[i]<-max(listvars[[i]]$ 20 gamma) 21 #Obtener variogramas omni. a varias } distancias limY<-c(0,max(maxtemp)*1.1) 22 secDists<-c(5,10,15,20) 23 colores<-c("red","green","blue") maxtempA<-array(0,nvars) 24 postscript("imgs/05/omni %1d.ps",onefile= maxtempB<-array(0,nvars) FALSE) for (i in 1:nvars){ maxtempA[i]<-min(listvars[[i]]$ 25 for (i in 1:length(secDists)){ 26 si<-secDists[i]; dist); letra<-as.character(si); maxtempB[i]<-max(listvars[[i]]$ 27 28 varOmni<-variogram(gstatV,width=si, dist) cutoff=101) } pintavars(list(varOmni),"black",as. limX<-c(min(maxtempA),max(maxtempB)* 29 expression(substitute(t==sii, 1.1) list(sii=si))),c(3,4000),"",c (40,2e4)) #Generar la matriz de separaciones } seps2<-matrix(seps,ncol=2,byrow=TRUE 30 31 dev.off() )
32
38 39 40
#Pintar el primer variograma 33 #Obtener el mapa superficial plotvar(listvars[[1]],titulo,colores 34 varOmni.map<-variogram(gstatV,width=10, [1],seps2[1,1],seps2[1,2],limY, cutoff=150,map=TRUE) limX,nps) 35 varOmni.map.df<-as.data.frame(varOmni. map) if (nvars>=2){ 36 postscript("imgs/05/omnimap1.ps") #A nadir el resto 37 grises<-grey(seq(0.99,0.3,l=30))
37
38 lp<-levelplot(map.Vmap.dx*map.dy,
39 40 41 42
43 44 45 46 #Crear familias de variogramas por 47 48 49 50 51 52 53 54 55 56
varOmni.map.df,col.regions=grises, aspect=mapasp(varOmni.df)) print(lp) dev.off() postscript("imgs/05/omnimap2.ps") cp<-contourplot(map.Vmap.dx*map.dy, varOmni.map.df,labels=FALSE,cuts=7, aspect=mapasp(varOmni.df)) print(cp) dev.off()
84 85 86 87
88 89
57 58 59 60 61 62 63 64 65 66 67
68
69 70 71 #Buscar direcciones de los ejes de 72 73 74 75 76 77 78 79 80 81 82 83
direcciones sec<-seq(0,180,5) separaciones<-rep(c(4,2000),length(sec)) posleyenda<-c(60,60000) 95 varsDirs<-list() 96 varsDirsP1<-list();sec1<-list() 97 varsDirsP2<-list();sec2<-list() 98 for (j in 1:length(sec)){ 99 sj<-sec[j] 100 101 letra<-as.character(sj) varsDirs[[letra]]<-variogram(gstatV, width=10,cutoff=100,alpha=sj,tol 102 103 .hor=40) 104 if (!sj % % 30) { varsDirsP1[[letra]]<-varsDirs[[ 105 106 letra]]; sec1[[letra]]<-letra 107 } if (!sj % % 15 && sj % % 30) { varsDirsP2[[letra]]<-varsDirs[[ 108 letra]]; 109 sec2[[letra]]<-letra } 110 } postscript("imgs/05/dirs %1d.ps",onefile= FALSE) pintavars(varsDirsP1,rainbow(length(sec1 111 )),as.character(sec1),separaciones," 112 113 ",posleyenda) pintavars(varsDirsP2,rainbow(length(sec2 114 )),as.character(sec2),separaciones," 115 ",posleyenda) 116 dev.off()
117
90 91 92 93 94 #Crear familias de variogramas por
} npMat<-cbind(grad,dist4) postscript("imgs/05/ejes.ps") plot(npMat,type="b",pch=21,bg="lightgray ",xlab="Direcci on",ylab="Distancia en lag = 4",axes=FALSE) axis(1,at=seq(0,180,10),labels=FALSE); axis(1,at=seq(0,180,20)) axis(3,at=seq(0,180,10),labels=FALSE); axis(3,at=seq(0,180,20)) axis(2);box() dev.off()
tolerancias sec<-c(90,160) lsec<-as.character(sec) tols<-seq(15,50,5) colores<-c("red","green") separaciones<-rep(c(4,2000),length(sec)) posleyenda<-c(60,40000) postscript("imgs/05/tols %1d.ps",onefile= FALSE) for (i in 1:length(tols)){ vars<-list() #Sengudo bucle por direcciones for (j in 1:length(sec)){ sj<-sec[j];letra<-as.character( sj) vars[[letra]]<-variogram(gstatV, width=10,cutoff=100,alpha=sj ,tol.hor=tols[i]) } letrat<-as.expression(substitute(tol ==t,list(t=tols[i]))) pintavars(vars,colores,lsec, separaciones,letrat,posleyenda, TRUE) } dev.off() #Variogramas cruzados #Para la U hay que eliminar localizaciones sin dato wlm2<-subset(wlm,!is.na(U)) gVU<-gstat(id="V",formula=V1,locations= X+Y,data=wlm) gVU<-gstat(gVU,id="U",formula=U1, locations=X+Y,data=wlm2) #Ver los variogramas separaciones<-rep(c(4,2000),3) posleyenda<-c(60,150000) colores<-c("black","red","green") etiqs<-c("Cruzado","85","160") varsCruz<-list() varsCruz[[1]]<-subset(variogram(gVU, width=10),id=="V.U") varsCruz[[2]]<-subset(variogram(gVU, width=10,alpha=85,tol.h=40),id=="V.U ") varsCruz[[3]]<-subset(variogram(gVU, width=10,alpha=160,tol.h=40),id=="V.
118 anisotrop a nVars<-length(varsDirs) grad<-sec 119 dist4<-rep(0,nVars) 120 yc<-90000 121 for (i in 1:nVars){ 122 j<-1; 123 while(varsDirs[[i]]$gamma[j]<yc){j<- 124 j+1}; 125 y2<-varsDirs[[i]]$gamma[j]; y1<-varsDirs[[i]]$gamma[j-1]; 126 x2<-varsDirs[[i]]$dist[j]; x1<-varsDirs[[i]]$dist[j-1]; dist4[i]<-(yc-y1)*(x2-x1)/(y2-y1) + 127 x1
38
U")
128 postscript("imgs/05/cross.ps") 129 pintavars(varsCruz,colores,etiqs,
130 dev.off()
separaciones,"",posleyenda)
39
del variograma experimental Tema 6 Modelizacion

Con el paquete gstat se puede modelizar el variograma a partir de los modelos mas comunes (efecto pepita puro, exponencial, esf erico, gaussiano) y muchos mas. Para verlos, basta con ejecutar la orden show.vgms()
0.0 0.5 1.0 1.5 2.0 2.5 3.0 0.0 0.5 1.0 1.5 2.0 2.5 3.0
vgm(1,"Nug",0)
vgm(1,"Exp",1)
vgm(1,"Sph",1)
vgm(1,"Gau",1)
vgm(1,"Exc",1)
3 2 1 0
vgm(1,"Mat",1)
3 2 1
vgm(1,"Cir",1)
vgm(1,"Lin",0)
vgm(1,"Bes",1)
vgm(1,"Pen",1)
semivariance
vgm(1,"Per",1)
vgm(1,"Hol",1)
vgm(1,"Log",1)
vgm(1,"Pow",1)
vgm(1,"Spl",1)
3 2 1 0
vgm(1,"Err",0)
3 2 1 0 0.0 0.5 1.0 1.5 2.0 2.5 3.0
vgm(1,"Int",0)
distance
Figura 32: Modelos de variogramas disponibles Para crear un modelo basta con utilizar la orden v.model. Una vez creado un modelo se puede anadir al objeto de tipo gstat para incluirlo junto con el variograma experimental. Por ejemplo si deseamos un modelo esf erico bastar a con indicar la meseta y el rango. Este modelo y el variograma se presentan gracamente en la gura 33 en la pagina siguiente.
> vgm(92000,"Sph",30) model psill range 1 Sph 92000 30
Pero los modelos pueden combinarse y formar ((estructuras imbricadas)) simplemente anadiendo su denicion unos sobre otros. As , podemos anadir crear un modelo o para el variograma direccional de 160 de la Variable V segun la siguiente denicion: (h) = 22,000 + 40,000 Esf30 (h) + 45,000 Esf150 (h) Crea el modelo de la gura 34 y su denicion en R es la siguiente:
> vgm(40e3,"Sph",30,add.to=vgm(22000,"Nug",add.to=vgm(45e3,"Sph",150))) model psill range 1 Sph 45000 150 2 Nug 22000 0 3 Sph 40000 30
(6)
40
80000
Semivariograma
60000
40000
20000
20
40
60
80
100
Distancia
Figura 33: Modelos esf erico de rango 30 y meseta parcial de 92000ppm
80000
Semivariograma
60000
40000
20000
20
40
60
80
100
Distancia
Figura 34: Modelos de variograma combinados
41
automatizada del modelo 6.1. Estimacion

6.1.1. gstat Tanto el paquete gstat como el paquete geoR proporcionan funciones para la estimacion automatizada del modelo de variograma. En ambos casos deberemos ajustar manualmente unos parametros iniciales para pasarlos como parametros de entrada, junto con el variograma. De este modo, si queremos modelar el variograma omnidireccional de V mediante un modelo exponencial con efecto pepita, basta con proceder del siguiente modo (g. 35):
> m<-vgm(72e3,"Exp",15,22e3) > mfit<-fit.variogram(variogram(g["V"],width=10),m) > g<-gstat(g,id="V",model=mfit) > print(g) data: V : formula = V1 ; locations = X + Y ; data dim = 470 x 6 variograms: model psill range V[1] Nug 186.3867 0.00000 V[2] Exp 93809.2943 12.01007
80000
Semivariograma
60000
40000
20000
20
40
60
80
100
Distancia
Figura 35: Modelo de variograma de V ajustado
6.1.2. geoR El paquete geoR proporciona una herramienta ciertamente interesante. Se trata de estimar de forma manual el modelo pero haciendo uso de una interfaz graca de usuario escrita en el lenguaje TclTk. Se pasa a la funcion eyefit el variograma a modelar y aparece un cuadro de dialogo donde es posible elegir el tipo de modelo y sus parametros mediante barras deslizantes.
> library(geoR) ------------------------------------------------------------Functions for geostatistical data analysis For an Introduction to geoR go to http://www.est.ufpr.br/geoR geoR version 1.5-7 (built on 2005/06/07) is now loaded ------------------------------------------------------------> geoV<-as.geodata(wlm,2:3,4)
42
> vario.b<-variog(geoV,max.dist=130,breaks=seq(0,130,10)) variog: computing omnidirectional variogram > vario.m<-eyefit(vario.b) Loading required package: tcltk
Figura 36: Modelado interactivo del variograma GeoR ademas genera un graco de resumen de los datos bastante descriptivo donde muestra la localizacion de los mismos, un histograma de la variable y dos gracos de dispersion de la variable contra las coordenadas X e Y respectivamente (g. 37). Por ultimo, geoR tambi en ajusta por m nimos cuadrados el modelo de variograma de forma analoga a gstat aunque algo mas simple, pasando unicamente un valor de meseta, un rango y un m etodo (g. 38 en la pagina siguiente). En este sentido gstat es mas potente al pasar a su funcion directamente un objeto de tipo variogram.model.
> vario.m2<-variofit(vario.b,c(93809,12),cov.model="exponential") variofit: weights used: npairs variofit: minimisation function used: optim > vario.m2 variofit: model parameters estimated by WLS (weighted least squares): covariance model is: exponential parameter estimates: tausq sigmasq phi 14873.1088 79431.7551 13.4973 variofit: minimised weighted sum of squares = 612531182238
43
300
250
Y Coord 150 200
100
50
50
100 150 X Coord
200
250
0 0
50
100
Coord Y 150 200
250
300
500 data
1000
1500
1500
1000
50
100 150 Coord X
200
250
0 0
20
Frequency 40 60 80
data
500
100
500 data
1000
1500
Figura 37: Descripcion de geoR del conjunto de datos
Semivariograma
0e+00 0
2e+04
4e+04
6e+04
8e+04
1e+05
20
40
60 Distancia
80
100
Figura 38: Modelo ajustado por geoR
44
6.2. Scripts
N DEL VARIOGRAMA 1 # TEMA 6 MODELIZACIO 2 ps.options(family="Bookman",pointsize =15) 3 rm(list=ls())
4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 31 var<-variogram(g["V"],alpha=160,tol.h
=40,width=10)
32 print(plot(var,model=m,xlab="Distancia",
ylab="Semivariograma",pch=3,bg="grey ")) 33 dev.off()

34 35 #Ajustar al variograma omnidireccional
library(lattice) library(gstat) #A nadir las funciones programadas source("Variogramas.R")
de V, un modelo ajustado "a ojo"

36 m<-vgm(72e3,"Exp",15,22e3) 37 mfit<-fit.variogram(variogram(g["V"], 38 39 40 41
26 27 28 29 30
width=10),m) g<-gstat(g,id="V",model=mfit) #Cargar los datos wlm print(g) data(walker) postscript("imgs/06/mfitV.ps") wlm<-walker;rm(walker) print(plot(variogram(g["V"]),model=g$ attach(wlm) model$V,xlab="Distancia",ylab=" # Crear el objeto gstat Semivariograma",pch=3,bg="grey")) g<-gstat(id="V",formula=V1,locations=X 42 dev.off() +Y,data=wlm) 43 44 # #Usando geoR para estimar y ajustar el #Ver modelos de variograma variograma postscript("imgs/06/modelos.ps") 45 library(geoR) print(show.vgms()) 46 geoV<-as.geodata(wlm,2:3,4) dev.off() 47 postscript("imgs/06/geoV.ps") 48 plot(geoV) #Modelo esf erico 49 dev.off() postscript("imgs/06/esferico.ps") 50 vario.b<-variog(geoV,max.dist=100,breaks print(plot(variogram(g["V"]),model=vgm =seq(0,130,10)) (92000,"Sph",30),xlab="Distancia", 51 # vario.m<-eyefit(vario.b) ylab="Semivariograma",pch=3,bg="grey 52 vario.m2<-variofit(vario.b,c(93809,12), ")) cov.model="exponential") dev.off() 53 postscript("imgs/06/geoVfit.ps") 54 plot(vario.b,pch=21,bg="grey",xlab=" #A nadir diferentes modelos. Distancia",ylab="Semivariograma"); postscript("imgs/06/imbr.ps") lines(vario.m2) m<-vgm(40e3,"Sph",30,add.to=vgm(22000," 55 dev.off() Nug",add.to=vgm(45e3,"Sph",150)))
45
Tema 8 Kriging
Denotado en la literatura en espanol como ((krigeado)), este conjunto de m etodos desarrollado por Denis Krige en los anos 50 del sigo XX es el m etodo geoestad stico mas ampliamente aceptado. El n ultimo del estudio y modelizacion del variograma es su aplicacion en este m etodo para el calculo del valor de la variable en puntos arbitrarios (generalmente una malla regular) mediante esta familia de m etodos.
8.1. wlc
Antes de empezar a modelar los datos de wlm, se presentan las distribuciones de U y V en wlc para poder compararlas con los modelos siguientes.
1600 5000
1400 250 250 4000
1200
200 1000
200
3000
150
800
150
600 100 100
2000
400 1000 50 200 50
0 50 100 150 200 50 100 150 200
(a) V
(b) U
Figura 39: Conjunto de datos wlc
8.2. Krigeado Ordinario (KO)

Para generar el krigeado antes necesitamos un conjunto de puntos sobre el que calcular el krigeado. Para realizar las pruebas se ha empleado una malla regular de unas 5000 celdas y para obtener los resultados nales se han empleado las localizaciones del conjunto de datos wlc (78000).
wlmGrid<-makegrid(wlm$X,wlm$Y,10000) names(wlmGrid)<-c("X","Y")
Ya solo falta ejecutar el comando krige con la formula necesaria para el KO V1, pasando el conjunto de datos de muestra, el conjunto de datos destino y el modelo del variograma. El comando genera un data.frame con las coordenadas de las localizaciones, la prediccion y su error (desviacion t pica). En las guras 45(a) y 46(a) se muestran ambos conjuntos de datos.
KO<-krige(V1,X+Y,wlm,wlc,g$model$V)
Si comparamos los datos con el conjunto wlc restando localizacion a localizacion comprobamos que las diferencias en general son pequenas y sim etricas ( 40 en la pagina siguiente).
46
30000
29491
Frecuencia
20000
25000
15000
15587 12274 8762
5000
10000
4409 3190 1900 30 150 610 1128 326 99 29 8 6 2
500
0 Error
500
1000
Figura 40: Error del Krigeado ordinario
8.3. Krigeado Universal (KU)

El krigeado Universal introduce un modelo de tendencia al considerar que la media local para cada localizacion no es constante. Basta con cambiar la formula de la orden krige a VX+Y para que se ejecute esta variedad del krigeado.
KU<-krige(VX+Y,X+Y,wlm,wlc,g$model$V)
8.4. Krigeado por bloques (KUB)

Este sistema calcula la media local discretizando una zona (bloque) en puntos individuales en lugar de utilizar la localizacion. Se puede utilizar tanto el el krigeado ordinario como en el krigeado universal y se indica unicamente pasando un nuevo parametro con el tamano del bloque en un vector de una a tres dimensiones (dependiendo del numero de variables espaciales de la muestra) o un data.frame con una a tres columnas indicando los puntos que deniran un bloque de forma irregular. Se ha realizado el krigeado por bloques universal de la variable V, con bloques de 10x10 metros y se han obtenido igualmente la prediccion y la varianza (guras 45(c) y 46(c)). Las diferencias son similares a las de los otros dos m etodos (42).
KBU<-krige(VX+Y,X+Y,wlm,wlc,g$model$V,block=c(10,10))
8.5. Krigeado Local (KUL)

Este krigeado utiliza un numero determinado de puntos para la estimacion de la localizacion: Se puede especicar numero maximo de puntos (nmax). O una distancia maxima (maxdist). Si se especica tambi en un numero m nimo (nmin) y el numero de puntos localizados por maxdist es menor, se genera un valor nulo.
47
30000
30154
Frecuencia
20000
25000
15000
15624
10000
11348 8874
5000
4405 3206 1919 618 22 142 1160 383 101 29 8 6 2
500
0 Error
500
1000
Figura 41: Error del Krigeado universal
30000
29782
Frecuencia
15000
20000
25000
15307 11473 9198
5000
10000
4673 3318 1869 16 111 584 1207 335 89 26 6 5 2
500
0 Error
500
1000
Figura 42: Error del Krigeado universal por bloques
48
Si se utilizan los parametros nmax y distmax operan ambos criterios Este m etodo tambi en se puede combinar con el krigeado ordinario o universal y el krigeado puntual o por bloques. En este trabajo se ha realizado el krigeado local universal de la variable V indicando una distancia maxima de 15 metros y un numero m nimo de 4 puntos (gs 45(d) y 46(d)).
KUL<-krige(VX+Y,X+Y,wlm,wlc,g$model$V,nmin=5,maxdist=30)
30000
29782
Frecuencia
15000
20000
25000
15307 11473 9198
5000
10000
4673 3318 1869 16 111 584 1207 335 89 26 6 5 2
500
0 Error
500
1000
Figura 43: Error del Krigeado local universal
8.6. Cokrigeado (CKO)

Tal y como se vio en el tema 3, las variables U y V estan correladas (aunque menos en wlm que en wle), Esto signica que podemos obtener informacion de una de ellas teniendo en cuenta tambi en los valores de la otra. Este principio es el que utiliza el cokrigeado, en el cual para estimar una variable se utilizan tanto los variogramas de ambas como el cruzado. gstat utiliza otro comando para realizar este m etodo. Ademas impone algunas restricciones para dar por bueno el modelo de corregionalizacion: Las variables deben encontrarse en todas las localizaciones. El rango en los modelos de los variogramas debe ser el mismo. El proceso es un poco diferente al de los anteriores, primero se crea el objeto gstat y se anaden las localizaciones, a continuacion se anaden los tres modelos iniciales. Despu es, con el comando fit.lmc se genera el modelo de corregionalizacion lineal ajustado. Se obtiene tambi en el variograma experimental del objeto gstat (el cual incluye las dos variables y el cruzado). Por ultimo ya se puede ejecutar el comando predict que devuelve el mismo tipo de data.frame que el comando krige pero en lugar de contener la prediccion y la varianza para una variable contiene las dos predicciones y varianzas y la covarianza entre U y V. Las guras 48(a) y 49(a) presentan la prediccion y la varianza de la variable U respectivamente.
49
#Crear el objeto gstat y a nadir las localizaciones g<-gstat(id="U",formula=U1,locations=X+Y,data=wlm2) g<-gstat(g,id="V",formula=V1,locations=X+Y,data=wlm2) #Crear los modelos mv<-vgm(72e3,"Exp",15,22e3) mu<-vgm(26e4,"Exp",15,26e4) muv<-vgm(.6e5,"Exp",15,3e5) #A nadir los modelos al objeto g<-gstat(g,id="U",model=mu) g<-gstat(g,id="V",model=mv) g<-gstat(g,id=c("U","V"),model=muv) #Calcular el variograma y el modelo ajustado x<-variogram(g,cutoff=100) g.fit=fit.lmc(x,g) #Calcular el cokrigeado ordinario CKO<-predict(g.fit, newdata = wlc)
Se ha comprobado de igual forma a los ejemplos anteriores el resultado con wlc y se ha obtenido el krigeado ordinario de U (UKO) para compararlo tambi en (gs 48(b) y 49(b)).
50
Frecuencia
20000
30000
40000
50000
50962
14892
10000
7493 2600 1131 492 195 96 62 27 12 6
22
2000
4000 Error
6000
8000
(a) Error del CKO
42942 40000 Frecuencia 30000
10000
20000
22081
8044 2657 1077 435 19995 50 23 7 4 1 3 0 2 1 1 0 0 1
3 22 59 294
2000
2000
4000 Error
6000
8000
10000
(b) Error del UKO
Figura 44: Diferencias con wlc de la modelizacion de U
51
8.7. Resultados
de V 8.7.1. Modelizacion A continuacion se muestran los resultados de los m etodos de krigeado empleados para modelar la variable V, as como sus predicciones de error. Por ultimo se comparan se muestran los estad sticos de las diferencias entre los valores estimados y los valores reales (wlc). En el diagrama de cajas (g. 47 en la pagina 53) se aprecia como el krigeado por bloques parece que ofrece el resultado con menos outliers al contrario que el krigeado local. Ademas el krigeado local no consigue dar valores a todos los puntos, pese a indicar una distancia maxima de 30 metros. El krigeado por bloques ofrece diferencias tanto maximas como m nimas mas pequenas que el resto (cuadro 2 en la pagina siguiente).
1600 1600
1400 250 250
1400
1200
1200
200 1000
200 1000
150
800
150
800
600 100 100
600
400
400
50 200
50 200
0 50 100 150 200 50 100 150 200
(a) KO
1600
(b) KU
1600
1400 250 250
1400
1200
1200
200 1000
200 1000
150
800
150
800
600 100 100
600
400
400
50 200
50 200
0 50 100 150 200 50 100 150 200
(c) KUB
(d) KUL
Figura 45: Prediccion en la modelizacion de V (wlm)
de U 8.7.2. Modelizacion Viendo la prediccion del cokrigeado, se observa que el resultado no es demasiado satisfactorio, seguramente porque la modelizacion del variograma, al restringir al
52
350
350
300 250 250
300
250 200 200
250
200
200
150 150
150 150
100 100
100 100
50 50
50 50
0 50 100 150 200 50 100 150 200
(a) KO
350
(b) KU
350
300 250 250
300
250 200 200
250
200
200
150 150
150 150
100 100
100 100
50 50
50 50
0 50 100 150 200 50 100 150 200
(c) KUB
(d) KUL
Figura 46: Desviacion t pica en la modelizacion de V (wlm)
Cuadro 2: Estad sticos de los errores en los m etodos de krigeado de V Estad sticos KU KUB 4.35 4.47 -684.74 -645.48 -72.55 -76.87 25.94 28.71 84.83 86.51 967.42 913.93 21344.87 21325.76 146.10 146.03 -0.34 -0.32 1.47 1.23 33.56 32.68
Media M nimo Q1 Q2 Q3 Max. Var. Desv. CS K CV
KO 5.77 -686.11 -71.55 28.6 89.31 967.39 21394.29 146.27 -0.39 1.45 25.34
KUL -4.13 -928.96 -80.85 12.572 77.19 969.61 22621.25 150.40 -0.32 1.86 -36.42
53
KO 1000
KU
KUB
KUL
500
500
1000
Figura 47: Diagramas de caja y bigote de las diferencias
54
mismo rango a los modelos tanto de U, como de V y el cruzado, a generado una peor modelizacion nal. As y todo, el cokrigeado a dado valores maximos y m nimos menores al ordinario y el resto de parametros estad sticos no son muy diferentes (cuadro 3 en la pagina siguiente).
5000 5000
250 4000
250 4000
200
200
3000
3000
150
Y
2000
150
2000 100
100
1000 50 50
1000
0 50 100 150 200 50 100 150 200
(a) CKO
(b) UKO
Figura 48: Prediccion en la modelizacion de U (wlm)
800
800
700 250 250
700
600
600
200 500
200 500
150
400
150
400
300 100 100
300
200
200
50 100
50 100
0 50 100 150 200 50 100 150 200
(a) CKO
(b) UKO
Figura 49: Desviacion t pica en la modelizacion de U (wlm)
55
Cuadro 3: Estad sticos de los errores en los m etodos de krigeado de U Estad sticos CKO UKO Media -250.44 -255.68 M nimo -1366.57 -2625.11 Q1 -487.89 -506.54 Q2 -404.07 -418.417 Q3 -196.54 -158.85 Max. 8973.018 9034.47 Var. 197185.3 197649.9 Desv. 444.05 444.58 CS 3.78 3.30 K 25.17 21.42 CV -1.77 -1.74
56
8.8. Scripts
1 # TEMA 8 KRIGING 2 ps.options(family="Bookman",pointsize 46 47 #Generar el objeto gstat 48 g<-gstat(id="V",formula=V1,locations=X
+Y,data=wlm)
49 #Establecer un modelo inicial de V 50 mv<-vgm(72e3,"Exp",15,22e3) 51 #Ajustar el modelo con el variograma
=15) 3 # rm(list=ls())
4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43
library(lattice) library(gstat) #Funci on que devuelve estad sticas calcestads<-function(datosO){ Desv<-function(datos){sqrt(var(datos ))} CS<-function(V){sum((V-mean(V))3)/( length(V)*Desv(V)3)} K<-function(V){sum((V-mean(V))4/ length(V))/Desv(V)4-3} CV<-function(V){Desv(V)/mean(V)} datos<-datosO[!is.na(datosO)] data.frame(media=mean(datos),Q0=min( datos), Q1=as.double(quantile(datos,0.25,na. rm=TRUE)), Q2=median(datos), Q3=as.double(quantile(datos,0.75,na. rm=TRUE)), Q4=max(datos), var=var(datos),desv=Desv(datos), CS=CS(datos),K=K(datos),CV=CV(datos) ) } #Niveles y secuencias niv<-seq(0,1600,200) niv2<-seq(0,350,50) nivH<-seq(-700,1000,100) nivH2<-seq(-1000,1000,100) gris1<-gray(seq(0.9,0.3,l=30)) gris2<-gray(seq(0.99,0.5,l=7)) gris3<-gray(seq(0.9,0.3,l=8)) #Cargar los datos wlm print("Cargando datos...") data(walker) wlm<-walker;rm(walker) attach(wlm) #Importar e imprimir los valores reales wlc wlc<-read.csv("wlc.csv",sep=";",dec=".", header=TRUE) postscript("imgs/08/wlc %01d.ps",onefile= FALSE) print(levelplot(VX+Y,wlc,at=niv,aspect= mapasp(wlc),col.regions=gris1)) print(levelplot(UX+Y,wlc,at=seq (0,5200,400),aspect=mapasp(wlc),col. regions=gris1)) dev.off()
omnidireccional
52 mvfit<-fit.variogram(variogram(g["V"],
width=10,cutoff=100),mv)
53 #A nadir el ajuste al objeto gstat 54 g<-gstat(g,id="V",model=mvfit) 55 56 #Crear una malla regular de unas 5000 57 58 59 60 61 62 63 64 65 66 67
celdas (aprox) wlmGrid<-makegrid(X,Y,5000) names(wlmGrid)<-c("X","Y") #KRIGEADO ORDINARIO print("Calculando KO...") KO<-read.csv("KO.csv") # KO<-krige(V1,X+Y,wlm,wlc,g$model$V) # write.table(round(KO,3),file="KO.csv", sep=",",row.names=FALSE) difKO<-KO$var1.pred-wlc$V postscript("imgs/08/KO %01d.ps",onefile= FALSE) graf1<-levelplot(var1.predX+Y,KO,at=niv ,aspect=mapasp(KKULO),contour=TRUE, col.regions=gris1,labels=FALSE) graf2<-levelplot(sqrt(var1.var)X+Y,KO, at=niv2,aspect=mapasp(KO),col. regions=gris2) print(graf1);print(graf2) hist(difKO,breaks=nivH,labels=TRUE,main= "",col="lightgray",xlab="Error",ylab ="Frecuencia") dev.off()
68
69 70
71 72 73 74 75 76 77 78 79 80 81
82
83 84
44 45
#KRIGEADO UNIVERSAL print("Calculando KU...") KU<-read.csv("KU.csv") # KU<-krige(VX+Y,X+Y,wlm,wlc,g$model$V ) # write.table(round(KU,3),file="KU.csv", sep=",",row.names=FALSE) difKU<-KU$var1.pred-wlc$V postscript("imgs/08/KU %01d.ps",onefile= FALSE) graf1<-levelplot(var1.predX+Y,KU,at=niv ,aspect=mapasp(KU),contour=TRUE,col. regions=gris1,labels=FALSE) graf2<-levelplot(sqrt(var1.var)X+Y,KU, at=niv2,aspect=mapasp(KU),col. regions=gris2) print(graf1);print(graf2) hist(difKU,breaks=nivH,labels=TRUE,main= "",col="lightgray",xlab="Error",ylab ="Frecuencia")
57
85 86 87 88 89 90 91 92 93 94 95
dev.off()
96
97 98
99 100 101 102 103 104 105 106 107 108
#KRIGEADO UNIVERSAL POR BLOQUES print("Calculando KUB...") KUB<-read.csv("KUB.csv") 129 # KUB<-krige(VX+Y,X+Y,wlm,wlc,g$model$ 130 V,block=c(10,10)) 131 # write.table(round(KUB,3),file="KUB.csv 132 ",sep=",",row.names=FALSE) difKUB<-KUB$var1.pred-wlc$V 133 postscript("imgs/08/KUB %01d.ps",onefile= FALSE) 134 graf1<-levelplot(var1.predX+Y,KUB,at= 135 niv,aspect=mapasp(KUB),contour=TRUE, 136 137 col.regions=gris1,labels=FALSE) graf2<-levelplot(sqrt(var1.var)X+Y,KUB, 138 at=niv2,aspect=mapasp(KUB),col. 139 140 regions=gris2) 141 print(graf1);print(graf2) hist(difKUB,breaks=nivH,labels=TRUE,main 142 ="",col="lightgray",xlab="Error", 143 144 ylab="Frecuencia") 145 dev.off()
125 126 # COKRIGEADO 127 print("Iniciando CKO...") 128 #Rehacer el objeto gstat s olo con las
localizaciones donde U existe wlm2<-subset(wlm,!is.na(wlm$U)) g<-"null" rm(g) g<-gstat(id="U",formula=U1,locations=X +Y,data=wlm2) g<-gstat(g,id="V",formula=V1,locations= X+Y,data=wlm2) #Crear los modelos mv<-vgm(72e3,"Exp",15,22e3) mu<-vgm(26e4,"Exp",15,26e4) muv<-vgm(.6e5,"Exp",15,3e5) #A nadir los modelos al objeto g<-gstat(g,id="U",model=mu) g<-gstat(g,id="V",model=mv) g<-gstat(g,id=c("U","V"),model=muv) #Calcular el modelo ajustado x<-variogram(g,cutoff=100) g.fit=fit.lmc(x,g)#,fit.ranges=T) #Ver variogramas y modelos postscript("imgs/08/UyV.ps") print(plot(variogram(g),model=g.fit,xlab ="Distancia",ylab="Semivariograma")) dev.off()
146 147 #KRIGEADO UNIVERSAL LOCAL 148 print("Calculando KUL...") KUL<-read.csv("KUL.csv") # KUL<-krige(VX+Y,X+Y,wlm,wlc,g$model$ 149 V,nmin=5,maxdist=30) 150 # write.table(round(KUL,3),file="KUL.csv 151 ",sep=",",row.names=FALSE) 152 difKUL<-KUL$var1.pred-wlc$V 153 postscript("imgs/08/KUL %01d.ps",onefile= 154 155 FALSE) graf1<-levelplot(var1.predX+Y,KUL,at= niv,aspect=mapasp(KUL),contour=TRUE, 156 157 col.regions=gris3,labels=FALSE) graf2<-levelplot(sqrt(var1.var)X+Y,KUL, at=niv2,aspect=mapasp(KUB),col. 158 regions=gris2) print(graf1);print(graf2) hist(difKUL,breaks=nivH2,labels=TRUE, main="",col="lightgray",xlab="Error" 159 ,ylab="Frecuencia") dev.off()
109
110 111
112 113 114 #RESUMEN DE DIFERENCIAS 115 print("Calculando Resumen...") 116 diff.l<-list(KO=difKO,KU=difKU,KUB=
difKUB,KUL=difKUL)
#Ejecutar el cokrigeado print("Calculando CKO...") CKO<-read.csv("CKO.csv") # CKO<-predict(g.fit, newdata = wlc) # write.table(round(CKO,3),file="CKO.csv ",sep=",",row.names=FALSE) difCKO<-wlc$U-CKO$U.pred postscript("imgs/08/CKO %01d.ps",onefile= FALSE) graf1<-levelplot(U.predX+Y,CKO,at=seq (0,5200,400),aspect=mapasp(CKO), contour=TRUE,col.regions=gray(seq (0.99,0.3,l=30)),labels=F) graf2<-levelplot(sqrt(U.var)X+Y,CKO,at= seq(0,800,100),aspect=mapasp(CKO), col.regions=gris1) 160 print(graf1);print(graf2) 161 hist(difCKO,labels=TRUE,main="",col=" lightgray",xlab="Error",ylab=" Frecuencia") 162 dev.off()
163 117 diff.R<-lapply(diff.l,calcestads) 118 diff.Rdf<-data.frame(rbind(diff.R[[1]], 164 #Calcular el krigeaado ordinario de U
para comparar diff.R[[2]],diff.R[[3]],diff.R[[4]]) 165 print("Calculando UKO...") ) 119 row.names(diff.Rdf)<-c("KO","KU","KUB"," 166 mufit<-fit.variogram(variogram(g["U"], KUL") width=10,cutoff=100),vgm(26e4,"Exp" 120 ,15,26e4)) 121 postscript("imgs/08/boxplot.ps") 167 UKO<-read.csv("UKO.csv") 122 boxplot(diff.l,horizontal=TRUE,col=" 168 # UKO<-krige(U1,X+Y,wlm2,wlc,mufit) lightgray") 169 # write.table(round(UKO,3),file="UKO.csv 123 dev.off() ",sep=",",row.names=FALSE) 124 print(diff.Rdf) 170 difUKO<-wlc$U-UKO$var1.pred
58
,col.regions=gris1) FALSE) 174 print(graf1);print(graf2) 172 graf1<-levelplot(var1.predX+Y,UKO,at= 175 hist(difUKO,labels=TRUE,main="",col=" seq(0,5200,400),aspect=mapasp(UKO), lightgray",xlab="Error",ylab=" contour=TRUE,col.regions=gray(seq Frecuencia") (0.99,0.3,l=30)),labels=F) 176 dev.off() 173 graf2<-levelplot(sqrt(var1.var)X+Y,UKO, at=seq(0,800,100),aspect=mapasp(UKO)
171 postscript("imgs/08/UKO %01d.ps",onefile=
59
Referencias
C ARLOS y C OLL A LIAGA , E LOINA, Apuntes de geoestad [1] M AR T I NEZ L LARIO , J OS E stica basica. Aplicaciones mediante Sistemas de Informacion Geograca , UPV, Valencia, 2005. [2] P EBESMA , E DZER J., gstats user manual, Dept. of Physical Geography, Utrecht University, Utretch, Holanda, 2001. URL http://www.gstat.org [3] P EBESMA , E DZER J., Multivariable geostatistics in S: the gstat package, Computers & Geosciences, tomo 30, pags. 683691, 2004. [4] R D EVELOPMENT C ORE T EAM, Introduccion a R, R Foundation for Statistical Computing, Vienna, Austria, 2000. [5] R D EVELOPMENT C ORE T EAM, R: A language and environment for statistical computing, R Foundation for Statistical Computing, Vienna, Austria, 2004, ISBN 3900051-07-0. URL http://www.R-project.org [6] R IBEIRO , PAULO J. J R y D IGGLE , P ETER J., geoR: a package for geostatistical analysis, R-NEWS, tomo 1(2), pags. 1418, 2001, iSSN 1609-3631. URL CRAN://doc/Rnews/ [7] S ARKAR , D EEPAYAN, lattice: Lattice Graphics, 2004, r package version 0.10-16.

Geoestadistica Con R

Cargado por

Copyright:

Formatos disponibles

Geoestadistica Con R

Cargado por

Información del documento

Derechos de autor

Formatos disponibles

Compartir este documento

Compartir o incrustar documentos

Opciones para compartir

¿Le pareció útil este documento?

¿Este contenido es inapropiado?

Copyright:

Formatos disponibles

Geoestadistica Con R

Cargado por

Copyright:

Formatos disponibles

Geoestad stica con R

Jorge Gaspar Sanz Salinas Septiembre de 2005

Modelizacion geoestad stica con R

Modelizacion geoestad stica con R

10. 11. 12. 13.

Modelizacion geoestad stica con R

Modelizacion geoestad stica con R

0.4. Datos de trabajo

Modelizacion geoestad stica con R

univariada Tema 1 Descripcion

24 103 34 110 22 97 27 103 20 94 16 86 15 85 15 83 11 74 4 47

27 112 36 121 24 105 27 111 27 110 23 101 16 90 15 87 29 108 32 111 14

30 123 29 119 25 112 32 122 29 116 24 109 17 97 16 94 37 121 38 124

0 19 7 77 10 91 4 64 19 108 25 113 18 101 17 99 55 143 20 109 16 X

18 111 18 111 19 115 15 105 16 107 15 102 6 72 2 48 3 52 14 98 18

univariada 1.2. M etodos gr acos para la descripcion

Modelizacion geoestad stica con R

Figura 3: Histograma acumulado de V

Figura 4: Graco de probabilidad uniforme de V

Modelizacion geoestad stica con R

Figura 5: Graco de probabilidad normal

4.0 Sample Quantiles

Figura 6: Graco de probabilidad lognormal

Modelizacion geoestad stica con R

1.3. M etodos num ericos

Figura 7: Graco de caja y bigotes de V

Modelizacion geoestad stica con R

[1] 695.3409 [1] 26.36932 [1] 34.5

El coeciente de curtosis o apuntalamiento se calcula como K=

=15) rm(list=ls()) #Cargar wle y redondearlo wle<-read.delim2("walker10.asc") wle<-round(wle,0)

#Adjuntar los datos de wle para acceder directamente 10 attach(wle)

#Gr afico de probabilidad normal y lognormal postscript("imgs/01/probNormal.ps") qqnorm(V,datax=TRUE,pch=3,main=""); qqline(V,datax=TRUE) 47 dev.off()

Modelizacion geoestad stica con R

48 49 logV<-log(V)[is.finite(log(V))] 50 postscript("imgs/01/probLogNormal.ps") 51 qqnorm(logV,datax=TRUE,pch=3,main="");

Modelizacion geoestad stica con R

bivariada Tema 2 Descripcion

30 75% 20 25% 10 Cuartiles 0 0% 50%

Modelizacion geoestad stica con R

pueden ser susceptibles de ser eliminados.

Figura 10: Graco de dispersion

2.2. M etodos num ericos

mRx ) (Ryi mRy ) Rx Ry

> cov(U,V);cor(U,V,method="pearson");cor(U,V,method="spearman") [1] 218.2778 [1] 0.8393395 [1] 0.8576604

Modelizacion geoestad stica con R

=15) rm(list=ls()) #Cargar wle y redondearlo wle<-read.delim2("walker10.asc") wle<-round(wle,0)

#Adjuntar los datos de wle para acceder directamente 10 attach(wle)

Modelizacion geoestad stica con R

espacial Tema 3 Descripcion

77 61 74 70 88 82 80 80 84 100 Mximos Mnimos 12

103 110 97 103 94 86 85 83 74 47

112 121 105 111 110 101 90 87 108 111

123 119 112 122 116 109 97 94 121 124

19 77 91 64 108 113 101 99 143 109

111 111 115 105 107 102 72 48 52 98

Figura 11: Distribucion de V

Modelizacion geoestad stica con R

Figura 12: Mapa graduado de color de V

0 81.75 100.5 116.25 145

Figura 13: Mapa graduado de tamano de V

Modelizacion geoestad stica con R