##Abrir o aRT
require (aRT)
##Conectar com o MYSQL
con=openConn()
##Apagar banco já existente
if(any(showDbs(con)=="geomedicina")) deleteDb(con, "geomedicina", force=T)
##Criar novo banco
db=createDb(con, "geomedicina")
##Criar novo layer
proj="+proj=latlong +ellps=GRS67 +towgs84=-66.8700,-4.3700,-38.5200"
l=createLayer(db, "Parana", proj=proj)
##Adicionar o contorno do Paraná
addShape(l, tab="Parana", file="parana_pol.shp", id="CODIGO", length=10)
##Ler tabela de dados químicos
agua<-read.table("agua2.csv", header=T, sep=";")
grade<-read.table("grade.csv", header=T, sep=";")
##Plotar os dados juntos
poly=getPolygons(l)
aRTplot(poly)
points(grade$LONGITUDE, grade$LATITUDE, pch=20, cex=0.3, col="blue")
points(agua$LONGITUDE,  agua$LATITUDE, pch=16, cex=0.3)


##Unir municípios em micro-regiões
require(foreign)
atr=read.dbf('41mu2500gc.dbf')
names(atr)[1]='CODIGO'
microgeom=list()
micro=unique(atr$MICROREGIA)

for(i in 1:length(micro))
{
  # selecionando municipios dentro da micro-regiao i
  id=as.character(atr$CODIGO[which(atr$MICROREGIA == micro[i])])
  cat(paste("Processando micro-regiao",micro[i], "com",length(id), "municipios\n"))
  # unindo os poligonos dos municipios para formar pol. da regiao
  set=getSetOperation(l,"union",id=id)
  set@polygons[[1]]@ID = paste(micro[i])
  # adicionando ao objeto (list) das regioes
  microgeom[i] = set@polygons
  aRTplot(set, col=terrain.colors(length(micro)+1)[i+1], add=T, lwd=2)
}
aRTplot(poly, add=T)
##convertendo para classe SpatialPolygons
res = SpatialPolygons(microgeom, 1:length(microgeom))

##insere os dados das micro-regioes, e cria temas para visualizacao
l2=createLayer(db,l="microreg", proj=proj)
addPolygons(l2, res)
aRTplot(res)
points(agua$LONGITUDE,  agua$LATITUDE, pch=16, cex=0.3)


#idf <- which(sapply(agua, class)=='factor')
#idf
#sapply(agua[, idf], levels)


aguadbf=read.dbf('lgqbd_agua.dbf')
#id.numer <- which(sapply(aguadbf, is.numeric))
#id.numer

require(splancs)
medias <- t(sapply(1:39, function(i) {
  id <- which(inout(cbind(agua$LONGITUDE, agua$LATITUDE),
                    res@polygons[[i]]@Polygons[[1]]@coords))
  sapply(aguadbf[id, id.numer], mean)
}))
dim(medias)
colnames(medias)