####################################################################
####Script para auxilio à resolução dos exercícios de geoestatística
####
####Exercício 1 - Teor de argila (%) na Terra Grande da Tapada da Ajuda
####

load("argila20.RData")exercício
attach(argila)

#a) Análise exploratória dos dados
summary(argila)
var(z); sd(z)

#histograma e boxplot
par(mfrow=c(1,2))
hist(z)
boxplot(z)

#b) Visualizar os dados
dev.off()
plot(x,y,pch=20)
text(x,y+3,round(z,0),cex=0.6)


#c) Estimação em pontos não amostrados por métodos determinísticos
A<-c(75,40)
B<-c(120,60)  
C<-c(100,120)

#visualização dos pontos no gráfico anterior
text(A[1],A[2],"A",col="red")
text(B[1],B[2],"B",col="green")
text(C[1],C[2],"C",col="blue")

#determinação da distância
dist.A<-sqrt((x-A[1])^2+(y-A[2])^2)
dist.B<-sqrt((x-B[1])^2+(y-B[2])^2)
dist.C<-sqrt((x-C[1])^2+(y-C[2])^2)

#i. Vizinho mais próximo
z[dist.A==min(dist.A)]
z[dist.B==min(dist.B)]
z[dist.C==min(dist.C)]

#ii. Média simples com raio 50,  75
mean(z[dist.A<=50])
mean(z[dist.A<=75])

mean(z[dist.B<=50])
mean(z[dist.B<=75])

mean(z[dist.C<=50])
mean(z[dist.C<=75])

#iii. Inverso da distância (k=1) com raio 50, 75
sum(1/dist.A[dist.A<=50]*z[dist.A<=50])/sum(1/dist.A[dist.A<=50])
sum(1/dist.A[dist.A<=75]*z[dist.A<=75])/sum(1/dist.A[dist.A<=75])

sum(1/dist.B[dist.B<=50]*z[dist.B<=50])/sum(1/dist.B[dist.B<=50])
sum(1/dist.B[dist.B<=75]*z[dist.B<=75])/sum(1/dist.B[dist.B<=75])

sum(1/dist.C[dist.C<=50]*z[dist.C<=50])/sum(1/dist.C[dist.C<=50])
sum(1/dist.C[dist.C<=75]*z[dist.C<=75])/sum(1/dist.C[dist.C<=75])


#iii. Inverso do quadrado da distância (k=2) com raio 50, 75
sum(1/dist.A[dist.A<=50]^2*z[dist.A<=50])/sum(1/dist.A[dist.A<=50]^2)
sum(1/dist.A[dist.A<=75]^2*z[dist.A<=75])/sum(1/dist.A[dist.A<=75]^2)

sum(1/dist.B[dist.B<=50]^2*z[dist.B<=50])/sum(1/dist.B[dist.B<=50]^2)
sum(1/dist.B[dist.B<=75]^2*z[dist.B<=75])/sum(1/dist.B[dist.B<=75]^2)

sum(1/dist.C[dist.C<=50]^2*z[dist.C<=50])/sum(1/dist.C[dist.C<=50]^2)
sum(1/dist.C[dist.C<=75]^2*z[dist.C<=75])/sum(1/dist.C[dist.C<=75]^2)


########################################################################
### Exercício 4 
### a) análise exploratória dos dados
install.packages("geoR")
load("GeoGrelha.RData")
attach(exGE)

summary(exGE)
var(z); sd(z)

# histograma e boxplot
par(mfrow=c(1,2)) 
hist(z); boxplot(z)

# visualizar os dados
dev.off()
plot(x,y,pch=20,xlab="X Coord",ylab="Y Coord")
text(x,y+0.2,z,cex=0.6)

### transformar os dados em tipo geodata com o comando as.geodata
library(geoR)
exGE.geo<-as.geodata(exGE)
exGE.geo
points(exGE.geo)
library(scatterplot3d)
plot(exGE.geo,scatter3d=T)

#b) variograma experimental
exGE.variog<-variog(exGE.geo,max.dist=4)
plot(exGE.variog)

#c) iii.
eyefit(exGE.variog)

# ajustar os parametros "Sill (patamar)", "Range (amplitude)" 
# e "Nugget (efeito pepita)"
# modelo esférico com sigmasq=10.2486,  phi=3.4418, tausq=1.3480

#d) 
distA<-sqrt((5-4.5)^2+(1-0.5)^2)
distB<-sqrt((5-6.5)^2+(1-5.5)^2)
distA; distB