extend point data extraction to also work vor mbtile layers

Author: hambsch

server/environment.yml

@@ -27,4 +27,5 @@ dependencies:
   - rio-tiler
   - prometheus-fastapi-instrumentator
   - pygeofilter
+  - mapbox-vector-tile
 prefix: /opt/mapbase3/server/conda-env

server/guppy/db/schemas.py

@@ -36,7 +36,7 @@ class Config:
 class PointResponse(CamelModel):
     type: str
     layer_name: str
-    value: Optional[float] = None
+    value: Optional[float | dict] = None
 
 
 class StatsResponse(CamelModel):

server/guppy/endpoints/endpoints.py

@@ -26,6 +26,7 @@
 from guppy.db import schemas as s
 from guppy.endpoints.endpoint_utils import get_overview_factor, create_stats_response, _extract_area_from_dataset, _extract_shape_mask_from_dataset, _decode, sample_coordinates_window, \
     create_quantile_response
+from guppy.endpoints.tile_utils import latlon_to_tilexy, get_tile_data, pbf_to_geodataframe
 
 logger = logging.getLogger(__name__)
 
@@ -292,10 +293,10 @@ def get_multi_line_data_list_for_wkt(db: Session, body: s.MultiLineGeometryListB
     return Response(status_code=status.HTTP_404_NOT_FOUND)
 
 
-def get_point_value_from_raster(db: Session, layer_name: str, x: float, y: float):
+def get_point_value_from_layer(db: Session, layer_name: str, x: float, y: float):
     t = time.time()
     layer_model = db.query(m.LayerMetadata).filter_by(layer_name=layer_name).first()
-    if layer_model:
+    if layer_model and not layer_model.is_mbtile:
         path = layer_model.file_path
         if os.path.exists(path) and x and y:
             with rasterio.open(path) as src:
@@ -307,6 +308,19 @@ def get_point_value_from_raster(db: Session, layer_name: str, x: float, y: float
                     logger.info(f'get_point_value_from_raster 200 {time.time() - t}')
                     return s.PointResponse(type='point value', layer_name=layer_name, value=None if math.isclose(float(v[0]), nodata) else float(v[0]))
         logger.warning(f'file not found {path}')
+    elif layer_model and layer_model.is_mbtile:
+        path = layer_model.file_path
+        if os.path.exists(path) and x and y:
+            tile_z, tile_x, tile_y = latlon_to_tilexy(x, y, 14)
+            tile = get_tile_data(layer_name=layer_name, mb_file=path, z=tile_z, x=tile_x, y=tile_y)
+            tile_df = pbf_to_geodataframe(tile, tile_x, tile_y, tile_z)
+            # get the value of the point
+            point = wkt.loads(f'POINT ({x} {y})')
+            values = tile_df[tile_df.intersects(point)].drop(columns=['geometry'])
+            result = {'type': 'point value', 'layer_name': layer_name, 'value': values.to_dict(orient='records')[0]}
+            logger.info(f'get_point_value_from_raster 200 {time.time() - t}')
+            return result
+        logger.warning(f'file not found {path}')
     logger.info(f'get_point_value_from_raster 204 {time.time() - t}')
     return Response(status_code=status.HTTP_204_NO_CONTENT)
 

server/guppy/endpoints/tile_utils.py

@@ -3,12 +3,21 @@
 import time
 from threading import Lock
 
+import geopandas as gpd
+import mapbox_vector_tile
+import mercantile
 import numpy as np
+from shapely.geometry import shape
+from shapely.ops import transform
 
 from guppy.db.db_session import SessionLocal
 from guppy.db.models import TileStatistics
+from guppy.error import create_error
 
 logger = logging.getLogger(__name__)
+from typing import Optional
+import gzip
+import sqlite3
 
 
 def tile2lonlat(x, y, z):
@@ -33,6 +42,14 @@ def tile2lonlat(x, y, z):
     return (lon_left, -lat_bottom, lon_right, -lat_top)
 
 
+def latlon_to_tilexy(lon, lat, z):
+    lat_rad = math.radians(lat)
+    n = 2.0 ** z
+    xtile = int((lon + 180.0) / 360.0 * n)
+    ytile = int((1.0 - math.asinh(math.tan(lat_rad)) / math.pi) / 2.0 * n)
+    return z, xtile, ytile
+
+
 # In-memory request counter, structured by layer_name -> z/x/y -> count
 request_counter = {}
 request_counter_lock = Lock()
@@ -163,3 +180,70 @@ def get_field_mapping(conn):
     cursor.execute("PRAGMA table_info(tiles)")
     columns = cursor.fetchall()
     return {col[1]: col[1] for col in columns}  # Simple mapping of name to name
+
+
+def get_tile_data(layer_name: str, mb_file: str, z: int, x: int, y: int) -> Optional[bytes]:
+    """
+    Args:
+        layer_name: The name of the layer for which the tile data is being retrieved.
+        mb_file: The path to the MBTiles file from which the tile data is being retrieved.
+        z: The zoom level of the tile.
+        x: The X coordinate of the tile.
+        y: The Y coordinate of the tile.
+
+    Returns:
+        Optional[bytes]: The tile data as bytes if found, or None if no tile data exists for the given parameters.
+
+    Raises:
+        HTTPException: If there is an error retrieving the tile data from the MBTiles file.
+
+    """
+    # Flip Y coordinate because MBTiles grid is TMS (bottom-left origin)
+    y = (1 << z) - 1 - y
+    logger.info(f"Getting tile for layer {layer_name} at zoom {z}, x {x}, y {y}")
+    try:
+        uri = f'file:{mb_file}?mode=ro'
+        with sqlite3.connect(uri, uri=True) as conn:
+            cursor = conn.cursor()
+            cursor.execute("SELECT tile_data FROM tiles WHERE zoom_level=? AND tile_column=? AND tile_row=?", (z, x, y))
+            tile = cursor.fetchone()
+            if tile:
+                return gzip.decompress(bytes(tile[0]))
+            else:
+                return None
+    except Exception as e:
+        create_error(code=404, message=str(e))
+
+
+def pbf_to_geodataframe(pbf_data, x, y, z):
+    """
+    Converts PBF data to a GeoDataFrame.
+
+    :param pbf_data: The PBF data to be decoded.
+    :param x: The x-coordinate of the tile.
+    :param y: The y-coordinate of the tile.
+    :param z: The zoom level of the tile.
+    :return: A GeoDataFrame containing the decoded PBF data in GeoJSON format.
+    """
+    # Decode PBF data
+    decoded_data = mapbox_vector_tile.decode(pbf_data)
+    tile_bounds = mercantile.bounds(x, y, z)
+    # Collect features and convert them to GeoJSON format
+    features = []
+    for layer_name, layer in decoded_data.items():
+        for feature in layer['features']:
+            geom = shape(feature['geometry'])
+
+            def scale_translate(x, y, bounds=tile_bounds, tile_dim=4096):
+                # Adjust for the flipped tiles by inverting the y-axis calculation
+                lon = (x / tile_dim) * (bounds.east - bounds.west) + bounds.west
+                lat = (y / tile_dim) * (bounds.north - bounds.south) + bounds.south
+                return lon, lat
+
+            geom_transformed = transform(scale_translate, geom)
+            properties = feature['properties']
+            properties["featureId"] = feature['id']
+            features.append({'type': 'Feature', 'geometry': geom_transformed, 'properties': properties})
+
+    gdf = gpd.GeoDataFrame.from_features(features, crs='EPSG:4326')
+    return gdf

server/guppy/routes/data_router.py

@@ -2,10 +2,10 @@
 from fastapi.responses import ORJSONResponse
 from sqlalchemy.orm import Session
 
-from guppy.endpoints import endpoints
 from guppy.config import config as cfg
 from guppy.db import schemas as s
 from guppy.db.dependencies import get_db
+from guppy.endpoints import endpoints
 
 router = APIRouter(
     prefix=f"{cfg.deploy.path}/layers",
@@ -38,9 +38,9 @@ def get_multi_line_data_list_for_wkt(body: s.MultiLineGeometryListBody, db: Sess
     return endpoints.get_multi_line_data_list_for_wkt(db=db, body=body)
 
 
-@router.get("/{layer_name}/point", response_model=s.PointResponse, description="Get point value for a given coordinate from raster within a layer.")
-def get_point_value_from_raster(layer_name: str, x: float, y: float, db: Session = Depends(get_db)):
-    return endpoints.get_point_value_from_raster(db=db, layer_name=layer_name, x=x, y=y)
+@router.get("/{layer_name}/point", response_model=s.PointResponse, description="Get point value for a given coordinate (in 4326) from a layer.")
+def get_point_value_from_layer(layer_name: str, x: float, y: float, db: Session = Depends(get_db)):
+    return endpoints.get_point_value_from_layer(db=db, layer_name=layer_name, x=x, y=y)
 
 
 @router.post("/{layer_name}/object", response_class=ORJSONResponse, description="Get object list for a given line wkt geometry within a layer.")