Desenvolvimento

Planilha com coordenadas e valores do atributo


library(openxlsx)
library(gstat)
data(oxford)
XCOORD <- c(oxford[,2])
YCOORD <- c(oxford[,3])
ELEV <- c(oxford[,4])
df <- data.frame(XCOORD, YCOORD, ELEV)
write.xlsx(df, 'coord_atributos.xlsx')



Mapa base bubble plot

coordinates(df)<-c("XCOORD","YCOORD")
bubble(df,"ELEV", maxsize = 1.1, col="dark green", xlab="XCOORD", ylab="YCOORD")


Mapa base sp plot (extra)

coordinates(oxford)<-c("XCOORD","YCOORD")
spplot(oxford, "ELEV",xlab="XCOORD",ylab="YCOORD",
main="Elevação", col.regions=bpy.colors(5))




Sumário da variável


summary(df$ELEV)

 Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  540.0   558.0   573.0   573.6   584.5   632.0 


Histograma de densidade


h<-hist(df$ELEV)
xhist<-c(min(h$breaks),h$breaks)
yhist<-c(0,h$density,0)
xfit<-seq(540,640, by=1.0)
yfit<-dnorm(xfit,mean=mean(df$ELEV),sd=sd(df$ELEV))
plot(xhist,yhist,type='s',ylim=c(0,max(yhist,yfit)),
     main="Histograma de elevação e curva normal ajustada",
     xlab="Elevação",ylab="Densidade")
lines(xfit,yfit,col="dark green")





Box-plot com outliers


boxplot(df$ELEV,horizontal=TRUE,range=1.5,main='Boxplot de ELEV', xlab='ELEV')





Variograma experimental 


data(oxford)
g <- gstat(id="ELEV", formula=ELEV~1, locations=~XCOORD+YCOORD,data= oxford)
graf<-variogram(g)
 plot(graf, , xlab='ELEV', ylab='Semivariância')




Estimação visual do alcance, efeito de pepita e patamar






🟥 => Efeito Pepita(Co)   🟨 => Alcance(A)    🟩 => Patamar(C)




Ajuste de um modelo de variograma



vgm1=variogram(ELEV~1, locations=~XCOORD+YCOORD, data=oxford)
f<-fit.variogram(vgm1, vgm(350,"Gau", 610, 44)) #Foi escolhido o modelo gaussiano por apresentar o melhor resultado.
ff<-variogramLine(f,maxdist=650,n=500,min=1.0e-6)
plot(ff,col='green')
points(vgm1[,2],vgm1[,3], col='red')



print(f)

  model     psill    range
1   Nug  19.78417   0.0000
2   Gau 367.66942 356.2184

Grid com máscara


Para a montagem do grid com máscara, primeiramente foi escolhido um polígono de seis lados, sendo possível visualizar as coordenadas escolhidas seguindo o link abaixo: 

library(sp)
library(gstat)
library(lattice)
library(latticeExtra)
data(oxford)
pts = read.table("C:/Users/Elson/Desktop/geoestatistica/coordenadas.txt", header=TRUE, sep="")
pts = rbind(pts, pts[1,])
spdf = SpatialPolygons( list( Polygons(list(Polygon(pts)), 1))) 
plot(spdf)


grd<-makegrid(spdf,n=10000,c(10,10))
colnames(grd) <- c('x','y')
grd_pts <- SpatialPoints(coords = grd)
grd_pts_in <- grd_pts[spdf, ]
novogrid<- as.data.frame(coordinates(grd_pts_in))
plot(novogrid)

Mapa de densidade do atributo

data(oxford)
coordinates(oxford)=~XCOORD+YCOORD
gridded(novogrid) = ~x+y
m <- vgm(367.66942, "Gau", 356.2184, 19.78417)
pred <- krige(ELEV~1, oxford, novogrid, model = m)
l2 = list("SpatialPolygonsRescale", layout.north.arrow(),offset = c(50,1800),scale = 250 )
scale <- list("SpatialPolygonsRescale",layout.scale.bar(height=0.05),offset = c(100,100),scale = 300,fill=c("transparent","black"))
text1 <- list("sp.text", c(70,100), "0")
text2 <- list("sp.text", c(500,100), "300m")
spplot(pred["var1.pred"], sp.layout=list(l2,scale,text1, text2), col.regions=bpy.colors(20), key.space=list(x=0.1,y=.95,corner=c(0,1)), main = "MAPA DENSIDADE - ELEVAÇÃO", xlab="XCOORD(m)", ylab= "YCOORD(m)" , xlim=c(0, 800), ylim=c(0,2100))






Diagrama de bloco

s.grid <- GridTopology(c(100,100),c(10,40),c(50,50))
s.grid <- SpatialPoints(s.grid)
data(oxford)
m <- vgm(367.66942, "Gau", 356.2184, 19.78417)
xx <- krige(ELEV~1, ~XCOORD+YCOORD, model = m, data = oxford, newd = s.grid)
dfxx <- as.data.frame(xx)
mz <- matrix(dfxx[,3], nrow=50, ncol=50, byrow=FALSE)
nmz <- matrix(nrow=50, ncol=50)
for (i in 1:50)
  for (j in 1:50)
  {nmz[i,j]=mz[i,51-j]}
persp(x =seq(100,599,by= 10), y=seq(100,2099,by=40),nmz,xlab="XCOORD",ylab="YCOORD",
      main="Diagrama de Bloco ELEV", theta=30, phi=30, r=100, d=20, scale=FALSE, col="green")




Validação cruzada

data(oxford)
coordinates(oxford) <- ~XCOORD+YCOORD
m <- vgm(367.66942, "Gau", 356.2184, 19.78417)
x <- krige.cv(ELEV~1, oxford, m, nmax = 40, nfold=5)
bubble(x, "residual", main = "ELEV: 5-fold CV
residuals")


Simulação


s.grid<-GridTopology(c(0,0),c(6,21),c(100,100))
s.grid<-SpatialPoints(s.grid)
gridded(s.grid)<-TRUE
data(oxford)
m <- vgm(367.66942, "Gau", 356.2184, 19.78417)
xx <- krige(ELEV~1, ~XCOORD+YCOORD, model = m,
            data = oxford, newd = s.grid, nsim=4, nmax=10)
spplot(xx["sim1"],xlab="XCOORD(0-4)",ylab="YCOORD(0-4)",main="Valores
Simulados-sim1")





Comentários

Postar um comentário