[GeoSTA] Semivariogramas teóricos con R

Empezando a aprender conceptos básicos de geoestadística, estoy lidiando con algunos de ellos duros inicialmente, pero interesantes e importantes para entender bien la modelización espacio-temporal. Así que en estas tareas, creo que he asimilado el concepto de variograma y semivariograma (pronto publicare unos pequeños apuntes sobre conceptos básicos); y he tenido que realizar un resumen de los principales modelos teóricos de semivariogramas. La cuestión es que me apetecía hacer unas gráficas básicas sobre los modelos teóricos y en primer lugar pense en hacerlos con power point, pero luego, pense que mejor los hacía con R, y así de camino aprendía algo más.

Para ello he utilizado el paquete geoR sobre geoestadística y he conseguido representar mediante una simulación tres modelos teóricos del semivariograma: exponencial, gausiano y esférico). A continuación os dejo la gráfica que he producido y el código que he utilizado para ello.

La figura muestra la semivarianza (eje y) frente a la distancia (eje x). En el siguiente código se explica como se ha realizado tal gráfico:

plot(0:2, 0:2, type="n",axes=F, xlab="|h|", ylab=expression(gamma(h)), font.lab=2)
axis(1, at=c(0,6), labels=F)
axis(2, at=c(-.5,6), labels=F, pos=0)
lines.variomodel(cov.m="sph", cov.p =c(1, .8), nug=.06, max.dist=3, lty=1, scaled=T, col="black", lwd=2)
lines.variomodel(cov.m="exp", cov.p =c(1, .8), nug=.06, max.dist=3, lty=1, scaled=T, col="red", lwd=2)
lines.variomodel(cov.m="gau", cov.p =c(1, .8), nug=.06, max.dist=3, lty=1, scaled=T, col="blue", lwd=2)
legend("top",
expression(list("Esférico"),list("Exponencial"), list("Gaussiano")),
lty=c(1,1,1), lwd=c(2,2,2), col=c("black","red","blue"))

Los parámetros mas importantes de lines.variomodel son:

• cov.m que nos indica el modelo teórico, como por ejemplo: esférico ("sph"), exponencial ("exp"), gaussiano ("gau"), etc.
• cov.p nos indica los parámetros meseta o sill y el rango
• nug nos permite especificar el efecto "pepita" (nugget) (vamos, la separación que observarmos en el eje de ordenadas respecto del origen)

Pues nada mas por ahora, pronto vendrán los apuntes de los conceptos básicos en geoestadística.

Escribir código R en tu blog... y que quede bonico.

Pues sí, mucha gente utiliza R para llevar a cabo sus análisis, y casi todos coinciden en la potencialidad de esta herramienta. Son muchos los investigadores que lo han incorporado como una herramienta diaria en su trabajo, y otros muchos, los que además de hacer increíbles cosas en R, lo publican en sus websites o en sus blogs (ej.: L. Cayuela, B. Benito, por citar gente cercana).

Como neo-neofito en esto de R, me ha sorprendido una simple aplicación que te permite poner bonico el código que publicas en tu blog, es decir, realizas tu script para explicar cualquier análisis y luego lo publicas en tu blog, y además, queda bonito. Pues si, es posible. Esta herramienta se llama Pretty R y transforma tu código R en código html que puedes fácilmente publicar en tu blog.

Así que ya tenemos otra excusa menos para seguir poniendo a disposición de la sociedad nuestros avances y las metodologías utilizadas para conseguirlos.

Baste un ejemplo tonto:

Imaginemos que tenemos que realizar un análisis cluster, lo hacemos con R, y queremos publicar nuestro código. La diferencia entre usar o no esta herramienta la podemos ver a continuación:

A. Código si usar Pretty R:

data <- read.table("datos_variables.txt", header=T)

d1 <- dis(data, method="euclidean")

cluster1 <- hclust(d, method = "complete")

B. Código usando Pretty R:

data <- read.table("datos_variables.txt", header=T)
d1 <- dis(data, method="euclidean")
cluster1 <- hclust(d, method = "complete") 

Pues nada, juzguen y decidan.