## ## Exemplos espaciais ## ## Exemplo 1 ## require(spatstat) data(redwoodfull) plot(redwoodfull) redwoodfull.extra$plot() data(swedishpines) X <- swedishpines) plot(X) X summary(X) plot(density(X, 10)) contour(density(X, 10), axes = FALSE) Q <- quadratcount(X, nx = 4, ny = 3) Q plot(X) plot(Q, add = TRUE, cex = 2) plot(envelope(X, fun=Gest, r=seq(0, sqrt(2)/6, by=0.005), nrank=2, nsim=99)) plot(envelope(X, fun=Fest, r=seq(0, sqrt(2)/6, by=0.005), nrank=2, nsim=99)) plot(envelope(X, fun=Kest, r=seq(0, sqrt(2)/6, by=0.005), nrank=2, nsim=99)) K <- Kest(X) plot(K) plot(dirichlet(X)) X <- runifpoint(100) plot(X) Z <- dirichlet(runifpoint(16)) plot(Z, add=T) plot(cut(X, Z)) plot(split(X, Z)) ## ## 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") ## BE local ## example(auckland) res <- EBlocal(auckland$M77_85, 9*auckland$Und5_81, auckland.nb) res <- EBlocal(auckland$Deaths.1977.85, 9*auckland$Under.5.1981, auckland.nb) brks <- c(-Inf,2,2.5,3,3.5,Inf) cols <- grey(6:2/7) plot(auckland, col=cols[findInterval(res$est*1000, brks, all.inside=TRUE)]) legend("bottomleft", fill=cols, legend=leglabs(brks), bty="n") title(main="Local moment estimator of infant mortality per 1000 per year") ## ## 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)))