✎ software:art:scripts:scriptgm-pj

[ Por favor, preencha este local (:wiki:site_notice) ]

Você está aqui: start » software » art » scripts » scriptgm-pj

Você carregou uma revisão antiga desse documento! Se você salvá-la, irá criar uma nova versão com esses dados.

Script by PJ on 16/11/2007

##
## 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)
## para arquivo do IBGE (Vila Alta == Alto Paraíso)

figado$municipio <- as.character(figado$municipio)
figado[figado$municipio == "Alto Paraíso", "municipio"] <- "Vila Alta"

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 campo a mais que o do DATASUS
figado$id <- as.character(figado$id)

cbind(figado$id, municDt$CODIGO)
figado$id <- paste(figado$id,substr(municDt$CODIGO, 7, 7),sep="")
all.equal(municDt$CODIGO,figado$id)

## Q: suponha agora que eu tenha varios atributos no R, provenientes de outra
##    fonte e associados aos polígonos como acima.
##    Qual a melhor maneira de adicionar à tabela já existente (Paraná)?
##    (tem que ser coluna a coluna com updatecolumns ou melhor abrir outra tabela ?
## 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:
## (o que tentei a seguir nao funcionou, tem que ser coluna a coluna...???)

## 1.1 tentando adcionar todas as colunas de uma só vez:

## CODIGO foi usado como ID para munic no addShape() portanto uando aqui tb para figado
all.equal(municDt$CODIGO,figado$id)
municTab
updateColumns(municTab, figado) ## nao funcionou mas pelo vignette aRTtable deveria...
## Erro em if (poskey != 1) { : argumento tem comprimento zero

## 1.2 adicionando as colunas uma a uma... (arrrgh)
names(figado)
for(cc in names(figado)[3:4]){
  tmp <- data.frame(id=figado$id, data = figado[,cc])
  names(tmp)[2] <- cc
  print(head(tmp))
  createColumn(municTab, cc, type="int")
  updateColumns(municTab, tmp)
}
## mesmo erro que acima:
##Erro em if (poskey != 1) { : argumento tem comprimento zero

## 2. adicionando como uma nova tabela

## 2.1 usando creatTable()
figadoTab <- createTable(l.pr, "CancerFigado")
figadoTab
l.pr
updateColumns(figadoTab, figado)
figadoTab

## usando importTable()
figadoTab1 <- importTable(l.pr, "FigadoCancer", id="id", figado)

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?

##
## 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)
coordinates(teores) <- c("LONGITUDE","LATITUDE")
teores$ID <- as.character(1:nrow(teores))
## 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")
## 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=";")
teores$ID <- as.character(1:nrow(teores))
#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")
#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")
teores$ID <- as.character(1:nrow(teores))
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...)
##
## alternativa a importSpData() seria fazer passo a passo
teores <- read.table("aguafix.csv", head=T, sep=";")
teores$ID <- as.character(1:nrow(teores))
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 etando em outro layer !!!!!!
tb4.geoq <- importTable(l4.geoq, "Elem", id="ID", lldf)
l4.geoq
tb4.geoq
</code>

Por favor digite todas as letras dentro da caixa para provar que você é humano. W A R Y D Por favor deixe esse campo em branco:

Resumo da edição

Observe: editando esta página você aceita disponibilizar o seu conteúdo sob a seguinte licença: GNU Free Documentation License 1.3