Implement SQLAlchemy native MVT query

Fix typo

Rename base tile provider

Fix Tile HTML highlighting

Fix missing tile rename

Update CRS conversion

Change fallback for vector tile id

Undo rename tile provider

Reduce diff

Cleanup mvt-postgres provider

Re-add ST_CurveToLine

97ba024

Author: webb-ben

docs/source/data-publishing/ogcapi-tiles.rst

@@ -140,9 +140,6 @@ MVT-postgresql
 .. note::
    Must have PostGIS installed with protobuf-c support 
 
-.. note::
-   Geometry must be using EPSG:4326
-
 This provider gives support to serving tiles generated using `PostgreSQL <https://www.postgresql.org/>`_ with `PostGIS <https://postgis.net/>`_.
 The tiles are rendered on-the-fly using `ST_AsMVT <https://postgis.net/docs/ST_AsMVT.html>`_ and related methods.
 

pygeoapi/api/tiles.py

@@ -246,7 +246,7 @@ def get_collection_tiles_data(
         if content is None:
             msg = 'identifier not found'
             return api.get_exception(
-                HTTPStatus.NO_CONTENT, headers, format_, 'NocContent', msg)
+                HTTPStatus.NO_CONTENT, headers, format_, 'NoContent', msg)
         else:
             return headers, HTTPStatus.OK, content
 

pygeoapi/provider/mvt_postgresql.py

@@ -3,10 +3,12 @@
 # Authors: Prajwal Amaravati <prajwal.s@satsure.co>
 #          Tanvi Prasad <tanvi.prasad@cdpg.org.in>
 #          Bryan Robert <bryan.robert@cdpg.org.in>
+#          Benjamin Webb <bwebb@lincolninst.edu>
 #
 # Copyright (c) 2025 Prajwal Amaravati
 # Copyright (c) 2025 Tanvi Prasad
 # Copyright (c) 2025 Bryan Robert
+# Copyright (c) 2025 Benjamin Webb
 #
 # Permission is hereby granted, free of charge, to any person
 # obtaining a copy of this software and associated documentation
@@ -34,20 +36,20 @@
 from copy import deepcopy
 import logging
 
-from sqlalchemy.sql import text
-
+from sqlalchemy.sql import func, select
+from sqlalchemy.orm import Session
 from pygeoapi.models.provider.base import (
     TileSetMetadata, TileMatrixSetEnum, LinkType)
 from pygeoapi.provider.base import ProviderConnectionError
 from pygeoapi.provider.base_mvt import BaseMVTProvider
 from pygeoapi.provider.postgresql import PostgreSQLProvider
 from pygeoapi.provider.tile import ProviderTileNotFoundError
-from pygeoapi.util import url_join
+from pygeoapi.util import url_join, get_crs_from_uri
 
 LOGGER = logging.getLogger(__name__)
 
 
-class MVTPostgreSQLProvider(BaseMVTProvider):
+class MVTPostgreSQLProvider(PostgreSQLProvider, BaseMVTProvider):
     """
     MVT PostgreSQL Provider
     Provider for serving tiles rendered on-the-fly from
@@ -62,48 +64,37 @@ def __init__(self, provider_def):
 
         :returns: pygeoapi.provider.MVT.MVTPostgreSQLProvider
         """
-
-        super().__init__(provider_def)
-
         pg_def = deepcopy(provider_def)
         # delete the zoom option before initializing the PostgreSQL provider
         # that provider breaks otherwise
-        del pg_def["options"]["zoom"]
-        self.postgres = PostgreSQLProvider(pg_def)
-
-        self.layer_name = provider_def["table"]
-        self.table = provider_def['table']
-        self.id_field = provider_def['id_field']
-        self.geom = provider_def.get('geom_field', 'geom')
-
-        LOGGER.debug(f'DB connection: {repr(self.postgres._engine.url)}')
-
-    def __repr__(self):
-        return f'<MVTPostgreSQLProvider> {self.data}'
+        del pg_def['options']['zoom']
+        PostgreSQLProvider.__init__(self, pg_def)
+        BaseMVTProvider.__init__(self, provider_def)
 
-    @property
-    def service_url(self):
-        return self._service_url
+    def get_fields(self):
+        """
+        Get Postgres columns
 
-    @property
-    def service_metadata_url(self):
-        return self._service_metadata_url
+        :returns: `list` of columns
+        """
+        return [
+            column for column in self.table_model.__table__.columns
+            if column.name != self.geom
+        ]
 
     def get_layer(self):
         """
-        Extracts layer name from url
+        Use table name as layer name
 
-        :returns: layer name
+        :returns: `str` of layer name
         """
-
-        return self.layer_name
+        return self.table
 
     def get_tiling_schemes(self):
-
         return [
-                TileMatrixSetEnum.WEBMERCATORQUAD.value,
-                TileMatrixSetEnum.WORLDCRS84QUAD.value
-            ]
+            TileMatrixSetEnum.WEBMERCATORQUAD.value,
+            TileMatrixSetEnum.WORLDCRS84QUAD.value
+        ]
 
     def get_tiles_service(self, baseurl=None, servicepath=None,
                           dirpath=None, tile_type=None):
@@ -118,13 +109,14 @@ def get_tiles_service(self, baseurl=None, servicepath=None,
         :returns: `dict` of item tile service
         """
 
-        super().get_tiles_service(baseurl, servicepath,
-                                  dirpath, tile_type)
+        BaseMVTProvider.get_tiles_service(self,
+                                          baseurl, servicepath,
+                                          dirpath, tile_type)
 
         self._service_url = servicepath
         return self.get_tms_links()
 
-    def get_tiles(self, layer=None, tileset=None,
+    def get_tiles(self, layer='default', tileset=None,
                   z=None, y=None, x=None, format_=None):
         """
         Gets tile
@@ -141,64 +133,55 @@ def get_tiles(self, layer=None, tileset=None,
         if format_ == 'mvt':
             format_ = self.format_type
 
-        fields_arr = self.postgres.get_fields().keys()
-        fields = ', '.join(['"' + f + '"' for f in fields_arr])
-        if len(fields) != 0:
-            fields = ',' + fields
-
-        query = ''
-        if tileset == TileMatrixSetEnum.WEBMERCATORQUAD.value.tileMatrixSet:
-            if not self.is_in_limits(TileMatrixSetEnum.WEBMERCATORQUAD.value, z, x, y): # noqa
-                raise ProviderTileNotFoundError
-
-            query = text("""
-                WITH
-                    bounds AS (
-                        SELECT ST_TileEnvelope(:z, :x, :y) AS boundgeom
-                    ),
-                    mvtgeom AS (
-                        SELECT ST_AsMVTGeom(ST_Transform(ST_CurveToLine({geom}), 3857), bounds.boundgeom) AS geom {fields}
-                        FROM "{table}", bounds
-                        WHERE ST_Intersects({geom}, ST_Transform(bounds.boundgeom, 4326))
-                    )
-                SELECT ST_AsMVT(mvtgeom, 'default') FROM mvtgeom;
-            """.format(geom=self.geom, table=self.table, fields=fields)) # noqa
-
-        if tileset == TileMatrixSetEnum.WORLDCRS84QUAD.value.tileMatrixSet:
-            if not self.is_in_limits(TileMatrixSetEnum.WORLDCRS84QUAD.value, z, x, y): # noqa
-                raise ProviderTileNotFoundError
-
-            query = text("""
-                WITH
-                    bounds AS (
-                        SELECT ST_TileEnvelope(:z, :x, :y,
-                                'SRID=4326;POLYGON((-180 -90,-180 90,180 90,180 -90,-180 -90))'::geometry) AS boundgeom
-                    ),
-                    mvtgeom AS (
-                        SELECT ST_AsMVTGeom(ST_CurveToLine({geom}), bounds.boundgeom) AS geom {fields}
-                        FROM "{table}", bounds
-                        WHERE ST_Intersects({geom}, bounds.boundgeom)
-                    )
-                SELECT ST_AsMVT(mvtgeom, 'default') FROM mvtgeom;
-            """.format(geom=self.geom, table=self.table, fields=fields)) # noqa
-
-        with self.postgres._engine.connect() as session:
-            result = session.execute(query, {
-                    'z': z,
-                    'y': y,
-                    'x': x
-            }).fetchone()
-
-            if len(bytes(result[0])) == 0:
-                return None
-            return bytes(result[0])
+        [tileset_schema] = [
+            schema for schema in self.get_tiling_schemes()
+            if tileset == schema.tileMatrixSet
+        ]
+        if not self.is_in_limits(tileset_schema, z, x, y):
+            LOGGER.warning(f'Tile {z}/{x}/{y} not found')
+            return ProviderTileNotFoundError
+
+        storage_srid = get_crs_from_uri(self.storage_crs).to_string()
+        out_srid = get_crs_from_uri(tileset_schema.crs).to_string()
+
+        tile_envelope = func.ST_TileEnvelope(z, x, y)
+        envelope = (
+            select(tile_envelope.label('bounds'))
+            .cte('envelope')
+        )
+
+        geom_column = getattr(self.table_model, self.geom)
+        mvtgeom = (
+            func.ST_AsMVTGeom(
+                func.ST_Transform(geom_column, out_srid),
+                func.ST_Transform(envelope.c.bounds, out_srid))
+            .label('mvtgeom')
+        )
+
+        geom_filter = geom_column.intersects(
+            func.ST_Transform(tile_envelope, storage_srid))
+        mvtrow = (
+            select(*self.get_fields(), mvtgeom)
+            .filter(geom_filter)
+            .select_from(self.table_model)
+            .cte('mvtrow')
+            .table_valued()
+        )
+
+        mvt_query = select(func.ST_AsMVT(mvtrow, layer))
+
+        with Session(self._engine) as session:
+            result = session.execute(mvt_query).scalar()
+
+        return bytes(result) or None
 
     def get_html_metadata(self, dataset, server_url, layer, tileset,
                           title, description, keywords, **kwargs):
 
         service_url = url_join(
             server_url,
-            f'collections/{dataset}/tiles/{tileset}/{{tileMatrix}}/{{tileRow}}/{{tileCol}}?f=mvt')  # noqa
+            f'collections/{dataset}/tiles/{tileset}'
+            '{tileMatrix}/{tileRow}/{tileCol}?f=mvt')
         metadata_url = url_join(
             server_url,
             f'collections/{dataset}/tiles/{tileset}/metadata')
@@ -217,9 +200,10 @@ def get_default_metadata(self, dataset, server_url, layer, tileset,
 
         service_url = url_join(
             server_url,
-            f'collections/{dataset}/tiles/{tileset}/{{tileMatrix}}/{{tileRow}}/{{tileCol}}?f=mvt')  # noqa
+            f'collections/{dataset}/tiles/{tileset}',
+            '{tileMatrix}/{tileRow}/{tileCol}?f=mvt'
+            )
 
-        content = {}
         tiling_schemes = self.get_tiling_schemes()
         # Default values
         tileMatrixSetURI = tiling_schemes[0].tileMatrixSetURI
@@ -231,17 +215,18 @@ def get_default_metadata(self, dataset, server_url, layer, tileset,
                 tileMatrixSetURI = schema.tileMatrixSetURI
 
                 tiling_scheme_url = url_join(
-                    server_url, f'/TileMatrixSets/{schema.tileMatrixSet}')
-                tiling_scheme_url_type = "application/json"
+                    server_url, 'TileMatrixSets', schema.tileMatrixSet)
+                tiling_scheme_url_type = 'application/json'
                 tiling_scheme_url_title = f'{schema.tileMatrixSet} tile matrix set definition' # noqa
 
-                tiling_scheme = LinkType(href=tiling_scheme_url,
-                                         rel="http://www.opengis.net/def/rel/ogc/1.0/tiling-scheme", # noqa
-                                         type_=tiling_scheme_url_type,
-                                         title=tiling_scheme_url_title)
+                tiling_scheme = LinkType(
+                    href=tiling_scheme_url,
+                    rel='http://www.opengis.net/def/rel/ogc/1.0/tiling-scheme',
+                    type_=tiling_scheme_url_type,
+                    title=tiling_scheme_url_title)
 
         if tiling_scheme is None:
-            msg = f'Could not identify a valid tiling schema'  # noqa
+            msg = 'Could not identify a valid tiling schema'
             LOGGER.error(msg)
             raise ProviderConnectionError(msg)
 
@@ -250,9 +235,9 @@ def get_default_metadata(self, dataset, server_url, layer, tileset,
                                   tileMatrixSetURI=tileMatrixSetURI)
 
         links = []
-        service_url_link_type = "application/vnd.mapbox-vector-tile"
+        service_url_link_type = 'application/vnd.mapbox-vector-tile'
         service_url_link_title = f'{tileset} vector tiles for {layer}'
-        service_url_link = LinkType(href=service_url, rel="item",
+        service_url_link = LinkType(href=service_url, rel='item',
                                     type_=service_url_link_type,
                                     title=service_url_link_title)
 
@@ -261,4 +246,7 @@ def get_default_metadata(self, dataset, server_url, layer, tileset,
 
         content.links = links
 
-        return content.dict(exclude_none=True)
+        return content.model_dump(exclude_none=True)
+
+    def __repr__(self):
+        return f'<MVTPostgreSQLProvider> {self.data}'

pygeoapi/provider/tile.py

@@ -56,7 +56,6 @@ def __init__(self, provider_def):
         self.mimetype = provider_def['format']['mimetype']
         self.options = provider_def.get('options')
         self.tile_type = None
-        self.fields = {}
 
     def get_layer(self):
         """

pygeoapi/templates/collections/tiles/index.html

@@ -108,7 +108,7 @@ <h3>Tiles</h3>
             maxZoom: {{ data['maxzoom'] }},
             indexMaxZoom: 5,
             getFeatureId: function(feat) {
-                return feat.properties["id"]
+                return feat.properties.id || feat.properties.fid || feat.properties.uri;
             }
         };
 
@@ -129,7 +129,7 @@ <h3>Tiles</h3>
                     .openOn(map);
 
                 clearHighlight();
-                highlight = e.layer.properties.id;
+                highlight = e.layer.properties.id || e.layer.properties.fid || e.layer.properties.uri;
                 tilesPbfLayer.setFeatureStyle(highlight, {
                     weight: 2,
                     color: 'red',