Merge pull request #43 from OpenGeoscience/postgis-vector-tiles

Generate vector tiles dynamically from extracted features

pyproject.toml

@@ -19,3 +19,7 @@ target-version = "py310"
 
 [tool.ruff.format]
 quote-style = "single"
+
+[tool.ruff.lint.isort]
+force-sort-within-sections = true # Sort by name, don't cluster "from" vs "import"
+combine-as-imports = true # Combines "as" imports on the same line

uvdat/core/admin.py

@@ -11,8 +11,8 @@
     RasterMapLayer,
     SimulationResult,
     SourceRegion,
+    VectorFeature,
     VectorMapLayer,
-    VectorTile,
 )
 
 
@@ -53,14 +53,8 @@ def get_dataset_name(self, obj):
         return obj.file_item.dataset.name
 
 
-class VectorTileAdmin(admin.ModelAdmin):
-    list_display = ['id', 'get_fileitem_name', 'get_map_layer_index', 'x', 'y', 'z']
-
-    def get_fileitem_name(self, obj):
-        return obj.map_layer.file_item.name
-
-    def get_map_layer_index(self, obj):
-        return obj.map_layer.index
+class VectorFeatureAdmin(admin.ModelAdmin):
+    list_display = ['id', 'map_layer']
 
 
 class SourceRegionAdmin(admin.ModelAdmin):
@@ -107,7 +101,7 @@ class SimulationResultAdmin(admin.ModelAdmin):
 admin.site.register(Chart, ChartAdmin)
 admin.site.register(RasterMapLayer, RasterMapLayerAdmin)
 admin.site.register(VectorMapLayer, VectorMapLayerAdmin)
-admin.site.register(VectorTile, VectorTileAdmin)
+admin.site.register(VectorFeature, VectorFeatureAdmin)
 admin.site.register(SourceRegion, SourceRegionAdmin)
 admin.site.register(DerivedRegion, DerivedRegionAdmin)
 admin.site.register(NetworkNode, NetworkNodeAdmin)

uvdat/core/migrations/0002_vector_features.py

@@ -0,0 +1,36 @@
+# Generated by Django 4.1 on 2024-06-24 18:24
+
+import django.contrib.gis.db.models.fields
+from django.db import migrations, models
+import django.db.models.deletion
+
+
+class Migration(migrations.Migration):
+    dependencies = [
+        ('core', '0001_models_redesign'),
+    ]
+
+    operations = [
+        migrations.CreateModel(
+            name='VectorFeature',
+            fields=[
+                (
+                    'id',
+                    models.BigAutoField(
+                        auto_created=True, primary_key=True, serialize=False, verbose_name='ID'
+                    ),
+                ),
+                ('geometry', django.contrib.gis.db.models.fields.GeometryField(srid=4326)),
+                ('properties', models.JSONField()),
+                (
+                    'map_layer',
+                    models.ForeignKey(
+                        on_delete=django.db.models.deletion.CASCADE, to='core.vectormaplayer'
+                    ),
+                ),
+            ],
+        ),
+        migrations.DeleteModel(
+            name='VectorTile',
+        ),
+    ]

uvdat/core/models/__init__.py

@@ -2,7 +2,7 @@
 from .context import Context
 from .dataset import Dataset
 from .file_item import FileItem
-from .map_layers import RasterMapLayer, VectorMapLayer, VectorTile
+from .map_layers import RasterMapLayer, VectorFeature, VectorMapLayer
 from .networks import NetworkEdge, NetworkNode
 from .regions import DerivedRegion, SourceRegion
 from .simulations import SimulationResult
@@ -14,7 +14,7 @@
     FileItem,
     RasterMapLayer,
     VectorMapLayer,
-    VectorTile,
+    VectorFeature,
     SourceRegion,
     DerivedRegion,
     NetworkEdge,

uvdat/core/models/dataset.py

@@ -89,43 +89,3 @@ def get_map_layers(self):
             return VectorMapLayer.objects.filter(file_item__dataset=self)
 
         raise NotImplementedError(f'Dataset Type {self.dataset_type}')
-
-    def get_map_layer_tile_extents(self):
-        """
-        Return the extents of all vector map layers of this dataset.
-
-        Returns `None` if the dataset is not a vector dataset.
-        """
-        if self.dataset_type != self.DatasetType.VECTOR:
-            return None
-
-        from uvdat.core.models import VectorMapLayer, VectorTile
-
-        # Retrieve all layers
-        layer_ids = VectorMapLayer.objects.filter(file_item__dataset=self).values_list(
-            'id', flat=True
-        )
-
-        # Return x/y extents by layer id and z depth
-        vals = (
-            VectorTile.objects.filter(map_layer_id__in=layer_ids)
-            .values('map_layer_id', 'z')
-            .annotate(
-                min_x=models.Min('x'),
-                min_y=models.Min('y'),
-                max_x=models.Max('x'),
-                max_y=models.Max('y'),
-            )
-            .order_by('map_layer_id')
-        )
-
-        # Deconstruct query into response format
-        layers = {}
-        for entry in vals:
-            map_layer_id = entry.pop('map_layer_id')
-            if map_layer_id not in layers:
-                layers[map_layer_id] = {}
-
-            layers[map_layer_id][entry.pop('z')] = entry
-
-        return layers

uvdat/core/models/map_layers.py

@@ -2,6 +2,7 @@
 from pathlib import Path
 import tempfile
 
+from django.contrib.gis.db import models as geomodels
 from django.core.files.base import ContentFile
 from django.db import models
 from django_extensions.db.models import TimeStampedModel
@@ -58,41 +59,8 @@ def read_geojson_data(self) -> dict:
         """Read and load the data from geojson_file into a dict."""
         return json.load(self.geojson_file.open())
 
-    def get_tile_extents(self):
-        """Return a dict that maps z tile values to the x/y extent at that depth."""
-        return {
-            entry.pop('z'): entry
-            for entry in (
-                VectorTile.objects.filter(map_layer=self)
-                .values('z')
-                .annotate(
-                    min_x=models.Min('x'),
-                    min_y=models.Min('y'),
-                    max_x=models.Max('x'),
-                    max_y=models.Max('y'),
-                )
-                .order_by()
-            )
-        }
-
-
-class VectorTile(models.Model):
-    EMPTY_TILE_DATA = {
-        'type': 'FeatureCollection',
-        'features': [],
-    }
 
+class VectorFeature(models.Model):
     map_layer = models.ForeignKey(VectorMapLayer, on_delete=models.CASCADE)
-    geojson_data = models.JSONField(blank=True, null=True)
-    x = models.IntegerField(default=0)
-    y = models.IntegerField(default=0)
-    z = models.IntegerField(default=0)
-
-    class Meta:
-        constraints = [
-            # Ensure that a full index only ever resolves to one record
-            models.UniqueConstraint(
-                name='unique-map-layer-index', fields=['map_layer', 'z', 'x', 'y']
-            )
-        ]
-        indexes = [models.Index(fields=('z', 'x', 'y'), name='vectortile-coordinates-index')]
+    geometry = geomodels.GeometryField()
+    properties = models.JSONField()

uvdat/core/rest/dataset.py

@@ -29,13 +29,8 @@ def map_layers(self, request, **kwargs):
         if dataset.dataset_type == Dataset.DatasetType.RASTER:
             serializer = uvdat_serializers.RasterMapLayerSerializer(map_layers, many=True)
         elif dataset.dataset_type == Dataset.DatasetType.VECTOR:
-            # Inject tile extents
-            extents = dataset.get_map_layer_tile_extents()
-            for layer in map_layers:
-                layer.tile_extents = extents.pop(layer.id)
-
             # Set serializer
-            serializer = uvdat_serializers.ExtendedVectorMapLayerSerializer(map_layers, many=True)
+            serializer = uvdat_serializers.VectorMapLayerSerializer(map_layers, many=True)
         else:
             raise NotImplementedError(f'Dataset Type {dataset.dataset_type}')
 

uvdat/core/rest/map_layers.py

@@ -1,19 +1,77 @@
 import json
 
+from django.db import connection
 from django.http import HttpResponse
 from django_large_image.rest import LargeImageFileDetailMixin
 from rest_framework.decorators import action
 from rest_framework.response import Response
 from rest_framework.viewsets import ModelViewSet
 
 from uvdat.core.models import RasterMapLayer, VectorMapLayer
-from uvdat.core.models.map_layers import VectorTile
 from uvdat.core.rest.serializers import (
     RasterMapLayerSerializer,
     VectorMapLayerDetailSerializer,
     VectorMapLayerSerializer,
 )
 
+VECTOR_TILE_SQL = """
+WITH
+tilenv as (
+    SELECT ST_TRANSFORM(ST_TileEnvelope(%(z)s, %(x)s, %(y)s), %(srid)s) as te
+),
+tilenvbounds as (
+    SELECT
+        ST_XMin(te) as xmin,
+        ST_YMin(te) as ymin,
+        ST_XMax(te) as xmax,
+        ST_YMax(te) as ymax,
+        (ST_XMax(te) - ST_XMin(te)) / 4 as segsize
+    FROM tilenv
+),
+env as (
+    SELECT ST_Segmentize(
+        ST_MakeEnvelope(
+            xmin,
+            ymin,
+            xmax,
+            ymax,
+            %(srid)s
+        ),
+        segsize
+    ) as seg
+    FROM tilenvbounds
+),
+bounds as (
+    SELECT
+        seg as geom,
+        seg::box2d as b2d
+    FROM env
+),
+mvtgeom as (
+    SELECT
+        ST_AsMVTGeom(
+            ST_Transform(t.geometry, %(srid)s),
+            bounds.b2d
+        ) AS geom,
+    t.properties as properties
+    FROM
+        core_vectorfeature t,
+        bounds
+    WHERE
+        t.map_layer_id = %(map_layer_id)s
+        AND ST_Intersects(
+            ST_Transform(t.geometry, %(srid)s),
+            ST_Transform(bounds.geom, %(srid)s)
+        )
+        AND (
+            ST_GeometryType(ST_AsText(t.geometry)) != 'ST_Point'
+            OR %(z)s >= 16
+        )
+)
+SELECT ST_AsMVT(mvtgeom.*) FROM mvtgeom
+;
+"""
+
 
 class RasterMapLayerViewSet(ModelViewSet, LargeImageFileDetailMixin):
     queryset = RasterMapLayer.objects.select_related('file_item__dataset').all()
@@ -48,10 +106,22 @@ def retrieve(self, request, *args, **kwargs):
         url_name='tiles',
     )
     def get_vector_tile(self, request, x: str, y: str, z: str, pk: str):
-        # Return vector tile or empty tile
-        try:
-            tile = VectorTile.objects.get(map_layer_id=pk, x=x, y=y, z=z)
-        except VectorTile.DoesNotExist:
-            return Response(VectorTile.EMPTY_TILE_DATA, status=200)
+        with connection.cursor() as cursor:
+            cursor.execute(
+                VECTOR_TILE_SQL,
+                {
+                    'z': z,
+                    'x': x,
+                    'y': y,
+                    'srid': 3857,
+                    'map_layer_id': pk,
+                },
+            )
+            row = cursor.fetchone()
 
-        return Response(tile.geojson_data, status=200)
+        tile = row[0]
+        return HttpResponse(
+            tile,
+            content_type='application/octet-stream',
+            status=200 if tile else 204,
+        )

uvdat/core/rest/serializers.py

@@ -85,14 +85,6 @@ class Meta:
         fields = '__all__'
 
 
-class ExtendedVectorMapLayerSerializer(serializers.ModelSerializer, AbstractMapLayerSerializer):
-    tile_extents = serializers.JSONField()
-
-    class Meta:
-        model = VectorMapLayer
-        exclude = ['geojson_file']
-
-
 class VectorMapLayerSerializer(serializers.ModelSerializer, AbstractMapLayerSerializer):
     class Meta:
         model = VectorMapLayer
@@ -101,17 +93,13 @@ class Meta:
 
 class VectorMapLayerDetailSerializer(serializers.ModelSerializer, AbstractMapLayerSerializer):
     derived_region_id = serializers.SerializerMethodField('get_derived_region_id')
-    tile_extents = serializers.SerializerMethodField('get_tile_extents')
 
     def get_derived_region_id(self, obj):
         dr = obj.derivedregion_set.first()
         if dr is None:
             return None
         return dr.id
 
-    def get_tile_extents(self, obj: VectorMapLayer):
-        return obj.get_tile_extents()
-
     class Meta:
         model = VectorMapLayer
         exclude = ['geojson_file']

uvdat/core/tasks/dataset.py

@@ -9,8 +9,9 @@
     SourceRegion,
     VectorMapLayer,
 )
+from uvdat.core.tasks.map_layers import save_vector_features
 
-from .map_layers import create_raster_map_layer, create_vector_map_layer, save_vector_tiles
+from .map_layers import create_raster_map_layer, create_vector_map_layer
 from .networks import create_network
 from .regions import create_source_regions
 
@@ -58,7 +59,7 @@ def convert_dataset(
 
             # Create vector tiles after geojson_data may have
             # been altered by create_network or create_source_regions
-            save_vector_tiles(vector_map_layer)
+            save_vector_features(vector_map_layer=vector_map_layer)
 
     dataset.processing = False
     dataset.save()