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)

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:

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:

Level 13:

Level 12:

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.

Building on @Spacedmans's answer and using the tpk file that he created, I show a range of zoom levels using the terra package

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)

Here using an extent to zoom in a bit to be able to better see the differences at the highest resolutions

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

And that looks the same in ArcMap with the original tpk file.