0

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.

5
  • 1
    you can just request the tiles that you need directly in a web browser or with curl
    – Ian Turton
    Commented Apr 25, 2023 at 7:25
  • Do you want to analyze the tiles in .pbf format or as converted into some GIS format?
    – user30184
    Commented Apr 25, 2023 at 8:56
  • I want to convert them into a GIS format. So they are still georeferenced. This is why just downloading the tiles using a browser doesn't work as they lose the location.
    – RobShaw_UK
    Commented Apr 25, 2023 at 9:13
  • Vector tiles are not georeferenced. Or, they are georeferenced only indirectly by the tile number- if tile number is this and tiling schema is that then the software knows where to render the data of the tile. MVT layer is not a real vector layer for QGIS qgis.org/pyqgis/3.22/core/QgsVectorTileLayer.html and that's why export option in QGIS is greyed out.
    – user30184
    Commented Apr 25, 2023 at 9:49
  • Thanks. Do you know if there is a way to identify how QGIS generates the x and y values to pass? The layer CRS is ESPG 3857.
    – RobShaw_UK
    Commented Apr 25, 2023 at 10:43

1 Answer 1

0

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)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.