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

Posted 2023-02-21

【中文标题】读取/导入 .tpk 地图到 R 或 QGIS 并用作 shapefile

【英文标题】：Reading/importing .tpk maps into R or QGIS and use as shapefile

【发布时间】：2022-01-02 04:02:54

【问题描述】：

是否可以在 R 或 中导入 .tpk 地图文件>QGIS 并将其用作 shapefile？我需要 .tpk 地图 中特定位置的坐标，可以使用 shapefile 地图提取。我无权访问 ArcGIS。

感谢任何指导！

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

【问题讨论】：

有一个可以下载试用的吗？我说“R 不能这样做（还）”，因为 GDAL 不理解这里的文件：github.com/consbio/tpkutils/tree/main/tests/data - 你可以尝试使用该 python 代码转换为 GDAL 确实理解的东西......

另请注意该存储库中有关 ESRI 更改格式（10.1 与 10.3）的评论。他们总是在移动球门柱。

感谢 python 包信息。这里有一个 .tpk 文件。如果 R 可以将其读取为 shapefile 那就太好了，这绝对是一个超出我 R 技能的项目：drive.google.com/file/d/10jfFw4ovfdfnw3qnLVsJjMSZQ8HvjTDG/…

这些是光栅数据而不是矢量数据，因此如果您的意思是线和多边形，您将无法将其“作为 shapefile”读取。它的网格图像，具有不同分辨率的不同地图。 python 代码可以以任何给定的分辨率创建地图，但 R 将所有分辨率解释为全局地图，因此在缩放 17 以查看您的研究区域时，栅格为数百万 x 数百万。考虑到边界，这可能是可裁剪的，但它仍然不是“shapefile”矢量数据。

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

【参考方案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 维度。某处星星中可能存在的错误。无论如何，这让您无需在任何地方使用 gdal_translate 即可获得缩放的栅格。

【讨论】：

library(terra); r = rast("04010588800801.mbtiles", opts=c("USE_BOUNDS=NO","ZOOM_LEVEL=17")) 怎么样？

@RobertHijmans “未使用的参数”错误与 opts。 1.3-22 ...从 CRAN (1.4.22) 获得了最新的并且可以运行，但似乎每次都没有获得正确的缩放级别 - 级别 15、16、17 看起来它们都是 15 级。他们具有预期的更高分辨率，但似乎重新缩放缩放 = 15 块。错误报告？

可以分享mbtiles文件吗？

dropbox.com/s/weaiqu7fkia3fbp/04010588800801.mbtiles?dl=0 我们去别处讨论？

【参考方案2】：

基于@Spacedmans 的回答并使用他创建的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)

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