## ## Exemplos de análises espaciais ## ## Exemplo 1 ## require(spatstat) data(redwoodfull) plot(redwoodfull) redwoodfull.extra$plot() ## ## Exemplo 2 ## library(maptools) library(spdep) data(auckland) brks <- c(-Inf, 2, 2.5, 3, 3.5, Inf) cols <- c("#F7F7F7", "#CCCCCC", "#969696", "#636363", "#252525") res <- probmap(auckland$Deaths.1977.85, 9*auckland$Under.5.1981) plot(auckpolys, col=cols[findInterval(res$raw*1000, brks)], forcefill=FALSE) legend(c(70,90), c(70,95), fill=cols, legend=leglabs(brks), bty="n") ## ## Exemplo 3: package spatial - processo pontual ## require(spatial) towns <- ppinit("towns.dat") 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") ## ## Exemplo 4: package spatial - geoestatistica ## require(spatial) require(MASS);require(geoR) data(topo, package="MASS") points(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) eqscplot(sesurf, type = "n") contour(sesurf, levels = c(22, 25), add = TRUE) ## ## Exemplo 5: splancs - Função K ## require(splancs) data(cardiff) plot(cardiff, asp=1) UL.khat <- Kenv.csr(length(cardiff$x), cardiff$poly, nsim=29, 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 6: splancs - kernel2d ## data(bodmin) plot(bodmin, asp=1) plot(bodmin$poly, asp=1, type="n") image(kernel2d(as.points(bodmin), bodmin$poly, h0=2, nx=100, ny=100), add=TRUE, col=terrain.colors(20)) pointmap(as.points(bodmin), add=TRUE, pch=19) polymap(bodmin$poly, add=TRUE, lwd=2) ## ## Exemplo 7: spatstats - simulações de processos ## X <- rThomas(15, 0.2, 5) plot(X) plot(Gest(X)) pp <- rSSI(0.05, 200) plot(pp) plot(Gest(pp)) ## ## Exemplo 8: spdep - Empirical Bayes ## library(maptools) library(spdep) data(auckland) summary(auckland$Under.5.1981) summary(auckland$Deaths.1977.85) res <- probmap(auckland$Deaths.1977.85, 9*auckland$Under.5.1981) # taxas brutas anuais X11() 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)], forcefill=FALSE) legend(c(70,90), c(70,95), fill=cols, legend=leglabs(brks), bty="n") # taxas padronizadas brks <- c(-Inf,47,83,118,154,190,Inf) cols <- c("#F7F7F7", "#CCCCCC", "#969696", "#636363", "#252525") plot(auckpolys, col=cols[findInterval(res$relRisk, brks)], forcefill=FALSE) legend(c(70,90), c(70,95), fill=cols, legend=leglabs(brks), bty="n") # taxas modelo Poisson brks <- c(0,0.05,0.1,0.2,0.8,0.9,0.95,1) cols <- c("#2166AC", "#67A9CF", "#D1E5F0", "#F7F7F7", "#FDDBC7", "#EF8A62", "#B2182B") plot(auckpolys, col=cols[findInterval(res$pmap, brks)], forcefill=FALSE) legend(c(70,90), c(70,95), fill=cols, legend=leglabs(brks), bty="n") # taxas segundo Bayes Empírico brks <- c(-Inf,2,2.5,3,3.5,Inf) cols <- c("#F7F7F7", "#CCCCCC", "#969696", "#636363", "#252525") res <- EBest(auckland$Deaths.1977.85, 9*auckland$Under.5.1981) plot(auckpolys, col=cols[findInterval(res$estmm*1000, brks)], forcefill=FALSE) legend(c(70,90), c(70,95), fill=cols, legend=leglabs(brks), bty="n") ## ## Exemplo 9: geoR - simulação ## require(geoR) ## ## Comparando simulações com diferentes valores de $\phi$ ## for(i in 0:9){ jpeg(paste("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() } ## ## Exemplo 10: RandomFields - simulação ## require(RandomFields) # processo Gaussiano x <- y <- seq(0,20,0.1) f <- GaussRF(x=x, y=x, model="stable", grid=TRUE,param=c(0,4,1,10,1)) image(x, x, f, col=gray(seq(1,0,l=15))) # Processos Max-Stable com marginais Fréchet unitárias n <- 100 x <- y <- 1:n ms <- MaxStableRF(x, y, grid=TRUE, model="exponen", param=c(0,1,0,40), maxstable="extr") image(x,y,ms, col=gray(seq(1,0,l=15))) ## ## Exemplo 11: geoR - análise Bayesiana ## 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) ## ## Exemplo 12: Programando métodos de estatística espacial em R ## set.seed(234) source("~/Projetos/Rcitrus/ar/InDisp.r") source("~/Projetos/Rcitrus/ar/conv.R") source("~/Projetos/Rcitrus/ar/MinDist.r") ## Lendo os dados arara2 <- read.csv2('~/Projetos/Rcitrus/ar/dcast2.csv') arara2 <- arara2[, 2:dim(arara2)[2]] da2 <- conv(arara2, rd=4, cd=7, ret="ge") # Convertendo os dados de Castanharo 2 ## Visualizando plot(da2$coords, asp=1) points(da2$coords[da2$data==1,], col=2, pch=20) ## aplicando o teste de minima distancia ad.a2 <- RDist.test(da2, 1, 99) ad.a2 hist(ad.a2$rd$V1) abline(v=ad.a2$obs) ### Analise por quadrats - talhoes de Castanharo e Itajobi RA2 <- InDisp1(arara2, q = matrix(rep(2:13, 2), 12, 2)) RA2[,4:8] <- round(RA2[,4:8], dig=4) RA2