Jorge Kazuo Yamamoto
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
Portuguese