QGIS expression function for intersecting points with underlying MapBox vector tile layer to capture attributes

Question

How can I use the QGIS 3.34.0 field calculator to intersect points with an underlying MapBox vector layer and add certain attributes of the intersected features as virtual point layer fields?

If this is actually not possible, how can I do this using a Processing Tool, etc. (i.e. GDAL)?

Here's the vector tile layer I'm using:

URL: https://bodenkarte.at/data/bodenkarte-tiles/{z}/{x}/{y}.pbf
Min Zoom: 0
Max Zoom: 21
Style URL: https://bodenkarte.at/styles/bodentyp.json

• christoph: @Taras I just wanted to make sure I hadn't overlooked anything. Actually, I just want to replicate the QGIS internal "Identity Features" function that returns all attributes. It looks like we have to get the vector tile first and then do the intersection with the point.
• christoph: Even more interesting is how we can get the entire feature geometry back. This is something that even QGIS with the "Identify Features" function is not yet able to do. But Openlayers does this perfectly: openlayers.org/en/latest/examples/vector-tile-selection.html
• Babel: @Taras I'm not sure why you say it's not possible - what about overlay functions? With this it's easy to get an attribute from another layer. Do I overlook some caveat or did I missunderstand the question?
• christoph: @Babel it looks like overlay functions don't work with vector tile layers: overlay_intersects( layer:='Bodenkarte',expression:='kurzbezeichnung') results in Eval Error: Layer 'Bodenkarte' could not be loaded.
• Babel: Oh sorry, seems I indde overlooked the tiles part... So yes, Taras' solution is the way to go!

Answer 1

As far as I know, there is no such function in the Field Calculator that can intersect point features with tiles and provide data from the vector tile layer.

I can suggest the following approach to get a set of attributes and their corresponding values as well as geometry (by means of the attributeMap() and geometry() of the QgsFeature class) for a specific location defined by latitude and longitude.

Here I have modified a great solution provided by @BenW in this thread: Selecting feature by coordinates from the Vector Tile Layer using PyQGIS. It utilizes the selectByGeometry() method of the QgsVectorTileLayer class.

Proceed with Plugins > Python Console > Show Editor and paste the script below:

# imports
from qgis.utils import iface
from qgis.core import (QgsVectorTileLayer, QgsGeometry, QgsPoint, QgsPointXY,
 QgsCoordinateReferenceSystem, QgsCoordinateTransform,
 QgsProject, QgsSelectionContext)
# vector tile layer settings
name = "bodenkarte"
url = "https://bodenkarte.at/data/bodenkarte-tiles/{z}/{x}/{y}.pbf"
style_url = "https://bodenkarte.at/styles/bodentyp.json"
min_zoom = 0
max_zoom = 21
uri = f"styleUrl={style_url}&type=xyz&url={url}&zmax={min_zoom}&zmax={max_zoom}"
# loading vector tile layer with corresponding styles
vector_tile_layer = QgsVectorTileLayer(path=uri, baseName=name)
vector_tile_layer.loadDefaultStyle()
# referring to the current QGIS project
project = QgsProject.instance()
# adding vector tile layer to the project
if not len(project.mapLayersByName(name)) > 0 and vector_tile_layer.isValid():
 project.addMapLayer(vector_tile_layer)
# a random location in Vienna
lat, lon = 48.240855, 16.555004
# creating QgsGeometry from that point
vienna = QgsGeometry(QgsPoint(lon, lat))
# reprojecting the point to vector tile layer crs
origin_crs = QgsCoordinateReferenceSystem('EPSG:4326')
target_crs = vector_tile_layer.crs()
trans = QgsCoordinateTransform(origin_crs, target_crs, project)
vienna.transform(trans)
# creating a tiny rectangular around that point
vienna_bb = QgsGeometry().fromRect(vienna.buffer(0.00001, 1).boundingBox())
# setting up proper canvas center and zoom
canvas = iface.mapCanvas()
# coords from the Austrian bbox center
austria_bbox_center = QgsPoint(13.334746, 47.693278)
austria_bbox_center.transform(trans)
canvas.setCenter(QgsPointXY(austria_bbox_center))
# number from zoom to the Austria bbox extent
canvas.zoomScale(3149952)
# setting up selection context
context = QgsSelectionContext()
context.setScale(3149952)
# performing selection in the vector tile layer
vector_tile_layer.selectByGeometry(vienna_bb, context)
# accessing attributes in the selection result
if vector_tile_layer.selectedFeatureCount() > 0:
 for feat in vector_tile_layer.selectedFeatures():
 print(feat.attributeMap())

Press Run script run script and get the output that will look like this:

76055280879434
{'ackerwert': 5, 'ausgangsmaterial': 5000, 'bodenart': 'lZ', 'bodenreaktion': 2, 'bodentyp': 'TS', 'bofo_id': 7019, 'durchlaessigkeit': 5, 'gruendigkeit': 5, 'gruenlandwert': None, 'humusart': 1, 'humusbilanz': 6, 'humuswert': 3, 'kalkgehalt': 9, 'kurzbezeichnung': 'TS', 'nitratrueckhalt': 4, 'nutzbare_feldkapazitaet': 4, 'potentielle_verdichtungsempfindlichkeit': 3, 'tile_layer': 'bodenform_mpoly', 'tile_zoom': 7, 'typengruppe': 'S', 'wasserverhaeltnisse': 5}
MultiPolygon (((1856502.54299035854637623 6146760.06658071931451559, 1856502.54299035854637623 6144314.08167559374123812, 1854056.55808523274026811 6144314.08167559374123812, 1854056.55808523274026811 6139422.11186534259468317, 1854056.55808523274026811 6136976.1269602170214057, 1851610.57318010716699064 6136976.1269602170214057, 1851610.57318010716699064 6139422.11186534259468317, 1851610.57318010716699064 6144314.08167559374123812, 1851610.57318010716699064 6146760.06658071931451559, 1849164.58827498136088252 6146760.06658071931451559, 1849164.58827498136088252 6149206.05148584488779306, 1851610.57318010716699064 6149206.05148584488779306, 1851610.57318010716699064 6151652.03639097046107054, 1854056.55808523274026811 6151652.03639097046107054, 1854056.55808523274026811 6149206.05148584488779306, 1856502.54299035854637623 6149206.05148584488779306, 1856502.54299035854637623 6146760.06658071931451559)),((1863840.49770573526620865 6139422.11186534259468317, 1866286.48261086083948612 6139422.11186534259468317, 1866286.48261086083948612 6136976.1269602170214057, 1871178.45242111221887171 6136976.1269602170214057, 1871178.45242111221887171 6134530.14205509144812822, 1871178.45242111221887171 6132084.15714996494352818, 1868732.46751598664559424 6132084.15714996494352818, 1868732.46751598664559424 6134530.14205509144812822, 1863840.49770573526620865 6134530.14205509144812822, 1863840.49770573526620865 6136976.1269602170214057, 1863840.49770573526620865 6139422.11186534259468317)),((1861394.51280060969293118 6144314.08167559374123812, 1863840.49770573526620865 6144314.08167559374123812, 1866286.48261086083948612 6144314.08167559374123812, 1866286.48261086083948612 6141868.09677046816796064, 1863840.49770573526620865 6141868.09677046816796064, 1863840.49770573526620865 6139422.11186534259468317, 1861394.51280060969293118 6139422.11186534259468317, 1861394.51280060969293118 6141868.09677046816796064, 1861394.51280060969293118 6144314.08167559374123812)),((1854056.55808523274026811 6136976.1269602170214057, 1856502.54299035854637623 6136976.1269602170214057, 1856502.54299035854637623 6134530.14205509144812822, 1856502.54299035854637623 6132084.15714996494352818, 1854056.55808523274026811 6132084.15714996494352818, 1854056.55808523274026811 6134530.14205509144812822, 1854056.55808523274026811 6136976.1269602170214057)),((1861394.51280060969293118 6139422.11186534259468317, 1861394.51280060969293118 6136976.1269602170214057, 1858948.5278954841196537 6136976.1269602170214057, 1856502.54299035854637623 6136976.1269602170214057, 1856502.54299035854637623 6139422.11186534259468317, 1858948.5278954841196537 6139422.11186534259468317, 1861394.51280060969293118 6139422.11186534259468317)),((1849164.58827498136088252 6141868.09677046816796064, 1849164.58827498136088252 6139422.11186534259468317, 1846718.60336985578760505 6139422.11186534259468317, 1846718.60336985578760505 6141868.09677046816796064, 1849164.58827498136088252 6141868.09677046816796064)),((1844272.61846473021432757 6146760.06658071931451559, 1841826.6335596046410501 6146760.06658071931451559, 1841826.6335596046410501 6149206.05148584488779306, 1844272.61846473021432757 6149206.05148584488779306, 1844272.61846473021432757 6146760.06658071931451559)),((1846718.60336985578760505 6154098.02129609603434801, 1846718.60336985578760505 6151652.03639097046107054, 1844272.61846473021432757 6151652.03639097046107054, 1844272.61846473021432757 6154098.02129609603434801, 1846718.60336985578760505 6154098.02129609603434801)),((1856502.54299035854637623 6144314.08167559374123812, 1858948.5278954841196537 6144314.08167559374123812, 1858948.5278954841196537 6141868.09677046816796064, 1856502.54299035854637623 6141868.09677046816796064, 1856502.54299035854637623 6144314.08167559374123812)),((1861394.51280060969293118 6144314.08167559374123812, 1858948.5278954841196537 6144314.08167559374123812, 1858948.5278954841196537 6146760.06658071931451559, 1861394.51280060969293118 6146760.06658071931451559, 1861394.51280060969293118 6144314.08167559374123812)),((1844272.61846473021432757 6146760.06658071931451559, 1846718.60336985578760505 6146760.06658071931451559, 1846718.60336985578760505 6144314.08167559374123812, 1844272.61846473021432757 6144314.08167559374123812, 1844272.61846473021432757 6146760.06658071931451559)))

Why reprojection of the point feature into EPSG:3857 is needed? See, @Mike's brilliant answer to this question Displaying Vector Tile Layer in different Coordinate System with other layers using OpenLayers:

Vector tiles can only be used in their native projection, but you can reproject other layers on the map to that projection. Also vector tiles in projections other than EPSG:3857 will use a custom tile grid which must be defined.

And it is true when you check the metadata with vector_tile_layer.htmlMetadata():

<h1>Coordinate Reference System (CRS)</h1>
<hr>
<table class="list-view">
<tr><td class="highlight">Name</td><td>EPSG:3857 - WGS 84 / Pseudo-Mercator</td></tr>
<tr><td class="highlight">Units</td><td>meters</td></tr>
<tr><td class="highlight">Method</td><td>Mercator</td></tr>
<tr><td class="highlight">Celestial body</td><td>Earth</td></tr>
<tr><td class="highlight">Accuracy</td><td>Based on <i>World Geodetic System 1984 ensemble</i> (EPSG:6326), which has a limited accuracy of <b>at best 2 meters</b>.</td></tr>
<tr><td class="highlight">Reference</td><td>Dynamic (relies on a datum which is not plate-fixed)</td></tr>
</table>

The above script can be also transformed into a custom function via the Function Editor.

---

References:

• Adding vector tile layer (OGC API - Tiles) does not add styling from styleUrl with PyQGIS
• TypeError: QgsGeometry.fromMultiPolygonXY(): argument 1 has unexpected type 'list'
• Getting dictionary for each feature with their attribute values using PyQGIS
• Saving shapes from vector tiles to a shapefile
• Changing the scale of the screen display in QGIS using PyQGIS
• Getting layer extent in PyQGIS?
• Zooming to selected feature using PyQGIS

• christoph: many thanks for your awesome reply! I'm going to put everything into a Python expression function. Unfortunately, it's more work to retrieve the uncut geometry. Therefore I opened a new QGIS issue: github.com/qgis/QGIS/issues/55364

Answer 2

With the helpful support of @Taras, I was able to put together the following Python expression function, which returns all attributes as JSON string. We can fetch individual attributes with the expression function 'map_get()'.

from qgis.utils import iface
from qgis.core import (QgsProject, QgsGeometry, QgsCoordinateReferenceSystem,
 QgsCoordinateTransform, QgsSelectionContext)
@qgsfunction(args='auto', group='Custom', usesGeometry=False, referencedColumns=[])
def getMVTFeatureAttributes(geom,epsg,mvt_layer_name,sscale):
 """
 <style>
 span { color: red }
 </style>
 <h2>returns attributes of a vector tile (MVT) feature intersected by a point geometry</h2>
 getMVTFeatureAttributes(<span>geom</span>,<span>epsg</span>,<span>mvt_layer</span>,<span>sscale</span>)<br/>
 
 <table>
 <tr><td>geom</td><td>... point geometry</td></tr>
 <tr><td>epsg</td><td>... EPSG Code of the point geometry (i.e. 'EPSG:31255')</td></tr>
 <tr><td>mvt_layer</td><td>... Vector Tile Layer Name</td></tr>
 <tr><td>sscale</td><td>... Scale of Selection Context</td></tr>
 </table>
 <br/><br/>
 Examples:
 <ul>
 <li>getMVTFeatureAttributes(@geometry,'EPSG:31255','Bodenkarte',50000) -> JSON String (i.e. { 'ackerwert': NULL, 'ausgangsmaterial': 5820, 'bod... })</li>
 <li>map_get(getMVTFeatureAttributes(@geometry,'EPSG:31255','Bodenkarte',50000),'kurzbezeichnung') -> 'wsEG'</li>
 </ul>
 
 """
 mvt_layer = QgsProject.instance().mapLayersByName(mvt_layer_name)[0]
 origin_crs = QgsCoordinateReferenceSystem(epsg)
 target_crs = mvt_layer.dataProvider().crs()
 trans = QgsCoordinateTransform(origin_crs, target_crs, QgsProject.instance())
 geom.transform(trans)
 print(geom)
 bbox = QgsGeometry().fromRect(geom.buffer(0.00001, 1).boundingBox())
 # setting up selection context
 context = QgsSelectionContext()
 context.setScale(sscale)
 mvt_layer.selectByGeometry(bbox, context)
 if mvt_layer.selectedFeatureCount() > 0:
 return mvt_layer.selectedFeatures()[0].attributeMap()

• Taras: So, happy that I could help you ^_^ Freut mich :)