将 .tpk 地图读取/导入到 R 或 QGIS，并将其用作 shapefile。

Question

有没有可能在R或QGIS中导入.tpk地图文件，并将其用作shapefile？我需要从.tpk地图中提取特定位置的坐标，这些坐标可通过shapefile地图提取。我无法访问ArcGIS。

非常感谢任何指导！

有关.tpk的其他信息，请参见：https://pro.arcgis.com/en/pro-app/latest/help/sharing/overview/tile-package.htm

- edo

---

有可以下载和尝试的吗？我说“R现在还不能做到这一点”，因为GDAL不理解来自此处的文件：https://github.com/consbio/tpkutils/tree/main/tests/data - 你可以尝试使用那个Python代码将其转换为GDAL所理解的格式... - Spacedman

还要注意该存储库中关于ESRI更改格式（10.1 vs 10.3）的评论。他们总是在变动目标。 - Spacedman

感谢提供 Python 包信息。这里是一个 .tpk 文件。如果 R 能够将其识别为 shapefile，那就太好了，这绝对超出了我的 R 技能范围：https://drive.google.com/file/d/10jfFw4ovfdfnw3qnLVsJjMSZQ8HvjTDG/view?usp=sharing - edo

这些是栅格数据而不是矢量数据，因此您不能将其作为“shapefile”阅读，如果您的意思是作为线条和多边形。它是由不同分辨率的网格图像组成的地图。Python代码可以在任何给定分辨率下创建地图，但是R将所有分辨率解释为全局地图，因此在缩放17以查看您的研究区域时，光栅是数百万x数百万。这可能可以通过范围进行裁剪，但仍然不是“shapefile”矢量数据。 - Spacedman

感谢您的见解和时间。您认为使用您提到的Python包，我能否获得这些地图的质心，或更一般地说，这些地图中特定点的坐标？您介意分享我可以尝试自己的潜在方法吗？ - edo

Answer 1

以下是我完成的步骤：

首先使用tpk转换工具将tpk文件转换为mbtiles格式，该工具可从https://github.com/consbio/tpkutils获取。

tpk export mbtiles 04010588800801.tpk 04010588800801.mbtiles

gdalinfo可以显示有关mbtiles文件的一些信息：

$ gdalinfo 04010588800801.mbtiles 
Driver: MBTiles/MBTiles
Files: 04010588800801.mbtiles
Size is 23418224, 20166662
...
Metadata:
  ZOOM_LEVEL=17
...
  minzoom=0
  maxzoom=17
...

我们可以将其加载到R中，但是这个文件太大了，我找不到一种方法来选择给定的缩放级别和基于有效瓷砖的给定范围。使用gdal_translate命令行或通过gdalUtils软件包创建一个给定缩放级别的GeoTIFF，使用USE_BOUNDS=NO将输出限制在只有瓷砖存在的区域内：

命令行:

gdal_translate -oo ZOOM_LEVEL=17 -oo USE_BOUNDS=NO 04010588800801.mbtiles zoom17.tiff

gdalUtils 包：

gdal_translate("04010588800801.mbtiles","z17.tiff",
     oo=c("USE_BOUNDS=NO","ZOOM_LEVEL=17"))

阅读并绘制这个17级RGB图像可以这样完成：

> z17 = raster::stack("z17.tiff")
[ignore CRS warnings...]
> plotRGB(z17)

请注意，这是一张高分辨率的图片，因此您无法阅读标签，但如果您放大（或加载到QGIS中并在那里进行交互），则可以阅读标签。以下是在QGIS中对17级别进行极端缩放，显示分辨率限制的示例：

请记住，这只是栅格图像数据，因此如果您想要这些点的坐标，您需要在QGIS中创建一个新图层，并在图像上手动创建一个点数据集。如果这是您想要的，则强烈建议您尝试从供应商获取矢量数据，而不必做所有这些工作！

其他缩放级别可能对您有用，因此请使用上面的gdal_translate过程将它们转换。随着缩小，您会失去细节，在13级以下，您只能获得概览地图。

16级：

等级13：

等级12：

更新

从.mbtiles文件中读取的一种更直接的方法。使用stars包，它允许您传递GDAL选项：

s = stars::read_stars("04010588800801.mbtiles",
    options=c("USE_BOUNDS=NO","ZOOM_LEVEL=17"), proxy=FALSE)
rs = as(s,"Raster")
raster::plotRGB(rs)

需要使用proxy=FALSE，否则在使用as(..,"Raster")进行转换时，输出会恢复到全球级别17缩放栅格的20166662x23418224维度。可能是stars中的错误。无论如何，这样可以让你获得缩放后的栅格图像，而不必在任何地方使用gdal_translate。

- Spacedman

---

如何使用 library(terra); r = rast("04010588800801.mbtiles", opts=c("USE_BOUNDS=NO","ZOOM_LEVEL=17")) ？ - Robert Hijmans

@RobertHijmans 在使用 opts.1.3-22 时遇到了 "未使用的参数" 错误......尝试从 CRAN 获取了最新版本 (1.4.22)，虽然可以运行，但似乎每次无法获取正确的缩放级别 - 级别 15、16、17 看起来都像是级别 15。这些瓷砖的分辨率比预期的要高，但似乎被重新缩放为缩放等级 15 的瓷砖。这是否是一个 bug 报告？ - Spacedman

你能分享一下 mbtiles 文件吗？ - Robert Hijmans

请下载此链接中的文件：https://www.dropbox.com/s/weaiqu7fkia3fbp/04010588800801.mbtiles?dl=0，我们可以在其他地方继续讨论。 - Spacedman

Answer 2

在@Spacedman的回答基础上，使用他创建的tpk文件，我使用terra包展示了一系列缩放级别。

library(terra)
plotm <- function(zoom, e=NULL) {
    # read file with GDAL options
    r = rast("04010588800801.mbtiles", opts=c("USE_BOUNDS=NO", paste0("ZOOM_LEVEL=",zoom))); 
    # declare the RGB channels
    RGB(r) <- 1:4
    plot(r, ext=e)
}
par(mfrow=c(3,3))
x <- lapply(6:14, plotm)

这里使用了一个扩展来稍微放大一点，以便更好地看到最高分辨率下的差异。

par(mfrow=c(2,2))
e <- ext(3935206, 3935662, 1085835, 1086294)
x <- lapply(14:17, plotm, e=e)

原始的tpk文件在ArcMap中看起来与此相同。

- Robert Hijmans