I read the references you provided and also this article: FAQ: What is a root tile and how are they used to make a vector tile package with a local coordinate system? https://support.esri.com/en/technical-article/000022396
Using all the info, I tried to decode the a vector tile myself, using World Topographic Map (https://basemaps.arcgis.com/arcgis/rest/services/World_Basemap_v2/VectorTileServer). I downloaded this particular tile in pbf, https://basemaps.arcgis.com/arcgis/rest/services/World_Basemap_v2/VectorTileServer/tile/6/24/17, and start decoding in R language (sorry it's the only language I know).
However, as it turned out, for some unknown reason, the origin of the selected tile seems to be the bottom left corner (x axis positive to the right, y axis positive upward), instead of top left corner. It is shown in the following codes:
```{r}
library(sf)
librart(ggplot2)
#the st_read() function automatically uses `MVT' driver, and convert the geometry data into tile coordinates
test <- st_read("17.pbf", layer = "City small scale")
#function that convert tile coordinates into the real mapping coordinates (EPSG: 3857)
shift_tile_coord <- function(sf, level, row, column, crs = 3857) {
if(crs != 3857) stop("currently only supports EPSG:3857 Web Mercator projection")
origin_x <- -20037508.342787
origin_y <- 20037508.342787
maxscale <- 2.958287637957775E8
tile_size <- 512 #pixels of a row/ column
dpi <- 96
inch2mapUnit <- 39.3700787 #1 meter = 39.3700787 inch
tile_extent <- 4096
root_tile_width <- maxscale * (tile_size/(dpi * inch2mapUnit))
selected_tile_width <- root_tile_width/2^level
sf_shifted <- sf
#x = tile_x_coor*selected_tile_width/tile_extent + origin_x + selected_tile_width*column
#y = tile_y_coor*selected_tile_width/tile_extent + origin_y - selected_tile_width*(row + 1)
sf_shifted$geometry <- sf$geometry * matrix(c(selected_tile_width/tile_extent, selected_tile_width/tile_extent), ncol = 2)
sf_shifted$geometry <- sf_shifted$geometry + matrix(c(origin_x + selected_tile_width * column, origin_y - selected_tile_width * (row + 1)), ncol = 2)
sf_shifted <- st_set_crs(sf_shifted, crs)
return(sf_shifted)
}
#visulaize the data
test_shifted <- shift_tile_coord(test, level = 6, row = 24, column = 17)
ggplot(data = test_shifted) +
geom_sf_label(aes(label = X_name_en)) +
geom_sf() +
coord_sf(datum = st_crs(3857))
```
Test Data
If you use map viewer to check (https://www.arcgis.com/apps/mapviewer/index.html?layers=7dc6cea0b1764a1f9af2e679f642f0f5#), the coordinates (in Web Mercator, EPSG: 3857) are correct. However, the conversion of a y coordinate = tile_y_coor*selected_tile_width/tile_extent + origin_y - selected_tile_width*(row + 1). This means a y coordinate first starts at the bottom of the tile (thus (row + 1)), and moves upward (positive sign in tile_y_coor*selected_tile_width/tile_extent). This seems to be against the documentation you provided. Is it designed to be this way, or is there something wrong?
@mgeorge