################# Computación estadistica 2018-3 ####################
######## Julián David Gómez Ballén
####### Estudiante Ing. Agronomica
####################### Dependencia espacial(isotropia) ##################
library(MASS)
library(gstat) #geostadistica
## Warning: package 'gstat' was built under R version 3.5.2
library(deldir) #Tessellation
## deldir 0.1-15
library(ape) # Para correr moran
## Warning: package 'ape' was built under R version 3.5.2
library(sp)
library(GISTools)
## Warning: package 'GISTools' was built under R version 3.5.2
## Loading required package: maptools
## Warning: package 'maptools' was built under R version 3.5.2
## Checking rgeos availability: TRUE
## Loading required package: RColorBrewer
## Loading required package: rgeos
## Warning: package 'rgeos' was built under R version 3.5.2
## rgeos version: 0.4-2, (SVN revision 581)
## GEOS runtime version: 3.6.1-CAPI-1.10.1
## Linking to sp version: 1.3-1
## Polygon checking: TRUE
library(gstat)
library(geoR)
## Warning: package 'geoR' was built under R version 3.5.2
## --------------------------------------------------------------
## Analysis of Geostatistical Data
## For an Introduction to geoR go to http://www.leg.ufpr.br/geoR
## geoR version 1.7-5.2.1 (built on 2016-05-02) is now loaded
## --------------------------------------------------------------
capas <- array(NA, dim=c(50,2,9)) #Creación de arreglo vacio
# creación de datos dentro del arreglo anterior con distribución normal multivariada
# y difernetes niveles de correlación
for(r in 1:9){
z <- (mvrnorm(50,mu= c(11,15),
Sigma = matrix(c(1,r/10,r/10,1),ncol = 2),
empirical = TRUE))
capas[,,r]=z
}
# Data frame con datos asigandos
df <- data.frame(capas[,,2])
# Creación de funcion para analisis bivariado de las dos capas asiganadas
CVM <- function(datos){
CEa=as.matrix(datos,ncol=2)
medias=colMeans(CEa)
covar=cov(CEa)
CV=sqrt((t(medias)%*%covar%*%(medias))/((t(medias)%*%(medias))**2))
return(CV)
}
#correlación datos asignados. Baja correlación
CVM(capas[,,2])
## [,1]
## [1,] 0.05866411
# analisis de correlación graficamente, no se observa tendencia que indique correlación
plot(df,xlab=list("CEa 10 cm",cex=0.8),ylab=list("CEa 20 cm",cex=0.8), pch=19, col="red",main=list("CV multivariado CEa",cex=1))
grid(5,5)
#install.packages("gstat")
#Crea datos para coordenadas de puntos muestreales
xy <- expand.grid(x = seq(0, 450, 50), y= seq(0, 225, 50))
#Actualización de datos con coordenadas
df<- data.frame(df,xy)
#Grafico de datos y cambio de tipo de punto deacuerdo a valor de CEa
par(mfrow=c(1,2))
plot(xy, col=heat.colors(df$X1),pch=17,cex=2,main="CEa 10 cm")
grid(10,10)
##Mapa de calor de los puntos muestreales de acuerdo a valores de CEa
plot(xy, col=heat.colors(df$X2),pch=17,cex=2,main="CEa 20 cm")
grid(10,10)
par(mfrow=c(1,1))
# Creación 20 datos al azar de humedad dentro de las coordenads de CEa
humedad <- rnorm(20,40,2)
# Establecimiento de semilla y creación de las coordenadas para los 20 datos de humedad
set.seed(3542)
xh<- runif(20,0,500)
yh<- runif(20,0,200)
################### Moran ################
###########################################
cea.d=as.matrix(dist(xy, diag=T, upper=T)) #matriz de distancias, calcula la distancia entro los puntos
cea.d.inv <-as.matrix(1/cea.d) #Inversa de la distancia de la matriz
diag(cea.d.inv) <- 0
W=as.matrix(cea.d.inv) #inversa de los datos
dim(W)
## [1] 50 50
SUMAS=apply(W,1,sum); SUMAS
## 1 2 3 4 5 6 7
## 0.2510632 0.2945859 0.3186374 0.3323501 0.3386700 0.3386700 0.3323501
## 8 9 10 11 12 13 14
## 0.3186374 0.2945859 0.2510632 0.2855502 0.3381863 0.3664653 0.3821548
## 15 16 17 18 19 20 21
## 0.3892754 0.3892754 0.3821548 0.3664653 0.3381863 0.2855502 0.2950639
## 22 23 24 25 26 27 28
## 0.3502586 0.3799772 0.3963786 0.4037904 0.4037904 0.3963786 0.3799772
## 29 30 31 32 33 34 35
## 0.3502586 0.2950639 0.2855502 0.3381863 0.3664653 0.3821548 0.3892754
## 36 37 38 39 40 41 42
## 0.3892754 0.3821548 0.3664653 0.3381863 0.2855502 0.2510632 0.2945859
## 43 44 45 46 47 48 49
## 0.3186374 0.3323501 0.3386700 0.3386700 0.3323501 0.3186374 0.2945859
## 50
## 0.2510632
EstD=W/SUMAS #Estandariza todos los datos
# Inversos
mi1<-Moran.I(capas[,1,1], cea.d.inv) #Capa 1
mi2<-Moran.I(capas[,2,1], cea.d.inv) #Capa 2
# Estandarizados
me1<-Moran.I(capas[,1,1], EstD) #Capa 1
me2<-Moran.I(capas[,2,1], EstD) #Capa 2
#P-value de los valores de la prueba de moran
data.frame("Datos"=c("CEa 10 cm","CEa 20 cm"),"Inversos"=c(mi1$p.value,mi2$p.value),"Estandarizados"=c(me1$p.value,me2$p.value))
## Datos Inversos Estandarizados
## 1 CEa 10 cm 0.4749788 0.4749788
## 2 CEa 20 cm 0.9216578 0.9216578
########## Tessellation ######
#######################
#ubicación espacial de datos de humedad
plot(xh,yh,pch=1,col="darkgreen",xlab="Longitud",ylab="Latitud",main=list("CEa y Humedad del suelo",cex=1))
vt <- deldir(xh,yh,sort=T,digits=4)
##
## PLEASE NOTE: The components "delsgs" and "summary" of the
## object returned by deldir() are now DATA FRAMES rather than
## matrices (as they were prior to release 0.0-18).
## See help("deldir").
##
## PLEASE NOTE: The process that deldir() uses for determining
## duplicated points has changed from that used in version
## 0.0-9 of this package (and previously). See help("deldir").
listaT <- tile.list(vt)
#Plot resultado de funcion de teselación
plot(listaT,close=T, border=T, showpoints=T,number=T, cex=1, add=T)
points(xy,pch=18,col="darkblue",cex=2)
points(xh,yh,pch=19,col="darkgreen")
#Poligonos Creación de poligonos
df_xy <- data.frame(xh,yh,humedad) #data frame con coordenas y tesella
coordinates(df_xy) <- ~xh+yh #Convertir a coordenadas
z <- deldir(xh,yh,plot=TRUE,wl='triang',wp="none")
#Generar lista con información geografica por tesella
w <- tile.list(z)
#lista vacio de poligonos
polys <- vector(mode="list", length=length(w))
for (i in seq(along=polys)){
pcrds <- cbind(w[[i]]$x, w[[i]]$y)
pcrds <- rbind(pcrds, pcrds[1,])
polys[[i]] <- Polygons(list(Polygon(pcrds)), ID=as.character(i))
}
class(polys)
## [1] "list"
#Establecer clase SpatialPolygons
SP<- SpatialPolygons(polys)
#SP <- SpatialPolygonsDataFrame(SP,df,match.ID = T)
plot(SP)
class(SP)
## [1] "SpatialPolygons"
## attr(,"package")
## [1] "sp"
#Generación de poligonos y extracción de información de CEa por tesella
df50 <- data.frame(PN=1:50,C1=capas[,1,2],C2=capas[,2,2])
#Establecer clase SpatialPointsDataFrame
SS <- SpatialPointsDataFrame(xy,df50)
class(SS)
## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"
#Escalas de numeros para gener color traslucido
Color<- matrix((floor(runif(60,0,240)))/240, ncol = 3)
#plot inicial con ubicación de datos de CEa
plot(xy,main="Zonificación con teselas",xlab="Longitud",ylab="Latitud")
#For para insertar los 20 poligonos de tesellas
for (i in 1:20) {
AOI.XY <- gIntersection(SP[i,],SS,byid = T)
plot(SP[i,],add=T,col=rgb(Color[i,1],Color[i,2],Color[i,3],.15))
points(AOI.XY, add=T, pch=i, col=i)
text(xh[i], yh[i], c(i), cex = 1)
}
## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a
## graphical parameter
## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a
## graphical parameter
## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a
## graphical parameter
## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a
## graphical parameter
## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a
## graphical parameter
## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a
## graphical parameter
## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a
## graphical parameter
## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a
## graphical parameter
## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a
## graphical parameter
## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a
## graphical parameter
## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a
## graphical parameter
## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a
## graphical parameter
## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a
## graphical parameter
## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a
## graphical parameter
## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a
## graphical parameter
## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a
## graphical parameter
## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a
## graphical parameter
## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a
## graphical parameter
## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a
## graphical parameter
## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a
## graphical parameter
#Creación de vectores para construcción de tabla de datos y generar valores de
# CEa por cada tesella en este caso se utilizo dos for para extraer los valore
CEa1<-c();CEa2<-c();CEa.H1<-c();CEa.H2<-c()
for (i in 1:20) {
AOI.XY.DF<- as.data.frame(gIntersection(SP[i,],SS,byid = T))
for (j in 1:length(AOI.XY.DF$x)) {
CEa1[j]<-df$X1[which(AOI.XY.DF[j,1]==df$x & AOI.XY.DF[j,2]==df$y)]
CEa2[j]<-df$X2[which(AOI.XY.DF[j,1]==df$x & AOI.XY.DF[j,2]==df$y)]
}
CEa.H1[i]<- mean(CEa1)
CEa.H2[i]<- mean(CEa1)
}
#Construcción de base de datos de valor de CEa a 10 y 20 cm promedio por tesella
hCe<-data.frame(Tesella=c(1:20),xh,yh,humedad,"CEa 10 cm"=CEa.H1,"CEa 20 cm"=CEa.H2)
hCe
## Tesella xh yh humedad CEa.10.cm CEa.20.cm
## 1 1 9.922637 155.099092 41.00373 11.12957 11.12957
## 2 2 339.258301 189.825618 40.92120 10.15740 10.15740
## 3 3 341.407724 143.653932 41.78654 11.04267 11.04267
## 4 4 161.338173 181.459025 38.67353 10.50294 10.50294
## 5 5 5.858352 86.273620 37.84633 11.45867 11.45867
## 6 6 491.881970 196.625404 40.90012 12.62305 12.62305
## 7 7 368.613237 170.289416 39.23790 11.08550 11.08550
## 8 8 184.241928 43.918768 39.36721 10.94053 10.94053
## 9 9 85.331983 4.295758 43.51806 11.03696 11.03696
## 10 10 302.711655 109.723986 41.77441 11.57913 11.57913
## 11 11 280.387327 113.553232 38.94842 11.23977 11.23977
## 12 12 61.457638 185.056855 44.22395 11.57287 11.57287
## 13 13 103.363560 113.385836 44.09239 11.32855 11.32855
## 14 14 310.909901 157.856323 39.95973 11.60452 11.60452
## 15 15 22.439754 164.725537 41.46503 11.64954 11.64954
## 16 16 127.814394 141.906642 40.13757 11.45986 11.45986
## 17 17 494.337955 131.673948 37.69927 11.47126 11.47126
## 18 18 122.264813 105.540411 41.78268 11.33681 11.33681
## 19 19 82.785694 26.409680 43.35692 11.70243 11.70243
## 20 20 479.115306 38.455016 46.22240 11.31079 11.31079
############ Moran con nueva base de datos de CEa por tesella########
######################################################################
cea.d=as.matrix(dist(hCe[,c(2,3)], diag=T, upper=T)) #matriz de distancias, calcula la distancia entro los puntos
cea.d.inv <-as.matrix(1/cea.d) #Inversa de la distancia de la matriz
diag(cea.d.inv) <- 0
W=as.matrix(cea.d.inv) #inversa de los datos
dim(W)
## [1] 20 20
SUMAS=apply(W,1,sum); SUMAS
## 1 2 3 4 5 6
## 0.17036012 0.15492114 0.17402812 0.13437719 0.11637965 0.08383741
## 7 8 9 10 11 12
## 0.15243136 0.12214473 0.13975949 0.17428420 0.16976773 0.13962020
## 13 14 15 16 17 18
## 0.19335571 0.17394719 0.17766692 0.17384297 0.08860286 0.19275871
## 19 20
## 0.15059642 0.07608352
EstD=W/SUMAS #Estandariza todos los datos
# Inversos
mh<-Moran.I(hCe$humedad, cea.d.inv) #humedad
mcea1<-Moran.I(hCe$CEa.10.cm, cea.d.inv) #Capa 1
mcea2<-Moran.I(hCe$CEa.20.cm, cea.d.inv) #Capa 1
# Estandarizados
meh<-Moran.I(hCe$humedad, EstD) #humedad
mecea1<-Moran.I(hCe$CEa.10.cm, EstD) #Capa 1
mecea2<-Moran.I(hCe$CEa.20.cm, EstD) #Capa 1
#P-value de los valores de la prueba de moran
##Valores de p-valor prueba de moran para cada caso
data.frame("Datos"=c("Humedad","CEa 10 cm","CEa 20 cm"),"Inversos"=c(mh$p.value,mcea1$p.value,mcea2$p.value),"Estandarizados"=c(meh$p.value,mecea1$p.value,mecea2$p.value))
## Datos Inversos Estandarizados
## 1 Humedad 0.9454998 0.9454998
## 2 CEa 10 cm 0.8660538 0.8660538
## 3 CEa 20 cm 0.8660538 0.8660538
###########Funcion de auto correlacion##########3
##############################################
Prp <- sort(rnorm(120,200,1))
acf(Prp)
#Autocorrelación de datos de CEa 10 y 20 cm
acf(capas[,1,2])
acf(capas[,1,2])
#Generación de la grilla
xy <- expand.grid(x = seq(10, 200, 10), y= seq(10, 100, 10))
plot(xy,col="darkgreen",ylab="Latitud",xlab="Longitud")
set.seed(1966)
names(xy) <- c('x','y')
max(dist(xy))
## [1] 210.238
capa1 <- gstat(formula=z~1, locations=~x+y,
dummy=T, beta=1,
model=vgm(psill=0.03, range=100, model='Exp'), nmax=20) #range: dependencia espacial maxima (0 minima)
class(capa1)
## [1] "gstat" "list"
#Asociación datos COO con coordenas
COO <- predict(capa1, newdata=xy, nsim=1)
## [using unconditional Gaussian simulation]
COOc=data.frame(x=COO$x,y=COO$y,COO=8*(COO$sim1))
#plot de calor de los datos deacuerdo a COO
plot(COOc$x,COOc$y,pch=19,cex=3,col=heat.colors(COOc$COO),xlab="Longitud",ylab="Longitud",main="mapa de calor de CO")
#Calculo de angulos deacuerdo a coordenads
angulos <- (180*atan(COOc$y/COOc$x))/pi
#Asociación de estas coordenadas
COOc=data.frame(x=COO$x,y=COO$y,COO=8*(COO$sim1),angulos)
#Categorización
cate<-ifelse(angulos<=30,1,ifelse(angulos<=60,2,3))
plot(COOc$x,COOc$y,col=cate,pch=cate+15,cex=3,xlab="Longitud",ylab="Longitud")
#Base de datos de acuerdo a categorización
men30<-COOc[which(cate==1),]
bet3060 <- COOc[which(cate==2),]
mas60<- COOc[which(cate==3),]
#Moran dependencia espacial#
#General
cea.d=as.matrix(dist(COOc[,1:2], diag=T, upper=T)) #matriz de distancias, calcula la distancia entro los puntos
cea.d.inv <-as.matrix(1/cea.d) #Inversa de la distancia de la matriz
diag(cea.d.inv) <- 0
W=as.matrix(cea.d.inv) #inversa de los datos
dim(W)
## [1] 200 200
SUMAS=apply(W,1,sum); SUMAS
## 1 2 3 4 5 6 7 8
## 2.531343 2.825535 3.023424 3.170851 3.283888 3.370660 3.435951 3.482775
## 9 10 11 12 13 14 15 16
## 3.513048 3.527915 3.527915 3.513048 3.482775 3.435951 3.370660 3.283888
## 17 18 19 20 21 22 23 24
## 3.170851 3.023424 2.825535 2.531343 2.790818 3.144178 3.375363 3.543199
## 25 26 27 28 29 30 31 32
## 3.669546 3.765281 3.836651 3.887501 3.920235 3.936271 3.936271 3.920235
## 33 34 35 36 37 38 39 40
## 3.887501 3.836651 3.765281 3.669546 3.543199 3.375363 3.144178 2.790818
## 41 42 43 44 45 46 47 48
## 2.944497 3.329791 3.583760 3.767108 3.904023 4.007002 4.083314 4.137441
## 49 50 51 52 53 54 55 56
## 4.172177 4.189163 4.189163 4.172177 4.137441 4.083314 4.007002 3.904023
## 57 58 59 60 61 62 63 64
## 3.767108 3.583760 3.329791 2.944497 3.036263 3.438777 3.706443 3.899877
## 65 66 67 68 69 70 71 72
## 4.043972 4.151980 4.231756 4.288191 4.324337 4.341992 4.341992 4.324337
## 73 74 75 76 77 78 79 80
## 4.288191 4.231756 4.151980 4.043972 3.899877 3.706443 3.438777 3.036263
## 81 82 83 84 85 86 87 88
## 3.079725 3.489922 3.763985 3.962332 4.110026 4.220596 4.302152 4.359775
## 89 90 91 92 93 94 95 96
## 4.396650 4.414651 4.414651 4.396650 4.359775 4.302152 4.220596 4.110026
## 97 98 99 100 101 102 103 104
## 3.962332 3.763985 3.489922 3.079725 3.079725 3.489922 3.763985 3.962332
## 105 106 107 108 109 110 111 112
## 4.110026 4.220596 4.302152 4.359775 4.396650 4.414651 4.414651 4.396650
## 113 114 115 116 117 118 119 120
## 4.359775 4.302152 4.220596 4.110026 3.962332 3.763985 3.489922 3.079725
## 121 122 123 124 125 126 127 128
## 3.036263 3.438777 3.706443 3.899877 4.043972 4.151980 4.231756 4.288191
## 129 130 131 132 133 134 135 136
## 4.324337 4.341992 4.341992 4.324337 4.288191 4.231756 4.151980 4.043972
## 137 138 139 140 141 142 143 144
## 3.899877 3.706443 3.438777 3.036263 2.944497 3.329791 3.583760 3.767108
## 145 146 147 148 149 150 151 152
## 3.904023 4.007002 4.083314 4.137441 4.172177 4.189163 4.189163 4.172177
## 153 154 155 156 157 158 159 160
## 4.137441 4.083314 4.007002 3.904023 3.767108 3.583760 3.329791 2.944497
## 161 162 163 164 165 166 167 168
## 2.790818 3.144178 3.375363 3.543199 3.669546 3.765281 3.836651 3.887501
## 169 170 171 172 173 174 175 176
## 3.920235 3.936271 3.936271 3.920235 3.887501 3.836651 3.765281 3.669546
## 177 178 179 180 181 182 183 184
## 3.543199 3.375363 3.144178 2.790818 2.531343 2.825535 3.023424 3.170851
## 185 186 187 188 189 190 191 192
## 3.283888 3.370660 3.435951 3.482775 3.513048 3.527915 3.527915 3.513048
## 193 194 195 196 197 198 199 200
## 3.482775 3.435951 3.370660 3.283888 3.170851 3.023424 2.825535 2.531343
EstD=W/SUMAS #Estandariza todos los datos
# Inversos
cooi<-Moran.I(COOc$COO, cea.d.inv)
# Estandarizados
cooe<-Moran.I(COOc$COO, EstD)
#30°
cea.d=as.matrix(dist(men30[,1:2], diag=T, upper=T)) #matriz de distancias, calcula la distancia entro los puntos
cea.d.inv <-as.matrix(1/cea.d) #Inversa de la distancia de la matriz
diag(cea.d.inv) <- 0
W=as.matrix(cea.d.inv) #inversa de los datos
dim(W)
## [1] 110 110
SUMAS=apply(W,1,sum); SUMAS
## 2 3 4 5 6 7 8 9
## 1.189167 1.407219 1.617500 1.797452 1.954607 2.092408 2.210007 2.308907
## 10 11 12 13 14 15 16 17
## 2.390541 2.455465 2.503457 2.533642 2.544508 2.533814 2.498312 2.433085
## 18 19 20 24 25 26 27 28
## 2.329980 2.172981 1.916368 1.625763 1.871118 2.102035 2.292678 2.445357
## 29 30 31 32 33 34 35 36
## 2.570568 2.673471 2.755926 2.818264 2.859750 2.878718 2.872521 2.837167
## 37 38 39 40 46 47 48 49
## 2.766399 2.649442 2.464633 2.153406 1.999688 2.278575 2.480234 2.641092
## 50 51 52 53 54 55 56 57
## 2.773108 2.879367 2.961437 3.018937 3.050208 3.052354 3.020931 2.948899
## 58 59 60 67 68 69 70 71
## 2.824115 2.622808 2.284742 2.072497 2.328381 2.552943 2.730101 2.870962
## 72 73 74 75 76 77 78 79
## 2.981753 3.062579 3.112387 3.128238 3.105732 3.037801 2.912007 2.704238
## 80 89 90 91 92 93 94 95
## 2.354649 2.287563 2.522557 2.728649 2.886050 3.002868 3.081541 3.118693
## 96 97 98 99 100 111 112 113
## 3.111393 3.053759 2.933925 2.728157 2.377260 2.423436 2.642602 2.830206
## 114 115 116 117 118 119 120 133
## 2.958957 3.028331 3.044164 3.004150 2.897710 2.701879 2.358133 2.518119
## 134 135 136 137 138 139 140 154
## 2.733197 2.849943 2.899938 2.887572 2.803512 2.626300 2.298294 2.374148
## 155 156 157 158 159 160 176 177
## 2.544850 2.658769 2.692095 2.643425 2.496263 2.194472 2.282469 2.376217
## 178 179 180 198 199 200
## 2.395257 2.296616 2.035529 2.020807 1.980076 1.780684
EstD=W/SUMAS #Estandariza todos los datos
# Inversos
coo30i<-Moran.I(men30$COO, cea.d.inv)
# Estandarizados
coo30e<-Moran.I(men30$COO, EstD)
#30-60°
cea.d=as.matrix(dist(bet3060[,1:2], diag=T, upper=T)) #matriz de distancias, calcula la distancia entro los puntos
cea.d.inv <-as.matrix(1/cea.d) #Inversa de la distancia de la matriz
diag(cea.d.inv) <- 0
W=as.matrix(cea.d.inv) #inversa de los datos
dim(W)
## [1] 63 63
SUMAS=apply(W,1,sum); SUMAS
## 1 22 23 42 43 44 45
## 0.8694751 1.1886835 1.2905561 1.2868381 1.5123741 1.5753901 1.5339067
## 63 64 65 66 83 84 85
## 1.5700022 1.7613963 1.8083965 1.7187586 1.5209660 1.8001699 1.9408871
## 86 87 88 104 105 106 107
## 1.9684012 1.9002859 1.7743509 1.6984967 1.9549629 2.0784974 2.1052222
## 108 109 110 125 126 127 128
## 2.0548626 1.9311576 1.7611745 1.8660878 2.0812712 2.1702896 2.1799517
## 129 130 131 132 145 146 147
## 2.1285352 2.0202534 1.8503282 1.6388483 1.7098227 1.9911479 2.1325541
## 148 149 150 151 152 153 166
## 2.1820227 2.1695100 2.1081087 2.0025928 1.8456014 1.6035931 1.8046754
## 167 168 169 170 171 172 173
## 1.9941071 2.0668745 2.0775080 2.0439521 1.9729088 1.8639349 1.7094060
## 174 175 186 187 188 189 190
## 1.5088866 1.2895024 1.5391412 1.7166467 1.7913362 1.8121764 1.7952398
## 191 192 193 194 195 196 197
## 1.7473884 1.6713131 1.5680502 1.4401240 1.2852552 1.0967499 0.8977884
EstD=W/SUMAS #Estandariza todos los datos
# Inversos
coo3060i<-Moran.I(bet3060$COO, cea.d.inv)
# Estandarizados
coo3060e<-Moran.I(bet3060$COO, EstD)
#60°
cea.d=as.matrix(dist(mas60[,1:2], diag=T, upper=T)) #matriz de distancias, calcula la distancia entro los puntos
cea.d.inv <-as.matrix(1/cea.d) #Inversa de la distancia de la matriz
diag(cea.d.inv) <- 0
W=as.matrix(cea.d.inv) #inversa de los datos
dim(W)
## [1] 27 27
SUMAS=apply(W,1,sum); SUMAS
## 21 41 61 62 81 82 101
## 0.5980605 0.7704590 0.9273039 0.9209747 1.0437894 1.0983985 1.1241852
## 102 103 121 122 123 124 141
## 1.2460011 1.1268226 1.1669314 1.3315239 1.2941418 1.0794056 1.1628220
## 142 143 144 161 162 163 164
## 1.3461452 1.3475898 1.1842213 1.0970169 1.2780657 1.2988388 1.1950112
## 165 181 182 183 184 185
## 0.9491436 0.9248329 1.0729817 1.0978855 1.0330588 0.8579061
EstD=W/SUMAS #Estandariza todos los datos
# Inversos
coo60i<-Moran.I(mas60$COO, cea.d.inv) #
# Estandarizados
coo60e<-Moran.I(mas60$COO, EstD) #
#P-value de los valores de la prueba de moran 30 30-60 60
data.frame("Datos"=c("General","30°","30-60°","60°"),"Inversos"=c(cooi$p.value,coo30i$p.value,coo3060i$p.value,coo60i$p.value),"Estandarizados"=c(cooe$p.value,coo30e$p.value,coo3060e$p.value,coo60e$p.value))
## Datos Inversos Estandarizados
## 1 General 0.000000e+00 0.000000e+00
## 2 30° 0.000000e+00 0.000000e+00
## 3 30-60° 3.352874e-14 3.352874e-14
## 4 60° 6.492305e-03 6.492305e-03
####################Variogramas#####################
#################################################
#semivarigrama general
COOc
## x y COO angulos
## 1 10 10 7.599606 45.000000
## 2 20 10 8.095129 26.565051
## 3 30 10 8.742248 18.434949
## 4 40 10 9.019807 14.036243
## 5 50 10 8.704293 11.309932
## 6 60 10 9.059224 9.462322
## 7 70 10 8.128114 8.130102
## 8 80 10 8.459307 7.125016
## 9 90 10 8.236819 6.340192
## 10 100 10 8.121514 5.710593
## 11 110 10 7.192771 5.194429
## 12 120 10 7.375476 4.763642
## 13 130 10 7.231232 4.398705
## 14 140 10 7.391270 4.085617
## 15 150 10 7.026750 3.814075
## 16 160 10 7.241300 3.576334
## 17 170 10 8.803706 3.366461
## 18 180 10 7.823838 3.179830
## 19 190 10 7.439845 3.012788
## 20 200 10 6.134178 2.862405
## 21 10 20 7.335439 63.434949
## 22 20 20 8.275702 45.000000
## 23 30 20 7.969560 33.690068
## 24 40 20 8.646882 26.565051
## 25 50 20 8.491668 21.801409
## 26 60 20 7.906565 18.434949
## 27 70 20 7.925427 15.945396
## 28 80 20 7.894993 14.036243
## 29 90 20 8.465109 12.528808
## 30 100 20 8.871341 11.309932
## 31 110 20 7.897435 10.304846
## 32 120 20 7.883980 9.462322
## 33 130 20 8.443875 8.746162
## 34 140 20 7.322126 8.130102
## 35 150 20 7.355387 7.594643
## 36 160 20 7.083032 7.125016
## 37 170 20 8.279132 6.709837
## 38 180 20 7.941939 6.340192
## 39 190 20 7.082733 6.009006
## 40 200 20 6.332145 5.710593
## 41 10 30 8.190212 71.565051
## 42 20 30 7.841665 56.309932
## 43 30 30 8.038933 45.000000
## 44 40 30 7.878435 36.869898
## 45 50 30 7.592060 30.963757
## 46 60 30 8.013796 26.565051
## 47 70 30 7.924535 23.198591
## 48 80 30 7.319558 20.556045
## 49 90 30 8.011812 18.434949
## 50 100 30 8.385343 16.699244
## 51 110 30 8.202766 15.255119
## 52 120 30 7.728922 14.036243
## 53 130 30 8.159523 12.994617
## 54 140 30 7.963435 12.094757
## 55 150 30 7.540375 11.309932
## 56 160 30 7.316779 10.619655
## 57 170 30 7.874661 10.007980
## 58 180 30 7.150154 9.462322
## 59 190 30 6.821391 8.972627
## 60 200 30 6.129259 8.530766
## 61 10 40 8.635670 75.963757
## 62 20 40 8.240382 63.434949
## 63 30 40 8.399942 53.130102
## 64 40 40 8.254744 45.000000
## 65 50 40 7.401454 38.659808
## 66 60 40 7.619006 33.690068
## 67 70 40 7.861113 29.744881
## 68 80 40 6.491736 26.565051
## 69 90 40 7.680521 23.962489
## 70 100 40 7.304784 21.801409
## 71 110 40 7.283616 19.983107
## 72 120 40 7.059701 18.434949
## 73 130 40 7.893995 17.102729
## 74 140 40 8.312197 15.945396
## 75 150 40 8.221010 14.931417
## 76 160 40 7.906387 14.036243
## 77 170 40 8.416477 13.240520
## 78 180 40 7.314311 12.528808
## 79 190 40 6.532924 11.888658
## 80 200 40 6.864060 11.309932
## 81 10 50 7.966257 78.690068
## 82 20 50 8.795851 68.198591
## 83 30 50 8.612638 59.036243
## 84 40 50 7.733980 51.340192
## 85 50 50 7.203474 45.000000
## 86 60 50 8.080185 39.805571
## 87 70 50 7.892415 35.537678
## 88 80 50 7.578161 32.005383
## 89 90 50 8.103988 29.054604
## 90 100 50 7.761403 26.565051
## 91 110 50 7.133402 24.443955
## 92 120 50 7.603435 22.619865
## 93 130 50 8.112029 21.037511
## 94 140 50 7.758917 19.653824
## 95 150 50 7.763667 18.434949
## 96 160 50 7.968409 17.354025
## 97 170 50 7.674114 16.389540
## 98 180 50 8.079547 15.524111
## 99 190 50 8.170059 14.743563
## 100 200 50 8.211835 14.036243
## 101 10 60 7.295954 80.537678
## 102 20 60 7.346365 71.565051
## 103 30 60 7.930076 63.434949
## 104 40 60 8.824402 56.309932
## 105 50 60 8.556037 50.194429
## 106 60 60 7.598719 45.000000
## 107 70 60 7.351245 40.601295
## 108 80 60 7.580835 36.869898
## 109 90 60 7.017160 33.690068
## 110 100 60 7.306309 30.963757
## 111 110 60 8.217058 28.610460
## 112 120 60 7.908980 26.565051
## 113 130 60 7.474753 24.775141
## 114 140 60 7.127772 23.198591
## 115 150 60 6.572388 21.801409
## 116 160 60 7.888066 20.556045
## 117 170 60 8.179105 19.440035
## 118 180 60 8.447729 18.434949
## 119 190 60 8.793820 17.525568
## 120 200 60 7.557568 16.699244
## 121 10 70 7.084392 81.869898
## 122 20 70 6.863608 74.054604
## 123 30 70 7.271312 66.801409
## 124 40 70 7.952480 60.255119
## 125 50 70 7.889137 54.462322
## 126 60 70 6.954197 49.398705
## 127 70 70 7.650751 45.000000
## 128 80 70 7.417158 41.185925
## 129 90 70 7.592549 37.874984
## 130 100 70 8.284936 34.992020
## 131 110 70 8.606913 32.471192
## 132 120 70 8.283666 30.256437
## 133 130 70 7.553545 28.300756
## 134 140 70 7.197820 26.565051
## 135 150 70 6.895631 25.016893
## 136 160 70 7.909714 23.629378
## 137 170 70 7.683515 22.380135
## 138 180 70 8.104276 21.250506
## 139 190 70 8.749349 20.224859
## 140 200 70 7.099870 19.290046
## 141 10 80 8.085214 82.874984
## 142 20 80 7.733637 75.963757
## 143 30 80 7.164853 69.443955
## 144 40 80 7.271214 63.434949
## 145 50 80 6.760005 57.994617
## 146 60 80 6.829668 53.130102
## 147 70 80 6.715952 48.814075
## 148 80 80 6.828267 45.000000
## 149 90 80 7.305087 41.633539
## 150 100 80 8.236507 38.659808
## 151 110 80 8.448192 36.027373
## 152 120 80 7.699719 33.690068
## 153 130 80 7.681319 31.607502
## 154 140 80 7.404845 29.744881
## 155 150 80 7.862743 28.072487
## 156 160 80 7.721190 26.565051
## 157 170 80 7.174370 25.201124
## 158 180 80 7.313687 23.962489
## 159 190 80 8.898555 22.833654
## 160 200 80 7.392265 21.801409
## 161 10 90 8.038233 83.659808
## 162 20 90 7.805701 77.471192
## 163 30 90 7.785034 71.565051
## 164 40 90 7.100927 66.037511
## 165 50 90 7.462243 60.945396
## 166 60 90 6.676769 56.309932
## 167 70 90 6.479807 52.125016
## 168 80 90 6.613892 48.366461
## 169 90 90 6.449153 45.000000
## 170 100 90 7.235559 41.987212
## 171 110 90 8.259346 39.289407
## 172 120 90 7.279548 36.869898
## 173 130 90 6.762874 34.695154
## 174 140 90 7.593700 32.735226
## 175 150 90 7.698771 30.963757
## 176 160 90 7.411485 29.357754
## 177 170 90 6.501895 27.897271
## 178 180 90 7.611600 26.565051
## 179 190 90 8.661821 25.346176
## 180 200 90 8.426834 24.227745
## 181 10 100 9.075759 84.289407
## 182 20 100 9.009368 78.690068
## 183 30 100 7.880638 73.300756
## 184 40 100 8.126968 68.198591
## 185 50 100 7.392014 63.434949
## 186 60 100 6.810459 59.036243
## 187 70 100 7.236910 55.007980
## 188 80 100 7.454287 51.340192
## 189 90 100 7.001478 48.012788
## 190 100 100 7.620176 45.000000
## 191 110 100 7.946354 42.273689
## 192 120 100 7.343051 39.805571
## 193 130 100 7.136161 37.568592
## 194 140 100 7.850639 35.537678
## 195 150 100 7.834189 33.690068
## 196 160 100 7.230361 32.005383
## 197 170 100 6.853681 30.465545
## 198 180 100 7.527069 29.054604
## 199 190 100 8.301206 27.758541
## 200 200 100 8.464602 26.565051
dfv=data.frame(x=COOc$x,y=COOc$y,z=COOc$COO/8)
coords<-dfv[,1:2]
varg <- variog(coords=coords, data=dfv$z,estimator.type='classical')
## variog: computing omnidirectional variogram
plot(varg,pch=19,col="brown",cex=1.5,main="Semivarigrama general")
#Semivariograma menor a 30
dfv=data.frame(x=men30$x,y=men30$y,z=men30$COO/8)
coords<-dfv[,1:2]
varg <- variog(coords=coords, data=dfv$z,estimator.type='classical')
## variog: computing omnidirectional variogram
plot(varg,pch=19,col="red",cex=1.5,main="Semivarigrama <30°")
#Semivariograma 60-60
dfv=data.frame(x=bet3060$x,y=bet3060$y,z=bet3060$COO/8)
coords<-dfv[,1:2]
varg <- variog(coords=coords, data=dfv$z,estimator.type='classical')
## variog: computing omnidirectional variogram
plot(varg,pch=19,col="blue",cex=1.5,main="Semivarigrama 30-60°")
#Semivariograma mayor a 60
dfv=data.frame(x=mas60$x,y=mas60$y,z=mas60$COO/8)
coords<-dfv[,1:2]
varg <- variog(coords=coords, data=dfv$z,estimator.type='classical')
## variog: computing omnidirectional variogram
plot(varg,pch=19,col="darkgreen",cex=1.5,main="Semivarigrama >60°")
