Exporting vector tiles or saving them locally using QGIS

Question

I'm using QGIS 3.30 and have added a vector tile layer in the format https://privateurl.com/tiles/{z}/{x}/{y}.pbf. I would like to be able to export this at zoom level 14 to local files so I can do additional analysis on the polygons. Does anyone know if this is possible?

The usual right-mouse-click > export > save as... is greyed out.

Answer 1

vector tiles can't be downloaded directly in QGIS, but with some python programming they can be requested tile by tile and then georeferenced. All vector tile server work differently and have a different scope so you will definitely have to adjust the code below.

import pandas as pd
import requests
import mapbox_vector_tile
from shapely.geometry import shape, Point
import geopandas as gpd
import math
from shapely import geometry

def fetch_tile(url, z, x, y):
    tile_url = url.format(z=z, x=x, y=y)
    response = requests.get(tile_url)
    if response.status_code == 200:
        return response.content
    else:
        #print("Error fetching tile:", response.status_code)
        return None

def coordinate_traslation(xtile, ytile, xcoord, ycoord, zoom):
    #adjustments to server
    xcoord = xcoord * 32
    ycoord =  ycoord * -32
    ytile = ytile +1

    x, y = xtile + xcoord / (2 ** zoom), ytile + ycoord/ (2 ** zoom)

    n = 2.0 ** zoom
    lon_deg = x / n * 360.0 - 180.0
    lat_rad = math.atan(math.sinh(math.pi * (1 - 2 * y / n)))
    lat_deg = math.degrees(lat_rad)
    return lon_deg, lat_deg

def to_gdf(tile_data):
    tile = mapbox_vector_tile.decode(tile_data)
    geometries = []
    attributes = []

    for layer_name, layer in tile.items():
        # Iterate through features
        for feature in layer['features']:
            # Extract geometry and attributes
            geometries.append(shape(feature['geometry']))
            attributes.append(feature['properties'])

    # Create GeoDataFrame
    gdf = gpd.GeoDataFrame(attributes, geometry=geometries, crs='EPSG:4326')

    # multipolygons to single
    gdf = gdf.explode(index_parts=True)
    return gdf


def translate(gdf, x, y, z):  # translate to EPSG 4326
    new_geometry = []
    for feture in gdf['geometry']:
        coordinates = []
        for coord in list(feture.exterior.coords):
            coordinates.append(list(coordinate_traslation(x, y, coord[0], coord[1], z)))
        new_geometry.append(geometry.Polygon(coordinates))

    # Create GeoDataFrame
    gdf['geometry'] = new_geometry
    return gdf

#main
# URL of the vector tile server
tile_url =  "https://bodenkarte.at/data/bodenkarte-tiles/{z}/{x}/{y}.pbf"

# Example usage
z = 14  # zoom level
# tile coordinates
minX = 8624
maxX = 8974
minY = 5624
maxY = 5

totalNumberOfTiles = (maxX - minX + 1) * (maxY - minY + 1)
numberOfTilesProcessed = 0
errors = 0
print("That's %d tiles total." % (totalNumberOfTiles))

gdf_list =[]

# Loop through all tiles
for x in range(minX, maxX + 1):
    for y in range(minY, maxY + 1):
        try:
            #fetch data from server
            tile_data = fetch_tile(tile_url, z, x, y)

            #convert to gdf
            gdf = to_gdf(tile_data)

            #transalte coordinates to epsg 4326
            gdf = translate(gdf, x, y, z)

            gdf_list.append(gdf)
            numberOfTilesProcessed += 1
            if numberOfTilesProcessed % 100 == 0:
                print(numberOfTilesProcessed)
        except:
            errors += 1
            if errors%100 == 0:
                print(errors, 'e')
            pass

print(errors, 'errors')
print('all ', totalNumberOfTiles, 'are processed' )

gdf = pd.concat(gdf_list)

#save as gpkg
gdf.to_file('/path/to(file.gpkg', driver='GPKG', layer='name')
print('saved')

the tile coordinates can be easily calculate with following script:

import math

def lat_lon_to_tile(lat, lon, zoom):
    n = 2.0 ** zoom
    x_tile = int((lon + 180.0) / 360.0 * n)
    lat_rad = math.radians(lat)
    y_tile = int((1.0 - math.log(math.tan(lat_rad) + (1 / math.cos(lat_rad))) / math.pi) / 2.0 * n)
    return x_tile, y_tile

# Example usage
latitude = 49.0208# Example latitude (e.g., New York City)
longitude = 17.170574# Example longitude (e.g., New York City)
zoom_level = 14  # Example zoom level
x, y = lat_lon_to_tile(latitude, longitude, zoom_level)
print("Tile Coordinates (x, y):", x, y)