em

Script: Como calcular variogramas experimentais para dados com distribuição irregular (Parte 2)

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.

Escrito por Jorge Kazuo Yamamoto

Prof. Dr. Jorge Kazuo Yamamoto, fundador da Geokrigagem, é geólogo, foi pesquisador do IPT e docente do Instituto de Geociências da USP, onde se aposentou como Professor Titular do Departamento de Geologia Sedimentar e Ambiental. Atualmente, atua como Professor Sênior do Departamento de Engenharia de Minas e de Petróleo – Escola Politécnica – USP. É responsável pela disciplina “Métodos geoestatísticos” na Pós-Graduação do IPT – Investigação do subsolo: Geotecnia e Meio Ambiente. Dedica-se ao ensino de geoestatística, com ênfase no desenvolvimento de algoritmos e pesquisa de novas aplicações, tais como: variância de interpolação, cálculo da variância global de depósitos minerais e correção do efeito de suavização da krigagem. Ultimamente, seu interesse está voltado para o ensino e divulgação da linguagem R.

Deixe uma resposta

O seu endereço de e-mail não será publicado.

Cálculo de variogramas experimentais para dados com distribuição irregular: parte 1

Distribuição normal: O que é e sua grande importância na estatística