mapas_datasus

# http://r-br.2285057.n4.nabble.com/R-br-Ler-MAP-no-R-Shapefile-td4661880.html
# http://www2.datasus.gov.br/DATASUS/tabwin/DocTabWin.htm
# http://datasus.saude.gov.br/versao-3
# http://www2.datasus.gov.br/DATASUS/tabwin/rx/Importa_mapa.htm
# http://datasus.saude.gov.br/cadastros-nacionais/294-download-mapas-tabwin


read.map = function(filename){
  zz=file(filename,"rb")
  #
  # header do .map
  #
  versao = readBin(zz,"integer",1,size=2)  # 100 = versao 1.00
  #Bounding Box
  Leste = readBin(zz,"numeric",1,size=4)
  Norte = readBin(zz,"numeric",1,size=4)
  Oeste = readBin(zz,"numeric",1,size=4)
  Sul   = readBin(zz,"numeric",1,size=4)
  
  geocodigo = ""
  nome = ""
  xleg = 0
  yleg = 0
  sede = FALSE
  poli = list()
  i = 0
  
  #
  # repete para cada objeto no arquivo
  #
  repeat{  
    tipoobj = readBin(zz,"integer",1,size=1) # 0=Poligono, 1=PoligonoComSede, 2=Linha, 3=Ponto
    
    if (length(tipoobj) == 0) break
    i = i + 1
    
    Len = readBin(zz,"integer",1,size=1)  # length byte da string Pascal
    geocodigo[i] = readChar(zz,10)
    Len = readBin(zz,"integer",1,size=1)  # length byte da string Pascal
    nome[i] = substr(readChar(zz,25),1,Len)
    xleg[i] = readBin(zz,"numeric",1,size=4)
    yleg[i] = readBin(zz,"numeric",1,size=4)
    numpontos = readBin(zz,"integer",1,size=2)
    
    sede = sede || (tipoobj = 1)
    
    x=0
    y=0   
    for (j in 1:numpontos){
      x[j] = readBin(zz,"numeric",1,size=4)
      y[j] = readBin(zz,"numeric",1,size=4)
    }
    
    
    # NAs separam v?rios pol?gonos no mesmo objeto
    # BUG a corrigir: Assim como est? o primeiro pol?gono n?o fecha e, em multiplos poligonos, h? um NA a mais no final 
    xInic = x[1]
    yInic = y[1]  
    for (j in 2:numpontos){
      if (x[j] == xInic & y[j] == yInic) {x[j]=NA; y[j] = NA}
    }
    
    poli[[i]] = c(x,y)
    dim(poli[[i]]) = c(numpontos,2)
  }
  
  class(poli) = "polylist"
  attr(poli,"region.id") = geocodigo
  attr(poli,"region.name") = nome
  attr(poli,"centroid") = list(x=xleg,y=yleg)
  attr(poli,"sede") = sede
  attr(poli,"maplim") = list(x=c(Oeste,Leste),y=c(Sul,Norte))
  
  close(zz)
  return(poli)
}


datasus_as_st <- function(map_file){
  
  map_datasus <- read.map(map_file)
  
  lst_polygon <- purrr::map(map_datasus, ~sp::Polygon(na.omit(.)))
  # function to convert Polygon-class to Polygons-class
  foo <- function(x, y) sp::Polygons(list(x), ID = y)
  
  id_cidades <- paste0(attributes(map_datasus)$region.id,
                       ";",
                       attributes(map_datasus)$region.name)
  
  x <- purrr::map2(lst_polygon,
                   id_cidades,
                   foo)
  
  geo <- sf::st_as_sf(sp::SpatialPolygons(x))
  geo[["cidade"]] <- id_cidades
  
  
  geo
}

dir("data/mapbr/")

mapa_br_municipios <- datasus_as_st("data/mapbr/BR.MAP")