rm(list=ls())
##
## Exemplos de análises espaciais
##
## Exemplo 1
##
require(spatstat)
data(redwoodfull)
plot(redwoodfull)
redwoodfull.extra$plot()
##
rm(list=ls())
##
## Exemplo 2
##
library(spdep)
##data(auckland)

auckland <- sf::st_read(system.file("shapes/auckland.gpkg", package="spData")[1], quiet=TRUE)
help(auckland)
head(auckland)
dim(auckland)
str(auckland)
plot(auckland)
plot(auckland[c("M77_85", "Und5_81")], col=gray(seq(1,0,l=21)), border="navyblue", lwd=2)

auckpolys <- sf::st_geometry(auckland)
plot(auckpolys)
##plot(auckpolys, col="gray80", border="darkgray", lwd=0.5)

brks <- c(-Inf, 5, 10, 20, 30, Inf)
cols <- gray(seq(1,0, len=5))
plot(auckpolys, col=cols[findInterval(auckland$M77_85, brks)])
legend(c(60,90), c(70,95), fill=cols, legend=SpatialEpi:::leglabs(brks), bty="n")

res <- with(auckland, probmap(M77_85, 9*Und5_81))
head(res, n = 10)
head(round(cbind(auckland[,3:4, drop = TRUE],res), dig=4))
summary(res)

## alternativa para facilitar: combinar dados e resultados em um único objeto
res.sf <- cbind(auckland, res)
head(res.sf)
head(res.sf[,-c(1,2,9), drop = T])

x11()
plot(res[,"expCount"], auckland[,"M77_85", drop = TRUE], asp = 1)
abline(0,1)
dev.off()


brks <- c(-Inf, 2, 2.5, 3, 3.5, Inf)
cols <- c("#F7F7F7", "#CCCCCC", "#969696", "#636363", "#252525")
plot(auckpolys, col=cols[findInterval(res$raw*1000, brks)])
legend(c(60,90), c(70,95), fill=cols, legend=SpatialEpi:::leglabs(brks), bty="n")

brks <- c(-Inf, 50, 100, 200, 500, Inf)
cols <- c("#F7F7F7", "#CCCCCC", "#969696", "#636363", "#252525")
plot(auckpolys, col=cols[findInterval(res$relRisk, brks)])
legend(c(60,90), c(70,95), fill=cols, legend=SpatialEpi:::leglabs(brks), bty="n")

brks <- c(-Inf, 0.1, 0.9, Inf)
cols <- gray(seq(1,0, len=3))
plot(auckpolys, col=cols[findInterval(res$pmap, brks)])
legend(c(60,90), c(70,95), fill=cols, legend=SpatialEpi:::leglabs(brks), bty="n")

# taxas Bayesianas Empíricas (global)
EBglobal <- with(auckland, EBest(M77_85, 9*Und5_81))
head(EBglobal)

brks <- c(-Inf,2,2.5,3,3.5,Inf)
cols <- c("#F7F7F7", "#CCCCCC", "#969696", "#636363", "#252525")
plot(auckpolys, col=cols[findInterval(EBglobal$estmm*1000, brks)])
legend(c(60,90), c(70,95), fill=cols, legend=SpatialEpi:::leglabs(brks), bty="n")
##
par(mfrow=c(1,2), mar=c(0,0,0,0))
brks <- c(-Inf,2,2.5,3,3.5,Inf)
cols <- c("#F7F7F7", "#CCCCCC", "#969696", "#636363", "#252525")
plot(auckpolys, col=cols[findInterval(res$raw*1000, brks)])
legend(c(60,90), c(70,95), fill=cols, legend=SpatialEpi:::leglabs(brks), bty="n")
plot(auckpolys, col=cols[findInterval(EBglobal$estmm*1000, brks)])
legend(c(60,90), c(70,95), fill=cols, legend=SpatialEpi:::leglabs(brks), bty="n")
dev.off()

x11()
plot(res$raw, EBglobal$estmm, asp = 1, pch = 19, cex = 0.5)
abline(0, 1)
dev.off()

# taxas Bayesianas Empíricas (local)
auckland.nb <- spdep::poly2nb(auckland)
head(auckland.nb)
length(auckland.nb)
EBlocal <- with(auckland, EBlocal(M77_85,  9*Und5_81, auckland.nb))

brks <- c(-Inf,2,2.5,3,3.5,Inf)
cols <- c("#F7F7F7", "#CCCCCC", "#969696", "#636363", "#252525")
plot(auckpolys, col=cols[findInterval(EBlocal$est*1000, brks)])
legend(c(60,90), c(70,95), fill=cols, legend=SpatialEpi:::leglabs(brks), bty="n")

par(mfrow=c(1,3), mar=c(0,0,0,0))
brks <- c(-Inf,2,2.5,3,3.5,Inf)
cols <- c("#F7F7F7", "#CCCCCC", "#969696", "#636363", "#252525")
plot(auckpolys, col=cols[findInterval(res$raw*1000, brks)])
legend(c(60,90), c(70,95), fill=cols, legend=SpatialEpi:::leglabs(brks), bty="n")
plot(auckpolys, col=cols[findInterval(EBglobal$estmm*1000, brks)], forcefill=FALSE)
legend(c(60,90), c(70,95), fill=cols, legend=SpatialEpi:::leglabs(brks), bty="n")
plot(auckpolys, col=cols[findInterval(EBlocal$est*1000, brks)])
legend(c(60,90), c(70,95), fill=cols, legend=SpatialEpi:::leglabs(brks), bty="n")

x11()
par(mfrow=c(1,2), mar = c(3,3,0,0), mgp = c(2,1,0))
plot(res$raw, EBglobal$estmm, asp = 1, pch = 19, cex = 0.5)
abline(0, 1)
plot(res$raw, EBlocal$est, asp = 1, pch = 19, cex = 0.5)
abline(0, 1)
dev.off()

##
rm(list=ls())
##
## Exemplo 3
##
require(spatial)
towns <- ppinit("towns.dat")
names(towns)
plot(towns, asp=1)
par(pty="s")
plot(Kfn(towns, 40), type="b")
plot(Kfn(towns, 10), type="b", xlab="distance", ylab="L(t)")
for(i in 1:10) lines(Kfn(Psim(69), 10))
lims <- Kenvl(10,100,Psim(69))
lines(lims$x,lims$l, lty=2, col="green")
lines(lims$x,lims$u, lty=2, col="green")
lines(Kaver(10,25,Strauss(69,0.5,3.5)), col="red")
##
rm(list=ls())
##
## Exemplo 4 
##
require(spatial)
data(topo, package="MASS")
points(geoR::as.geodata(topo))

topo.kr <- surf.gls(2, expcov, topo, d=0.7)
trsurf <- trmat(topo.kr, 0, 6.5, 0, 6.5, 50)
contour(trsurf, add = TRUE)

prsurf <- prmat(topo.kr, 0, 6.5, 0, 6.5, 50)
contour(prsurf, levels=seq(700, 925, 25))
sesurf <- semat(topo.kr, 0, 6.5, 0, 6.5, 30)
MASS::eqscplot(sesurf, type = "n")
contour(sesurf, levels = c(22, 25), add = TRUE)

require(geoR)
topo.gd <- as.geodata(topo)
plot(variog(topo.gd, max.dist=4, trend="2nd"))
ml <- likfit(topo.gd, ini=c(1000, 0.7))
gr <- expand.grid(seq(0,6.5, l=50),seq(0,6.5, l=50))
kr <- krige.conv(topo.gd, loc=gr, krige=krige.control(obj.m=ml))
contour(kr)
image(kr, col=gray(seq(1,0,l=21)), asp=1)
##
rm(list=ls())
##
## Exemplo 5
##
require(splancs)
data(bodmin)
plot(bodmin, asp=1)
plot(bodmin$poly, asp=1, type="n")
bodmin.kernel <- kernel2d(as.points(bodmin), bodmin$poly, h0=2, nx=100, ny=100)
image(bodmin.kernel, add=TRUE, col=terrain.colors(20))
pointmap(as.points(bodmin), add=TRUE, pch=19)
polymap(bodmin$poly, add=TRUE, lwd=2)
##
## Exemplo 6
##
require(splancs)
data(cardiff)
plot(cardiff, asp=1)
UL.khat <- Kenv.csr(length(cardiff$x), cardiff$poly, nsim=99, seq(2,30,2))
plot(seq(2,30,2), sqrt(khat(as.points(cardiff), cardiff$poly, seq(2,30,2))/pi)-seq(2,30,2), type="l", xlab="Splancs - polygon boundary", ylab="Estimated L", ylim=c(-1,1.5))
lines(seq(2,30,2), sqrt(UL.khat$upper/pi)-seq(2,30,2), lty=2)
lines(seq(2,30,2), sqrt(UL.khat$lower/pi)-seq(2,30,2), lty=2)
##
## Exemplo 7
##
x11()
plot(c(0,1), c(0,1), asp=1, ty="n")
polygon(rbind(c(0,0), c(1,0), c(1,1), c(0,1)))
pp <- locator(50, type="p")
## clique em 50 pontos completamente ao acaso dentro do quadrado
borda <- rbind(c(0,0), c(1,0), c(1,1), c(0,1), c(0,0))
Kseq <- seq(0.1, 0.25, l=20)
pp.K <- khat(as.points(pp), borda, Kseq) 
pp.Kenv <- Kenv.csr(length(pp$x), borda, nsim=29, Kseq)
pp.low <- sqrt(pp.Kenv$lower/pi)-Kseq
pp.upp <- sqrt(pp.Kenv$upper/pi)-Kseq
plot(Kseq, sqrt(pp.K/pi) - Kseq, ty="l", ylim = 1.2*c(min(pp.low), max(pp.upp)))
lines(Kseq, pp.low, lty=2)
lines(Kseq, pp.upp, lty=2)
plot(Kseq, sqrt(pp.K/pi) - Kseq, ty="l", ylim = 1.2*c(min(pp.low), max(pp.upp)))
lines(Kseq, pp.upp, lty=2)
lines(Kseq, pp.low, lty=2)
##
dev.off()
rm(list=ls())
##
## Exemplo 8 
##
require(spatstat)
X <- rThomas(15, 0.2, 5)
plot(X)
plot(Gest(X))

pp <- rSSI(0.05, 200)
plot(pp)
plot(Gest(pp))
##
rm(list=ls())
##
## Exemplo 9 
##
require(geoR)
##
## Comparando simulações com diferentes valores de $\phi$
##
for(i in 0:9){
    system("mkdir anima_phi")
    jpeg(paste("anima_phi/phi",i, ".jpg",  sep=""), wid=600, hei=600)
    par(mfrow=c(2,2), mar=c(1.5,.5,1.5,0), mgp=c(1, .5, 0))
    set.seed(234+i)
    ap1 <- grf(961, grid="reg", cov.pars=c(1, 0))
    set.seed(234+i)
    ap2 <- grf(961, grid="reg", cov.pars=c(1, .1))
    set.seed(234+i)
    ap3 <- grf(961, grid="reg", cov.pars=c(1, .25))
    set.seed(234+i)
    ap4 <- grf(961, grid="reg", cov.pars=c(1, .75))
    iis <- range(c(ap1$data, ap2$data, ap3$data, ap4$data))
    image(ap1, xlab="", ylab="", col=gray(seq(1,0,l=21)), zlim=iis)
    mtext(expression(phi==0), cex=1.5)
    image(ap2, xlab="", ylab="", col=gray(seq(1,0,l=21)), zlim=iis)
    mtext(expression(phi==0.10), cex=1.5)
    image(ap3, xlab="", ylab="", col=gray(seq(1,0,l=21)), zlim=iis)
    mtext(expression(phi==0.25), cex=1.5)
    image(ap4, xlab="", ylab="", col=gray(seq(1,0,l=21)), zlim=iis)
    mtext(expression(phi==0.75), cex=1.5)
    dev.off()
}
library(magick)
frames <- list.files("anima_phi", full.names=TRUE, pattern="\\.jpg$")
foo <- image_read(frames)
foo <- image_animate(foo, fps = 1)
image_write(foo, "phis.gif")
system("eog phis.gif &")
##
rm(list=ls())
##
## Exemplo 10
##
require(geoR)
ex.data <- grf(70, cov.pars=c(10, .25))
plot.geodata(ex.data)
ex.grid <- as.matrix(expand.grid(seq(0,1,l=11), seq(0,1,l=11)))

PC <- prior.control(phi.discrete=seq(0, 2, l=21))
OC <- output.control(n.post=500)                     
ex.bayes <- krige.bayes(ex.data, loc=ex.grid, prior = PC, output=OC)
names(ex.bayes)

plot(ex.bayes)

op <- par(no.readonly = TRUE)
par(mfrow=c(2,2))
par(mar=c(3,4,1,1))
par(mgp = c(2,1,0))
image(ex.bayes, main="predicted values")
image(ex.bayes, val="variance", main="prediction variance")
image(ex.bayes, val= "simulation", number.col=1,
      main="a simulation from the \npredictive distribution")
image(ex.bayes, val= "simulation", number.col=2,
      main="another simulation from \nthe predictive distribution")
par(op)
##
rm(list=ls())


## Exemplo 11
##
## Bairros de Curitiba
##
## lendo dados tipo shapefiles
require(sf)
require(RColorBrewer) ## usada aqui para criar palhetas de cores nas visualizações em mapas
require(classInt)     ## rotinas para faclitar a divisão de dados em classes por vários critérios
require(spdep)        ## funções análises de dados de áreas

##
## Deve baixar dados/arquivos da página do curso
##
## bairros.dbf  bairros.shp  bairros.shx  tabe.csv

ctba <- sf::st_read("bairros.shp", quiet=TRUE, stringsAsFactors=FALSE)
plot(ctba)
plot(sf::st_geometry(ctba))

head(ctba)
ctba$NOME <- as.character(ctba$NOME)
## mudança de codificação de caracteres (só use se necessário, 
## depende do sistema operacional e do enconding do sistema)
## Veja os nomes dos unicípios e se os acentos aparecem corretamente
Encoding(ctba$NOME)
Encoding(ctba$NOME) <- "latin1"
ctba$NOME <- enc2native(ctba$NOME)
head(ctba)

tb <- read.table("tabe.csv", sep=";", head=T, enc="latin1")
head(tb)
str(tb)

## e agora ordenando os dados corretamente, para compatibilizar a ordem
## dos minicípios no shape e na tabela de atributos
head(ctba)
head(tb)
ind <- match(ctba$CODE, tb$CODE)
ind

tb <- tb[ind,]
row.names(tb) <- ctba$CODE

ctba <- cbind(ctba, tb)
head(ctba)
names(ctba)

## Visualizando um mapa do atributo: Segurança
#INT <- classIntervals(ctba$Segurança, n=5, style="quantile")
INT <- classIntervals(ctba$Segurança, style="fixed", fixedBreaks=c(0, 20, 40, 60, 80, 100))
CORES.5 <- c(rev(brewer.pal(3, "Blues")), brewer.pal(3, "Reds"))[-1]  ## tirando uma cor no final
COL <- findColours(INT, CORES.5)
plot(ctba["Segurança"], col=CORES.5)
TB <- attr(COL, "table")
legtext <- paste(names(TB))  #, " (", TB, ")", sep="")
legend("bottomright", fill=attr(COL, "palette"), legend=legtext, bty="n")

## experimente outros estilos de "fatiamento" na função classIntervals
## e note quemo o aspecto visual se altera

## calculando vizinhanças (segundo um critério de vizinhança -- note que há outros!!!)
ctba.nb1 <- poly2nb(ctba)
args(poly2nb)
class(ctba.nb1)
ctba.nb1[[1]]
ctba[6,]
ctba.nb1[[6]]
ctba[ctba.nb1[[6]],"NOME"]

## OBS: outros tipos de vizinhanças podem ser calculados:
##      - mudando argumentos de poly2nb
##      - usando outras funções de vizinhança tal como:  dnearneigh() and knearneigh()

## Moran global
moran.test(ctba$Segurança, listw=nb2listw(ctba.nb1))
args(moran.test)

moran.mc(ctba$Segurança, listw=nb2listw(ctba.nb1), nsim=99)
args(moran.mc)

## Uma alternativa ao Moran: estatística de Geary
geary.test(ctba$Segurança, listw=nb2listw(ctba.nb1))