Script by PJ on 18/11/2007 (23:00 BRT)

##
## Importação dos dados de arquivos
## para o R e para um banco Terralib/MySQL
##    - teores de elementos medidos em água e localizações das amostras
##    - poligonos com municípios do Paraná
##
## carregando os pacotes que serão usados (aRT)
require(aRT)
 
##Conectar com o SGBD MYSQL (local)
con=openConn()
con
 
##Apagar banco já existente
if(any(showDbs(con)=="geomedicina")) deleteDb(con, "geomedicina", force=T)
##Criando um novo banco
db=createDb(con, "geomedicina")
db
## string definindo a projeção
proj="+proj=latlong +ellps=GRS67 +towgs84=-66.8700,-4.3700,-38.5200"
##
## Layer 1: municípios
##
## Criar no BD layer para polígonos dos municípios
l.pr <- createLayer(db, "Parana", proj=proj)
l.pr
 
##Adicionar ao BD o contorno do Paraná disponível em shapefile (fonte: IBGE)
addShape(l.pr, tab="Parana", file="parana_pol.shp", id="CODIGO", length=10)
## verificando status do banco e layer
db
l.pr
 
## visualizando polygonos
plot(l.pr)
 
## importando polygonos do banco para o R
munic <- getPolygons(l.pr)
class(munic)
 
## Q: qual a diferença em usar os dois comandos abaixo?
aRTplot(munic)
plot(munic)
 
## importando atributos dos polygonos do BD para o R
municTab <- openTable(l.pr, "Parana")
municTab
municDt <- getData(municTab)
dim(municDt)
names(municDt)
head(municDt)
 
## centroids e áreas de cada município
center <- getOperation(l.pr, "centroid")
center
## Q: como converter para UTM (ou obter um UTM) ???
 
## Q: area deveria usar UTM para o cálculo?
## Q: área calculada com LatLong é correta? 
area <- getMetric(l.pr, "area")
area
 
##
## Lendo atributos dos poligonos (munícipios) de outra fonte (DATASUS)
## (serão necessárias operações para compatibilizar com dados do IBGE)
figado <- read.csv("cancerfigadopop2005.csv")
head(figado)
 
names(figado)
names(municDt)
 
## checando e corrigindo **nomes** dos municípios:
 
## um município trocou o nome.....
## fix (temporário até obter mapa atualizado do IBGE.... se tiver...
## no site do IBGE consta Alto Paraíso -- o nome atual--, mas no arquivo apra download ainda está
## Alta Vista, mesmo para 2005)
## para arquivo do IBGE (Vila Alta == Alto Paraíso)
 
figado$municipio <- as.character(figado$municipio)
figado[figado$municipio == "Alto Paraíso", "municipio"] <- "Vila Alta"
 
## verificando igualdade entre nomes de Municípios
all.equal(municDt$NOME,figado$municipio)
##estão em ordem diferente...
all.equal(sort(municDt$NOME),sort(figado$municipio))
## ainda 1 incongruencia:
foo1 <- which(sort(municDt$NOME) != sort(figado$municipio))
cbind(sort(municDt$NOME),sort(figado$municipio))[foo1,]
foo2 <- cbind(sort(municDt$NOME),sort(figado$municipio))[foo1,]
foo2
municDt[municDt$NOME == foo2[1], "NOME"] <- foo2[2]
## nomes agora OK!
all.equal(sort(municDt$NOME),sort(figado$municipio))
 
## reordenando segundo a ordem do IBGE 
figado <- figado[match(municDt$NOME, figado$municipio),]
all.equal(municDt$NOME,figado$municipio)
 
## checando e corrigindo **CODIGOS** dos municípios:
## código do IBGE tem 1 dígito a mais que o do DATASUS ...
## vamos extrair do IBGE e anexar ao do DATASUS
figado$id <- as.character(figado$id)
 
head(cbind(figado$id, municDt$CODIGO))
figado$id <- paste(figado$id,substr(municDt$CODIGO, 7, 7),sep="")
all.equal(municDt$CODIGO,figado$id)
 
##
## Adicionando atributos de outra fonte de dados aos polígonos
##
## No exemplo abaixo quero adicionar o data frame "figado" ao banco 
 
## tentei pegar is ID's da tabela que veio do banco mas não deu...
#> getID(municDt)
#Erro em getID(municDt) : Not implemented yet
#> getID(municTab)
#Erro em getID(municTab) : Not implemented yet
 
## 1. adicionando à tabela já existente:
## 1.1 Adicionar/atualizar todas as colunas de uma só vez: updateColumns()
## Erros aconteciam abaixo qdo usava-se o nome errado da variável id:
## a msg poderia ser informativa?
## msg de erro era:
## Erro em if (poskey != 1) { : argumento tem comprimento zero
 
all.equal(municDt$CODIGO,figado$id)
municTab
names(figado)
 
## note que para updateColumns() a doc diz que o id deve estar na primeira coluna do data-frame!
## o ID deve ter o mesmo nome da tabela original (CODIGO neste exemplo)
## CODIGO foi usado como ID para munic no addShape() portanto usando aqui tb para figado
names(figado)[1] <- "CODIGO"
updateColumns(municTab, figado)
 
## 1.2 alternativa adicionando as colunas uma a uma... 
## note que o ID dev e ser chamado de CODIGO que foi definido anteriormente
#names(figado)
#for(cc in names(figado)[3:4]){
#  tmp <- data.frame(CODIGO=figado$id, data = figado[,cc])
#  names(tmp)[2] <- cc
#  print(head(tmp))
#  createColumn(municTab, cc, type="int")
#  updateColumns(municTab, tmp)
#}
 
## 2. adicionando como uma nova tabela
## 2.1 usando importTable()
figadoTab1 <- importTable(l.pr, "FigadoCancer", id="id", figado)
l.pr
 
## 2.2 usando creatTable()
figadoTab <- createTable(l.pr, "CancerFigado")
figadoTab
l.pr
updateColumns(figadoTab, figado)
figadoTab
l.pr
 
## Q: Pedro. qual a diferença entre usar createTable() e importTable() ?
##    A segunda simplesmente encapsula os comandos da primeira em uma unica chamada?
 
## removendo tabela(s):
## o seguinte não funcionou:
deleteTable(l.pr, "FigadoCancer")
#Erro em function (classes, fdef, mtable)  : 
#  unable to find an inherited method for function "deleteTable", for signature "aRTlayer"
 
## e o seguinte deu um crash!
#deleteTable(db, "FigadoCancer")
 
##
## Layer 2: dados geoquímicos
##
## Q:
## Pedro, aqui estou tentando usar 3 formas diferentes de colocar
## o dado de pontos no banco, mas as baseadas im importSpData()
## estao dando problemas
 
##
## Forma 1: (não funcionou, veja msg de erro abaixo)
##
teores <- read.table("aguafix.csv", head=T, sep=";")
head(teores)
summary(teores)
plot(munic)
points(teores[,c("LONGITUDE", "LATITUDE")], pch=20, cex=0.3, col="blue")
 
## convertendo teores para SpatialPointsDataFrame
names(teores)
teores$ID <- as.character(1:nrow(teores))
coordinates(teores) <- c("LONGITUDE","LATITUDE")
## Q: como atribuir projeção ao Layer criado com importSpData()? O seguinte é válido?
attr(teores, "proj4string") <- CRS(proj)
l1.geoq <- importSpData(db, teores, "GQ1", "geoq1")
## MSG DE ERRO:
#Erro em .aRTcall(object, "cppUpdateColumns", colnames = colnames(data),  : 
#  Could not update the table
l1.geoq
 
##
## Forma 2: (não funcionou, veja msg de erro abaixo)
##
teores <- read.table("aguafix.csv", head=T, sep=";")
names(teores)[1] <- "ID"
teores$ID <- as.character(teores$ID)
##ll <- SpatialPoints(teores[,c("LONGITUDE", "LATITUDE")], proj=CRS(proj))
lldf <- SpatialPointsDataFrame(coords=teores[,c("LONGITUDE", "LATITUDE")],
                               data=teores[,-(4:5)],
                               proj = CRS(proj))
l12.geoq <- importSpData(db, lldf, "GQ2", "geoq2")
#Erro em .aRTcall(object, "cppUpdateColumns", colnames = colnames(data),  : 
#  Could not update the table
 
##
## Forma 3: (Funcionou)
##
## alternativa a importSpData() seria fazer passo a passo
teores <- read.table("aguafix.csv", head=T, sep=";")
coordinates(teores) <- c("LONGITUDE","LATITUDE")
names(teores)[1] <- "ID"
teores$ID <- as.character(teores$ID)
attr(teores, "proj4string") <- CRS(proj)
 
l3.geoq <- createLayer(db, "GQ3", proj=proj)
l3.geoq
addPoints(l3.geoq, teores)
l3.geoq
tb3.geoq <- importTable(l3.geoq, "Elementos", id="ID", teores)
l3.geoq
tb3.geoq
 
##
## Forma 4: (Funcionou, mas...)
##
teores <- read.table("aguafix.csv", head=T, sep=";")
names(teores)[1] <- "ID"
teores$ID <- as.character(teores$ID)
lldf <- SpatialPointsDataFrame(coords=teores[,c("LONGITUDE", "LATITUDE")],
                               data=teores[,-(4:5)],
                               proj = CRS(proj))
 
l4.geoq <- createLayer(db, "GQ4", proj=proj)
l4.geoq
addPoints(l4.geoq, lldf)
l4.geoq
## ele nao aceitou usar o nome da tabela "Elementos" já usado acima
##  mesmo estando em outro layer !!!!!!
tb4.geoq <- importTable(l4.geoq, "Elem", id="ID", lldf)
l4.geoq
tb4.geoq