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.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$residualA2))
MSPEco=mean(c.v$residualA2)
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.