• Author / Uploaded
• Byron Lima Morillo

• Categories
• Áreas de informática
• Informática
• Tecnología
• Informática y tecnología de la información
• Matemáticas

Citation preview

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 superficiales mediante sistemas de informacion ´ geografica ´ 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 ´ Geografica, ´ 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

3.5. Diagramas de dispersion ´ cruzados . . . . . . . . . . . 20 3.6. Scripts . . . . . . . . . . . . . 20 . . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

4 4 4 5 5

´ univariada 1. Descripcion 6 1.1. Carga y visualizacion ´ de los datos . . . . . . . . . . . . . . 6 1.2. M´etodos graficos ´ para la descripcion ´ univariada . . . 6 1.3. M´etodos num´ericos . . . . . 9 1.4. Scripts . . . . . . . . . . . . . 10 ´ bivariada 2. Descripcion 12 2.1. M´etodos graficos ´ . . . . . . . 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

´ 4. Estimacion. M´etodos deterministas 24 4.1. Scripts . . . . . . . . . . . . . 27 5. Continuidad espacial de V 5.1. Variograma omnidireccional 5.2. Variograma superficial . . . 5.3. Variogramas direccionales . 5.4. Variogramas cruzados . . . . 5.5. Scripts . . . . . . . . . . . . .

30 30 30 30 33 36

´ 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

2

´Indice de figuras 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 Grafico ´ de probabilidad uniforme de V . . . . . . . . Grafico ´ de probabilidad normal . . . . . . . . . . . . . Grafico ´ de probabilidad lognormal . . . . . . . . . . . . . Grafico ´ de caja y bigotes de V Grafico ´ de caja y bigotes de VyU. . . . . . . . . . . . . . Grafico ´ de cuantiles de V y U Grafico ´ de dispersion ´ . . . . Distribucion ´ de V . . . . . . Mapa graduado de color de V Mapa graduado de tamano ˜ de V . . . . . . . . . . . . . . Mapas de indicadores . . . . Mapa de superficie interpolada . . . . . . . . . . . . . . Media y varianza en ventana de 3x3 . . . . . . . . . . . Grafico ´ 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 superficial . . . . . . . . . . . . Isol´ıneas del variograma superficial . . . . . . . . . . .

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

´Indice de cuadros 1.

Resumen de estad´ısticos de datos de validacion ´ y estimados . . . . . . . . . . . . . 26

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.

3.

32 33 34 35 35 39

40 40 41 42 43 43 45 46 47 47 48 50 51 52 53 54 54

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

Modelizacion ´ geoestad´ıstica con R

3

´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 . . . . . . . . .

8. 10 14 22 23 27 27 28

9.

10. 11. 12. 13.

Script para ps.map del m´etodo RST . . . . . . . . . . Script para ps.map del m´etodo Pol´ıgonos de influencia . . . . . . . . . . . . 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

Modelizacion ´ geoestad´ıstica con R

4

´ 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 graficos ´ 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 graficas ´ 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 grafica ´ 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 graficos ´ lattice [7]. En definitiva, se dispone de un sistema ampliable que se maneja como una consola de entrada de comandos que permite adquirir datos desde ficheros, manipularlos, crear nuevos datos y por ultimo ´ o bien ver las graficas ´ por pantalla o mandarlas a ficheros 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 ficheros de parametros) ´ como interactiva mostrando los resultados utilizando el programa para presentacion ´ de graficos ´ 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

Modelizacion ´ geoestad´ıstica con R

5

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 graficas ´ 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 ´ geografica ´ 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 influencia (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 ficheros de texto separados por comas (CSV) para poder operar con ellos.

Modelizacion ´ geoestad´ıstica con R

6

´ univariada Tema 1 Descripcion ´ de los datos 1.1. Carga y visualizacion

242

244

Y

246

248

250

En este cap´ıtulo se van a usar los datos wle redondeados a valores enteros. El primer paso sera´ por tanto cargar el fichero walker10.asc, que es un fichero de texto separado por tabuladores que se puede importar directamente con la orden read.delim2 para a continuacion ´ redondearlo. En la figura 1 se muestra la distribucion ´ de los datos, as´ı como los valores que toman las variables U y V.

15 81

12 77

24 103

27 112

30 123

0 19

2 40

18 111

18 114

18 120

16 82

7 61

34 110

36 121

29 119

7 77

4 52

18 111

18 117

20 124

16 82

9 74

22 97

24 105

25 112

10 91

7 73

19 115

19 118

22 129

21 88

8 70

27 103

27 111

32 122

4 64

10 84

15 105

17 113

19 123

21 89

18 88

20 94

27 110

29 116

19 108

7 73

16 107

19 118

22 127

15 77

16 82

16 86

23 101

24 109

25 113

7 79

15 102

21 120

20 121

14 74

15 80

15 85

16 90

17 97

18 101

14 96

6 72

28 128

25 130

14 75

15 80

15 83

15 87

16 94

17 99

13 95

2 48

40 139

38 145

16 77

17 84

11 74

29 108

37 121

55 143

11 91

3 52

34 136

35 144

22 87

28 100

4 47

32 111

38 124

20 109

0 0

14 98

31 134

34 144

12

14

16

18

20

X

Figura 1: Distribucion ´ de U y V

´ univariada 1.2. M´etodos gr´aficos para la descripcion El m´etodo grafico ´ 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 (fig. 2 en la pagina ´ siguiente). Otro grafico ´ interesante es el histograma acumulado en el que a partir de las mismas clases del histograma anterior se muestra la suma acumulada (fig 3 en la pagina ´ siguiente). El grafico ´ de probabilidad acumulada muestra la proporcion ´ de datos para cada punto que son menores que e´ l (fig 4 en la pagina ´ siguiente) Las figuras 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.

Modelizacion ´ geoestad´ıstica con R

7

16

17 15 14

14

14

10 8 6

Frequency

12

12 11

4

4 3

3 2

2 1

1

2

1

0

0

0

50

100

150

V

60 40 0

20

Frecuencia acumulada

80

100

Figura 2: Histograma de V

0

10

20

30

40

50

60

70

80

90 100

120

140

Variable

0.6 0.4 0.2 0.0

Probabilidad acumulada

0.8

1.0

Figura 3: Histograma acumulado de V

0

50

100

150

V

Figura 4: Grafico ´ de probabilidad uniforme de V

8

1 0 −1 −2

Theoretical Quantiles

2

Modelizacion ´ geoestad´ıstica con R

0

50

100

150

Sample Quantiles

1 0 −1 −2

Theoretical Quantiles

2

Figura 5: Grafico ´ de probabilidad normal

3.0

3.5

4.0

4.5

Sample Quantiles

Figura 6: Grafico ´ de probabilidad lognormal

5.0

Modelizacion ´ geoestad´ıstica con R

9

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 definido, pero a partir de la tabla definida del corte de V (tcutV) donde se almacenan las frecuencias relativas en las clases definidas 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 dt dtˆ2;dt;as.numeric(cuantiles["75 %"]-cuantiles["25 %"])

Modelizacion ´ geoestad´ıstica con R

10

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

1.3.3. Medidas de forma El coeficiente de sesgo o asimetr´ıa (skewness) se calcula a partir de la formula ´ Pn

1 n

CS =

i=1 (xi σ3

− m)3

(1)

El coeficiente de curtosis o apuntalamiento se calcula como K=

Pn

i=1

(xi −m)4 n σ4

−3

(2)

El coeficiente de variacion ´ no es mas ´ que el cociente entre la desviacion ´ t´ıpica y la media, siendo trivial su calculo. ´ En R estos tres coeficientes se calculan como: > media print(CS print(K print(CV