em

Script – Como fazer Interpolação de Dados Geoespaciais e sua Representação Em Mapa usando a Linguagem R

A explicação do script a seguir está no artigo Como fazer Interpolação de Dados Geoespaciais e sua Representação em Mapa usando a Linguagem R

#script emqGlobal.R
#escrito por Jorge Kazuo Yamamoto
setwd("C:\\IPT2021\\dados\\walker_data")
#script - leitura dos parametros e dados
#leitura do arquivo de parametros
par<- read.csv("walker_dat_par.csv",sep=";",header=TRUE)
xmin=par$xmin; xmax=par$xmax; dx=par$dx; nx=par$nx
ymin=par$ymin; ymax=par$ymax; dy=par$dy; ny=par$ny
code=par$code  #valor nao calculado
#leitura do arquivo de dados
dados <- read.csv("walker_dat.csv",sep=";",header=TRUE)
x=dados$Xlocation[dados$V!=-99.00]
y=dados$Ylocation[dados$V!=-99.00]
z=dados$V[dados$V!=-99.00]
n=length(z)
print(c(n,min(z),max(z)))

#script sistema de equacoes multiquadricas
C=1  #constante multiquadrica
Q=matrix(c(rep(0,n*n)),nrow=n)
#calculando a matriz dos coeficientes
for (i in 1:n){
  for (j in 1:n){
    Q[i,j]=sqrt((x[i]-x[j])^2+(y[i]-y[j])^2+C)}}
z=matrix(c(z),ncol=1)
c=solve(Q,z)

#Script interpolacao da malha regular
print(c(nx,ny))
x0=c(rep(0,nx*ny)); y0=c(rep(0,nx*ny)); z0=c(rep(0,nx*ny))
for (i in 1:ny){
  for (j in 1:nx){
    kb=j+(i-1)*nx
    x0[kb]=xmin+(j-1)*dx+dx/2
    y0[kb]=ymin+(i-1)*dy+dy/2
    z0[kb]=0
    for (k in 1:n){
      z0[kb]=z0[kb]+c[k]*sqrt((x[k]-x0[kb])^2+(y[k]-y0[kb])^2+C)}
  }
}

#Script gravacao arquivo csv da malha regular
saida=data.frame(x0,y0,z0)
names(saida)[1]=paste("X0")
names(saida)[2]=paste("Y0")
names(saida)[3]=paste("Z0")
setwd("C:\\IPT2021\\dados\\walker_data")
write.csv(saida,file="walker_dat_EMQ.csv",row.names=FALSE,quote=F)

#Script leitura do arquivo de cores
def.par=par(no.readonly=TRUE)
setwd("C:\\IPT2021\\dados\\cores")
cores <- read.csv("cores.csv",sep=";",header=TRUE)
r=cores$r; g=cores$g; b=cores$b
ncores=length(r)
#Script leitura dos arquivos: parametros e malha regular
setwd("C:\\IPT2021\\dados\\walker_data")
#leitura do arquivo de parametros
par<- read.csv("walker_dat_par.csv",sep=";",header=TRUE)
xmin=par$xmin; xmax=par$xmax; dx=par$dx; nx=par$nx
ymin=par$ymin; ymax=par$ymax; dy=par$dy; ny=par$ny
code=par$code  #valor nao calculado
#leitura do arquivo da malha regular *.csv
#o arquivo csv gerado pelo R tem a "," como separador
dados <- read.csv("walker_dat_EMQ.csv",sep=",",header=TRUE)
x=dados[,1]; y=dados[,2]; z=dados[,3]
n=length(z)
#Script - definicao da dimensao do mapa no dispositivo grafico novo
if ((xmax-xmin) > (ymax-ymin)){
  cx=5
  cy=(ymax-ymin)*cx/(xmax-xmin)
} else {
  cy=5
  cx=(xmax-xmin)*cy/(ymax-ymin)
}
dev.new(width=cx, height=cy, unit="in")
setwd("C:\\IPT2021\\figuras")
#jpeg("walker_dat_map.jpeg",width=cx,height=cy,
#units="in",res=300)
par(mar=c(5,5,10,10), xpd=TRUE,cex=0.75)
plot(NA,NA,xlim=c(xmin,xmax), ylim=c(ymin,ymax), type="n", 
frame=TRUE, xlab=colnames(dados[1]), ylab=colnames(dados[2]))
#Script - definicao dos limites (zmin,zmax)
zmin=min(z[z!=code]); zmax=max(z[z!=code])
print(c(zmin,zmax))
zmin=9.9e+20; zmax=-zmin
for (i in 1:n){
  if (z[i] != code)    {
      if (z[i] < zmin) {zmin=z[i]}
      if (z[i] > zmax) {zmax=z[i]}
    }
}
print(c(zmin,zmax))
#Script - plotagem dos retangulos
delta=(zmax-zmin)*1.0000001
#plotagem como retangulos
for (i in 1:n)  {
  if (z[i] != code) {
    x1=x[i]-dx/2
    y1=y[i]-dy/2
    x2=x[i]+dx/2
    y2=y[i]+dy/2
    cor=trunc((z[i]-zmin)*ncores/delta)+1
    rect(x1,y1,x2,y2,border=NA,col=rgb(r[cor],g[cor],b[cor])) 
  }
}
#Script - plotagem da legenda de cores
dx1=(xmax-xmin)*0.10; dx2=(xmax-xmin)*0.14
eixox=c(rep(xmax+dx2,ncores+1))
eixoy=c(rep(0,ncores+1))
deltay=(ymax-ymin)/ncores
for (i in 1:(ncores+1))
{eixoy[i]=ymin+(i-1)*deltay}
eixoz=c(rep(0,ncores+1))
deltaz=(zmax-zmin)/ncores
for (i in 1:(ncores+1))
{eixoz[i]=zmin+(i-1)*deltaz}
vetorz=c(rep(0,ncores+1))
vetorz=sprintf("%.3f",eixoz)
text(x=c(eixox),y=c(eixoy),c(vetorz),xpd=NA,cex=0.65, pos=4)
for (i in 1:ncores)
{rect(xmax+dx1,eixoy[i],xmax+dx2,eixoy[i+1],border=NA,col=rgb(r[i],g[i],b[i]))}
#dev.off()
#fim

Arquivos para utilização:

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.

Como Calcular a Redução de Prazo de Financiamento na Tabela Price Após Amortização Extra

Como fazer Interpolação de Dados Geoespaciais e sua Representação em Mapa usando a Linguagem R