Jorge Kazuo Yamamoto
Você encontra a explicação desse script clicando aqui.
Você pode fazer o download desse script e dos arquivos de parâmetros (simetrica49.csv e simetrica49_varpar.csv) clicando aqui.
def.par=par(no.readonly=TRUE)
varPares=function(x,y,n,dir,tol,largBanda){
vd=as.vector(c(sin(dir),cos(dir)))
altCone=largBanda/tan(tol)
largBanda=largBanda^2
library(pracma)
k=0
pares=matrix(c(rep(0,n*(n-1))),
nrow=(n*(n-1)/2),ncol=2)
if (tol < pi/2) {
for (i in 1:(n-1)){
for (j in (i+1):n){
vp=as.vector(c((x[i]-x[j]),(y[i]-y[j])))
dist=abs(dot(vd,vp)) #produto escalar
lb=dot(vp,vp)-dist*dist
if (dist > altCone){
if (lb < largBanda){
k=k+1;pares[k,1]=i;pares[k,2]=j}
} else{
lbCone=(dist*tan(tol))^2
if (lb < lbCone){
k=k+1;pares[k,1]=i;pares[k,2]=j}
}
}
}
} else {
for (i in 1:(n-1)){
for (j in (i+1):n){
k=k+1;pares[k,1]=i;pares[k,2]=j}
}
}
npares=k
outPares=pares[1:npares,]
return(outPares)
}
varClasses=function(osPares,x,y,z,npassos,passo,tolPasso){
pmin=passo-tolPasso
pmax=pmin+passo*npassos
delta=(pmax-pmin)*1.0000001
#processa os pares de valores do variograma
d=c(rep(0,npassos))
pairs=c(rep(0,npassos))
gama=c(rep(0,npassos))
npares=nrow(osPares)
for (k in 1:npares){
i=osPares[k,1]; j=osPares[k,2]
dist=sqrt((x[i]-x[j])^2+(y[i]-y[j])^2)
if (dist<=pmax){
l=trunc((dist-pmin)*npassos/delta)+1
pairs[l]=pairs[l]+1
d[l]=d[l]+dist
gama[l]=gama[l]+(z[i]-z[j])^2
}
}
for (i in 1:npassos){
if (pairs[i] > 0){
d[i]=d[i]/pairs[i]
gama[i]=0.5*gama[i]/pairs[i]}}
res=cbind(d,gama,pairs)
return(res)
}
print("+++++++++++++++++++++++++++++++++++++++++++++++++++++++++")
setwd("C:\\GK21\\dados\\simetrica49")
dados=read.csv("simetrica49.csv",sep=";",header=T)
x=dados$X; y=dados$Y; z=dados$Zgauss
n=nrow(dados)
deg2rad=pi/180
#parametros do variograma
param=read.csv("simetrica49_varpar.csv",sep=";",header=T)
ndir=nrow(param)
dir=param$dir*deg2rad; tol=param$tol*deg2rad; largBanda=param$largBanda
passo=param$passo; tolPasso=param$tolPasso; npassos=param$npassos
gamaDat=c(rep(0,npassos[1]*3*ndir))
hmax=c(rep(0,ndir)); gmax=c(rep(0,ndir))
for (k in 1:ndir){
osPares=varPares(x,y,n,dir[k],tol[k],largBanda[k])
print(osPares)
res=round(varClasses(osPares,x,y,z,npassos[k],passo[k],
tolPasso[k]),digits=3)
hmax[k]=max(res[,1]); gmax[k]=max(res[,2])
for (i in 1:npassos[k]){
for (j in 1:3){
kb=i+(j-1)*npassos[k]
kc=j+(i-1)*3+(k-1)*3*npassos[k]
gamaDat[kc]=res[kb]}
}
}
recKC=function(j,k,npasso,vetor){
vet=c(rep(0,npasso))
for (i in 1:npasso){
kc=j+(i-1)*3+(k-1)*3*npasso
vet[i]=vetor[kc]}
return(vet)
}
xmax=max(hmax); ymax=max(gmax)
print(c(xmax,ymax))
#plotagem do variograma
setwd("C:\\GK21\\figuras")
jpeg("simetrica49_var.jpeg",width=7.5,height=5,
units="in",res=300)
par(mar=c(7,5,1,1))
plot(NA,NA,xlim=c(-5,30), ylim=c(-1,20), type="n",
frame=FALSE, xlab="", ylab="",asp=1,
xaxt="n",yaxt="n")
nx=4; cx=25
segments(0,0,cx,0,col="black") #eixo X
eixox=seq(from=0,to=xmax,length=nx+1)
valorx=c(rep(0,nx+1))
for (j in 1:(nx+1))
{valor=round(eixox[j],digits=2)
labx=paste(valor)
text(eixox[j]*cx/xmax,0.25,pos=1,labx,cex=0.90)
segments(eixox[j]*cx/xmax,0,eixox[j]*cx/xmax,-0.25,col="blue")
}
ny=5; cy=15
segments(0,0,0,cy,col="black") #eixo Y
eixoy=seq(from=0,to=ymax,length=ny+1)
for (i in 1:(ny+1)){
valor=round(eixoy[i],digits=3)
laby=paste(valor)
text(0.25,eixoy[i]*cy/ymax,pos=2,laby,cex=0.9)
segments(0,eixoy[i]*cy/ymax,-0.25,eixoy[i]*cy/ymax,col="black")
}
#plotagem dos pontos de dados
cor=c("red","blue","darkgreen","maroon")
for (k in 1:ndir){
x=recKC(1,k,npassos[k],gamaDat)*cx/xmax
y=recKC(2,k,npassos[k],gamaDat)*cy/ymax
print(cbind(x,y))
points(x,y,col=cor[k],pch=21,bg=cor[k],type="o",cex=1)
}
#legenda dos variogramas
for (k in 1:ndir){
posx=(k-1)*cx/ndir
points(posx,cy+1.5,col=cor[k],pch=21,bg=cor[k],type="o",cex=1)
labx=paste(dir[k]/deg2rad)
text(posx,cy+1.5,labx,pos=4,cex=0.9)
}
text(cx,0,pos=3,"h",cex=0.9)
laby=expression(gamma(h))
text(0,cy,pos=1,laby,srt=90,cex=0.9)
print(dir/deg2rad)
par(def.par)
dev.off()
Você pode fazer o download desse script e dos arquivos de parâmetros (simetrica49.csv e simetrica49_varpar.csv) clicando aqui.
Portuguese