r - 如何在 R 中读取 .MAP 文件扩展名？

Question

有没有一种简单的方法可以在 R 中读取扩展名为 .MAP 的文件？我尝试了以下几个选项但没有成功。 Here is a .MAP file for a reproducible example .

背景:出于某种奇怪的原因，巴西卫生规划政策中使用的空间区域化仅适用于这种格式。我想将它转换为 geopackage，这样我们就可以将它添加到 geobr 包中。

# none of these options work
mp <- sf::st_read("./se_mapas_2013/se_regsaud.MAP")
mp <- rgdal::readGDAL("./se_mapas_2013/se_regsaud.MAP")
mp <- rgdal::readOGR("./se_mapas_2013/se_regsaud.MAP")
mp <- raster::raster("./se_mapas_2013/se_regsaud.MAP")
mp <- stars::read_stars("./se_mapas_2013/se_regsaud.MAP")

ps。 there is a similar question on SO focused on Python, unfortunately unanswered

更新

我们发现 a publication 使用自定义函数来读取 .MAP 文件。请参见下面的示例。但是，它返回一个 "polylist" 对象。有没有一种简单的方法可以将其转换为简单的功能？

原创自定义函数

read.map = function(filename){
  zz=file(filename,"rb")
  #
  # header of .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

  #
  # repeat of each object in file
  #
  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)
    }


    # separate polygons
    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)
}

使用原始自定义函数

mp <- read.map("./se_mapas_2013/se_regsaud.MAP")

class(mp)
>[1] "polylist"

# plot
plot(attributes(mp)$maplim, type='n', asp=1, xlab=NA, ylab=NA)
title('Map')
lapply(mp, polygon, asp=T, col=3)

Answer 1

问题是:使用带有尾随 nul 字节的 readChar - 更改为 readBin()； rawToChar() 不接受的 8 位字符(在我的 UTF-8 系统上)；一些需要删除的文件中有多个条子；和其他一些。我将上面已编辑的 read.map() 函数添加到 ma​​ptools，但名称不同且未导出。所以现在(使用 ma​​ptools rev 370 from https://r-forge.r-project.org/R/?group_id=943 当构建完成时):

library(maptools)
o <- maptools:::readMAP2polylist("se_regsaud.MAP")
oo <- maptools:::.makePolylistValid(o)
ooo <- maptools:::.polylist2SpP(oo, tol=.Machine$double.eps^(1/4))
rn <- row.names(ooo)
df <- data.frame(ID=rn, row.names=rn, stringsAsFactors=FALSE)
res <- SpatialPolygonsDataFrame(ooo, data=df)
library(sf)
res_sf <- st_as_sf(res)
res_sf
plot(st_geometry(res_sf))

这种方法重用了近 20 年前的 ma​​ptools 代码，稍作修改以处理读取二进制文件和修复碎片的后续更改。