Add DerivedRegions from Region Union/Intersection

sample_data/datasets.json

@@ -110,13 +110,13 @@
         }
     },
     {
-        "name": "Boston Census Blocks",
-        "description": "Boundaries of 2020 Census Blocks in the Boston Metropolitan area",
+        "name": "Boston Neighborhoods",
+        "description": "The Neighborhood boundaries data layer is a combination of zoning neighborhood boundaries, zip code boundaries and 2010 Census tract boundaries",
         "category": "region",
         "city": "Boston, MA",
         "raw_data_type": "geojson",
         "url": "https://data.kitware.com/api/v1/item/64caa3da77edef4e1ea8ef57/download",
-        "path": "boston/census2020.json",
+        "path": "boston/neighborhoods2020.json",
         "style": {
             "property_map": {
                 "name": "BlockGr202"
@@ -129,6 +129,43 @@
             }
         }
     },
+    {
+        "name": "Boston Census 2020 Block Groups",
+        "description": "Block groups (between 600 and 3,000 people per block) for Boston from the 2020 census",
+        "category": "region",
+        "city": "Boston, MA",
+        "raw_data_type": "shape_file_archive",
+        "url": "https://data.boston.gov/dataset/c478b600-3e3e-46fd-9f57-da89459e9928/resource/11282722-9386-4272-8a82-2fcec89e6d55/download/census2020_blockgroups.zip",
+        "path": "mass/blockgroups.json",
+        "style": {
+            "property_map": {
+                "name": "GEOID20"
+            },
+            "options": {
+                "outline": "white",
+                "palette": [
+                    "green"
+                ]
+            }
+        }
+    },
+    {
+        "name": "Boston Zip Codes",
+        "description": "Zip codes 01001-02791",
+        "category": "region",
+        "city": "Boston, MA",
+        "raw_data_type": "shape_file_archive",
+        "url": "https://data.kitware.com/api/v1/item/64fbb4c6e99e9e6006f00114/download",
+        "path": "boston/zipcodes.json",
+        "style": {
+            "options": {
+                "outline": "white",
+                "palette": [
+                    "grey"
+                ]
+            }
+        }
+    },
     {
         "name": "9-inch Sea Level Rise",
         "description": "From Analyze Boston; 9 inches of sea level rise is expected across emissions scenarios as likely to occur as early as the 2030s.",

uvdat/core/migrations/0009_derived_regions_and_region_constraints.py

@@ -0,0 +1,48 @@
+# Generated by Django 4.1 on 2023-09-27 19:45
+
+import django.contrib.gis.db.models.fields
+from django.db import migrations, models
+import django.db.models.deletion
+
+
+class Migration(migrations.Migration):
+
+    dependencies = [
+        ('core', '0008_simulation_results'),
+    ]
+
+    operations = [
+        migrations.CreateModel(
+            name='DerivedRegion',
+            fields=[
+                ('id', models.BigAutoField(auto_created=True, primary_key=True, serialize=False, verbose_name='ID')),
+                ('name', models.CharField(max_length=255)),
+                ('properties', models.JSONField(blank=True, null=True)),
+                ('boundary', django.contrib.gis.db.models.fields.MultiPolygonField(srid=4326)),
+                ('source_operation', models.CharField(choices=[('UNION', 'Union'), ('INTERSECTION', 'Intersection')], max_length=12)),
+            ],
+        ),
+        migrations.AlterField(
+            model_name='region',
+            name='name',
+            field=models.CharField(max_length=255),
+        ),
+        migrations.AddConstraint(
+            model_name='region',
+            constraint=models.UniqueConstraint(fields=('dataset', 'name'), name='unique-name-per-dataset'),
+        ),
+        migrations.AddField(
+            model_name='derivedregion',
+            name='city',
+            field=models.ForeignKey(on_delete=django.db.models.deletion.CASCADE, related_name='derived_regions', to='core.city'),
+        ),
+        migrations.AddField(
+            model_name='derivedregion',
+            name='source_regions',
+            field=models.ManyToManyField(related_name='derived_regions', to='core.region'),
+        ),
+        migrations.AddConstraint(
+            model_name='derivedregion',
+            constraint=models.UniqueConstraint(fields=('city', 'name'), name='unique-name-per-city'),
+        ),
+    ]

uvdat/core/models.py

@@ -46,12 +46,44 @@ class NetworkNode(models.Model):
 
 
 class Region(models.Model):
-    name = models.CharField(max_length=255, unique=True)
+    name = models.CharField(max_length=255)
     properties = models.JSONField(blank=True, null=True)
     boundary = geo_models.MultiPolygonField()
     city = models.ForeignKey(City, on_delete=models.CASCADE, related_name='regions')
     dataset = models.ForeignKey(Dataset, on_delete=models.CASCADE, related_name='regions')
 
+    class Meta:
+        constraints = [
+            models.UniqueConstraint(name='unique-name-per-dataset', fields=['dataset', 'name'])
+        ]
+
+
+class DerivedRegion(models.Model):
+    """A Region that's derived from other regions."""
+
+    class VectorOperation(models.TextChoices):
+        UNION = 'UNION', 'Union'
+        INTERSECTION = 'INTERSECTION', 'Intersection'
+
+    name = models.CharField(max_length=255)
+    city = models.ForeignKey(City, on_delete=models.CASCADE, related_name='derived_regions')
+    properties = models.JSONField(blank=True, null=True)
+    boundary = geo_models.MultiPolygonField()
+
+    # Data from the source regions
+    source_regions = models.ManyToManyField(Region, related_name='derived_regions')
+    source_operation = models.CharField(
+        max_length=max(len(choice[0]) for choice in VectorOperation.choices),
+        choices=VectorOperation.choices,
+    )
+
+    class Meta:
+        constraints = [
+            # We enforce name uniqueness across cities, since
+            # DerivedRegions can consist of regions from multiple datasets
+            models.UniqueConstraint(name='unique-name-per-city', fields=['city', 'name'])
+        ]
+
 
 class Chart(models.Model):
     name = models.CharField(max_length=255, unique=True)

uvdat/core/serializers.py

@@ -1,6 +1,9 @@
+import json
+
+from django.contrib.gis.serializers import geojson
 from rest_framework import serializers
 
-from uvdat.core.models import Chart, City, Dataset, NetworkNode, SimulationResult
+from uvdat.core.models import Chart, City, Dataset, DerivedRegion, NetworkNode, SimulationResult
 from uvdat.core.tasks.simulations import AVAILABLE_SIMULATIONS
 
 
@@ -59,4 +62,43 @@ def get_name(self, obj):
 
     class Meta:
         model = SimulationResult
+
+
+class RegionFeatureCollectionSerializer(geojson.Serializer):
+    # Override this method to ensure the pk field is a number instead of a string
+    def get_dump_object(self, obj):
+        val = super().get_dump_object(obj)
+        val["properties"]["id"] = int(val["properties"].pop("pk"))
+
+        return val
+
+
+class DerivedRegionListSerializer(serializers.ModelSerializer):
+    class Meta:
+        model = DerivedRegion
+        fields = ['id', 'name', 'city', 'properties', 'source_regions', 'source_operation']
+
+
+class DerivedRegionDetailSerializer(serializers.ModelSerializer):
+    class Meta:
+        model = DerivedRegion
         fields = '__all__'
+
+    boundary = serializers.SerializerMethodField()
+
+    def get_boundary(self, obj):
+        return json.loads(obj.boundary.geojson)
+
+
+class DerivedRegionCreationSerializer(serializers.ModelSerializer):
+    class Meta:
+        model = DerivedRegion
+        fields = [
+            'name',
+            'city',
+            'regions',
+            'operation',
+        ]
+
+    regions = serializers.ListField(child=serializers.IntegerField())
+    operation = serializers.ChoiceField(choices=DerivedRegion.VectorOperation.choices)

uvdat/core/tasks/conversion.py

@@ -5,7 +5,6 @@
 import zipfile
 
 from celery import shared_task
-from django.contrib.gis.geos import MultiPolygon, Polygon
 from django.core.files.base import ContentFile
 from geojson2vt import geojson2vt, vt2geojson
 import geopandas
@@ -14,8 +13,9 @@
 import rasterio
 import shapefile
 
-from uvdat.core.models import Dataset, Region
+from uvdat.core.models import Dataset
 from uvdat.core.tasks.networks import save_network_nodes
+from uvdat.core.tasks.regions import save_regions
 from uvdat.core.utils import add_styling
 
 
@@ -41,6 +41,7 @@ def convert_raw_data(dataset_id):
 
 
 def convert_cog(dataset):
+    """Saves the raster_file attribute of the dataset."""
     with tempfile.TemporaryDirectory() as temp_dir:
         raw_data_path = Path(temp_dir, 'raw_data.tiff')
         raster_path = Path(temp_dir, 'raster.tiff')
@@ -112,15 +113,17 @@ def convert_cog(dataset):
 
 
 def convert_geojson(dataset, geodata_path=None):
+    """Saves the vector_tiles_file and geodata_file attributes of the dataset."""
     with tempfile.TemporaryDirectory() as temp_dir:
         if geodata_path is None:
             geodata_path = Path(temp_dir, 'geo.json')
             with open(geodata_path, 'wb') as geodata_file:
                 original_data = dataset.raw_data_archive.open('rb').read()
                 original_data = json.loads(original_data)
                 original_data['features'] = add_styling(original_data['features'], dataset.style)
-                original_projection = original_data.get('crs').get('properties').get('name')
                 geodata_file.write(json.dumps(original_data).encode())
+
+            original_projection = original_data.get('crs', {}).get('properties', {}).get('name')
             if original_projection:
                 geodata = geopandas.read_file(geodata_path)
                 geodata = geodata.set_crs(original_projection, allow_override=True)
@@ -198,25 +201,3 @@ def convert_shape_file_archive(dataset):
         print(f'\t Shapefile to GeoJSON conversion complete for {dataset.name}.')
 
         convert_geojson(dataset, geodata_path=geodata_path)
-
-
-def save_regions(dataset):
-    dataset.regions.all().delete()
-    property_map = dataset.style.get('property_map')
-    name_property = property_map.get('name') if property_map else None
-
-    geodata = json.loads(dataset.geodata_file.read().decode())
-    for feature in geodata.get('features'):
-        geometry = feature.get('geometry')
-        properties = feature.get('properties')
-        polygons = [Polygon(p if len(p) != 1 else p[0]) for p in geometry.get('coordinates')]
-        region = Region(
-            name=properties.get(name_property) if name_property else "",
-            boundary=MultiPolygon(*polygons),
-            properties=properties,
-            dataset=dataset,
-            city=dataset.city,
-        )
-        region.save()
-
-    print(f"Saved regions for {dataset.name}")

uvdat/core/tasks/regions.py

@@ -0,0 +1,98 @@
+import json
+import secrets
+from typing import List
+
+from django.contrib.gis.db.models.aggregates import Union
+from django.contrib.gis.geos import GEOSGeometry
+from django.db import transaction
+
+from uvdat.core.models import DerivedRegion, Region
+
+
+class DerivedRegionCreationException(Exception):
+    pass
+
+
+def create_derived_region(name: str, city_id: int, region_ids: List[int], operation: str):
+    # Ensure at least two regions provided
+    source_regions = Region.objects.filter(pk__in=region_ids)
+    if source_regions.count() < 2:
+        raise DerivedRegionCreationException("Derived Regions must consist of multiple regions")
+
+    # Ensure all regions are from one city
+    source_cities = list((source_regions.values_list('city', flat=True).distinct()))
+    if len(source_cities) > 1:
+        raise DerivedRegionCreationException(
+            f"Multiple cities included in source regions: {source_cities}"
+        )
+
+    # Only handle union operations for now
+    if operation == DerivedRegion.VectorOperation.INTERSECTION:
+        raise DerivedRegionCreationException("Intersection Operation not yet supported")
+
+    # Simply include all multipolygons from all source regions
+    # Convert Polygon to MultiPolygon if necessary
+    geojson = json.loads(source_regions.aggregate(polys=Union('boundary'))['polys'].geojson)
+    if geojson['type'] == 'Polygon':
+        geojson['type'] = 'MultiPolygon'
+        geojson['coordinates'] = [geojson['coordinates']]
+
+    # Form proper Geometry object
+    new_boundary = GEOSGeometry(json.dumps((geojson)))
+
+    # Check for duplicate derived regions
+    existing = list(
+        DerivedRegion.objects.filter(city=city_id, boundary=GEOSGeometry(new_boundary)).values_list(
+            'id', flat=True
+        )
+    )
+    if existing:
+        raise DerivedRegionCreationException(
+            f"Derived Regions with identical boundary already exist: {existing}"
+        )
+
+    # Save and return
+    with transaction.atomic():
+        derived_region = DerivedRegion.objects.create(
+            name=name,
+            city=city_id,
+            properties={},
+            boundary=new_boundary,
+            source_operation=operation,
+        )
+        derived_region.source_regions.set(source_regions)
+
+    return derived_region
+
+
+def save_regions(dataset):
+    dataset.regions.all().delete()
+    property_map = dataset.style.get('property_map')
+    name_property = property_map.get('name') if property_map else None
+
+    geodata = json.loads(dataset.geodata_file.read().decode())
+    for feature in geodata['features']:
+        properties = feature['properties']
+        geometry = feature['geometry']
+
+        # Ensure a name field
+        name = secrets.token_hex(10)
+        if name_property and name_property in properties:
+            name = properties[name_property]
+
+        # Convert Polygon to MultiPolygon if necessary
+        if geometry['type'] == 'Polygon':
+            geometry['type'] = 'MultiPolygon'
+            geometry['coordinates'] = [geometry['coordinates']]
+
+        # Create region with properties and MultiPolygon
+        region = Region(
+            name=name,
+            boundary=GEOSGeometry(str(geometry)),
+            properties=properties,
+            dataset=dataset,
+            city=dataset.city,
+        )
+        region.save()
+
+    print(f"Saved regions for {dataset.name}")