From 6821d568564008e224a9ddc5abcae008e613dd86 Mon Sep 17 00:00:00 2001 From: Benjamin Webb Date: Thu, 29 May 2025 23:03:42 -0400 Subject: [PATCH 01/15] 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 https://github.com/geopython/pygeoapi/pull/1979/commits/97ba024f4aa8a253ff04ac718bfb3feba528080a --- docs/source/data-publishing/ogcapi-tiles.rst | 3 - pygeoapi/api/tiles.py | 2 +- pygeoapi/provider/mvt_postgresql.py | 188 ++++++++---------- pygeoapi/provider/tile.py | 1 - .../templates/collections/tiles/index.html | 4 +- 5 files changed, 91 insertions(+), 107 deletions(-) diff --git a/docs/source/data-publishing/ogcapi-tiles.rst b/docs/source/data-publishing/ogcapi-tiles.rst index 66ab76c91..c4941bd48 100644 --- a/docs/source/data-publishing/ogcapi-tiles.rst +++ b/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 `_ with `PostGIS `_. The tiles are rendered on-the-fly using `ST_AsMVT `_ and related methods. diff --git a/pygeoapi/api/tiles.py b/pygeoapi/api/tiles.py index 0d457e596..cb957eee8 100644 --- a/pygeoapi/api/tiles.py +++ b/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 diff --git a/pygeoapi/provider/mvt_postgresql.py b/pygeoapi/provider/mvt_postgresql.py index c8fa2fc1a..c27d52ce1 100644 --- a/pygeoapi/provider/mvt_postgresql.py +++ b/pygeoapi/provider/mvt_postgresql.py @@ -3,10 +3,12 @@ # Authors: Prajwal Amaravati # Tanvi Prasad # Bryan Robert +# Benjamin Webb # # 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.sql 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' {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' {self.data}' diff --git a/pygeoapi/provider/tile.py b/pygeoapi/provider/tile.py index 7e6352ae3..a57488db1 100644 --- a/pygeoapi/provider/tile.py +++ b/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): """ diff --git a/pygeoapi/templates/collections/tiles/index.html b/pygeoapi/templates/collections/tiles/index.html index a9bd56cd9..75e1f9d36 100644 --- a/pygeoapi/templates/collections/tiles/index.html +++ b/pygeoapi/templates/collections/tiles/index.html @@ -108,7 +108,7 @@

Tiles

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 @@

Tiles

.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', From 08ada270bd8aeaa44f38d2750c1187205e87b929 Mon Sep 17 00:00:00 2001 From: Benjamin Webb Date: Mon, 2 Jun 2025 19:30:45 -0400 Subject: [PATCH 02/15] Add ST_CurveToLine for output geom --- pygeoapi/provider/mvt_postgresql.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pygeoapi/provider/mvt_postgresql.py b/pygeoapi/provider/mvt_postgresql.py index c27d52ce1..7b96f88ed 100644 --- a/pygeoapi/provider/mvt_postgresql.py +++ b/pygeoapi/provider/mvt_postgresql.py @@ -153,7 +153,7 @@ def get_tiles(self, layer='default', tileset=None, geom_column = getattr(self.table_model, self.geom) mvtgeom = ( func.ST_AsMVTGeom( - func.ST_Transform(geom_column, out_srid), + func.ST_Transform(func.ST_CurveToLine(geom_column), out_srid), func.ST_Transform(envelope.c.bounds, out_srid)) .label('mvtgeom') ) From b3a9e9ae6545e10e35995fba538477eef03ab867 Mon Sep 17 00:00:00 2001 From: Benjamin Webb Date: Mon, 2 Jun 2025 19:38:37 -0400 Subject: [PATCH 03/15] Set id field for vector files --- pygeoapi/provider/mvt_postgresql.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/pygeoapi/provider/mvt_postgresql.py b/pygeoapi/provider/mvt_postgresql.py index 7b96f88ed..a27c2aa66 100644 --- a/pygeoapi/provider/mvt_postgresql.py +++ b/pygeoapi/provider/mvt_postgresql.py @@ -78,8 +78,9 @@ def get_fields(self): :returns: `list` of columns """ return [ - column for column in self.table_model.__table__.columns - if column.name != self.geom + c.label('id') if c.name == self.id_field else c + for c in self.table_model.__table__.columns + if c.name != self.geom ] def get_layer(self): From 2fb298d14b898f624d21484474c3f4cc46a999f5 Mon Sep 17 00:00:00 2001 From: Benjamin Webb Date: Tue, 10 Jun 2025 11:48:13 -0400 Subject: [PATCH 04/15] Add MVT Postgres Provider test Fix Flake8 --- pygeoapi/provider/mvt_postgresql.py | 26 +-- tests/test_postgresql_mvt_provider.py | 261 ++++++++++++++++++++++++++ 2 files changed, 274 insertions(+), 13 deletions(-) create mode 100644 tests/test_postgresql_mvt_provider.py diff --git a/pygeoapi/provider/mvt_postgresql.py b/pygeoapi/provider/mvt_postgresql.py index a27c2aa66..995428964 100644 --- a/pygeoapi/provider/mvt_postgresql.py +++ b/pygeoapi/provider/mvt_postgresql.py @@ -144,37 +144,37 @@ def get_tiles(self, layer='default', tileset=None, storage_srid = get_crs_from_uri(self.storage_crs).to_string() out_srid = get_crs_from_uri(tileset_schema.crs).to_string() + envelope = func.ST_TileEnvelope(z, x, y).label('bounds') - tile_envelope = func.ST_TileEnvelope(z, x, y) - envelope = ( - select(tile_envelope.label('bounds')) - .cte('envelope') + geom_column = getattr(self.table_model, self.geom) + geom_filter = geom_column.intersects( + func.ST_Transform(envelope, storage_srid) ) - geom_column = getattr(self.table_model, self.geom) mvtgeom = ( func.ST_AsMVTGeom( func.ST_Transform(func.ST_CurveToLine(geom_column), out_srid), - func.ST_Transform(envelope.c.bounds, out_srid)) + func.ST_Transform(envelope, out_srid)) .label('mvtgeom') ) - geom_filter = geom_column.intersects( - func.ST_Transform(tile_envelope, storage_srid)) mvtrow = ( - select(*self.get_fields(), mvtgeom) + select(mvtgeom, *self.fields) .filter(geom_filter) - .select_from(self.table_model) .cte('mvtrow') .table_valued() ) - mvt_query = select(func.ST_AsMVT(mvtrow, layer)) + mvtquery = select( + func.ST_AsMVT(mvtrow, layer) + ) with Session(self._engine) as session: - result = session.execute(mvt_query).scalar() + result = bytes( + session.execute(mvtquery).scalar() + ) or None - return bytes(result) or None + return result def get_html_metadata(self, dataset, server_url, layer, tileset, title, description, keywords, **kwargs): diff --git a/tests/test_postgresql_mvt_provider.py b/tests/test_postgresql_mvt_provider.py new file mode 100644 index 000000000..9317cdf0b --- /dev/null +++ b/tests/test_postgresql_mvt_provider.py @@ -0,0 +1,261 @@ +# ================================================================= +# +# Authors: Ben Webb +# +# Copyright (c) 2025 Ben Webb +# +# Permission is hereby granted, free of charge, to any person +# obtaining a copy of this software and associated documentation +# files (the "Software"), to deal in the Software without +# restriction, including without limitation the rights to use, +# copy, modify, merge, publish, distribute, sublicense, and/or sell +# copies of the Software, and to permit persons to whom the +# Software is furnished to do so, subject to the following +# conditions: +# +# The above copyright notice and this permission notice shall be +# included in all copies or substantial portions of the Software. +# +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +# EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES +# OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND +# NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT +# HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, +# WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +# FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR +# OTHER DEALINGS IN THE SOFTWARE. +# +# ================================================================= + +# Needs to be run like: python3 -m pytest +# See pygeoapi/provider/postgresql.py for instructions on setting up +# test database in Docker + +import os +import pytest + +from pygeoapi.provider.mvt_postgresql import MVTPostgreSQLProvider +from pygeoapi.provider.tile import ProviderTileNotFoundError + +PASSWORD = os.environ.get('POSTGRESQL_PASSWORD', 'postgres') +SERVER_URL = 'http://localhost' +DATASET = 'hotosm_bdi_waterways' + + +@pytest.fixture() +def config(): + return { + 'name': 'MVT-postgresql', + 'type': 'tile', + 'data': { + 'host': '127.0.0.1', + 'dbname': 'test', + 'user': 'postgres', + 'password': PASSWORD, + 'search_path': ['osm', 'public'] + }, + 'id_field': 'osm_id', + 'table': 'hotosm_bdi_waterways', + 'geom_field': 'foo_geom', + 'options': { + 'zoom': { + 'min': 0, + 'max': 15 + } + }, + 'format': { + 'name': 'pbf', + 'mimetype': 'application/vnd.mapbox-vector-tile' + }, + 'storage_crs': 'http://www.opengis.net/def/crs/EPSG/0/4326' + } + + +def test_metadata(config): + """Testing query for a valid JSON object with geometry""" + p = MVTPostgreSQLProvider(config) + + assert p.table == 'hotosm_bdi_waterways' + assert p.geom == 'foo_geom' + assert p.id_field == 'osm_id' + assert p.get_layer() == config['table'] + + +def test_get_tiling_schemes(config): + provider = MVTPostgreSQLProvider(config) + + schemes = provider.get_tiling_schemes() + assert any(s.tileMatrixSet == 'WebMercatorQuad' for s in schemes) + assert any(s.tileMatrixSet == 'WorldCRS84Quad' for s in schemes) + + +def test_get_default_metadata_WebMercatorQuad(config): + provider = MVTPostgreSQLProvider(config) + ts = 'WebMercatorQuad' + + md = provider.get_default_metadata( + dataset=DATASET, + server_url=SERVER_URL, + layer='layer1', + tileset=ts, + title='Waterways', + description='OpenStreetMap Waterways', + keywords=['osm', 'rivers'] + ) + + assert md['tileMatrixSetURI'] == \ + 'http://www.opengis.net/def/tilematrixset/OGC/1.0/WebMercatorQuad' + assert md['crs'] == \ + 'http://www.opengis.net/def/crs/EPSG/0/3857' + + assert 'links' in md + assert any(link['rel'].endswith('tiling-scheme') for link in md['links']) + + [tile_format,] = [link for link in md['links'] if link['rel'] == 'item'] + assert tile_format['href'].startswith(SERVER_URL) + assert DATASET in tile_format['href'] + assert ts in tile_format['href'] + + +def test_get_default_metadata_WorldCRS84Quad(config): + provider = MVTPostgreSQLProvider(config) + ts = 'WorldCRS84Quad' + + md = provider.get_default_metadata( + dataset=DATASET, + server_url=SERVER_URL, + layer='layer1', + tileset=ts, + title='Waterways', + description='OpenStreetMap Waterways', + keywords=['osm', 'rivers'] + ) + + assert md['tileMatrixSetURI'] == \ + 'http://www.opengis.net/def/tilematrixset/OGC/1.0/WorldCRS84Quad' + assert md['crs'] == \ + 'http://www.opengis.net/def/crs/OGC/1.3/CRS84' + + assert 'links' in md + assert any(link['rel'].endswith('tiling-scheme') for link in md['links']) + + [tile_format,] = [link for link in md['links'] if link['rel'] == 'item'] + assert tile_format['href'].startswith(SERVER_URL) + assert DATASET in tile_format['href'] + assert ts in tile_format['href'] + + +def test_get_html_metadata_WebMercatorQuad(config): + provider = MVTPostgreSQLProvider(config) + title = 'Waterways' + ts = 'WebMercatorQuad' + + md = provider.get_html_metadata( + dataset=DATASET, + server_url=SERVER_URL, + layer='layer1', + tileset=ts, + title=title, + description='OpenStreetMap Waterways', + keywords=['osm', 'rivers'] + ) + + assert md['id'] == DATASET + assert md['title'] == title + assert md['tileset'] == ts + + assert md['collections_path'].startswith(SERVER_URL) + assert DATASET in md['collections_path'] + assert ts in md['collections_path'] + + assert md['json_url'].startswith(SERVER_URL) + assert DATASET in md['json_url'] + assert ts in md['json_url'] + + +def test_get_html_metadata_WorldCRS84Quad(config): + provider = MVTPostgreSQLProvider(config) + title = 'Waterways' + ts = 'WorldCRS84Quad' + + md = provider.get_html_metadata( + dataset=DATASET, + server_url=SERVER_URL, + layer='layer1', + tileset=ts, + title=title, + description='OpenStreetMap Waterways', + keywords=['osm', 'rivers'] + ) + + assert md['id'] == DATASET + assert md['title'] == title + assert md['tileset'] == ts + + assert md['collections_path'].startswith(SERVER_URL) + assert DATASET in md['collections_path'] + assert ts in md['collections_path'] + + assert md['json_url'].startswith(SERVER_URL) + assert DATASET in md['json_url'] + assert ts in md['json_url'] + + +def test_get_tiles_WebMercatorQuad(config): + p = MVTPostgreSQLProvider(config) + tileset = 'WebMercatorQuad' + + # Valid tile, no content + z, x, y = 14, 10200, 10300 + tile = p.get_tiles( + tileset=tileset, + z=z, x=x, y=y, + ) + assert tile is None + + # Valid tile, content + z, x, y = 10, 595, 521 + tile = p.get_tiles( + tileset=tileset, + z=z, x=x, y=y, + ) + assert isinstance(tile, bytes) + assert len(tile) > 0 + + # Tile does not exist in matrixset + z, x, y = 1, 1000000, 1000000 + result = p.get_tiles( + tileset=tileset, + z=z, x=x, y=y + ) + assert result == ProviderTileNotFoundError + + +def test_get_tiles_WorldCRS84Quad(config): + p = MVTPostgreSQLProvider(config) + tileset = 'WorldCRS84Quad' + + # Valid tile, no content + z, x, y = 14, 10200, 10300 + tile = p.get_tiles( + tileset=tileset, + z=z, x=x, y=y, + ) + assert tile is None + + # Valid tile, content + z, x, y = 10, 595, 521 + tile = p.get_tiles( + tileset=tileset, + z=z, x=x, y=y, + ) + assert isinstance(tile, bytes) + assert len(tile) > 0 + + # Tile does not exist in matrixset + z, x, y = 1, 1000000, 1000000 + result = p.get_tiles( + tileset=tileset, + z=z, x=x, y=y + ) + assert result == ProviderTileNotFoundError From 6b7b3c8f5dbad036b556c37373fe443ff408680e Mon Sep 17 00:00:00 2001 From: Benjamin Webb Date: Tue, 10 Jun 2025 11:59:49 -0400 Subject: [PATCH 05/15] Fix indentation in MVT-Postgres docs --- docs/source/data-publishing/ogcapi-tiles.rst | 36 ++++++++++---------- 1 file changed, 18 insertions(+), 18 deletions(-) diff --git a/docs/source/data-publishing/ogcapi-tiles.rst b/docs/source/data-publishing/ogcapi-tiles.rst index c4941bd48..871a52cf1 100644 --- a/docs/source/data-publishing/ogcapi-tiles.rst +++ b/docs/source/data-publishing/ogcapi-tiles.rst @@ -149,24 +149,24 @@ This code block shows how to configure pygeoapi to render Mapbox vector tiles fr providers: - type: tile - name: MVT-postgresql - data: - host: 127.0.0.1 - port: 3010 # Default 5432 if not provided - dbname: test - user: postgres - password: postgres - search_path: [osm, public] - id_field: osm_id - table: hotosm_bdi_waterways - geom_field: foo_geom - options: - zoom: - min: 0 - max: 15 - format: - name: pbf - mimetype: application/vnd.mapbox-vector-tile + name: MVT-postgresql + data: + host: 127.0.0.1 + port: 3010 # Default 5432 if not provided + dbname: test + user: postgres + password: postgres + search_path: [osm, public] + id_field: osm_id + table: hotosm_bdi_waterways + geom_field: foo_geom + options: + zoom: + min: 0 + max: 15 + format: + name: pbf + mimetype: application/vnd.mapbox-vector-tile PostgreSQL-related connection options can also be added to `options`. Please refer to the :ref:`PostgreSQL OGC Features Provider` documentation for more information. From 23f35a7026ee7b0a5f075ab75ebe4c5871918e3b Mon Sep 17 00:00:00 2001 From: Benjamin Webb Date: Sun, 15 Jun 2025 18:40:19 -0400 Subject: [PATCH 06/15] Fix get_metadata() - Fix get_metadata function - Add unit test for get_metadata --- pygeoapi/provider/mvt_postgresql.py | 2 +- tests/test_postgresql_mvt_provider.py | 21 +++++++++++++++++++++ 2 files changed, 22 insertions(+), 1 deletion(-) diff --git a/pygeoapi/provider/mvt_postgresql.py b/pygeoapi/provider/mvt_postgresql.py index 995428964..0edaaf99b 100644 --- a/pygeoapi/provider/mvt_postgresql.py +++ b/pygeoapi/provider/mvt_postgresql.py @@ -49,7 +49,7 @@ LOGGER = logging.getLogger(__name__) -class MVTPostgreSQLProvider(PostgreSQLProvider, BaseMVTProvider): +class MVTPostgreSQLProvider(BaseMVTProvider, PostgreSQLProvider): """ MVT PostgreSQL Provider Provider for serving tiles rendered on-the-fly from diff --git a/tests/test_postgresql_mvt_provider.py b/tests/test_postgresql_mvt_provider.py index 9317cdf0b..82ed63d9f 100644 --- a/tests/test_postgresql_mvt_provider.py +++ b/tests/test_postgresql_mvt_provider.py @@ -74,12 +74,33 @@ def config(): def test_metadata(config): """Testing query for a valid JSON object with geometry""" p = MVTPostgreSQLProvider(config) + ts = 'WebMercatorQuad' assert p.table == 'hotosm_bdi_waterways' assert p.geom == 'foo_geom' assert p.id_field == 'osm_id' assert p.get_layer() == config['table'] + md = p.get_metadata( + dataset=DATASET, + server_url=SERVER_URL, + layer='layer1', + tileset=ts, + metadata_format='json', + title='Waterways', + description='OpenStreetMap Waterways', + ) + assert md['crs'] == \ + 'http://www.opengis.net/def/crs/EPSG/0/3857' + + assert 'links' in md + assert any(link['rel'].endswith('tiling-scheme') for link in md['links']) + + [tile_format,] = [link for link in md['links'] if link['rel'] == 'item'] + assert tile_format['href'].startswith(SERVER_URL) + assert DATASET in tile_format['href'] + assert ts in tile_format['href'] + def test_get_tiling_schemes(config): provider = MVTPostgreSQLProvider(config) From 8e0b7f3af2adde2628e245ef6c2be06d8740340a Mon Sep 17 00:00:00 2001 From: Benjamin Webb Date: Sun, 15 Jun 2025 18:59:23 -0400 Subject: [PATCH 07/15] Fix MVT feature properties - Fix MVT field selection - Add test for vector tile properties --- pygeoapi/provider/mvt_postgresql.py | 20 ++++++++++++++------ tests/test_postgresql_mvt_provider.py | 21 +++++++++++++++++++++ 2 files changed, 35 insertions(+), 6 deletions(-) diff --git a/pygeoapi/provider/mvt_postgresql.py b/pygeoapi/provider/mvt_postgresql.py index 0edaaf99b..3d18f44f2 100644 --- a/pygeoapi/provider/mvt_postgresql.py +++ b/pygeoapi/provider/mvt_postgresql.py @@ -77,11 +77,19 @@ def get_fields(self): :returns: `list` of columns """ - return [ - c.label('id') if c.name == self.id_field else c - for c in self.table_model.__table__.columns - if c.name != self.geom - ] + if not self._fields: + for column in self.table_model.__table__.columns: + LOGGER.debug(f'Testing {column.name}') + if column.name == self.geom: + continue + + self._fields[str(column.name)] = ( + column.label('id') + if column.name == self.id_field else + column + ) + + return self._fields def get_layer(self): """ @@ -159,7 +167,7 @@ def get_tiles(self, layer='default', tileset=None, ) mvtrow = ( - select(mvtgeom, *self.fields) + select(mvtgeom, *self.fields.values()) .filter(geom_filter) .cte('mvtrow') .table_valued() diff --git a/tests/test_postgresql_mvt_provider.py b/tests/test_postgresql_mvt_provider.py index 82ed63d9f..ee30caeab 100644 --- a/tests/test_postgresql_mvt_provider.py +++ b/tests/test_postgresql_mvt_provider.py @@ -239,10 +239,21 @@ def test_get_tiles_WebMercatorQuad(config): tile = p.get_tiles( tileset=tileset, z=z, x=x, y=y, + layer=p.get_layer() ) assert isinstance(tile, bytes) assert len(tile) > 0 + # Layer name + assert b'waterways' in tile + assert b'default' not in tile + + # Feature properties + assert b'id' in tile + assert b'waterway' in tile + assert b'name' in tile + assert b'z_index' in tile + # Tile does not exist in matrixset z, x, y = 1, 1000000, 1000000 result = p.get_tiles( @@ -273,6 +284,16 @@ def test_get_tiles_WorldCRS84Quad(config): assert isinstance(tile, bytes) assert len(tile) > 0 + # Layer name + assert b'waterways' not in tile + assert b'default' in tile + + # Feature properties + assert b'id' in tile + assert b'waterway' in tile + assert b'name' in tile + assert b'z_index' in tile + # Tile does not exist in matrixset z, x, y = 1, 1000000, 1000000 result = p.get_tiles( From 5eb616c64aac5422dcf6221d1e6882e45814c931 Mon Sep 17 00:00:00 2001 From: Benjamin Webb Date: Sun, 15 Jun 2025 19:02:51 -0400 Subject: [PATCH 08/15] Fix flake8 --- pygeoapi/provider/mvt_postgresql.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pygeoapi/provider/mvt_postgresql.py b/pygeoapi/provider/mvt_postgresql.py index 3d18f44f2..a6831ab8c 100644 --- a/pygeoapi/provider/mvt_postgresql.py +++ b/pygeoapi/provider/mvt_postgresql.py @@ -84,7 +84,7 @@ def get_fields(self): continue self._fields[str(column.name)] = ( - column.label('id') + column.label('id') if column.name == self.id_field else column ) From 740180c6ff27040eb8f381038d8e372dae91462d Mon Sep 17 00:00:00 2001 From: Benjamin Webb Date: Sun, 15 Jun 2025 19:03:07 -0400 Subject: [PATCH 09/15] Add storage_crs information to docs --- docs/source/data-publishing/ogcapi-tiles.rst | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/docs/source/data-publishing/ogcapi-tiles.rst b/docs/source/data-publishing/ogcapi-tiles.rst index 871a52cf1..173d08057 100644 --- a/docs/source/data-publishing/ogcapi-tiles.rst +++ b/docs/source/data-publishing/ogcapi-tiles.rst @@ -160,6 +160,7 @@ This code block shows how to configure pygeoapi to render Mapbox vector tiles fr id_field: osm_id table: hotosm_bdi_waterways geom_field: foo_geom + storage_crs: http://www.opengis.net/def/crs/EPSG/0/4326 options: zoom: min: 0 @@ -168,6 +169,9 @@ This code block shows how to configure pygeoapi to render Mapbox vector tiles fr name: pbf mimetype: application/vnd.mapbox-vector-tile +.. tip:: + Geometry must have correctly defined :ref:`storage_crs` + PostgreSQL-related connection options can also be added to `options`. Please refer to the :ref:`PostgreSQL OGC Features Provider` documentation for more information. WMTSFacade From 615288d4d5597b0d8042229a756725c5c0b3dc70 Mon Sep 17 00:00:00 2001 From: Benjamin Webb Date: Sun, 15 Jun 2025 19:11:22 -0400 Subject: [PATCH 10/15] Add test_postgresql_mvt_provider to GitHub CI --- .github/workflows/main.yml | 1 + 1 file changed, 1 insertion(+) diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml index 91b8d30af..3c556b84b 100644 --- a/.github/workflows/main.yml +++ b/.github/workflows/main.yml @@ -164,6 +164,7 @@ jobs: pytest tests/test_oracle_provider.py pytest tests/test_parquet_provider.py pytest tests/test_postgresql_provider.py + pytest tests/test_postgresql_mvt_provider.py pytest tests/test_mysql_provider.py pytest tests/test_rasterio_provider.py pytest tests/test_sensorthings_edr_provider.py From ddb1ba299c952ec7da3d374de01637d9e30a1c99 Mon Sep 17 00:00:00 2001 From: Benjamin Webb Date: Sun, 15 Jun 2025 20:14:45 -0400 Subject: [PATCH 11/15] Bump GitHub Action Postgis version Bump GitHub Postgis engine to >3 --- .github/workflows/main.yml | 16 +++++++++------- 1 file changed, 9 insertions(+), 7 deletions(-) diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml index 3c556b84b..373414f6c 100644 --- a/.github/workflows/main.yml +++ b/.github/workflows/main.yml @@ -24,6 +24,15 @@ jobs: - python-version: '3.10' env: PYGEOAPI_CONFIG: "$(pwd)/pygeoapi-config.yml" + + services: + postgres: + image: postgis/postgis:14-3.2 + ports: + - 5432:5432 + env: + POSTGRES_DB: test + POSTGRES_PASSWORD: ${{ secrets.DatabasePassword || 'postgres' }} steps: - name: Pre-pull Docker Images @@ -32,7 +41,6 @@ jobs: docker pull appropriate/curl:latest & docker pull elasticsearch:8.17.0 & docker pull opensearchproject/opensearch:2.18.0 & - docker pull mdillon/postgis:latest & docker pull mongo:8.0.4 & docker pull ghcr.io/cgs-earth/sensorthings-action:0.1.0 & docker pull postgis/postgis:14-3.2 & @@ -58,12 +66,6 @@ jobs: sudo sysctl -w vm.swappiness=1 sudo sysctl -w fs.file-max=262144 sudo sysctl -w vm.max_map_count=262144 - - name: Install and run PostgreSQL/PostGIS 📦 - uses: huaxk/postgis-action@v1 - with: - postgresql password: ${{ secrets.DatabasePassword || 'postgres' }} - postgresql db: 'test' - - name: "Install and run MySQL 📦" uses: mirromutth/mysql-action@v1.1 with: From f71a63f28daa74ddc0012089273488b81f061847 Mon Sep 17 00:00:00 2001 From: Benjamin Webb Date: Mon, 16 Jun 2025 21:02:08 -0400 Subject: [PATCH 12/15] Filter SQL connection options Prevent passing dict to SQL connection which throws an error because it cannot be cached --- pygeoapi/provider/mvt_postgresql.py | 7 +------ pygeoapi/provider/sql.py | 21 ++++++++++----------- 2 files changed, 11 insertions(+), 17 deletions(-) diff --git a/pygeoapi/provider/mvt_postgresql.py b/pygeoapi/provider/mvt_postgresql.py index a6831ab8c..2c925e3b4 100644 --- a/pygeoapi/provider/mvt_postgresql.py +++ b/pygeoapi/provider/mvt_postgresql.py @@ -33,7 +33,6 @@ # # ================================================================= -from copy import deepcopy import logging from sqlalchemy.sql import func, select @@ -64,11 +63,7 @@ def __init__(self, provider_def): :returns: pygeoapi.provider.MVT.MVTPostgreSQLProvider """ - pg_def = deepcopy(provider_def) - # delete the zoom option before initializing the PostgreSQL provider - # that provider breaks otherwise - del pg_def['options']['zoom'] - PostgreSQLProvider.__init__(self, pg_def) + PostgreSQLProvider.__init__(self, provider_def) BaseMVTProvider.__init__(self, provider_def) def get_fields(self): diff --git a/pygeoapi/provider/sql.py b/pygeoapi/provider/sql.py index 2775e705b..cdc90d2c7 100644 --- a/pygeoapi/provider/sql.py +++ b/pygeoapi/provider/sql.py @@ -107,7 +107,7 @@ def __init__( self, provider_def: dict, driver_name: str, - extra_conn_args: Optional[dict] + extra_conn_args: Optional[dict] = {} ): """ GenericSQLProvider Class constructor @@ -143,9 +143,7 @@ def __init__( LOGGER.debug(f'Configured Storage CRS: {self.storage_crs}') # Read table information from database - options = None - if provider_def.get('options'): - options = provider_def['options'] + options = provider_def.get('options', {}) self._store_db_parameters(provider_def['data'], options) self._engine = get_engine( driver_name, @@ -154,7 +152,7 @@ def __init__( self.db_name, self.db_user, self._db_password, - **(self.db_options or {}) | (extra_conn_args or {}) + **self.db_options | extra_conn_args ) self.table_model = get_table_model( self.table, self.id_field, self.db_search_path, self._engine @@ -443,7 +441,11 @@ def _store_db_parameters(self, parameters, options): # reflecting the table definition from the DB self.db_search_path = tuple(parameters.get('search_path', ['public'])) self._db_password = parameters.get('password') - self.db_options = options + self.db_options = { + k: v + for k, v in options.items() + if not isinstance(v, dict) + } def _sqlalchemy_to_feature(self, item, crs_transform_out=None): feature = {'type': 'Feature'} @@ -608,7 +610,7 @@ def get_engine( database: str, user: str, password: str, - **connection_options + **connect_args ): """Create SQL Alchemy engine.""" conn_str = URL.create( @@ -619,11 +621,8 @@ def get_engine( port=int(port), database=database ) - conn_args = { - **connection_options - } engine = create_engine( - conn_str, connect_args=conn_args, pool_pre_ping=True + conn_str, connect_args=connect_args, pool_pre_ping=True ) return engine From 75cac2348e3b02c21d0de7af558b0e8ab9b120b0 Mon Sep 17 00:00:00 2001 From: Benjamin Webb Date: Mon, 16 Jun 2025 21:04:01 -0400 Subject: [PATCH 13/15] Use geoalchemy2 functions --- pygeoapi/provider/mvt_postgresql.py | 22 ++++++++++++---------- 1 file changed, 12 insertions(+), 10 deletions(-) diff --git a/pygeoapi/provider/mvt_postgresql.py b/pygeoapi/provider/mvt_postgresql.py index 2c925e3b4..7b08fc35b 100644 --- a/pygeoapi/provider/mvt_postgresql.py +++ b/pygeoapi/provider/mvt_postgresql.py @@ -35,8 +35,12 @@ import logging -from sqlalchemy.sql import func, select +from geoalchemy2.functions import (ST_TileEnvelope, ST_Transform, ST_AsMVTGeom, + ST_AsMVT, ST_CurveToLine, ST_MakeEnvelope) + +from sqlalchemy.sql import select from sqlalchemy.orm import Session + from pygeoapi.models.provider.base import ( TileSetMetadata, TileMatrixSetEnum, LinkType) from pygeoapi.provider.base import ProviderConnectionError @@ -134,8 +138,7 @@ def get_tiles(self, layer='default', tileset=None, :returns: an encoded mvt tile """ - if format_ == 'mvt': - format_ = self.format_type + z, y, x = map(int, [z, y, x]) [tileset_schema] = [ schema for schema in self.get_tiling_schemes() @@ -151,13 +154,13 @@ def get_tiles(self, layer='default', tileset=None, geom_column = getattr(self.table_model, self.geom) geom_filter = geom_column.intersects( - func.ST_Transform(envelope, storage_srid) + ST_Transform(envelope, storage_srid) ) mvtgeom = ( - func.ST_AsMVTGeom( - func.ST_Transform(func.ST_CurveToLine(geom_column), out_srid), - func.ST_Transform(envelope, out_srid)) + ST_AsMVTGeom( + ST_Transform(ST_CurveToLine(geom_column), out_srid), + ST_Transform(envelope, out_srid)) .label('mvtgeom') ) @@ -165,11 +168,10 @@ def get_tiles(self, layer='default', tileset=None, select(mvtgeom, *self.fields.values()) .filter(geom_filter) .cte('mvtrow') - .table_valued() ) mvtquery = select( - func.ST_AsMVT(mvtrow, layer) + ST_AsMVT(mvtrow.table_valued(), layer) ) with Session(self._engine) as session: @@ -184,7 +186,7 @@ def get_html_metadata(self, dataset, server_url, layer, tileset, service_url = url_join( server_url, - f'collections/{dataset}/tiles/{tileset}' + f'collections/{dataset}/tiles/{tileset}', '{tileMatrix}/{tileRow}/{tileCol}?f=mvt') metadata_url = url_join( server_url, From 8b7d44f4666ebf98e13083ab9626f0f9d597d580 Mon Sep 17 00:00:00 2001 From: Benjamin Webb Date: Mon, 16 Jun 2025 21:05:16 -0400 Subject: [PATCH 14/15] Amend WorldCRS84Quad Tile bbox Align with WorldCRS84Quad spec using https://github.com/geopython/pygeoapi/pull/2042 --- pygeoapi/provider/mvt_postgresql.py | 37 ++++++++++++++++++++++++++- tests/test_postgresql_mvt_provider.py | 2 +- 2 files changed, 37 insertions(+), 2 deletions(-) diff --git a/pygeoapi/provider/mvt_postgresql.py b/pygeoapi/provider/mvt_postgresql.py index 7b08fc35b..aa791ec45 100644 --- a/pygeoapi/provider/mvt_postgresql.py +++ b/pygeoapi/provider/mvt_postgresql.py @@ -150,7 +150,7 @@ def get_tiles(self, layer='default', tileset=None, storage_srid = get_crs_from_uri(self.storage_crs).to_string() out_srid = get_crs_from_uri(tileset_schema.crs).to_string() - envelope = func.ST_TileEnvelope(z, x, y).label('bounds') + envelope = self.get_envelope(z, y, x, tileset) geom_column = getattr(self.table_model, self.geom) geom_filter = geom_column.intersects( @@ -254,5 +254,40 @@ def get_default_metadata(self, dataset, server_url, layer, tileset, return content.model_dump(exclude_none=True) + @staticmethod + def get_envelope(z, y, x, tileset): + """ + Calculate the Tile bounding box of a tile at zoom z, y, x. + + WorldCRS84Quad tiles have: + - origin top-left (y=0 is north) + - full lon: -180 to 180 + + :param tileset: mvt tileset + :param z: z index + :param y: y index + :param x: x index + + :returns: SQL Alchemy Tile Envelope + """ + + if tileset == TileMatrixSetEnum.WORLDCRS84QUAD.value.tileMatrixSet: + + tile_size = 180 / 2 ** z + + xmin = tile_size * x - 180 + ymax = tile_size * -y + 90 + + # getting bottom-right coordinates of the tile + xmax = xmin + tile_size + ymin = ymax - tile_size + + envelope = ST_MakeEnvelope(xmin, ymin, xmax, ymax, 4326) + + else: + envelope = ST_TileEnvelope(z, x, y) + + return envelope.label('bounds') + def __repr__(self): return f' {self.data}' diff --git a/tests/test_postgresql_mvt_provider.py b/tests/test_postgresql_mvt_provider.py index ee30caeab..172b4b1e5 100644 --- a/tests/test_postgresql_mvt_provider.py +++ b/tests/test_postgresql_mvt_provider.py @@ -276,7 +276,7 @@ def test_get_tiles_WorldCRS84Quad(config): assert tile is None # Valid tile, content - z, x, y = 10, 595, 521 + z, x, y = 9, 595, 267 tile = p.get_tiles( tileset=tileset, z=z, x=x, y=y, From 3384f69c3af5cd344bdf257367fe734e4833d94c Mon Sep 17 00:00:00 2001 From: Benjamin Webb Date: Tue, 17 Jun 2025 13:07:03 -0400 Subject: [PATCH 15/15] Support WorldCRS84Quad in default tile html Add minimal support stretching EPSG3857 TileLayer to EPSG4326. Really should be a CRS84 based layer like: L.tileLayer.wms("https://gibs.earthdata.nasa.gov/wms/epsg4326/best/wms.cgi", { layers: "BlueMarble_ShadedRelief_Bathymetry", attribution: "NASA GIBS", }) .addTo(map); --- .../templates/collections/tiles/index.html | 22 +++++++++++++++++++ 1 file changed, 22 insertions(+) diff --git a/pygeoapi/templates/collections/tiles/index.html b/pygeoapi/templates/collections/tiles/index.html index 75e1f9d36..25e8248d7 100644 --- a/pygeoapi/templates/collections/tiles/index.html +++ b/pygeoapi/templates/collections/tiles/index.html @@ -69,6 +69,7 @@

Tiles

{% block extrafoot %}