Store subdivided areas.

This creates a table containing the same areas as the Geometry table,
but split using ST_Subdivide. This is updated on save, and is used for
postcode/point lookups, as well as geometry intersection.

mapit/migrations/0005_geometrysubdivided.py

@@ -0,0 +1,27 @@
+# Generated by Django 3.2.18 on 2023-05-23 09:21
+
+import django.contrib.gis.db.models.fields
+from django.db import migrations, models
+import django.db.models.deletion
+from django.conf import settings
+
+
+class Migration(migrations.Migration):
+
+    dependencies = [
+        ('mapit', '0004_add_country_m2m'),
+    ]
+
+    operations = [
+        migrations.CreateModel(
+            name='GeometrySubdivided',
+            fields=[
+                ('id', models.AutoField(auto_created=True, primary_key=True, serialize=False, verbose_name='ID')),
+                ('division', django.contrib.gis.db.models.fields.PolygonField(srid=settings.MAPIT_AREA_SRID)),
+                ('geometry', models.ForeignKey(on_delete=django.db.models.deletion.CASCADE, related_name='subdivided', to='mapit.geometry')),
+            ],
+            options={
+                'verbose_name_plural': 'subdivided geometries',
+            },
+        ),
+    ]

mapit/migrations/0006_auto_20230523_0940.py

@@ -0,0 +1,18 @@
+# Generated by Django 3.2.18 on 2023-05-23 09:40
+
+from django.db import migrations
+
+
+class Migration(migrations.Migration):
+
+    dependencies = [
+        ('mapit', '0005_geometrysubdivided'),
+    ]
+
+    operations = [
+        migrations.RunSQL(
+            sql="""INSERT INTO mapit_geometrysubdivided (geometry_id, division)
+                SELECT id,ST_Subdivide(polygon) FROM mapit_geometry""",
+            reverse_sql="TRUNCATE mapit_geometrysubdivided RESTART IDENTITY"
+        ),
+    ]

mapit/models.py

@@ -6,6 +6,7 @@
 from django.contrib.gis.db import models
 from django.conf import settings
 from django.db import connection
+from django.db.models import F, Func
 from django.db.models.query import RawQuerySet
 from django.utils.encoding import smart_str
 from django.utils.functional import cached_property
@@ -144,7 +145,7 @@ def by_location(self, location, generation=None):
         if not location:
             return []
         return Area.objects.filter(
-            polygons__polygon__contains=location,
+            polygons__subdivided__division__contains=location,
             generation_low__lte=generation, generation_high__gte=generation
         )
 
@@ -173,25 +174,25 @@ def intersect(self, query_type, area, types, generation):
         else:
             query_area_type = ''
 
-        query_geo = ' OR '.join(['ST_%s(geometry.polygon, target.polygon)' % type for type in query_type])
+        query_geo = ' OR '.join(['ST_%s(geometry_sd.division, target_sd.division)' % type for type in query_type])
 
         query = '''
-WITH
-    target AS ( SELECT ST_collect(polygon) polygon FROM mapit_geometry WHERE area_id=%%s ),
-    geometry AS %s (
-        SELECT mapit_geometry.*
-          FROM mapit_geometry, mapit_area, target
-         WHERE mapit_geometry.area_id = mapit_area.id
-               AND mapit_geometry.polygon && target.polygon
-               AND mapit_area.id != %%s
-               AND mapit_area.generation_low_id <= %%s
-               AND mapit_area.generation_high_id >= %%s
-               %s
-    )
 SELECT DISTINCT mapit_area.*
-  FROM mapit_area, geometry, target
- WHERE geometry.area_id = mapit_area.id AND (%s)
-''' % (materialized(), query_area_type, query_geo)
+FROM
+    mapit_area,
+    mapit_geometry geometry, mapit_geometrysubdivided geometry_sd,
+    mapit_geometry target, mapit_geometrysubdivided target_sd
+WHERE
+    geometry_sd.geometry_id = geometry.id
+    AND target_sd.geometry_id = target.id
+    AND geometry.area_id = mapit_area.id
+    AND target.area_id = %%s
+    AND geometry.area_id != %%s
+    AND mapit_area.generation_low_id <= %%s
+    AND mapit_area.generation_high_id >= %%s
+    %s
+    AND (%s)
+''' % (query_area_type, query_geo)
         return RawQuerySet(raw_query=query, model=self.model, params=params, using=self._db)
 
     def get_or_create_with_name(self, country=None, type=None, name_type='', name=''):
@@ -331,7 +332,24 @@ class Meta:
         verbose_name_plural = 'geometries'
 
     def __str__(self):
-        return '%s, polygon %d' % (smart_str(self.area), self.id)
+        return '%s, polygon %s' % (smart_str(self.area), self.id)
+
+    def save(self, *args, **kwargs):
+        super().save(*args, **kwargs)
+        divided = Geometry.objects.filter(id=self.id).annotate(division=Func(F('polygon'), function='ST_Subdivide'))
+        GeometrySubdivided.objects.filter(geometry=self).delete()
+        GeometrySubdivided.objects.bulk_create(GeometrySubdivided(geometry=self, division=p.division) for p in divided)
+
+
+class GeometrySubdivided(models.Model):
+    geometry = models.ForeignKey(Geometry, related_name='subdivided', on_delete=models.CASCADE)
+    division = models.PolygonField(srid=settings.MAPIT_AREA_SRID)
+
+    class Meta:
+        verbose_name_plural = 'subdivided geometries'
+
+    def __str__(self):
+        return '%s, subdivision %s' % (smart_str(self.geometry), self.id)
 
 
 class NameType(models.Model):
@@ -412,12 +430,15 @@ class PostcodeQuerySet(models.QuerySet):
     def filter_by_area(self, area, limit=''):
         if limit:
             limit = 'LIMIT %s' % limit
-        polygons = 'SELECT ST_Transform(polygon, 4326) AS polygon from mapit_geometry where area_id=%s'
+        polygons = '''SELECT ST_Transform(division, 4326) AS division
+            FROM mapit_geometrysubdivided
+            JOIN mapit_geometry ON geometry_id = mapit_geometry.id
+            WHERE area_id = %s'''
         query = '''
 WITH target AS %s ( %s )
 SELECT "mapit_postcode"."id", "mapit_postcode"."postcode", "mapit_postcode"."location"::bytea
   FROM mapit_postcode, target
- WHERE ST_CoveredBy(location, target.polygon)
+ WHERE ST_CoveredBy(location, target.division)
  %s
 ''' % (materialized(), polygons, limit)
         return self.raw(query, params=[area.id])

mapit/views/areas.py

@@ -405,11 +405,11 @@ def areas_by_point(request, srid, x, y, bb=False, format=''):
         raise ViewException(format, _('Bad "within" parameter specified'), 400)
 
     if method == 'box':
-        q &= Q(polygons__polygon__bbcontains=location)
+        q &= Q(polygons__subdivided__division__bbcontains=location)
     elif within:
-        q &= Q(polygons__polygon__dwithin=(location, within))
+        q &= Q(polygons__subdivided__division__dwithin=(location, within))
     else:
-        q &= Q(polygons__polygon__contains=location)
+        q &= Q(polygons__subdivided__division__contains=location)
     areas = Area.objects.filter(q).distinct()
 
     return output_areas(