1

I'm trying to extract features using python from a vector tile layer's single .pbf tile :

https://indianopenmaps.fly.dev/google-buildings/{z}/{x}/{y}.pbf

This tile layer comes from : https://github.com/ramSeraph/google_buildings_india, and when we add it on QGIS, it loads up properly, screenshot:

QGIS vector tile layer preview

From this vector layer's properties, seeing these CRS details:

Name    EPSG:3857 - WGS 84 / Pseudo-Mercator
Units   meters
Method  Mercator
Celestial body  Earth
Accuracy    Based on World Geodetic System 1984 ensemble (EPSG:6326), which has a limited accuracy of at best 2 meters.
Reference   Dynamic (relies on a datum which is not plate-fixed)

So I'm assuming the source SRID is 3857.

One vector tile file: https://indianopenmaps.fly.dev/google-buildings/12/2934/1684.pbf

I'm using "mapbox-vector-tile" library (homepage: https://github.com/tilezen/mapbox-vector-tile) to decode the binary pbf data into geojson. Functional python code:

import requests, mapbox_vector_tile
from pyproj import Transformer

url = "https://indianopenmaps.fly.dev/google-buildings/12/2934/1684.pbf"
response = requests.get(url)
pbf_data = response.content
# when printed, this shows: b'\x1a\xe2\x8f\x0cx\x02\n\x10google_building....

vector1 = mapbox_vector_tile.decode(tile=pbf_data)
print(vector1['google_buildings']['features'][0])

Output:

{'geometry': {'type': 'Polygon',
   'coordinates': [[[5, 3922], [7, 3918], [1, 3916], [-1, 3919], [5, 3922]]]},
  'properties': {'confidence': 0.7949},
  'id': 0,
  'type': 'Feature'}

This is obviously not in lat-longs. So the Readme at "mapbox-vector-tile" github provides a way to convert to lat-longs, so I applied that:

SRID_LNGLAT = 4326
SRID_SPHERICAL_MERCATOR = 3857
reverse_transformer = Transformer.from_crs(crs_from=SRID_SPHERICAL_MERCATOR, crs_to=SRID_LNGLAT, always_xy=True)

vector2 = mapbox_vector_tile.decode(tile=pbf_data, default_options={"transformer": reverse_transformer.transform})
print(vector2['google_buildings']['features'][0])

Output:

{'geometry': {'type': 'Polygon', 'coordinates': [[[4.491576420597607e-05, 0.035231923222862935], [6.28820698883665e-05, 0.03519599061828461], [8.983152841195214e-06, 0.03517802431599025], [-8.983152841195214e-06, 0.035204973769430485], [4.491576420597607e-05, 0.035231923222862935]]]}, 'properties': {'confidence': 0.7949}, 'id': 0, 'type': 'Feature'}

Unfortunately, even now I'm not getting this in lat-longs.
Those original co-ords ([5, 3922], [7, 3918]..) were looking suspect anyways - didn't seem like 3857 co-ords. I'm guessing they're relative to some place.

So, how do I get these co-ords in latitude-longitude? I know that this tile layer's data is proper, since QGIS is not having any problem in rendering it out-of-the-box.

Related questions:

  1. Convert Mapbox Vector Tiles to Longitude and Latitude : Here the top-voted answer merely refers to the "mapbox-vector-tile" library's homepage, there's no code sample that shows how to convert to lat-longs, and following the transform process given in that homepage has given the above result. So, we're not having a complete answer there.

  2. Extracting specific layers from the PBF file : answers here refer to the same protobuf to json conversion which is already successfully happening, but don't mention anything about transforming the output co-ords to lat-longs. So, we're not having a complete answer there.

1
  • 1
    Read the MVT specification mapbox.github.io/vector-tile-spec. The coordinates in the vector tiles are always in the local tile units. The top-left corner of every tile is (0,0). For connecting MVT coordinates with real world coordinates you must know the tiling schema, and by the schema resolve where the tile 12/2936/1684 is located.
    – user30184
    Commented Feb 3 at 15:32

1 Answer 1

3

Update: Came across this answer : https://gis.stackexchange.com/a/460173/44746 posted at similar question that I hadn't found yet: Decoding Mapbox Vector Tiles?

Updating my code to have a full working solution from tile url to shapefile:

import requests, mapbox_vector_tile, math, json

def pixel2deg(xtile, ytile, zoom, xpixel, ypixel, extent = 4096):
    # from https://gis.stackexchange.com/a/460173/44746
    n = 2.0 ** zoom
    xtile = xtile + (xpixel / extent)
    ytile = ytile + ((extent - ypixel) / extent)
    lon_deg = (xtile / n) * 360.0 - 180.0
    lat_rad = math.atan(math.sinh(math.pi * (1 - 2 * ytile / n)))
    lat_deg = math.degrees(lat_rad)
    return (lon_deg, lat_deg)

url = "https://indianopenmaps.fly.dev/google-buildings/12/2934/1684.pbf"
tile_x = 2934
tile_y = 1684
tile_zoom = 12

response = requests.get(url)
pbf_data = response.content
# when printed, this shows: b'\x1a\xe2\x8f\x0cx\x02\n\x10google_building....

vector3 = mapbox_vector_tile.decode(tile=pbf_data, 
  transformer = lambda x, y: pixel2deg(tile_x, tile_y, tile_zoom, x, y)
)
print(vector3['google_buildings']['features'][0])

with open('shape1.geojson','w') as f:
    json.dump(vector3['google_buildings'],f)

Output:

{'geometry': {'type': 'Polygon',
  'coordinates': [[[77.8712010383606, 30.445454921154706],
    [77.87124395370483, 30.445380925465102],
    [77.87111520767212, 30.445343927599236],
    [77.87107229232788, 30.445399424392768],
    [77.8712010383606, 30.445454921154706]]]},
 'properties': {'confidence': 0.7949},
 'id': 0,
 'type': 'Feature'}

Now this looks like proper lat-longs!

And the .geojson file when dropped over https://geojson.io/ : enter image description here

The shoe fits!

Props to Mark Egge for sharing the pixel2deg function, it was the missing link.

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.