Georreferenciando dados com comandos do R (polígonos)

Last modified by Gustavo Basso on 2016/07/24 01:46

    O processo de georreferenciamento consiste em associar uma informação à sua localização geográfica. Mais precisamente, a informação é associada a um elemento geográfico que pode ser:

  • polígono: uma sequência de pontos (coordenadas geográficas) e linhas que definem os limites de uma determinada área, que pode ser um município, UF, macrorregião, um país, etc.
  • linha: sequência de pontos e linhas, não necessariamente definindo uma região delimitada
  • ponto: uma coordenada geográfica, que pode ser a localização de uma ocorrência (processos pontuais) ou ter algum valor numérico associado (dados de superfície)

     Nesta página serão apresentados exemplos de georreferenciamento usando ferramentas do próprio R, disponíveis nos pacotes sp e maptools;. Mais especificamente, executaremos o georreferenciamento de dados a polígonos, que é o caso mais comum. Vale salientar, entretanto, que os pacotes citados trabalham também com os outros tipos de dados espaciais.
    Georreferenciar dados a polígonos consiste em associar a cada polígono uma ou mais estatísticas, que são referentes à região que este polígono representa. No Brasil, essa associação é feita principalmente através do Código Municipal do IBGE, no caso dos municípios. Uma vez que cada código IBGE identifica um polígono específico, a relação é possível quando a informação tabular, não georreferenciada, utiliza esse mesmo código para identificar os municípios.
    Para acelerar a execução e viabilizar o compartilhamento de arquivos, nesta página trabalharemos com dados de microrregiões. Os polígonos de microrregiões do Brasil podem ser obtidos nesse arquivo, criados através do processo documentado nesta página.

1. Preparando o dado tabular

    Neste exemplo iremos georreferenciar os dados de produção municipal de soja de 2013, disponíveis na base Produção Agrícola Municipal (TAB_0006) em formato csv. Mais especificamente, construiremos um novo arquivo shp que armazene não só os atributos geográficos dos polígonos (coordenadas, área, etc.), mas também os dados da PAM. Os comandos abaixo irão declarar como diretório de trabalho a pasta onde essa base está armazenada. Em seguida, irá carregá-la para a workspace do R.

### endereço padrão para a pasta "Bases Consolidadas"
pasta.ODR  <- "//MISRV54/cgma-dados/BASES_NOVO_ODR/Bases Consolidadas"
setwd(pasta.ODR)

### carregando base de pam municipal
pam <- read.csv2(file.path(pasta.ODR,"TAB_0005.csv"))

    A base contém dados de 64 produtos agrícolas em vários anos. É necessário filtrar essa base para 2013 e considerando apenas a produção de soja. Essas operações são executadas nos comandos a seguir.

## verificando os nomes existentes na coluna "chave.prod"
unique(pam$chave.prod)

## filtrando para 2013 e considerando apenas produção de soja
pam_13_soja <- subset(pam,ano==2013 & chave.prod=="soja")

    Agora é necessário agregar esses valores por microrregião. Antes, no entanto, é preciso buscar os códigos de microrregião na tabela de DTB (Divisão Territorial do Brasil, TAB_0005) e associar aos municípios da PAM. Os códigos abaixo carregam a DTB, associa seus dados à PAM (por ano e código do município) e agrega por microrregião as somas dos valores de ára plantada, de área colhida, da quantidade produzida e do valor de produção.

## carregando dtb (tab_0005)
dtb <- read.csv2("TAB_0005.csv")

## merge com base anterior
pam_13_soja <- merge(pam_13_soja,dtb[c("codigo_mun7","codigo_micro","ano")],by=c("codigo_mun7","ano"))

## agreação por microrregião
pam_13_soja_micro <- aggregate(cbind(area_plnt_dest,area_colh,qt_prod,valor_prod_1000) ~ codigo_micro,
                               data=pam_13_soja,FUN=sum)

2. Georreferenciando dados em polígonos: relação um para um

  Uma vez confeccionados os dados tabulares por microrregião, o objetivo agora é associar esses valores aos seus respectivos polígonos. Primeiramente, é necessário carregar a malha de microrregiões disponibilizada no início desta página. O comando readShapePoly do pacote maptools carrega o arquvo .shp no R. Essa operação é exemplificada no código abaixo. Antes de executá-lo, siga os passos abaixo:

  1. Faça o download da malha disponibilizada no início da página
  2. Descompacte o arquivo selecionando "extrair para malha_bra_mun13"
  3. No R, instale o pacote "maptoos". Se já tiver instalado, não precisa instalar novamente

    Ao executar a rotina, o R abrirá uma janela para que se escolha uma pasta. Selecione o diretório "malha_bra_mun13" para onde os arquivos foram descompactados

library(maptools)

## diretório dos shapefiles
dir.shapes <- choose.dir()
setwd(dir.shapes)

### malha municipal de 2013:
malha.bra.mun13 <- readShapePoly("malha_bra_mun13.shp")

   A malha carregada foi salva na variável malha.bra.mun13, que é um objeto da classe sp. Essa classe de objetos contem vários componentes armazenados em slots. O slot "data" armazena os atributos dos polígonos, enquanto o slot "polygons" contém as coordenadas que definem os limites de cada região, além de outros metadados como área, centróide, etc. O código abaixo mostra a classe e os diferentes slots de dados.

class(malha.bra.mun13)
slotNames(malha.bra.mun13)

   Os comandos acima devem retornar as seguintes saídas no console do R:

> class(malha.bra.mun13)
[1] "SpatialPolygonsDataFrame"
attr(,"package")
[1] "sp"
> slotNames(malha.bra.mun13)
[1] "data"        "polygons"    "plotOrder"   "bbox"        "proj4string"

   Uma vez que temos na workspace do R a malha municipal e a base de produção de soja de 2013, o georreferenciamento da última é feita através de um merge específico do pacote sp. Esse merge é chamado um para um porque cada código de microrregião aparece apenas uma vez na tabela, e portanto cada polígono é associado a apenas uma linha de dados. Para testar essa funcionalidade, execute o código abaixo.

## merge com a malha
pam_13_soja_sp <- merge(malha.bra.mun13,pam_13_soja_micro,by.x="SP_ID",by.y="codigo_mun7",all.y=T)

   Os comandos abaixo confirmam que trata-se de um dado espacial (classe "sp"). O comando View no slot "data" permite visualizar os atributos dos polígonos, que agora contém as informações da base "pam_13_soja_micro".

class(pam_13_soja_sp)
slotNames(pam_13_soja_sp)
View(pam_13_soja_sp@data)

   Repare que todos os polígonos originais de microrregião estão no slot "data" do nosso novo objeto sp, existindo alguns valores NA nas microrregiões que não produzem soja. É possível filtrar esse objeto para excluir os valores NA através do comando subset(), da mesma forma que se faz com um data frame usual. Nos comandos abaixo, um novo objeto sp é formado a partir desse filtro. Os dois plots finais comparam os dois objetos. Note que o segundo plot, referente ao filtro, não contém todos os polígonos.

## filtrando apenas valores não 'NA' da base
pam_13_soja_sp_com <- subset(pam_13_soja_sp,!is.na(valor_prod_1000))
View(pam_13_soja_sp_com@data)

## comparando os polígonos
plot(pam_13_soja_sp)
plot(pam_13_soja_sp_com)

3. Exportando o dado georreferenciado em um novo shape

   A base de produção de soja pode ser salva em arquivo .shp. O comando writePolyShape() do pacote maptools exporta os dados para esse formato, criando também os complementos em .dbf e .shx.

## Exportando teste
writePolyShape(pam_13_soja_sp_com,"pam_13_soja_sp.shp")

   O arquivo foi salvo na mesma pasta que você selecionou para descompactar a malha de microrregiões. É possível carregar o arquivo .shp no ArcMAP normalmente e checar os atributos dos polígonos.

4 - Georreferenciando dados: relação um para muitos

   O georreferenciamento um para muitos ocorre quando no dado tabular cada ID de polígono aparece várias vezes, como em uma série histórica de produção de soja em uma microrregião, por exemplo. No código abaixo, a série de produção anual de soja por microrregião é criada e um merge e logo após tenta-se um merge com a malha. Após o último comando, aparecerá a mensagem de erro "'y' has multiple records for one or more 'by.y' key(s)", mencionando exatamente que há múltiplos IDs de polígonos na base tabular.

## filtrando PAM para soja
pam_soja <- subset(pam,chave.prod=="soja")

## merge com a DTB
pam_soja <- merge(pam_soja,dtb[c("codigo_mun7","codigo_micro","ano")],by=c("codigo_mun7","ano"))

## agregação por microrregiao e ano
pam_soja_micro <- aggregate(cbind(area_plnt_dest,area_colh,qt_prod,valor_prod_1000) ~ codigo_micro + ano,
                            data=pam_soja,FUN=sum)

## tentano georreferenciamento "um para muitos"
pam_soja_sp <- merge(malha.bra.mun13,pam_soja_micro,by.x="CD_GEOCMI",by.y="codigo_micro")

    Há uma alternativa para essa operação dentro do R. A idéia é gerar "quebras" ("splits") do data frame por ano. Assim, em cada quebra cria-se uma nova sequência de id's das microrregiões na tabela e a mesma sequência é atribuída nos id's dos polígonos através da função spChFIDs(). Em seguida, um merge da quebra com os polígonos é executado e a mesma operação é executada no próximo split, mas de modo que o menor ID da nova quebra é sequência do último ID da quebra anterior.
    No final, cada split é um objeto sp que é concatenado de forma similar ao que foi feito na criação do shape de microrregião usado neste exemplo. O código abaixo executa esse procedimento e resulta no arquivo "pam_soja_sp.shp". Antes é necessário instalar o pacote "plyr"

library(plyr)

geo_um_muitos <- function(data,malha){
 
  # pegando "splits"
  data.split <- split(data,f=data[c("ano")],drop=T)  
 
  # looping nos sp's criados
  ti <- Sys.time()
  i <- 0
  maxid <- 0
  for(split in data.split){
   
    # mostrando iteração
    i <- i + 1
    print(i)
   
   
    # juntando substrato usando merge do 'sp' no 'SPLIT'
    spdf <- merge(malha,split,by.x="SP_ID",by.y="codigo_micro",all.y=T)
   
   
    # spdf: cópia dos IDs originais
    spdf$ID_orig <- spdf$SP_ID
   
    # spdf: criando nomes que depois vão virar novos IDs
    novo_ID <- (maxid + 1):(maxid +NROW(spdf$SP_ID))
    maxid = max(novo_ID)
   
    # spdf: DEFININDO novo_ID COMO NOVOS IDS DE POLÍGONOS
    spdf <- spChFIDs(spdf,as.character(novo_ID))
   
    if(exists("data.sp_teste2")){
      # spdf E data.sp_teste: juntando os dois caras
      data.sp_teste2 <- rbind(data.sp_teste2,spdf)    
    } else{
      # se data.sp_teste2 nao existe, ele será o spdf atual
      data.sp_teste2 <- spdf
    }
  }
  (tf <- difftime(Sys.time(),ti))
  return(data.sp_teste2)
 
}


pam_soja_sp <- geo_um_muitos(pam_soja_micro,malha.bra.mun13)
View(pam_soja_sp@data)

## filtrando apenas valores não 'NA' da base
pam_soja_sp_com <- subset(pam_soja_sp,!is.na(valor_prod_1000))

writePolyShape(pam_soja_sp_com,"pam_soja_sp.shp")

   O procedimento acima, entretanto, é bastante limitado. A começar pela complexidade de comandos, que exige múltiplas quebras e redefinições de ID's. Além disso, esse método não funciona para muitas quebras, já que pode atingir o limite da memória RAM. Pelo mesmo motivo, o procedimento também é limitado para dados municipais. A melhor alternativa para georreferenciar dados no modo um para muitos sem sair do R é considerar comandos em PostGIS usando conexão com banco PostgreSQL via R

Tags:
Created by Wesley Silva on 2015/09/11 17:40
    
This wiki is licensed under a Creative Commons 2.0 license
XWiki Enterprise 6.3 - Documentation