Search code examples
rarcgisqgisr-sf

Reading/importing .tpk maps into R or QGIS and use as shapefile

Is it possible to import a .tpk map file either in R or QGIS and use it as a shapefile? I need the coordinates of specific locations in the .tpk map that could be extracted with a shapefile map. I do not have access to ArcGIS.

Grateful for any guidance!

For additional information about .tpk, please see: https://pro.arcgis.com/en/pro-app/latest/help/sharing/overview/tile-package.htm

Solution

  • Here's what I've done:

    First convert the tpk file to mbtiles using the tpk conversion utility from https://github.com/consbio/tpkutils

    tpk export mbtiles 04010588800801.tpk 04010588800801.mbtiles

    gdalinfo shows a few things about the mbtiles file:

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

    We could load this into R, but that size is enormous and I can't find a way to select a given zoom level and a given range based on the valid tiles. This I can do with gdal_translate on the command line or via the gdalUtils package to create a GeoTIFF for a given zoom level, using USE_BOUNDS=NO to constrain the output to only where tiles exist:

    Command line:

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

    gdalUtils package:

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

    Reading and then plotting this level 17 RGB image can then be done:

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

    enter image description here

    Note this is quite a high resolution image so you can't read the labels, but if you zoom in (or load into QGIS and interact there) you can read the labels. Here's an extreme zoom into level 17 in QGIS showing the limits of the resolution:

    enter image description here

    Remember this is only raster image data so if you want the coordinates of those points you'll have to create a new layer in QGIS and manually create a point data set over the imagery. If this is what you want then I strongly suggest you try and get vector data from the supplier and not have to do all this!

    The other zoom levels might be useful to you so convert them using the gdal_translate procedure above. As you go out you lose detail, and below level 13 you just get the overview maps.

    Level 16:

    enter image description here

    Level 13:

    enter image description here

    Level 12:

    enter image description here

    Update

    A slightly more direct way to read from the .mbtiles file. Use the stars package which does allow you to pass GDAL options:

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

    The proxy=FALSE is needed otherwise when converting with as(..,"Raster") the output reverts back to the 20166662x23418224 dimension of the full global level 17 zoom raster. Possible bug in stars somewhere. Anyway, this gets you the zoomed rasters without having to use gdal_translate anywhere.