a15apendice2

Apéndice 2

"Script" para realizar un ejemplo de cómo calcular la riqueza de especies con los métodos de Kriging y Co-Kriging en R. Para mayor detalle se pueden consultar las librerías gstat (Pebesma, 2004) y rgdal (Keitt et al., 2011). Los comentarios en el Script están marcados con el símbolo de #. En "C" se debe de crear una carpeta que se llame "Ejemplo k_ck" y copiar el archivo "coor_sr.csv", para después poder ejecutar el Script.

library(gstat) #Cargar librería

###########Cargar archivo##############

#Se escribe la ruta y el nombre del archivo "coor_sr.csv"

rte=read.csv("C:/Ejemplo k_ck/coor_sr.csv")

#Se muestra la estructura del objeto "rte". Es una hoja de datos

#que tiene 5 columnas y 48 filas.

str(rte)

####################Análisis Visual######################

#Cálculo de correlación entre riqueza de especies y Asteraceae, la cual es de 0.89

cor(rte$tsp,rte$asp)

#Graficar riqueza total de especies vs riqueza de Asteraceae plot(rte$asp,rte$tsp,xlab="Número de Asteraceae",ylab="Número Total de Especies", pch=20)

abline(lm(rte$tsp~rte$asp)) text(50,3000,"r=0.89")

#Evaluación de la distribución de las valores de riqueza de especies

#Se realiza una análisis visual de la distribución de la riqueza total de especies y

#Asteraceae. El histograma muestra que no existe distribución normal para

#los valores de riqueza total de especies par(mfrow=c(2,1))

hist(rte$tsp,main="",xlab="Riqueza total de especies",ylab="Frecuencia")

hist(rte$asp,,main="",xlab="Riqueza de Asteraceae",ylab="Frecuencia")

#Se aplica logaritmo a los valores de riqueza para ajustar a una distribución normal.

#En análisis posteriores los valores de riqueza total de especies serán transformados con #logaritmo. Existen otras funciones para realizar la transformación de las variables como

#raíz cuadrada y cuarta, log (1 + x), etc.

#Se grafican los valores de riqueza de especies transformados hist(log(rte$tsp),main="",xlab="Riqueza total de especies", ylab="Frecuencia (log transformación)")

###############################Selección de variograma y Kriging##########

#La hoja de datos "rte" se tranforma a "SpatialPointsDataFrame",

#formato requerido por la librería "gstat" para aplicarle sus funciones rte.coor=rte

coordinates(rte.coor)=~x+y

str(rte.coor)

###Seleccion y ajuste de semivariograma###

#Primero se evalúa la estructura espacial de los datos con la ayuda de variograma

#que grafica la semivarianza y la distancia entre pares de puntos

Variograma.Var=variogram(log(rte$tsp)~x+y,rte.coor)

#Se grafica el semivariograma para obtener los parámetros de inicio para

#hacer su ajuste (sill, range y nugget). plot(Variograma.Var)

# En el punto 4 se observa la asintota del variograma (sill).

#El valor de 4 se sustituye en los parámetros de sill y range para obtener los parámetros

# de inicio para el ajuste de modelos permisibles

sill=Variograma.Var[4,3]

range=Variograma.Var[4,2]

nugget=min(Variograma.Var[,3])

#Se generan 7 modelos permisibles con los parámetros de inicio

Esférico=vgm(sill,"Sph",range,nugget,variogramModel=Variog rama.Var)

Exponencial=vgm(sill,"Exp",range,nugget,variogramModel=Va riograma.Var)

Gaussiano=vgm(sill,"Gau",range,nugget,variogramModel=Vari ograma.Var)

Lineal=vgm(sill,"Lin",range,nugget,variogramModel=Variogra ma.Var)

Matern=vgm(sill,"Mat",range,nugget,variogramModel=Variogr ama.Var)

Bessel=vgm(sill,"Bes",range,nugget,variogramModel=Variogra ma.Var)

Pentaesf=vgm(sill,"Pen",range,nugget,variogramModel=Variog rama.Var)

#Se ajustan los modelos permisibles

fit.Esférico=fit.variogram(Variograma.Var,Esférico)

fit.Exponencial=fit.variogram(Variograma.Var,Exponencial)

fit.Gaussiano=fit.variogram(Variograma.Var,Gaussiano)

fit.Lineal=fit.variogram(Variograma.Var,Lineal)

fit.Matern=fit.variogram(Variograma.Var,Matern)

fit.Bessel=fit.variogram(Variograma.Var,Bessel)

fit.Pentaesf=fit.variogram(Variograma.Var,Pentaesf)

#Se evalúan los modelos y se selecciona el que tenga menor error, en este caso, fue el"Lineal"

attr(fit.Esférico,"SSErr")

attr(fit.Exponencial,"SSErr")

attr(fit.Gaussiano,"SSErr")

attr(fit.Lineal,"SSErr")

attr(fit.Matern,"SSErr")

attr(fit.Bessel,"SSErr")

attr(fit.Pentaesf,"SSErr")

#Se obsevan los parámetros con los que se ajustó el variograma

#sill parcial (psill)=0.574, nugget=0.078 y range=1.719

fit.Lineal

#Se grafica el variograma ajustado que fue el lineal plot(Variograma.Var, pl=F,model=fit.Lineal,

,col='T3lue",pch=20,main="Variograma ajustado con modelo Linear",xlab="Distancia",ylab="Semivarianza")

###Análisis de de validacion cruzada para calcular la precisión de3 Kriging###

###Aviso checar que no haya valores con iguales coordenadas###

#La validación cruzada se realiza empleando 10-fold para medir los errores de ajuste de3 #kriging, universal, ordinario #y simple

#Universal

U.cross=krige.cv(log(tsp)~x+y,rte.coor,fit.Lineal,nfold=10)

#Ordinario

O.cross=krige.cv(log(tsp)~1,rte.coor,fit.Lineal,nfold=10)

#Simple

S.cross=krige.cv(log(tsp)~1,rte.coor,fit.Lineal,nfold=10,beta=5)

#Transformar a hoja de datos

Ures=as.data.frame(U.cross)$residual

Ores=as.data.frame(O.cross)$residual

Sres=as.data.frame(S.cross)$residual

#Mean Error(ME)

U.ME=mean(Ures)

O.ME=mean(Ores)

S.ME=mean(Sres)

ME=c(U.ME,O.ME,S.ME)

#Root Mean Squared Error (RMSE)

U.RMSE=sqrt(mean(Ures^2))

O.RMSE=sqrt(mean(Ores^2))

S.RMSE=sqrt(mean(Sres^2))

RMSE=c(U.RMSE,O.RMSE,S.RMSE)

#Mean Standardized Prediction Error (MSPE)

U.MSPE=mean(UresA2)

O.MSPE=mean(OresA2)

S.MSPE=mean(SresA2)

MSPE=c(U.MSPE,O.MSPE,S.MSPE)

datos=c(ME,RMSE,MSPE)

Evaluación=matrix(datos,nrow=3,ncol=3,byrow=TRUE,

dimnames = list(c("ME", "RMSE","MSPE"), c("Universal", "Ordinario", "Simple")))

#Se imprime una tabla donde se muestran los errores obtenidos por

#la validación cruzada

Evaluación#Kriging Universal el de menor error

#####################Selección de co-variograma y Co-Kriging#################

#Se modela la co-regionalización de los datos empleando la función gstat. Donde los modelos se #ajustan simultáneamente en forma directa y con variograma-cruzado

#Se crea el objeto gstat para especificar los variogramas experimentales

#Se llena el primer marco del objeto gstat con los valores de riqueza total

g2=gstat(id = "rte", fórmula = log(tsp)~1, data = rte.coor) #Posteriormente, se adicionan los valores de riqueza de Asteraceae, el objeto gstat,

#ahora tiene 2 marcos

g2=gstat(g2,id = "Asteraceae", formula = asp~1,data = rte. coor)

#Se adicionan parámetros de inicio para ajustar los modelos

g2=gstat(g2,id = "rte", model=vgm(psill=0.5744, "Lin",range=1.72, nugget=0.078))

g2=gstat(g2,id = "Asteraceae", model=vgm(psill=0.5744, "Lin",range=1.72, nugget=0.078)) v.cross2 <- variogram(g2)

g2=gstat(g2,id = "rte", model=vgm(psill=0.5744, "Lin",range=1.72, nugget=0.078),fill.all=T)

#Se realiza el ajuste de los modelos (g2 <- fit.lmc(v.cross2, g2))

#Se grafica variogramas ajustados directo y cruzado plot(variogram(g2), model=g2$model,col="blue",pch=20)

###Validación cruzcada de CoKriging### c.v=gstat.cv(g2,nfold=10)

#Se obtienen los valores de error del ajuste de los modelos

MEco=mean(c.v$residual)

RMSEco=sqrt(mean(c.v$residual^A2))

MSPEco=mean(c.v$residual^A2)

MEco

RMSEco

MSPEco

#############################Predicción espacial#############################

#Delimitar área de estudio con base en las coordenadas mínimas y máximas

#de los puntos de muestreo

xmin=min(rte$x)

xmax=max(rte$x)

ymin=min(rte$y)

ymax=max(rte$y)

#Resolución de 1km aprox

reso=0.008333

#Hacer Grid donde se almacenan datos interpolados

grid.xy <- expand.grid(x = seq(xmin,xmax,by=reso),y= seq(ymax,ymin,by=-reso))

coordinates(grid.xy) <- ~x + y

gridded(grid.xy)=T

###Interpolación Kriging###

#Se realiza la predicción espacial empleando el modelo teórico lineal y Kriging universal

Predicción=krige(log1p(tsp)~ 1,rte.coor,grid.xy ,fit.Lineal)

#Se elimina el logaritmo a valores de riqueza total, si es que se aplica, y se guarda el resultado en #nueva columna la "tsp"

Predicción$tsp=exp(Predicción$var1.pred)

#Se elimina el logaritmo a valores de varianza de la riqueza total, si es que se aplica y se guarda #los valores en nueva columna la "vtsp"

Predicción$vrte1g=exp(Predicción$var1.var)

###Interpolación Co-Kriging###

Ck=predict.gstat(g2,grid.xy)

# Se elimina el logaritmo a valores de riqueza total, si es que se aplica, y se guarda el resultado en #en nueva columna la "ck"

Ck$ck=exp(Ck$rte.pred)

# Se elimina el logaritmo a valores de varianza de la riqueza total, si es que se aplica y se guarda #los valores en nueva columna la "ckv"

Ck$ckv=exp(Ck$rte.var)

###Presentación de mapas###

#Se grafican los mapas de riqueza total y la varianza de los valores de riqueza total

krig=spplot(Predicción[''tsp"],mam=list("a)",cex=2),scales = list(draw = T),ylab=list("Coordenada Y",cex=1.3), col.regions = bpy.colors(100))

co.krig=spplot(Ck["ck"],main=list("b)",cex=2),scales=list(draw = T),xlab=list("Coordenada X",cex=1.3), col.regions = bpy.colors(100))

print(krig, position=c(0,0,0.5,1), more = TRUE)#a)Riqueza de especies con Kriging

print(co.krig, position=c(Q.5,Q,1,1), more = FALSE)#b)Riqueza de especies con Co-Kriging

###############################Exportación de cobertura##################

library(rgdal) #cargar librería

###Se exporta la capa de riqueza de especie generada por Kriging en formato Geotiff###

rts=Predicción["tsp"]#riqueza total de especies

writeGDAL(rts,"C:/Ejemplo k_ck/rtsKU. tif",drivername="GTiff",options=NULL)

### Se exporta la capa de riqueza de especie generada por Co-Kriging en formato Geotiff###

rts_ck=Ck["ck"]#riqueza total de especies Co-Kriging

writeGDAL(rts_ck,"C:/Ejemplo k_ck/rtsCk. tif",drivername="GTiff",options=NULL)

Nota: Se recomienda consultar las siguientes fuentes de información: http://spatial-analyst.net/ Fortin, M. J. y M. Dala. 2008. Spatial analysis, a guide for ecologists. Cambridge University Press. 392 p.